*Editorial* **Complexity Science in Human Change: Research, Models, Clinical Applications**

**Franco Orsucci 1,2,\* and Wolfgang Tschacher <sup>3</sup>**


Complexity and entropy prevail in human behavior and social interaction because the systems underlying behavior and interaction are, without a doubt, highly complex. The human brain, body, language, society, and culture consist of vast numbers of components, and the degrees of freedom in behavior, cognition, and experience are just as immense. So why do we usually experience the world around us as structured and well-organized instead of disorganized and random? How do the patterns emerge? We are witnessing selforganization and pattern-formation processes, which organize and modulate complexity.

Increasingly, such processes are acknowledged as essential for human affairs and are gradually coming to the fore in psychology and social science research. Research informed by dynamical systems theory, synergetics, and complexity theory has introduced concepts such as attractor, synchrony, and coupling to psychology. In psychotherapy research, empirical findings show that regular patterns of interaction arise in all therapeutic relationships. The therapist–patient alliance is a paradigmatic case to highlight further how interactions evolve and can be changed and how humans can change. Attractors describe the stable states of a process, e.g., the stability or instability of personalities and disorders. They can be detected and described based on empirical time-series.

In the first paper, Orsucci [1] examines certain theoretical implications of empirical studies developed over recent years by his research groups. These experiments have explored the biosemiotic nature of communication streams from emotional neuroscience and embodied mind perspectives. Information combinatorics analysis enabled a deeper understanding of the coupling and decoupling dynamics of biosemiotics streams. They investigated intraindividual and interpersonal relations as the coevolution dynamics of hybrid couplings, synchronizations, and desynchronizations. Cluster analysis and Markov chains produced evidence of chimera states and phase transitions. A probabilistic and nondeterministic approach clarified the properties of these hybrid dynamics. As a result, multidimensional theoretical models can better represent the hybrid nature of human interactions.

In the second contribution, a study by Tomashin et al. [2], the authors consider how fractal properties in time series of human behavior and physiology are ubiquitous, and several methods to capture such properties have been proposed in the past decades. The paper takes this suggestion as a point of departure to propose and test several approaches to quantifying fractal fluctuations in synthetic and empirical time-series data using recurrencebased analysis. They show that such measures can be extracted based on recurrence plots and contrast the different approaches in terms of their accuracy and range of applicability.

In the third paper, Altmann et al. [3] compare eight algorithms that quantify synchronization in time series. The authors use a benchmark dataset that describes body movement in 30 dyadic interviews on somatic complaints conducted by medical students with 15 depressed and 15 healthy interviewees. Twenty-one different synchrony measures are tested,

**Citation:** Orsucci, F.; Tschacher, W. Complexity Science in Human Change: Research, Models, Clinical Applications. *Entropy* **2022**, *24*, 1670. https://doi.org/10.3390/e24111670

Received: 8 November 2022 Accepted: 8 November 2022 Published: 17 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

derived from four classes of algorithms: (windowed) cross-correlations, local regressions with or without peak-picking, mutual information, and cross-recurrence quantification. The intercorrelations of the results show that synchrony estimations are highly divergent, and no convergent validity is manifest. However, measures from the same class tend to be correlated, and cross-correlation-based measures form a factor. In contrast, the mutual information and the peak-picking measures load on a different factor. Most measures do not support the assumption underlying predictive validity that depression should have lowered synchrony measures. The authors conclude that more analyses and sensitivity studies are needed to clarify the psychological meaning of synchrony.

In the fourth contribution, research by Stamovlasis et al. [4], the authors investigate and propose a nonlinear model that might explain empirical data better than ordinary linear ones and elucidate the role of depression in a financial capacity. Financial incapacity is one of the cognitive deficits observed in amnestic mild cognitive impairment and dementia, while the combined interference of depression remains unexplored. Cusp catastrophe analysis was applied to the data, which suggested that the nonlinear model was superior to the linear and logistic alternatives, demonstrating that depression contributes to a bifurcation effect. Depressive symptomatology induces nonlinear effects, and a sudden decline in financial capacity is observed beyond a threshold value. Implications for theory and practice are discussed.

In the fifth contribution, Zubek et al. [5] reflected on the pandemic that forced our daily interactions to move into the virtual world. People had to adapt to new communication media that afford different ways of interaction. Remote communication decreases the availability and salience of some cues but also may enable and highlight others. Importantly, basic movement dynamics, which are crucial for any interaction as they are responsible for the informational and affective coupling, are affected. It is therefore essential to discover exactly how these dynamics change. In this exploratory study of six interacting dyads, they used tradi-tional variability measures and cross recurrence quantification analysis to compare the movement coordination dynamics in quasi-natural dialogues in four situations: (1) re-mote video-mediated conversations with a self-view mirror image present, (2) remote video-mediated conversations without a self-view, (3) face-to-face conversations with a self-view, and (4) face-to-face conversations without a self-view. They discovered that in remote interactions movements pertaining to communicative gestures were exagger-ated, while the stability of interpersonal coordination was greatly decreased. The presence of the self-view image made the gestures less exaggerated but did not affect the coordination. The dynamical analyses clarified the interaction processes and may be useful in explaining phenomena connected with video-mediated communication, such as "Zoom fatigue".

In the sixth contribution, Lauda´nska et al. [6] clarify how infants' limb movements evolve from disorganized to more selectively coordinated during the first year of life as they learn to navigate and interact with an ever-changing environment more efficiently. However, how these coordination patterns change during the first year of life and across different contexts is unknown. Here, they used wearable motion trackers to study the developmental changes in the complexity of limb movements (arms and legs) at 4, 6, 9, and 12 months of age in two different tasks: rhythmic rattle-shaking and free play. They applied multidimensional recurrence quantification analysis (MdRQA) to capture the nonlinear changes in infants' limb complexity. They show that the MdRQA parameters (entropy, recurrence rate, and mean line) are task-dependent only at 9 and 12 months of age, with higher values in rattle-shaking than free play. Infants' motor system becomes more stable and flexible with age, allowing for the flexible adaptation of behaviors to task demands.

The seventh contribution by Ganesh and Gabora [7] take a human dynamical systems approach to modeling therapeutic change, using reflexively autocatalytic food setderived (RAF) networks. RAFs have been used to model the self-organization of adaptive networks associated with the origin and early evolution of both biological life and the development of the kind of cognitive structure necessary for cultural evolution. The RAF approach is applicable in these seemingly disparate cases because it provides a

theoretical framework for formally describing under what conditions systems com-posed of elements that interact and "catalyze" the formation of new elements collec-tively become integrated wholes. This contribution develops in line with the growing recognition of the role of embodiment, affect, and implicit processes in psychotherapy, and several recent studies have examined the role of physiological synchrony in the process and outcome of psychotherapy. This study aims to introduce partial directed coherence (PDC) as a novel approach to calculating psychophysiological synchrony and examines its potential to contribute to our understanding of the therapy process. The study adopts a single-case, mixed-method design and examines physiological syn-chrony in one-couple therapy in relation to the therapeutic alliance and a narrative analysis of meaning construction in the sessions. The findings of this study point to the complex interplay between explicit and implicit levels of interaction and the potential contribution of including physiological synchrony in the study of interactional pro-cesses in psychotherapy.

The paper by Avdi et al. [8] addresses physiological synchrony in one family therapy of fifteen sessions and two physiological measurement sessions in which cardiac measures were recorded. The sessions concern a couple and two female psychotherapists. Physiological data are transformed into an index of sympathetic activity, and synchrony is computed using partial directed coherence. This method detects the direction of influence ("pacing/leading") in each pair of participants. In addition to the quantitative findings on synchrony, rating scales depict therapeutic alliance, and qualitative coding separates the measurement sessions into topical episodes that are semantically similar. Finally, the therapy process is described by the percentage of time windows synchronized concerning the couple's sympathetic activity, which is found to be reduced in the second measurement session, where the patterns of pacing and leading have changed towards a more balanced embodied relatedness. The authors conclude that a mixed-methods approach allows the linking of the quantitative synchrony findings to the qualitative clinical process in this successful couple therapy.

In the ninth paper, a study by Nkomidio et al. [9], the authors investigate the response characteristics of a two-dimensional neuron model exposed to an externally applied extremely low frequency (ELF) sinusoidal electric field and the synchronization of neurons weakly coupled with gap junction. They find, by numerical simulations, that neurons can exhibit different spiking patterns, which are well observed in the structure of the recurrence plot (R.P.). Then they further study the synchronization between weakly coupled neurons in chaotic regimes under the influence of a weak ELF electric field. In general, detecting the phases of chaotic spiky signals is not easy when using standard methods. Recurrence analysis provides a reliable tool for defining phases, even for noncoherent regimes or spiky signals. Recurrence-based synchronization analysis reveals that, even in the range of weak coupling, the phase synchronization of the coupled neurons occurs. By adding an ELF electric field, this synchronization increases depending on the amplitude of the externally applied ELF electric field. Authors further suggest a novel measure for RP-based phase synchronization analysis, which better considers the probabilities of recurrences.

In the tenth contribution, Webber [10] clarifies how the recurrence analyses of dynamical systems can only process the short sections of signals that may be infinitely long. By necessity, the recurrence plot and its quantifications are constrained within a truncated triangle that clips the signals at its borders. Recurrence variables defined within these confining borders can be influenced by truncation effects depending on the system under evaluation. In this study, the question being asked is, if the boundary borders were tilted, what would be the effect on all recurrence variables? This question is examined by comparing recurrence variables computed with the triangular recurrence area versus the boxed recurrence area. Examples include the logistic equation (mathematical series), the Dow Jones Industrial Average over a decade (real-world data), and a square wave pulse (toy series). Good agreement among the variables in terms of timing and amplitude was found for most, but not all, variables. These significant results are discussed.

In the eleventh paper, Ciompi and Tschacher [11] develop an account of modeling schizophrenia based on four different but related complexity theories: affect-logic, 4E-cognition/embodiment, synergetics, and the free-energy principle. All theories have in common that they are built on loop dynamics, so-called circular causality. In affect-logic, the loop is given by circular interactions between emotion ('affect') and cognition ('logic'), where emotion is the energy source for cognitive dynamics. Such interactions occur at the individual level, the level of micro-social interaction, and the societal level, which are structurally coupled. In synergetics, emotions act as control parameters that drive the system toward pattern formation. The embodiment and the free-energy principles are likewise built on circular dynamics between mind and body, respectively, between a generative model and sensory evidence. The article uses these commonalities for insights into the dynamics of schizophrenia spectrum symptoms: overly strong emotional tension then forces the cognitive system into dysfunctional patterns, which, however, are functional insofar as free energy is reduced. Ideas for therapeutic guidelines, as in the Soteria model, are also derived.

In the twelfth paper, Prinz et al. [12] compute the synchrony of electrodermal activity in psychotherapy interventions, focusing on the technique of imagery rescripting. This therapeutic technique was developed to modulate traumatic memories in a positive and desired direction. The activation of such memories commonly also affects the therapist involved in the session. Therefore, client-therapist synchrony based on cross-correlations is explored in 50 clients. Client-led synchrony is differentiated from therapist-led synchrony by the sign of the lags of cross-correlations. It is found that therapist-led synchrony is significantly associated with clients' emotional experiences of greater contentment, lower anxiety, and lower depression. In contrast, client-led synchrony is linked to the clients' more significant anxiety. The authors interpret their findings as supportive of the therapists' role in regulating mood.

In their contribution, Gennaro et al. [13] introduce a novel lexical method called the affective saturation index (ASI) to assess affectivity based on interview transcripts. Affect is semiotically defined as a sign that makes sense of the world in terms of patterns of bodily activation. Affect saturation in the ASI is then defined based on a "phase space of meaning". The ASI was correlated with several measures of semantic complexity, students' emotion regulation, and heart-rate variability in a sample of 40 students who participated in semi-structured interviews on neutral issues. The study shows that affective saturation is significantly and inversely linked to the semantic entropy index and heart-rate variability, consistent with the expectation that the ASI can detect the lexical-syntactic complexity of the interview text as well as physiological signs of affective arousal. It is concluded that ASI thus has potential applicability in clinical and community interventions, social communication, marketing, and media monitoring. Most approaches to computing interpersonal coupling are dyadic in that they focus on bivariate synchrony, such as that between client and therapist.

Meier and Tschacher [14] developed an algorithm for multivariate surrogate synchrony (mv-SUSY), which is based on the eigen-decomposition of the correlation matrix of multiple time series, like principal component analysis (mv-SUSY variant lambdamax). A further variant labeled omega is derived from the determinant of the variance-covariance matrix as a measure of actual entropy and standardized by potential entropy, the product of all variances. Computation is carried out in time-series segments, and segment-shuffled surrogate datasets are used as a control condition. The authors apply mv-SUSY to the simulated multivariate time series that realize the various types of regularities (random data, autocorrelations, trends, oscillatory behavior, intercorrelated random data) and to empirical multivariate time series (motion capture data from persons dancing and from a group discussion). It was found that mv-SUSY correctly identifies whether regular patterns exist in the datasets. The results of the multivariate algorithms are additionally validated by conventional dyadic synchrony methods.

We believe that change is generally studied in phase transitions when the dynamics move between different attractors, as evident in behaviors, mental states, and neurobiology. Theoretical models can represent dynamical change maps in mathematical equations and topological structures. Mapping theory to empirical research and vice versa is challenging but heuristic. Nevertheless, it paves the road to a future discipline of a general complexity theory of human change.

One feature of complexity and self-organization is the presence of scaling and fractal dynamics with the emergence of higher-order organizations. Moreover, heterogeneous human networks present specific kinds of self-similarity in the embodied mind in individual and social dynamics. Finally, translational processes and procedures from research to applications and vice versa are particularly relevant as they frequently include interdisciplinary collaborations.

Based on these thoughts, the Special Issue "Complexity Science in Human Change" has addressed an interdisciplinary community of scientists and practitioners interested in dynamical systems theory, especially approaches considering complex systems and applications to psychology and psychotherapy. Most of the contributions in this Issue analyze empirical data, predominantly time series. Some contributions contain theoretical models or methodological topics of complexity science.

This Special Issue closely represents the current work of complexity researchers in human behavior and change. Psychotherapy and communication systems are the backgrounds of six articles, mental health and psychopathology of four articles, and one paper concerns developmental psychology. Six papers put forward methods for the computation of interpersonal synchrony and innovations of existing synchrony algorithms. Regarding methodology, recurrence quantification analysis is frequently applied in articles on this topic. Four put forward correlation-based analyses of synchrony. Finally, network modeling, catastrophe theory, and nonlinear regression are tackled in one article.

This Special Issue highlights achievements of complexity science in studying patterns of organization and change in human dynamics. It also highlights new challenges that lie ahead. First, it clarifies how complex systems present plural structural forms and varieties of organization and disorganization [15,16]. These varieties are frequently distributed even within any singular system, creating rugged dynamical landscapes. This applies more specifically to human systems, which are hybrid by default [17]. They present multiple scales and heterogeneous subsystems; synchronous and asynchronous interactions; stable, unstable, and metastable states; and localized and generalized dynamics. Therefore, considering the hyper-complexity of the human dynamical landscapes, empirical studies can use mixed methods with several different approaches (sometimes all at once). Accordingly, multiple and varied interventions can induce change or facilitate its natural evolution and emersion [18,19]. This can explain how multiple and different therapeutic techniques can produce similar (though not identical) outcomes in the clinical field. It is the good old equifinality principle of complex open systems still at work.

We hope that with this, new questions and new research might be incited, as Voltaire once suggested: "The most useful books are those in which the readers themselves supply half of the meaning".

**Author Contributions:** Writing—original draft preparation, F.O. and W.T.; writing—review and editing, F.O. and W.T. All authors have read and agreed to the published version of the manuscript.

**Acknowledgments:** We express our thanks to the authors of the above contributions, and the journal Entropy and MDPI for their support during this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Perspective* **Human Synchronization Maps—The Hybrid Consciousness of the Embodied Mind**

**Franco Orsucci 1,2**

<sup>1</sup> Department of Psychology, University College London, London WC1E 6BT, UK; f.orsucci@ucl.ac.uk

<sup>2</sup> Norfolk and Suffolk NHS Foundation Trust, Drayton High Road, Norwich NR6 5BE, UK

**Abstract:** We examine the theoretical implications of empirical studies developed over recent years. These experiments have explored the biosemiotic nature of communication streams from emotional neuroscience and embodied mind perspectives. Information combinatorics analysis enabled a deeper understanding of the coupling and decoupling dynamics of biosemiotics streams. We investigated intraindividual and interpersonal relations as coevolution dynamics of hybrid couplings, synchronizations, and desynchronizations. Cluster analysis and Markov chains produced evidence of chimaera states and phase transitions. A probabilistic and nondeterministic approach clarified the properties of these hybrid dynamics. Thus, multidimensional theoretical models can represent the hybrid nature of human interactions.

**Keywords:** synchronization; semiotics; information; cognitive neuroscience; psychotherapy; conversation; mapping; chimaera states; statistical dynamics; coupling

Science is built up with facts, as a house is with stones.

However, a collection of facts is no more a science than a heap of stones is a house.

*Henri Poincaré, Science and Hypothesis*

## **1. Introduction. Complexity, Noise, and Orders**

We will try to expand some theoretical outcomes of empirical and experimental research on human interactions published by our laboratories in recent years. We built an advanced multidimensional methodology for analyzing human dynamics, mainly focusing on synchronization in an embodied mind framework [1,2]. Patterns of synchronization form the foundations of the cognition [3,4] continuum between healthy and disease states [5]. Structural coupling and synchronization arise in human dynamics in many ways, including coordination in conversations: speech, movement, emotions, and physiology [6–12]. It is a partially self-contained setting and practices to observe and facilitate transformation in human conditions and relations. Psychotherapy has been described as "one of the most complex bio-psycho-social systems in which patterns of language, cognition, emotion, and behavior are formed and changed through the dynamics of therapist and patient interactions" [13,14]. Beyond clinical research, studying such an exceptional human dynamics environment can lead to a general model of human dynamics, comprehending the linguistic, behavioral, and physiological realms. The integration of communication, action, bodies, and environments highlights our embodied interactions' multimodality and parallel multiactivity [15,16].

We started by focusing our studies on language. Language study is scaled in complex structures: from informational systems to mesoscopic morphological patterns to semantic and narrative streams. In verbal interactions, voice tonality, volume, pitch, cadence, rhythm, and turn-taking are relevant. Shannon [17] built the foundations of the information theory of texts and speech. His less famous work on the prediction and entropy of printed English [18] is a resource for inspiring new research. It might be interesting to consider

**Citation:** Orsucci, F. Human Synchronization Maps—The Hybrid Consciousness of the Embodied Mind. *Entropy* **2021**, *23*, 1569. https://doi.org/10.3390/ e23121569

Academic Editor: Milan Paluš

Received: 11 October 2021 Accepted: 18 November 2021 Published: 25 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the distribution of information and organization in different living and nonliving systems in the same perspective. In this perspective, a graph proposed by Schreiber [19] mapped scattered areas of various forms of order, entropy and knowledge still interspersed with regions of the unknown, as in old charts. Following his mapping, we can find periodic and noisy oscillations, deterministic and stochastic areas of chaos, stochastic resonance, self-organized criticality, nonlinearity, or noise. Then, there are a few other islands where a connection between our knowledge models and real-world phenomena is yet to be well established. This kind of dynamical mapping might be synchronic and diachronic, in spatial distribution and time transitions.

The structure of different systems can be known and modified through the emergence of self-organization or by external actions, by casual or planned perturbations, including measurements. Some interactions can lead to coupling between systems, and if they repeat in time, they might produce forms of synchronization. Maturana and Varela [20] considered synchronization a form of structural coupling occurring when two systems repeatedly perturb each other. "Synchronization is a nonlinear phenomenon discovered at the beginning of the scientific revolution", and in its classical definition, synchronization refers to adjustment or entrainment in frequencies or phase of periodic oscillators due to weak interactions that lead to structural coordination between systems" [21]. This process can lead to the emergence of adaptive behavior between interacting systems. Pecora and Carroll [22], Ott, Grebogi and Yorke [23], and Pyragas [24] found that synchronization can be used to change the dynamic behavior of complex systems.

#### **2. Materials and Methods. Biosemiotics Pattern Analysis**

Our initial approach was different from most of the studies mentioned above. We chose a method, Recurrence Quantification Analysis–RQA [25,26], that does not generate any specific hypothesis on the form of data and does not need to consider time series produced by a dynamic system. Our primary aim was to build a statistical tool for reliable quantitative measures of the degree of organization (as expressed by the recurrence of patterns) of a flow of signs. We demonstrated how this could be performed with a relatively simple mathematical model. The analysis of the informational structure of a text (irrespective of its meaning) could unveil the hidden matching of patterns between two speaking persons. The hidden matching relates to the flow and forms of information linking partners in conversation. Through the phonetic configuration of speech, as represented in orthography, we can extract relevant patterns in the dynamic structures of human interactions.

We used Recurrence Quantification Analysis, a methodology that can reliably measure the recurrence of patterns, determinism, and entropy. Recurrence Plots (RP) were first pioneered in physics by Eckmann, Kamphorst and Ruelle [25,26]. Later, Webber, Zbilut, Giuliani and Marwan augmented this technique by identifying nonlinear variables for the quantitative assessment of RPs, thus creating RQA. Since then, RQA has been used in different areas, from molecular dynamics [27,28] to physiology [29] and bioinformatics [30–32]. In performing RQA, the original time series must be placed into an embedding matrix by converting the original *n* elements column vector correspondent to the symbol series into a *p*-dimensional matrix with columns as the original *X<sup>n</sup>* series plus its lagged copies *Xn+*1, *Xn+*2, . . . , *Xn+*1, while *p* is the embedding dimension. The quantification of recurrences is acquired by many different 'counts' of repetitions within the matrix.

While testing the robustness of this methodology on written language [33], we had to set to three (letters) for dimensional embedding, as this amplifies its sensitivity while avoiding noise from low-level statistical features (for example, asymmetrical distribution of couplets of letters). We might notice that a three-letter dimension represents a mesoscopic information level in natural language, just between single letters and whole words. We will later see the theoretical implications of this seemingly technical specification. Our time series analysis used RQA and CRQA (cross recurrence) to measure the synchronization in conversations as semiotic interactions. These informational patterns represent a preverbal and a-conscious communication channel revealed by the frequent emergence of patterns

of prosodic structures (such as the musicality of phoneme sequences, stereotyped words, pauses and phrases).

Other independent centers started developing research on social and clinical interactions using a similar methodology based on recurrence analysis. For example, they studied postural or verbal time series of interpersonal coordination during conversations [34–36]. These studies usually took one type of time series (i.e., movement, speech, or physiology) while not considering the mutual influence between different kinds of interaction. However, as human relationship dynamics are naturally hybrid, one type of interaction can influence the coupling or uncoupling of the other streams: motor, semiotic or physiological.

#### **3. Results. Hybrid Couplings and Synchronizations**

Human interactions constantly involve multiple streams (language, movement, emotions) which undergo coupling, decoupling and synchronizing. These multiscale and hybrid interactions are better comprehended within the biosemiotics, embodied mind framework that we defined as Mind Force [37,38]. We built the empirical paradigm of this approach as a multidimensional analysis of speech and emotions in patients and therapists in psychotherapy [39]. We chose Galvanic Skin Response—GSR and verbal prosody, as both variables reveal, in different flows, the expression of emotions [40,41]. Our new experiments studied four signals: the therapist's speech transcription, the patient's speech transcription, the therapist's GSR, and the patient's GSR. We focused on how those four variables modulated, coupled, synchronized, or desynchronized with each other. First, we considered the combinatorics and patterns of letters (phonemes) and morphemes (the minor portion of words that communicate significance). As mentioned, we had established this methodology in previous studies, which validated robust informational measures of entropy and determinism. In this new study, we initially considered the synchronization with standard correlation coefficients of Principal Components Analysis. Afterwards, we clustered all four signals using k-means resulting in a model representing this complex system's phase space and state transitions. Then, using a Markov Transition Matrix (see Figure 1), we disclosed phase transition probabilities between linguistic and physiological time series. *Entropy* **2021**, *23*, x FOR PEER REVIEW 4 of 8

**Figure 1.** The Markov Transition Matrix map of attractors and phase transitions (Orsucci et al., 2016). Greek creature made up of parts of different animals and introduced some theoretical **Figure 1.** The Markov Transition Matrix map of attractors and phase transitions [1].

that can change over time. The Japanese physicist Yoshiki Kuramoto (1984b, 1984a) proposed a paradigmatic mathematical model to describe synchronization dynamics in a large set of coupled oscillators. The most frequent form of the model has the following

൫ − ൯, ൌ 1 … . , (1)

ሻ sinሺሺ. ሻ − ሺᇱ

, ሻ + ሻ <sup>ᇱ</sup>

 

with *ω*(*x*) = *ω* for all *x*.

ൌ +

ሺ, ሻ ൌ ሺሻ −නሺ−ᇱ

 

ൌ1

where the system is formed of *N* limit-cycle oscillators with phase *θi* and coupling *K*.

Then, in November 2002, Yoshiki Kuramoto and Dorjsuren Battogtokh published the paper "Coexistence of Coherence and Incoherence in Nonlocally Coupled Phase Oscillators" [42,43]. They observed the coexistence of coherence and incoherence in a network of identical, nonlocally coupled, complex Ginzburg–Landau oscillators. While coupled nonidentical oscillators were known to exhibit mixed complex behavior (frequency locking, phase synchronization, partial synchronization, and incoherence), identical oscillators were supposed to either synchronize in phase or incoherently drift. They showed that oscillators that were identically coupled with similar natural frequencies could behave differently from one another for specific initial conditions. Some could synchronize while others remained incoherent in a stable state. They considered the following equation, which they called the nonlocally coupled complex Ginzburg–Landau equation:

Later, Abrams and Strogatz [44,45] named it a chimaera state, from the mythological

equation:

**4. Discussion. The Chimera States in Human Interactions** 

The complex dyadic system evolves between two attractors. In the first attractor, state four, the therapist strives to attune and entrain with the patient presenting low values of GSR recurrence and determinism. The therapist has high recurrence and determinism in prosody with repetitive semiotic patterns, perhaps to direct the patient's emotional expressions. The second attractor, at state five, is characterized by a medium level of GSR recurrence and determinism for both patient and therapist. We evidence semiotic medium recurrence and determinism for the therapist and low recurrence and determinism for the patient. Overall, this phase represents a state in which the patient's physiological anxiety becomes more manageable and linguistic expressions are more integrated. In short, while state four is an erratic phase of the interaction in which semiotics seems independent from passions, state five shows an integration. This sequence in human interactions is consistent with the literature in psychotherapy and neuroscience research on the embodied mind.

#### **4. Discussion. The Chimera States in Human Interactions**

This data analysis and mapping highlight dynamical landscapes of mixed states of coupling, with mixed zones of synchronization, noninteraction, and drift in uncoupling that can change over time. The Japanese physicist Yoshiki Kuramoto (1984b, 1984a) proposed a paradigmatic mathematical model to describe synchronization dynamics in a large set of coupled oscillators. The most frequent form of the model has the following equation:

$$\frac{d\theta\_i}{dt} = \omega\_i + \frac{K}{N} \sum\_{j=1}^{N} \sin\left(\theta\_j - \theta\_i\right), \ i = 1 \ldots N,\tag{1}$$

where the system is formed of *N* limit-cycle oscillators with phase *θ<sup>i</sup>* and coupling *K*.

Then, in November 2002, Yoshiki Kuramoto and Dorjsuren Battogtokh published the paper "Coexistence of Coherence and Incoherence in Nonlocally Coupled Phase Oscillators" [42,43]. They observed the coexistence of coherence and incoherence in a network of identical, nonlocally coupled, complex Ginzburg–Landau oscillators. While coupled nonidentical oscillators were known to exhibit mixed complex behavior (frequency locking, phase synchronization, partial synchronization, and incoherence), identical oscillators were supposed to either synchronize in phase or incoherently drift. They showed that oscillators that were identically coupled with similar natural frequencies could behave differently from one another for specific initial conditions. Some could synchronize while others remained incoherent in a stable state. They considered the following equation, which they called the nonlocally coupled complex Ginzburg–Landau equation:

$$\frac{\delta}{\delta t}\psi(\mathbf{x},t) = \omega(\mathbf{x}) - \int G(\mathbf{x}-\mathbf{x}')\sin\left(\psi(\mathbf{x},t) - \psi(\mathbf{x}',t) + a\right)d\mathbf{x}'$$

with *ω*(*x*) = *ω* for all *x*.

Later, Abrams and Strogatz [44,45] named it a chimaera state, from the mythological Greek creature made up of parts of different animals and introduced some theoretical clarifications for such behavior. Finally, they studied the most straightforward system presenting a chimaera state, a ring of phase oscillators governed by:

$$\frac{\partial \Phi}{\partial t} \omega - \int\_{-\pi}^{\pi} G(\mathbf{x} - \mathbf{x}') \sin \left[ \phi(\mathbf{x}, \mathbf{x}t) - \phi(\mathbf{x}', t) + \alpha \right] d\mathbf{x}' \alpha$$

Here, *φ*(*x* 0 , *t*) is the phase of the oscillator at position *x* at time *t.* The space variable *x* runs from −π to π with periodic boundary conditions. The frequency *ω* plays no role in the dynamics; one can set *ω* = 0 by redefining *φ* → *φ* + *ωt* without otherwise changing the form of the equation.

Chimaera states were later found in limit-cycle oscillators, chaotic oscillators, chaotic maps and in neuronal systems. In the beginning, chimaera patterns were observed in nonlocally coupled networks, but afterwards, these states were also found globally and locally (nearest neighbor) coupled networks and in modular networks [46,47]. The usage of Markov chains for mapping couplings and chimaera states was also explored ([48,49]. C.R. Laing studied chimaera state in heterogeneous networks, analyzing the influence of heterogeneous coupling strengths. Of further interest for human dynamics is the emergence of chimaera states in multiscale networks that result from the networking of different networks [50,51]. The ubiquity of chimaera mapping of synchronization and its different typologies extended its original definition to areas that might include nonidentical coupling oscillators in hybrid networks and multiscale networking of networks that were already known to present chimeralike dynamics before this definition started to be used.

Our studies' dynamic mapping of heterogeneous synchronization indicates that similar dynamics involving different brain areas related to emotional, motor, and verbal interactions co-occur. Cognitive tasks constantly require a balance between segregated and integrated neural processing with relevant consequences for cognitive performance. Segregation enables efficient computations in specialized brain regions, while integrated systems ensure coordinated, robust performance. Focused states tend to involve shorter, local connections, while integration largely relies on subcortical regions and cortical hubs with diverse connections to other brain regions [52,53]. "Recognizing chimaera dynamics can help to clarify the hybrid complexity of synchronization in critical cognitive states where a balance between integration and segregation is required for adaptive cognition and social interactions" [54]. Brain chimaera dynamics might also be related to different neuronal interactions mediated by different electrical or chemical synapses in the nervous system.

Further neural interactions involve neuromodulators and hormones, faster or slower action, and different time frames [55]. Various types of neural interaction are undoubtedly an additional factor in the emergence of chimaera dynamically states in human hybrid synchronization [51]. As separate regions interact to perform neurocognitive tasks, variable patterns of partial synchrony form chimaera states [3].

#### **5. Conclusions: From Determinism to Statistical Dynamics**

Human dynamics are so complex and prone to indeterminacy and randomness that even deterministic chaos might be considered, in many cases, as a reductionist simplification. Therefore, we might consider probabilistic models including elements of randomness. The previous study highlighted how biopsychosocial dynamics are hybrid, discontinuous, and have many degrees of freedom. We also highlighted how cluster analysis and Markov states could help to clarify the dynamics. However, our knowledge of the state of the systems is always incomplete, and some uncertainty is part of the game. While standard dynamics usually consider the behavior of a single state, statistical dynamics define the statistical ensemble as a probability cloud of the possible conditions in the system [56]. The ensemble probability can be interpreted in two main ways:


Following this perspective, we used a probabilistic approach in the study of empathy [57]. Empathy plays a significant role in human coordination, collaboration, and change, in human interactions. Most authors agree that forms of resonance in imitation, emotions, and communication are relevant factors of empathy. Following the expanding literature on relational physiology, we explored if empathy would present physiological evidence. We applied a Principal Component Analysis (PCA) on simultaneous GSR and HR signals from a patient-therapist dyad. PCA revealed a 'shared' component in signals, and two 'individual' components of independent correlation. Regression analysis showed that the shared component predicted a therapy outcome (R<sup>2</sup> = 0.28). We further examined the common component dynamics in a symbolic Markovian discrete model and cluster analysis.

Several studies on cognitive neuroscience [58,59] established statistical dynamics in biological systems focusing on the reciprocal correlations between system descriptors. This scientific position focuses on the mesoscopic level [60,61], potentially expanding correlations among system variables. This is the midpoint between pure "bottom-up" and "top-down" approaches. The crucial role of mesoscopic dynamics was validated in our physiological analysis and semiotics, as highlighted by the robust evidence for our mesoscopic embedding in RQA since the first experiments. In this way, we focused on the level of morphemes as word subcomponents. Morphemes have meaning and grammatical functions. They can be decomposed into smaller morphemes without losing these two crucial properties. Morphemes can be considered as semiotic quanta of information in natural language, as they are the basic lexical item in a language. They are usually composed of more than one phoneme and several letters or informational units [62–64]. Therefore, we can consider morphemes as the information quanta structuring coupling and synchronization in natural human language.

We explored informational patterns in human interactions. We investigated intraindividual and interpersonal relations as coevolution dynamics of hybrid couplings, synchronizations, and desynchronization. Cluster analyses and Markov chains produced evidence of chimaera states and phase transitions. A probabilistic and nondeterministic approach can clarify relevant properties of human dynamics, focusing on the mesoscopic scale and statistical dynamics. Theoretical models of human interactions should be founded on the hybrid nature of human structural couplings.

**Funding:** Parts of this research have been funded by: 2012–2014 Grant, MURST (Italian Ministry for University, Research & Technology) with ENEA, Rome, Italy. 2005 Invited Lecturer fund, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany. 1999–2001 Research Project, EU ICT 11696, European Commission, Brussels, Belgium.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board (or Ethics Committee) of the University of Siena and of ENEA, Rome, Italy.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Conflicts of Interest:** Parts of this research have been presented at the 8th Experimental Chaos Conference (Florence, June 2004), at the Tandem Workshop sponsored by the Max Planck Institute for Human Cognitive and Brain Sciences, University of Potsdam (Berlin, March 2005), at the International Conferences of the Society for Psychotherapy Research (Jerusalem 2016, Oxford 2017, Amsterdam 2018).

#### **References**


## *Article* **Four Methods to Distinguish between Fractal Dimensions in Time Series through Recurrence Quantification Analysis**

**Alon Tomashin <sup>1</sup> , Giuseppe Leonardi <sup>2</sup> and Sebastian Wallot 3,4,\***


**Abstract:** Fractal properties in time series of human behavior and physiology are quite ubiquitous, and several methods to capture such properties have been proposed in the past decades. Fractal properties are marked by similarities in statistical characteristics over time and space, and it has been suggested that such properties can be well-captured through recurrence quantification analysis. However, no methods to capture fractal fluctuations by means of recurrence-based methods have been developed yet. The present paper takes this suggestion as a point of departure to propose and test several approaches to quantifying fractal fluctuations in synthetic and empirical time-series data using recurrence-based analysis. We show that such measures can be extracted based on recurrence plots, and contrast the different approaches in terms of their accuracy and range of applicability.

**Keywords:** recurrence quantification analysis; fractals; monofractals; fractal time series

Wallot, S. Four Methods to Distinguish between Fractal Dimensions in Time Series through Recurrence Quantification Analysis. *Entropy* **2022**, *24*, 1314. https:// doi.org/10.3390/e24091314

**Citation:** Tomashin, A.; Leonardi, G.;

Academic Editors: Franco Orsucci and Wolfgang Tschacher

Received: 19 July 2022 Accepted: 13 September 2022 Published: 19 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

## **1. Introduction**

Since Gilden et al.'s [1] seminal paper, showing the presence of 1/*f* <sup>α</sup> -fluctuations in human time estimation performance, a huge interest in the presence and meaning of fractal fluctuations in human behavior has emerged. On the one hand, fractal patterns have been found in virtually all aspects of human physiology and behavior across recent studies [2–11]. On the other hand, their meaning has been intensely discussed [12–19].

Through the same period, the development and refinement of different time-series analysis techniques gained momentum, so that fractal properties could be quantified with a variety of methods, based on the power spectrum of a time series [20], their standard deviation [21] or residual fluctuations [22]—each of which has particular advantages and downsides, as well as requirements for preprocessing [21,23]. This was of central importance, because methods that are suitable for special fractals, such as box counting, are not equally applicable to time-series data [24].

In the current paper, we want to present another way of quantifying fractal fluctuations in time-series data using recurrence quantification analysis [25,26]. Our motivation for the present work is two-fold: firstly, to extend the use of recurrence plot-based methods to capture fractal properties. This is something that recurrence plot-based analyses have not been capable of. Further, to pave the way to provide an easy-to-use tool to compare fractal dimensions of time series that are well-applicable to binary data, and in the future also to multidimensional time series using multidimensional recurrence plot methods [27]. As has been suggested elsewhere [28], fractal properties in time-series data can be wellcaptured by the concept of (imperfectly) recurring patterns over time, and this is—as the name implies—what recurrence quantification analysis is about. Specifically, Webber [28] encouraged researchers to explore RQA as a bridge to further understand fractal systems

in various fields. However, Webber did not specify how to quantify fractal fluctuations by means of recurrence plots.

Hence, the aim of the present paper is to take this next step and propose, as well as compare, novel recurrence-based approaches that can be used to quantify fractal fluctuations. In the following, each approach is introduced, tested on synthetic data, and evaluated; in addition, a Matlab (The MathWorks, Inc., Natick, MA, USA) implementation of the approaches presented in this paper is available on GitHub: https://github.com/alontom/FARQA (Accessed on 19 July 2022, see Appendix A). Finally, we discuss the individual strengths and weaknesses of each approach and relate the results to those obtained from detrended fluctuation analysis (DFA; [22]), as DFA is one of the most widely used methods with accurate performance in capturing fractal fluctuations in time-series data [29–31].

#### **2. Methods and Results**

#### *2.1. Synthetic Data*

In this section, we show four new approaches to differentiate between the power-law scaling exponent (α; 1/*f* <sup>α</sup>) based on several RQA properties. Each evaluation approach was applied to synthetic data consisting of 1026 data points with different fractal dimensions ranging from α = −1 (antipersistent) to 2 (persistent) generated by 'power noise' function [32] using Matlab version 2021b (The MathWorks, Inc.). For every fractal dimension, 100 time series were generated under two conditions: idealized fractal time series and a noisy fractal time series (SNR 2:1). The noise component added was drawn from a normal distribution with 50% of the SD of the idealized fractal time series. We conducted RQA without embedding (delay and embedding parameters of 1, euclidean normalization of the phase space, and radius = 0.4) on the z-scored generated time series and utilized its properties to discriminate between signals with different 1/*f* values. As a benchmark to compare against, alongside the true predetermined α, we also subjected the data to detrended fluctuation analysis (DFA; [22]). In the next sections, we will describe each of the methods and present the results of their application. After that, we will apply the methods to empirical data of a time-estimation task. Finally, we will provide a summary of the strengths and weaknesses of each method and the intercorrelations of their results.

#### 2.1.1. Detrended Fluctuations Analysis (DFA)

First, we tested the fractal properties of the dataset by applying a detrended fluctuation analysis [22]. To do so, we used the following DFA parameters: a minimum bin size of 10, a maximum of 510, linear detrending. The results are presented in Figure 1 and show that the Hurst exponents *H* estimated via the DFA scale well with the true α-values of the time series. In the absence of random noise, DFA distinguishes scaling relations well down to antipersistent fluctuations with α = −1 (Figure 1, left panel, *R* <sup>2</sup> = 0.997). When noise is added, the capacity of DFA to distinguish among antipersistent was slightly compromised (Figure 1, right panel, *R* <sup>2</sup> = 0.965).

2.1.2. First Approach: Estimating Scaling Using the *SD* of *%REC* over a Range of Bin Sizes (*%REC SD*)

To capture the change in fluctuations with scale, the RP was split into bins of various sizes (powers of two). In each, we calculated the recurrence rate. Then, the *SD* of all the bins of the same size was computed, and we fitted a linear line to the log–log plot of the *SD* vs. the bin sizes. Figure 2 illustrates the approach.

The rationale behind this approach is that a time series of i.i.d. white noise will yield a recurrence plot that is statistically uniformly populated by mostly isolated recurrence points, while the correlation structure of persistent fluctuations will yield a more clustered, nonuniform distribution, and will hence lead to a slower increase in *SD* compared to the white noise case (Figure 3, α > 0). However, antipersistent fluctuations tend to systematically decluster recurrences, and the result is likewise a relatively uniformly distributed recurrence plot (Figure 3, α = −1).

**Figure 1.** Detrended fluctuation analysis (DFA) results: Left panel: Box plots of the true α-values on the *x*-axis and the estimated Hurst exponents *H* on the *y*-axis from DFA. As can be seen, the DFA *H* scales well with the true alpha values down to antipersistent fluctuations (α = −1). Right panel: Box plots of the true α-values on the *x*-axis and the estimated Hurst exponents *H* on the *y*-axis from DFA, when random noise is added (SNR = 2:1). DFA still scales well for persistent fluctuations with the true α-values, but is relatively less sensitive to distinguishing between different types of antipersistent fluctuations. **Figure 1.** Detrended fluctuation analysis (DFA) results: Left panel: Box plots of the true α-values on the *x*-axis and the estimated Hurst exponents *H* on the *y*-axis from DFA. As can be seen, the DFA *H* scales well with the true alpha values down to antipersistent fluctuations (α = −1). Right panel: Box plots of the true α-values on the *x*-axis and the estimated Hurst exponents *H* on the *y*-axis from DFA, when random noise is added (SNR = 2:1). DFA still scales well for persistent fluctuations with the true α-values, but is relatively less sensitive to distinguishing between different types of antipersistent fluctuations. *Entropy* **2022**, *24*, x FOR PEER REVIEW 4 of 16

**Figure 2.** Demonstration of approach 1 over a simple recurrence plot: (**A**) A hypothetic 8 × 8 recurrence plot (RP) where blue squares stand for recurrence points and blank squares for nonrecurrent ones. In (**B**,**C**) the RP is split into bins of 2 and 4 (respectively, marked in brown). With approach 1, one finds the *%REC* in every bin and computes the *SD* between the recurrence percentages. Afterward, a linear trend is fitted to the log–log scaling plot (**D**) and the slope represents the scaling. **Figure 2.** Demonstration of approach 1 over a simple recurrence plot: (**A**) A hypothetic 8 × 8 recurrence plot (RP) where blue squares stand for recurrence points and blank squares for nonrecurrent ones. In (**B**,**C**) the RP is split into bins of 2 and 4 (respectively, marked in brown). With approach 1, one finds the *%REC* in every bin and computes the *SD* between the recurrence percentages. Afterward, a linear trend is fitted to the log–log scaling plot (**D**) and the slope represents the scaling.

The rationale behind this approach is that a time series of i.i.d. white noise will yield a recurrence plot that is statistically uniformly populated by mostly isolated recurrence

to the white noise case (Figure 3, α > 0). However, antipersistent fluctuations tend to systematically decluster recurrences, and the result is likewise a relatively uniformly distrib-

uted recurrence plot (Figure 3, α = −1).

*Entropy* **2022**, *24*, x FOR PEER REVIEW 5 of 16

**Figure 3.** RP and scaling plots for different alpha values: Examples of univariate RP time series generated with different α-values, and scaling plots demonstrate the association of bin sizes and the *SD* of the recurrence rates between bins. As can be seen, from α = 0 the slope tends to decrease, suggesting a lower fractal dimension (i.e., higher α). **Figure 3.** RP and scaling plots for different alpha values: Examples of univariate RP time series generated with different α-values, and scaling plots demonstrate the association of bin sizes and the *SD* of the recurrence rates between bins. As can be seen, from α = 0 the slope tends to decrease, suggesting a lower fractal dimension (i.e., higher α). **Figure 3.** RP and scaling plots for different alpha values: Examples of univariate RP time series generated with different α-values, and scaling plots demonstrate the association of bin sizes and the *SD* of the recurrence rates between bins. As can be seen, from α = 0 the slope tends to decrease, suggesting a lower fractal dimension (i.e., higher α).

Figure 4 shows the model coefficients for each α (*n* = 100 simulations each). This method seems appropriate for distinguishing among persistent signals (α > 0) for the idealized, but also noise data. It does not work for antipersistent fluctuations (α < 0). Here, the method simply does not distinguish between α-values of 0 and −1. For the noisy time series, we fitted linear regressions between α values and the power-law coefficients separately for the α ≥ 0 (*R*2 = 0.9) and α ≤ 0 (*R*2 = 0.04), which support the above statement. Figure 4 shows the model coefficients for each α (*n* = 100 simulations each). This method seems appropriate for distinguishing among persistent signals (α > 0) for the idealized, but also noise data. It does not work for antipersistent fluctuations (α < 0). Here, the method simply does not distinguish between α-values of 0 and −1. For the noisy time series, we fitted linear regressions between α values and the power-law coefficients separately for the α ≥ 0 (*R* <sup>2</sup> = 0.9) and <sup>α</sup> <sup>≤</sup> 0 (*<sup>R</sup>* <sup>2</sup> = 0.04), which support the above statement. Figure 4 shows the model coefficients for each α (*n* = 100 simulations each). This method seems appropriate for distinguishing among persistent signals (α > 0) for the idealized, but also noise data. It does not work for antipersistent fluctuations (α < 0). Here, the method simply does not distinguish between α-values of 0 and −1. For the noisy time series, we fitted linear regressions between α values and the power-law coefficients separately for the α ≥ 0 (*R*2 = 0.9) and α ≤ 0 (*R*2 = 0.04), which support the above statement.

sizes on the *y*-axis. As observed, the coefficient scales well with the true alpha values for the persistent fluctuations (α > 0). Right panel: Box plots of the true alpha values on the *x*-axis and the power-law coefficients for the association of *SD* of *%REC* between the bins and bin sizes on the *y*-axis, when random noise (SNR 2:1) is added. Still, the resulted coefficients scale well for persistent fluctuations with the true alpha values but are relatively insensitive to distinguishing between different types of antipersistent fluctuations. **Figure 4.** Results of approach 1 (*SD %REC*): Left panel: Box plots of the true alpha values on the *x*axis and the power-law coefficients for the association of *SD* of *%REC* between the bins and bin sizes on the *y*-axis. As observed, the coefficient scales well with the true alpha values for the persistent fluctuations (α > 0). Right panel: Box plots of the true alpha values on the *x*-axis and the powerlaw coefficients for the association of *SD* of *%REC* between the bins and bin sizes on the *y*-axis, when random noise (SNR 2:1) is added. Still, the resulted coefficients scale well for persistent fluctuations with the true alpha values but are relatively insensitive to distinguishing between different types of

#### 2.1.3. Second Approach: Estimating Scaling Using Laminarity *(%LAM)*

*Entropy* **2022**, *24*, x FOR PEER REVIEW 6 of 16

For this approach, we simply calculated the percentage of recurrence points that have a vertical/horizontal neighbor (*%LAM*, laminarity; [33]) over the whole plot (Figure 5). The rationale behind this approach is somewhat similar to the first approach, which is that fractal fluctuations tend to be manifested by patches or squares in the recurrence plot (see Figure 3). Hence, *%LAM* would represent the persistence of the data well. 2.1.3. Second Approach: Estimating Scaling Using Laminarity *(%LAM)* For this approach, we simply calculated the percentage of recurrence points that have a vertical/horizontal neighbor (*%LAM*, laminarity; [33]) over the whole plot (Figure 5). The rationale behind this approach is somewhat similar to the first approach, which is that fractal fluctuations tend to be manifested by patches or squares in the recurrence plot (see Figure 3). Hence, *%LAM* would represent the persistence of the data well.

antipersistent fluctuations.

**Figure 5.** Quantifying laminarity: An 8 × 8 RP where colored squares represent recurrence points. The orange-filled recurrence points have a vertical/horizontal neighbor, while the blue squares do not. *%LAM* is the percentage of the recurrence points that have a vertical neighbor (orange) out of all the recurrence points (colored). **Figure 5.** Quantifying laminarity: An 8 × 8 RP where colored squares represent recurrence points. The orange-filled recurrence points have a vertical/horizontal neighbor, while the blue squares do not. *%LAM* is the percentage of the recurrence points that have a vertical neighbor (orange) out of all the recurrence points (colored).

ity with and without noise (Figure 6). In addition, there was a tight connection between the *%LAM* values and the true α-values, marked by a high R2 (0.96) quantifying correlation between α and *%LAM* for the noise condition. While there is a mathematical relation between *%LAM* and autocorrelations in a time series, the method has a downside in that it does not capture scaling relations within the The results corroborate this: persistent fractal fluctuations lead to increased laminarity with and without noise (Figure 6). In addition, there was a tight connection between the *%LAM* values and the true α-values, marked by a high R<sup>2</sup> (0.96) quantifying correlation between α and *%LAM* for the noise condition. *Entropy* **2022**, *24*, x FOR PEER REVIEW 7 of 16

The results corroborate this: persistent fractal fluctuations lead to increased laminar-

on the *x*-axis and the *%LAM* on the *y*-axis, when random noise (SNR 2:1) is added. Still, the resulted

2.1.4. Third Approach: Estimating Scaling Relations via Diagonal Recurrence Rates *(Diag* 

**Figure 7.** Approaches 3 and 4—diagonals in RP: An 8 × 8 RP presents the diagonal lines from 0 (main diagonal) to 7; due to the univariate RP's symmetrical characteristic, only the bottom triangle was

The third approach is based on diagonal recurrence profiles of a time series. The diagonal recurrence profile quantifies the number of recurrences at different lags, similar to the autocorrelation function [34]. To obtain the diagonal recurrence profile, one simply counts the proportion of recurrence points in the off-diagonals towards the lower-right or lower-left of the recurrence plot and plots them as a function of distance from the main diagonal; that is, lag [35]. Figure 7 illustrates the computation of the diagonal recurrence

coefficients scale well for every alpha value (−1 < α < 2).

*%REC)*

profile.

for both persistent and antipersistent fluctuations (α > 0). Right panel: box plots of the true α-values on the *x*-axis and the *%LAM* on the *y*-axis, when random noise (SNR 2:1) is added. Still, the resulted coefficients scale well for every alpha value (−1 < α < 2).

**Figure 6.** Results of approach 2 (laminarity). Left panel: box plots of the true alpha values on the *x*-

*Entropy* **2022**, *24*, x FOR PEER REVIEW 7 of 16

While there is a mathematical relation between *%LAM* and autocorrelations in a time series, the method has a downside in that it does not capture scaling relations within the data per se, and hence represents more of a correlate of fractal fluctuations, albeit a very useful one. axis and the *%LAM* on the *y*-axis. As can be seen, the coefficient scales well with the true α-values for both persistent and antipersistent fluctuations (α > 0). Right panel: box plots of the true α-values on the *x*-axis and the *%LAM* on the *y*-axis, when random noise (SNR 2:1) is added. Still, the resulted coefficients scale well for every alpha value (−1 < α < 2).

#### 2.1.4. Third Approach: Estimating Scaling Relations via Diagonal Recurrence Rates *(Diag %REC)* 2.1.4. Third Approach: Estimating Scaling Relations via Diagonal Recurrence Rates *(Diag %REC)*

The third approach is based on diagonal recurrence profiles of a time series. The diagonal recurrence profile quantifies the number of recurrences at different lags, similar to the autocorrelation function [34]. To obtain the diagonal recurrence profile, one simply counts the proportion of recurrence points in the off-diagonals towards the lower-right or lower-left of the recurrence plot and plots them as a function of distance from the main diagonal; that is, lag [35]. Figure 7 illustrates the computation of the diagonal recurrence profile. The third approach is based on diagonal recurrence profiles of a time series. The diagonal recurrence profile quantifies the number of recurrences at different lags, similar to the autocorrelation function [34]. To obtain the diagonal recurrence profile, one simply counts the proportion of recurrence points in the off-diagonals towards the lower-right or lower-left of the recurrence plot and plots them as a function of distance from the main diagonal; that is, lag [35]. Figure 7 illustrates the computation of the diagonal recurrence profile.

**Figure 7.** Approaches 3 and 4—diagonals in RP: An 8 × 8 RP presents the diagonal lines from 0 (main diagonal) to 7; due to the univariate RP's symmetrical characteristic, only the bottom triangle was **Figure 7.** Approaches 3 and 4—diagonals in RP: An 8 × 8 RP presents the diagonal lines from 0 (main diagonal) to 7; due to the univariate RP's symmetrical characteristic, only the bottom triangle was used. In approach 3, we counted the recurrence points in each diagonal and divided them by the diagonal's length. Additionally, approach 4 utilizes the ratio of recurrence percentage between every two subsequent diagonal lines. Both approaches focused on the middle diagonals to avoid the main diagonal's 100% recurrence points and the short diagonals towards the edges of the recurrence plot.

The rationale behind the approach is that the diagonal recurrence profile is a modelfree type of autocorrelation [33,36], and hence captures the magnitude of autocorrelation at different lags, which is related to fractal fluctuations in a time series [37]. Accordingly, a scaling relation between the logarithm of the recurrence rate and the logarithm of the diagonal number (reflecting the frequency spectra) should be related to fractal scaling. Here, a sharper negative slope indicates dominance of lower frequencies. Hence, contrasting the previous approaches, a lower power-law coefficient evidence a more persistent fluctuation. Correspondingly to spectral scaling analysis, this method yielded a scaling exponent of 0 to white noise (α = 0)—a benchmark to determine whether the time series is persistent, random, or antipersistent.

As can be seen in Figure 8, this approach distinguishes comparatively well between the different exponents for persistent fluctuations, with and without noise, but is less sensitive to the antipersistent fluctuations (however, the exponents are still increasing with decreasing negative alpha-values). Moreover, the relation to the true α-values appears strong for this range, even with noise (*R* <sup>2</sup> = 0.88).

**Figure 8.** Results of approach 3 (*Diag %REC*). (**A**) Left panel: box plots of the true alpha values on the *x*-axis and the power-law coefficients for the association of diagonal *%REC* and the diagonal index (distance from the main diagonal) on the *y*-axis. As can be seen, the coefficient scales well with the true α-values for both persistent and antipersistent fluctuations (−1 < α < 2). Right panel: box plots of the true alpha values on the *x*-axis and the power-law coefficients for the association of diagonal *%REC* and the diagonal index on the *y*-axis when random noise (SNR = 2:1) is added. Still, the resulted coefficients scale well with the true α-values for persistent and antipersistent fluctuations, but are somewhat less sensitive to distinguishing between different types of antipersistent fluctuations (α < 0). (**B**) Example of scaling plots demonstrating the association of diagonal *%REC* and diagonal index for different α values.

Another version of this approach is derived from Zbilut and Marwan's [38] proposal, which applied the Wiener–Khinchin theorem [38] to the analysis of diagonal recurrence profiles. They show that one can detect (nonlinear) periodicities by applying a Fourier transform to the diagonal recurrence profile of an RP (Figure 7). Just as with the raw diagonal recurrence profile, we fitted a linear trend line to the log–log plot power spectrum (obtained via the Fourier Transform) of the diagonal recurrence profile (Figure 9). The results were similar to what we observed for the raw diagonal recurrence profile in that the method distinguished between persistent (*R* <sup>2</sup> = 0.68, <sup>α</sup> <sup>≥</sup> 0) fluctuations. However, the standard errors were higher, and the method did not capture antipersistent fluctuations (*R* <sup>2</sup> = 0.002, <sup>α</sup> <sup>≤</sup> 0). *Entropy* **2022**, *24*, x FOR PEER REVIEW 10 of 16

persistent fluctuations (α > 0). Right panel: box plots of the true α-values on the *x*-axis and the power-law coefficients for the association of FT of the diagonal *%REC* and the diagonal index on the *y*-axis, when random noise (SNR = 2:1) is added. The relation of the resulting coefficients to the

index on the *y*-axis. As can be seen, the coefficient scales well with the true α-values for the persistent fluctuations (α > 0). Right panel: box plots of the true α-values on the *x*-axis and the power-law coefficients for the association of FT of the diagonal *%REC* and the diagonal index on the *y*-axis, when random noise (SNR = 2:1) is added. The relation of the resulting coefficients to the true α-values is not as good for persistent fluctuations (cannot differentiate α = 0 and 0.5) and is relatively insensitive to distinguishing between different types of antipersistent fluctuations. (**B**) Example of scaling plots demonstrating the association of FT of the diagonal *%REC* and diagonal index for different α values. *Entropy* **2022**, *24*, x FOR PEER REVIEW 11 of 16 Example of scaling plots demonstrating the association of FT of the diagonal *%REC* and diagonal index for different α values.

#### 2.1.5. Fourth Approach: Consecutive Diagonals Recurrence Ratio *(Diag ratio)* 2.1.5. Fourth Approach: Consecutive Diagonals Recurrence Ratio *(Diag ratio)*

Until this point, the analysis techniques were more effective for persistence signals and did not distinguish between antipersistent signals well. Approach number four solves this issue to some degree. Similar to the third approach, we utilized the recurrence percentage of the diagonal lines. Here, however, we calculate the ratio between each couple of consecutive diagonal lines (Figure 7). The rationale behind the approach is that antipersistent fluctuations will tend to yield oscillations at high frequencies, and the ratio of recurrence rate of adjacent diagonals in the recurrence plot will capture the magnitude of such oscillations. Just as with the laminarity measure, however, this method is more of a correlate of antipersistent fractal scaling, and does not capture scaling properties directly. Until this point, the analysis techniques were more effective for persistence signals and did not distinguish between antipersistent signals well. Approach number four solves this issue to some degree. Similar to the third approach, we utilized the recurrence percentage of the diagonal lines. Here, however, we calculate the ratio between each couple of consecutive diagonal lines (Figure 7). The rationale behind the approach is that antipersistent fluctuations will tend to yield oscillations at high frequencies, and the ratio of recurrence rate of adjacent diagonals in the recurrence plot will capture the magnitude of such oscillations. Just as with the laminarity measure, however, this method is more of a correlate of antipersistent fractal scaling, and does not capture scaling properties directly.

As seen in Figure 10, with this measure, we can differentiate negative α-values (antipersistent) from α = 0, both with and without external noise. However, the method does not distinguish between the different alpha values of the persistent fluctuations. As seen in Figure 10, with this measure, we can differentiate negative α-values (antipersistent) from α = 0, both with and without external noise. However, the method does not distinguish between the different alpha values of the persistent fluctuations.

**Figure 10.** Results of approach 4 (consecutive diagonals *%REC* ratio). Left panel: box plots of the true alpha values on the *x*-axis and the mean ratio between subsequent diagonals' *%REC* on the *y*axis. As can be seen, the coefficient scales well with the true α-values for antipersistent fluctuations (α < 0) and converges to 1 from α = 0. Right panel: box plots of the true α-values on the *x*-axis and the mean ratio between subsequent diagonals' *%REC* on the *y*-axis when random noise (SNR 2:1) is added. Still, the resulted coefficients scale well for negative alpha value (*R*2 = 0.79). **Figure 10.** Results of approach 4 (consecutive diagonals *%REC* ratio). Left panel: box plots of the true alpha values on the *x*-axis and the mean ratio between subsequent diagonals' *%REC* on the *y*-axis. As can be seen, the coefficient scales well with the true α-values for antipersistent fluctuations (α < 0) and converges to 1 from α = 0. Right panel: box plots of the true α-values on the *x*-axis and the mean ratio between subsequent diagonals' *%REC* on the *y*-axis when random noise (SNR 2:1) is added. Still, the resulted coefficients scale well for negative alpha value (*R* <sup>2</sup> = 0.79).

The approaches were tested on a dataset of a tapping experiment during which participants listened to a certain beat and were then instructed to tap according to the tempo they had heard. Under one of the two within-participant conditions, participants received visual feedback on every trial to help them align their tapping performance with the target

*2.2. Empirical Example* 

series.

#### *2.2. Empirical Example* feedback would reduce long-range dependencies in the data related to cognitive-motor processes of timing, and hence yield a more random ('whiter') noise manifested by a lower

The approaches were tested on a dataset of a tapping experiment during which participants listened to a certain beat and were then instructed to tap according to the tempo they had heard. Under one of the two within-participant conditions, participants received visual feedback on every trial to help them align their tapping performance with the target tempo, while in the other condition no such feedback was provided. The sample was comprised of 36 time series from 18 participants with at about 1000 tapping intervals per time series. α exponent [39]. Our findings, displayed in Figure 11, support these expectations in several ways. Firstly, a negative power-law coefficient in approach 3 along with a ~1 ratio between subsequent diagonals (approach 4) indicate a persistent fluctuation in both conditions and is supported by a Hurst exponent 1.0 > *H* > 0.5, suggesting a pinkish noise. Further, approaches 1–3, as well as Wiener–Khinchin theorem's results, imply a lower αexponent for the feedback condition (see Table 1). While *SD %REC* and *%LAM* exhibit it by presenting a higher clustering characteristic for the no-feedback condition, *Diag %REC*

Drawing on previous research on cognitive processes, we expected the time series to show persistent fractal fluctuation. Moreover, previous research showed that receiving

*Entropy* **2022**, *24*, x FOR PEER REVIEW 12 of 16

Drawing on previous research on cognitive processes, we expected the time series to show persistent fractal fluctuation. Moreover, previous research showed that receiving feedback would reduce long-range dependencies in the data related to cognitive-motor processes of timing, and hence yield a more random ('whiter') noise manifested by a lower α exponent [39]. Our findings, displayed in Figure 11, support these expectations in several ways. Firstly, a negative power-law coefficient in approach 3 along with a ~1 ratio between subsequent diagonals (approach 4) indicate a persistent fluctuation in both conditions and is supported by a Hurst exponent 1.0 > *H* > 0.5, suggesting a pinkish noise. Further, approaches 1–3, as well as Wiener–Khinchin theorem's results, imply a lower α-exponent for the feedback condition (see Table 1). While *SD %REC* and *%LAM* exhibit it by presenting a higher clustering characteristic for the no-feedback condition, *Diag %REC* and the Wiener–Khinchin theorem display it with a stronger lower frequency dominance when no feedback is given. and the Wiener–Khinchin theorem display it with a stronger lower frequency dominance when no feedback is given. **Table 1.** Paired *t*-test results comparing the outcomes of feedback and no-feedback conditions. **Approach** *t df p* 1—*SD %REC* −5.12 17 >0.001 2—*%LAM* −3.33 17 0.004 3—*Diag %REC* 3.43 17 0.003 4—*Diag ratio* 0.82 17 0.42 Wiener–Khinchin theorem 3.68 17 0.002 DFA −4.08 17 >0.001

**Figure 11.** Box plots illustrating the outcomes of feedback and no-feedback conditions: Six sets of box plots represent a comparison between the outcomes of each approach for feedback (orange) and no-feedback (blue) conditions. While the frame of the boxplot is defined by the interquartile range, the notch represents a 95% confidence interval and the whiskers show the maximum and the minimum of each distribution (except outliers). As expected, due to its persistent noise characteristics (α > 0), behavioral data would be appropriately analyzed by approaches 1–3 but not approach 4. Approaches 1 and 2, as well as DFA, yield higher results for the no-feedback condition, indicating a larger α, meaning a more persistent behavior. Likewise, approach 3 and Wiener–Khinchin theorem suggest a lower frequency dominancy in the no-feedback condition. **Figure 11.** Box plots illustrating the outcomes of feedback and no-feedback conditions: Six sets of box plots represent a comparison between the outcomes of each approach for feedback (orange) and no-feedback (blue) conditions. While the frame of the boxplot is defined by the interquartile range, the notch represents a 95% confidence interval and the whiskers show the maximum and the minimum of each distribution (except outliers). As expected, due to its persistent noise characteristics (α > 0), behavioral data would be appropriately analyzed by approaches 1–3 but not approach 4. Approaches 1 and 2, as well as DFA, yield higher results for the no-feedback condition, indicating a larger α, meaning a more persistent behavior. Likewise, approach 3 and Wiener–Khinchin theorem suggest a lower frequency dominancy in the no-feedback condition.

24


**Table 1.** Paired *t*-test results comparing the outcomes of feedback and no-feedback conditions.

#### *2.3. Comparison of the Approaches*

To evaluate the presented approaches in relation to the true alpha values of the generated time series, we focus on three main parameters: (a) fractal dimension range, (b) sensitivity to noise, and (c) summary of the quantitative relation to the true alpha values. Furthermore, we investigated their applicability to empirical behavioral data.

#### 2.3.1. Range

As presented above, approaches 1 (*SD %REC*) and the Wiener–Khinchin-based analysis are sensitive to persistent fluctuations. Conversely, approach 4 (*Diag ratio*) differentiates only antipersistent fluctuations, whereas approaches 2 and 3 (*%LAM*, *Diag %REC*) are applicable throughout the whole tested range (−1 < α < 2), like DFA. Hence, with no estimation of the time series' fractal dimension, one should conduct an analysis according to approaches 2 or 3, otherwise the researcher might prefer to pick the analysis technique that best fits his data's characteristics. On a similar note, one can try to detect whether there are persistent fluctuations using approach 4, which yields a ~1 ratio for α ≥ 0.

#### 2.3.2. Robustness to Noise

Most of the analysis techniques that were applied were robust to noise. Except for the Wiener–Khinchin theorem approach, the rest distinguished between α-values within their range comparably with and without noise. Nevertheless, antipersistent fluctuations were less distinguishable by both approach 3 and DFA when i.i.d. noise (SNR = 2:1) was introduced.

#### 2.3.3. Quantitative Relation to True Alpha Values

Table 2 provides a summary of the *R* 2 -values that capture the relation between the true α-values and the estimated parameters of the different approaches, separately for persistent and antipersistent fluctuations. As DFA is the gold standard for fractal analyses in time-series methods, the comparison of the recurrence-based approaches to DFA is of particular interest here. Comparing the likelihoods of the linear models of each of our four approaches to DFA, we found that the association between the true α values and Hurst exponent is significantly stronger than in almost every other method (α < 0.05). On the contrary, approaches 2 and 4 yielded significantly higher association (than DFA) with the true α values for antipersistent fluctuation when noise is introduced, but somewhat below DFA under the no-noise condition. Nevertheless, approaches 1–3 show similar *R* 2 to DFA when analyzing persistent noise. It has to be kept in mind that the sample sizes here are quite large, and tests of significance are of limited value in this case.

#### 2.3.4. Applicability

All approaches were found applicable to behavioral data and concluded conformally despite the small sample size. The utilized data were most likely to behave persistently and hence were out of the fourth approach range. Yet, we suggest using approach 4 to confirm whether the time series is persistent or not (persistent fluctuations are indicated by a 1:1 ratio between subsequent diagonals). Table 3 provides an overview of the *R 2* for the different approaches (including DFA) when comparing the two time-estimation groups (i.e., with and without feedback). In a model comparison, approaches 2–4 were less predictive than DFA (α < 0.05), while approach 1 was not significantly lower.

**Table 2.** Comparison of approaches on simulated data.


**Table 3.** Comparison of approaches on empirical data.


#### **3. Conclusions**

In the current paper, we presented and compared several recurrence-based approaches to quantify the strength of monofractal autocorrelations in time-series data. This is a major step forward for integrating the quantification of scaling properties into recurrence quantification analysis, as previous research has suggested that such analyses are theoretically possible (e.g., [28]), but did not point to concrete means for how to deduce such properties. The proposed methods differ in quality, as well as in the range of applicability to particular types of colored noise, as we have shown on synthetic and empirical data. Based on our results, we recommend using approaches 3 and 4 to determine whether the data are persistent, antipersistent, or white noise. Then, approaches 1, 2, and 3 would be suitable to compare the fractal dimensionality of persistent data, while approaches 2 and 4 would fit antipersistent time series.

Thus, the present work lays the foundations for integrating fractal analysis into an RQA framework, and defining appropriate recurrence-based quantifies. Moreover, these methods might be amenable to quantifying time-dependent fractal fluctuations of not only univariate time series, but also strange attractor profiles, which possess fractal properties and are readily analyzable within the framework of recurrence quantification analysis [28]. In the future, these methods could be extended to capturing fractal dimensions in multidimensional systems via multidimensional recurrence quantification analysis. In addition, an evaluation and adaptation of these approaches to multifractals would be valuable [40,41].

**Author Contributions:** Conceptualization, A.T. and S.W.; methodology, A.T.; formal analysis, A.T. and S.W.; investigation, A.T., S.W. and G.L.; writing—original draft preparation, A.T., S.W. and G.L.; writing—review and editing, A.T., S.W. and G.L.; All authors have read and agreed to the published version of the manuscript.

**Funding:** S.W. acknowledges funding from the German Science Foundation (DFG; grant number 442405852).

**Institutional Review Board Statement:** The collection of the timing data was approved by the institutional review board of Leuphana University of Lüneburg ("EB-Antrag\_202202-03-Wallot\_Timing Distraction").

**Data Availability Statement:** The data is available upon request from the corresponding author. Also, a Matlab function (The MathWorks, Inc.) to perform the four analyses is available on GitHub, see Appendix A.

**Acknowledgments:** We thank Stine Hollah and Gerke Feindt for collecting the timing data presented in this manuscript. S.W. acknowledges funding by the German Research Foundation (DFG; grant numbers 442405852 and 442405919).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A**

A Matlab (The MathWorks, Inc.) implementation of the four approaches is available on GitHub: https://github.com/alontom/FARQA. Accessed on 19 July 2022.

#### **References**

