*Article* **Studying Physiological Synchrony in Couple Therapy through Partial Directed Coherence: Associations with the Therapeutic Alliance and Meaning Construction**

**Evrinomy Avdi 1,\*, Evangelos Paraskevopoulos <sup>2</sup> , Christina Lagogianni <sup>1</sup> , Panagiotis Kartsidis <sup>3</sup> and Fotis Plaskasovitis <sup>1</sup>**


**Abstract:** In line with the growing recognition of the role of embodiment, affect and implicit processes in psychotherapy, several recent studies examine 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 examine its potential to contribute to our understanding of the therapy process. The study adopts a single-case, mixedmethod design and examines physiological synchrony in one-couple therapy in relation to the therapeutic alliance and a narrative analysis of meaning construction in the sessions. Interpersonal Physiological Synchrony (IPS) was calculated, via a windowed approach, through PDC of a Heart Rate Variability-derived physiological index, which was measured in the third and penultimate sessions. Our mixed-method analysis shows that PDC quantified significant moments of IPS within and across the sessions, modeling the characteristics of interpersonal interaction as well as the effects of therapy on the interactional dynamics. 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 processes in psychotherapy.

**Keywords:** physiological synchrony; heart rate; therapeutic alliance; psychotherapy process; couple therapy

#### **1. Introduction**

This study rests on the assumption that psychotherapy relies on both implicit and explicit processes and that both need to be taken into account when studying clinical process [1,2]. It focuses on one aspect of implicit interaction, interpersonal physiological synchrony (IPS), and introduces the use of Partial Directed Coherence as a metric for operationalizing IPS in psychotherapy sessions. Using a single-case, mixed-method design on one couple therapy, physiological synchrony is examined in relation to the therapeutic alliance and a qualitative analysis that draws upon narrative principles of meaning reconstruction in the sessions.

Synchrony is observed in many complex biological systems and is assumed to occur through nonlinear dynamic processes rather than simple causal links. In social interaction, synchrony concerns the temporal covariation of behavior or internal states in interacting partners and can be broadly defined as 'the social coupling of two (or more) individuals in the here-and-now of a communication context that emerges alongside, and in addition to, their verbal exchanges' [3] (p. 558).

**Citation:** Avdi, E.; Paraskevopoulos, E.; Lagogianni, C.; Kartsidis, P.; Plaskasovitis, F. Studying Physiological Synchrony in Couple Therapy through Partial Directed Coherence: Associations with the Therapeutic Alliance and Meaning Construction. *Entropy* **2022**, *24*, 517. https://doi.org/10.3390/e24040517

Academic Editors: Franco Orsucci and Wolfgang Tschacher

Received: 31 January 2022 Accepted: 2 April 2022 Published: 6 April 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/).

A key concept in the literature on interactional synchrony is interpersonal coordination, which refers to the degree to which the behaviors of interacting partners are nonrandom, patterned or synchronized in timing and form [4]. There is ample evidence that behavioral matching and interactional synchrony are ubiquitous features of human interaction, on both verbal (e.g., vocal tone, word choice, laughter, speech accent, syntax, intonation) and nonverbal (e.g., posture, gesture, facial expression, orientation, etc.) levels. Interpersonal coordination emerges early in life and is an automatic, non-conscious process that is associated with liking, affiliation, rapport, cooperation, self–other merging, perspective taking, empathy, smoothness of interaction, prosocial behaviors, compassion and increased performance in tasks that rely on joint actions [5–7].

In the context of psychotherapy, 'being in sync' has been examined primarily in relation to nonverbal behaviors and has been shown to be associated with important psychotherapy processes, such as rapport [8], therapist empathy, the therapeutic alliance [6,7], session quality and therapy outcome [9–11], as well as mental state in relation to attachment [12–14]. Drawing upon developmental research, several authors have proposed that synchronous behaviors between therapist and client are crucial for the formation of the therapeutic alliance, which in turn promotes affect regulation in the client and fosters therapeutic change [7]. Similarly, research on infant development suggests that repeated experiences of biobehavioral synchrony between infants and their parents are central to the development of affect regulation capacities in the infant and security of attachment [15–18]. There is some evidence that synchrony is associated with affect regulation in adulthood as well, as interacting partners in close relationships coregulate their arousal around a homeostatic optimal level [19,20].

#### *1.1. Interpersonal Physiological Synchrony*

In addition to studying synchrony in observable behavior, in recent years, there has been a growing interest in the role of synchrony in physiological arousal in psychotherapy. This is in line with the recognition that psychological and social processes cannot be isolated from embodiment and affect [21,22]. The inclusion of affective and embodied aspects of interaction is arguably particularly relevant to psychotherapy, given that affect is intimately linked with meaning construction and forms an integral part of the work of therapy [23]. The Autonomic Nervous System (ANS) plays a key role in cognition, emotion and behavior [24], and although ANS activation is not specific to affect, most emotions are associated with increased physiological arousal [25,26]. As such, several recent studies include psychophysiological measures in psychotherapy process research and treat physiological activation, and particularly its arousal component, as an index of affect [27,28]. In this literature, it is assumed that measures of psychophysiology enable the study of aspects of the therapy process that may not be accessible through self-report or observation, and can therefore add another layer of information on clinical process [29]. In other words, psychophysiological measures may reflect non-conscious, implicit affective processes and can then be used as correlates of implicit intra- and interpersonal processes in therapy [23,30]. In addition to these theoretical developments, technological advances make the continuous recording of physiological states in therapy relatively easy and unobtrusive.

Research on interpersonal physiology concerns the temporal coregulation of physiological activation in interacting partners, using continuous measures of physiological activity. The indices of ANS arousal most commonly used include electrodermal activity (EDA), considered to reflect sympathetic arousal, and variables associated with heart rate (e.g., heart rate variability), which are associated with both sympathetic and parasympathetic activity. Due to the sufficient time resolution of these variables [31], their outcome may be used to estimate the influence that one person's physiological indices exert over another's, through a model of physiological interactions or coupling [32,33]. In this context, interpersonal physiological synchrony (IPS) is defined as 'any interdependent or associated activity identified in the physiological processes of two or more individuals' [34] (p. 2). Recent reviews of studies of IPS in different interactional contexts suggest that physiological synchrony is a robust phenomenon identifiable through different methods [11,34].

In the context of psychotherapy, physiological synchrony between therapists and clients was first examined in a series of studies in the 1950s in relation to rapport and empathy [35]. More recently, the role of IPS in the psychotherapy process has been examined in several studies of psychotherapy sessions [3,23,30,36–40], as well as simulated sessions [13,38,41]. In a recent review of this literature, Kleinbub [14] concluded that physiological synchrony in psychotherapy is an established fact, although its clinical meaning is far from known.

Physiological synchrony in psychotherapy has primarily been associated with empathy [11,16,42,43]. However, research in interactional contexts other than psychotherapy suggest that physiological synchrony is not uniquely associated with empathy and is not necessarily positive for interactions. For example, research on infant development [17,18] shows that attachment security is associated with medium-range synchrony in parent– infant interaction and that 'too much' synchrony is predictive of attachment insecurity. Similarly, findings regarding the role of physiological linkage in the quality of adult romantic relationships are mixed, with several studies showing that increased physiological linkage in couples tends to be associated with poorer relationship satisfaction and the escalation of negative affect [44]. The evidence to date suggests that, in the context of negative interactions, IPS is associated with relationship dissatisfaction and conflict, whereas in positive interactions, it is primarily associated with empathy and rapport [34]. In addition to the *affective valence* of interactions, the *degree of emotional arousal* may also moderate physiological synchrony; for example, in studies of mother–infant interactions, higher maternal heart rate, thought to reflect increased affective arousal, has been associated with lower physiological synchrony with her infant [45]. Drawing upon these findings, it seems important for future research to take into account the characteristics of the relational context when studying the role of IPS in psychotherapy.

A related issue concerns the way IPS is conceptualized, operationally defined and calculated. The majority of studies to date of IPS in psychotherapy examine only positive correlations, i.e., in-phase synchrony, where the therapist's and client's arousal covary in the same direction, and assume that negative correlations, or anti-phase synchrony, reflect lack of synchrony. Other studies, however, suggest that anti-phase synchrony, where one partner's physiological arousal decreases as the other partner's increases, reflects processes of coregulation or complementarity [46,47]. For example, in one study implicating a storytelling task, it was found that the narrator's autonomic arousal decreased when the listener's increased and he or she displayed affiliation; this was interpreted as reflecting a process of 'sharing the emotional load', whereby the listener's engagement regulated the teller's physiological arousal [48]. Similarly, in a study of ANS activation in psychoanalytic therapy, the therapists' empathic displays were associated with increased arousal in the therapist and decreased arousal in the client, whereas sequences of the therapists' challenges were associated with increases in both participants' arousal [49]. In line with these findings, Butler & Randell [19] suggest that asynchrony may be associated with stress buffering, whereby one individual moderates the stress level of another. Based on the above, including both in-phase and anti-phase synchrony in studies of IPS in psychotherapy is likely to provide a more nuanced approach to understanding this multifaceted interactional phenomenon.

The metric employed to estimate interaction is also of importance. Most studies investigating IPS use correlation-derived estimates, which are sensitive to spurious correlations and do not address causality or directionality in the interaction [14]. In order to overcome this issue, approaches that employ specific causality tests, such as Granger causality, adjusted for estimating the information flow between multivariate time series can be used in the frequency domain [50]. Combined with surrogate testing of the parameters used to estimate interactions [11], such approaches may be combined with a windowed analysis to reach a stable and fine-grained temporal resolution that can also provide directionality.

Another important issue when examining physiological synchrony in psychotherapy relates to the timescales employed in the analysis. Most studies calculate IPS over whole sessions, despite the fact that IPS is likely to be a transient phenomenon that fluctuates through sessions [1,34]. Similarly, recent studies approach the therapeutic alliance as a dynamic phenomenon and show that therapy sessions contain several periods characterized by ruptures in the alliance, often followed by interactive repair [51]. Indeed, several authors suggest that it is precisely such repairs that are important for optimal development and therapeutic change [51–53]. Therefore, examining synchrony on a more micro-level of interaction can shed light on processes that may not be apparent at the session level.

In sum, research on physiological synchrony in psychotherapy suggests that it can add important information regarding the psychotherapy process; given that IPS may reflect different interactional processes—including empathy, affect coregulation and conflict caution is needed when interpreting findings. Moreover, the field is fragmented on both conceptual and methodological levels, as reflected in the prevalent lack of agreement on terminology, data collection methods, research designs and statistical analyses [11,34,54]. Recent reviews suggest that, given how little we know about the context-specific factors that affect IPS, it may be preferable to use idiographic designs and theoretically informed analyses of the therapy process. Since the publication of these reviews, a few such studies have been published that shed light on the different functions of physiological synchrony in psychotherapy [23,30,37,42,55–57].

Before turning to the current study, we briefly discuss the concept of the therapeutic alliance, with a focus on couple therapy, given that it is a key clinical concept that has been associated with physiological synchrony.

#### *1.2. The therapeutic Alliance in Couple Therapy*

Several contemporary approaches to psychotherapy adopt a discursive and narrative perspective and conceptualize the process of change in psychotherapy in terms of meaning reconstruction [58]. In this framework, psychotherapy is described as a semantic process that relies on the creation of a dialogical space, which facilitates the reconstruction of clients' life narratives so that they become more complex, polyphonic, emotionally salient, inclusive and flexible [59]. The therapist's receptive and relationally responsive attitude towards the clients' storytelling and expression of affect are considered crucial elements in this process [60]. There is ample evidence that different therapist actions associated with responsiveness play an important role for the process and outcome of psychotherapy [61], with the therapeutic alliance being a key relational aspect in this process.

The therapeutic alliance is a pan-theoretical concept that is associated with the collaborative aspects of the therapeutic relationship and has been extensively studied as an important process variable in psychotherapy. It is usually conceptualized as comprising three interlinked aspects: a strong emotional bond between clients and therapists, and agreement and collaboration on the goals and the tasks of therapy [62]. The quality of the therapeutic alliance has consistently been shown to be a predictor of outcome in individual psychotherapy across different modalities [63], as well as couple and family therapy (CFT) [64–69]

In conjoint treatments, such as couple therapy, the therapeutic alliance consists of a web of interlinked relationships between participants and the various subsystems thus formed [63,66]. Several factors—such as power dynamics, conflict, trust, loyalties and secrets in the couple or family—affect the formation of the alliance in CFT [69–71]. A strong overall alliance in couple therapy requires a balanced alliance between the therapist and each partner, as well as agreement in the couple on the problems, goals and values of therapy; as such, the therapist is encouraged to foster an alliance with each partner, avoiding 'split alliances', and to promote within-couple alliance [66].

The current study is a mixed-method, single-case study aiming to illustrate the potential of the PDC metric as a useful way of examining IPS in relation to the therapy process; it assumes a theoretically driven idiographic design and examines whether the therapeutic alliance maps onto IPS findings.

#### **2. Materials and Methods**

The research material in this study is drawn from one-couple therapy, conducted in an outpatient Family Therapy Department in Greece, in the context of a wider naturalistic, multisite research study [30,39,72]. The treatment in this service follows systemic principles and includes the use of reflective conversations with a co-therapist. In usual clinical practice, sessions are provided monthly; a second therapist watches the session between the primary therapist and the couple behind a one-way mirror and joins them for a reflective conversation towards the end of each session [73]. Participating couples were informed about the study by a graduate researcher at the end of their first session. Participation in the project was voluntary, and ethical approval was granted by the Family Therapy Department's Scientific Board. Both clients and therapists gave permission for the data to be used for research purposes.

#### *2.1. The Case*

This therapy consisted of 15 sessions spanning 14 months. The couple, Costas and Demetra, is a white heterosexual couple in their mid-thirties. Demetra is a law graduate with a successful professional career. Costas has no university education; he worked as a technician in the past and is currently unemployed. The couple had been in a long-term relationship of over 10 years when they came to therapy. They sought therapy because of increasing tension in their relationship following the birth of their baby 10 months earlier. Two experienced female clinical psychologists and systemic family therapists in their fifties participated in this therapy. The therapy centered on Demetra's distress in her role as a mother, the expression of anger and conflict between the spouses, and Costas' low self-esteem associated with periods of unemployment. At the end of treatment, the couple reported an improvement in their personal lives and their relationship.

#### *2.2. Procedure*

All sessions were video-recorded in split-screen mode with four web-cameras. In addition, in two sessions (sessions 3 and 14), physiological measures of the participants' heart rates were recorded for the duration of the session. Within 24 h of the measurement sessions, a graduate researcher conducted separate Stimulated Recall interviews [74,75] with each client and therapist, each lasting approximately 30 min.

#### *2.3. Measures*

#### 2.3.1. Autonomic Nervous System Responses

The participants' autonomic nervous system (ANS) responses were recorded via Firstbeat Bodyguard (Firstbeat Technologies, Jyväskylä, Finland) [76] mobile heart rate (HR) monitors. Ag/AgCl electrodes, connected to the Firstbeat Bodyguard, were attached on two sites on the skin of the chest before the start of each measurement session and were removed the next day, with the guidance of a graduate researcher. HR was continuously recorded during this period.

#### 2.3.2. Clinical Outcomes in Routine Evaluation–Outcome Measure (CORE-OM)

The outcome of therapy was examined using the CORE-OM, administered at the start and end of therapy. The CORE-OM is a widely used, 34-item self-report measure that examines psychological distress in four domains: wellbeing, problems, functioning and risk [77].

#### 2.3.3. Session Rating Scale (SRS)

The SRS is a four-item, ultra-brief visual analogue instrument to assess the global strength of the alliance, designed to be used in routine outcome monitoring [78]. The four items measure the therapist–client emotional bond, agreement on goals, agreement on tasks and overall rating of the alliance. It is scored by summing the marks measured to the nearest centimeter on each of the four lines. Based on a total possible score of 40, any score lower than 36 overall, or 9 on any scale, could be a source of concern.

#### 2.3.4. System for Observing Family Therapy Alliances (SOFTA-o)

The SOFTA-o is an observer-based measure developed to study the therapeutic alliance in couple and family therapy [79]. It examines the contribution of each participant to the alliance by coding specific behaviors in four dimensions: Emotional connection, Engagement in the therapeutic process, Safety within the therapeutic system and Shared sense of purpose. The first three dimensions concern the therapist(s)–clients relationship, whereas the fourth concerns the couple sub-system. Following the coding of specific items, global ratings are provided for each dimension on a 7-point ordinal scale, ranging from -3 (extremely problematic) to +3 (extremely strong), with 0 denoting an unremarkable or neutral alliance. These dimensions are conceptually interdependent and moderately correlated and can be combined in a composite score [66].

#### *2.4. Data Analysis*

#### Interpersonal Physiological Synchrony

Data from the ANS were analyzed using Firstbeat PRO Wellness Analysis Software® version 1.4.1. This software uses neural network modeling to calculate Heart Rate Variability (HRV) indices second-by-second. This is achieved using a short-time Fourier Transform method (STFT) combining data from HR- and HRV-derived variables that describe respiration rate and oxygen consumption (VO2). In addition, the absolute stress vector (ASV) is calculated from the HR, high-frequency power (HFP), low-frequency power (LFP) and HRV-derived respiratory variables, as an index of the activity of the sympathetic nervous system. The ASV grounds on detecting sympathetic reactivity that exceeds the momentary metabolic requirements of the ANS. Hence, the ASV is high when the heart rate is elevated, HRV is low and respiration rate is low relative to HR and HRV [80]. The ASV is calculated at a 1 Hz rate.

#### *2.5. Partial Directed Coherence within Sessions*

Within-session, directed, interpersonal physiological synchrony based on ASV was estimated using Partial Directed Coherence (PDC) [50]. PDC analysis transforms the ASV time series into the frequency domain and provides time-lagged associations between two participants' multivariate signals, assessing their statistical independence or predictability [50]. Specifically, grounded on instantaneous Granger causality, it implies that, knowing the previous states of the first signal (the leading signal), one may achieve a better prediction of the second signal (the pacing signal), than just knowing the previous states of the second signal. Hence, it describes the direction of information flow between isolated pairs of time series, in a frequency-domain representation of the notion of Granger causality. This approach has recently been proposed as a method of choice for estimating IPS in psychotherapy by Kleinbub [54], due to its ability to establish direction, and thus causality, in interactions. Due to the time-varying conditional variance of HRV signals [81], PDC as a frequency-domain method for identifying causal interactions between the signals was preferred over the classical Granger causality, which estimates interactions in the time domain. In addition, PDC has previously been used to successfully estimate the frequencydomain causality in cardiovascular time series with Instantaneous Interactions [82].

The second-by-second ASV data of the measurement sessions were imported into Matlab (MathWorks Inc., Natick, MA, USA) as time series. The ASV time series were segregated into time-windows of 50 s, and the PDC for each window was estimated independently for each pair of participants in each session via an in-house script based on the work of Baccalá and Sameshima [50]. The length of the time-window was empirically determined on the basis of a series of tests comparing the number of significant PDC time-windows within independent sets of surrogate data generated via Matlab, aiming to achieve the best possible balance between the resolution of the analysis (i.e., smallest time-window) and the absence of false-positive significant PDC time-windows. Hence, for each 50 s time-window of the session, we retrieved two PDC values for each pair of participants (one for each direction, i.e., one in which participant 1 leads and participant 2 paces, and one in which participant 2 leads and participant 1 paces). Additionally, a statistical test based on Monte Carlo iterations of the corresponding data was performed for each pair, in order to identify time-windows with a significant PDC. The threshold of significance was defined as *p* = 0.05/3, accounting for the total number of comparisons in which the same set of data participated, thereby effectively controlling for multiple comparisons. Only significant PDC values were taken into account.

#### *2.6. Partial Directed Coherence between Sessions*

The number of significant PDC time-windows for each pair of participants was compared between sessions 3 and 14 as an index of the overall effect of therapy on interpersonal physiological synchronization. The aim was to identify differences in the global characteristics of IPS between sessions at the start and end of therapy.

#### *2.7. Qualitative Analysis of the Therapy Process*

#### 2.7.1. Topical Episodes

The measurement sessions were segmented into topical episodes, i.e., periods of time during which a specific topic was discussed [83]. This coding was initially carried out by two graduate researchers and was checked by third researcher, and any discrepancies were resolved through discussion. This initial thematic coding provides a description of the main themes discussed in a session. Session 3 was segmented into 14 topical episodes, ranging from 2 to 15 min' duration, and session 14 was segmented into 12 topical episodes, ranging from approximately 1 to 9 min' duration.

#### 2.7.2. Therapeutic Alliance

Two graduate psychologists, trained in using the SOFTA-o, coded each session. The raters coded the sessions independently and then discussed any discrepancies until consensus was reached. Next, in order to gain a more fine-grained coding of the development of the alliance through the session, the strength of the alliance was coded for each topical episode.

#### **3. Findings and Discussion**

With regards to the outcome of therapy, the clients' CORE-OM scores decreased significantly over the course of therapy, suggesting a clinically significant reduction in psychological distress (Table 1). At the onset of therapy, both partners reported a medium level of distress, and, importantly, Costas scored on items concerning the risk of harming himself. At the end of therapy, Demetra's CORE-OM score decreased to the cut-off point for clinical distress (<10), and Costas' showed clinically significant change (>5 clinical score points) [77]. In terms of the therapeutic alliance, Costas' scores indicated a positive alliance in session 3, which further increased in the penultimate session, whereas Demetra's scores indicated a problematic alliance in session 3, which improved in the penultimate session.

**Table 1.** Clients' CORE-OM and SRS scores.


Next, we present the key quantitative findings regarding interpersonal physiological synchrony (IPS) within and across the two measurement sessions. Then, the potential of PDC analysis as a useful way of examining the process of therapy is explored through a mixed-method analysis of session 3. cal synchrony (IPS) within and across the two measurement sessions. Then, the potential of PDC analysis as a useful way of examining the process of therapy is explored through a mixed-method analysis of session 3. The physiological activity of the couple, as reflected in their ASV, in the two sessions

Next, we present the key quantitative findings regarding interpersonal physiologi-

Demetra 12 10 0 1.6 5.6 8.0 Costas 19 11 5 0 8.9 9.8

(>5 clinical score points) [77]. In terms of the therapeutic alliance, Costas' scores indicated a positive alliance in session 3, which further increased in the penultimate session, whereas Demetra's scores indicated a problematic alliance in session 3, which improved

> **CORE-OM CORE-OM RISK SRS Session 1 Session 15 Session 1 Session 15 Session 3 Session 14**

The physiological activity of the couple, as reflected in their ASV, in the two sessions is presented in Figure 1. In both sessions, Demetra's autonomic arousal decreased as the session progressed, whereas Costas' remained relatively constant through. It is worth noting that Demetra's mean ASV score in the penultimate session was significantly higher than in the third session, and her arousal shows higher variance. The clinical relevance of this observation would require further investigation and lies beyond the scope of this study. is presented in Figure 1. In both sessions, Demetra's autonomic arousal decreased as the session progressed, whereas Costas' remained relatively constant through. It is worth noting that Demetra's mean ASV score in the penultimate session was significantly higher than in the third session, and her arousal shows higher variance. The clinical relevance of this observation would require further investigation and lies beyond the scope of this study.

*Entropy* **2022**, *24*, x FOR PEER REVIEW 8 of 18

in the penultimate session.

**Table 1.** Clients' CORE-OM and SRS scores.

**Figure 1.** Demetra and Costa's ASV values in the sessions, plotted against the mean ASV value for the session (red line), 2 Standard Deviations below (yellow line) and 2 SD above (grey line) the mean. **Figure 1.** Demetra and Costa's ASV values in the sessions, plotted against the mean ASV value for the session (red line), 2 Standard Deviations below (yellow line) and 2 SD above (grey line) the mean.

#### *3.1. Interpersonal Physiological Synchrony*

#### *3.1. Interpersonal Physiological Synchrony* 3.1.1. IPS in Session 3

3.1.1. IPS in Session 3 The PDC analysis identified 29 time-windows in which the participants' ASV were synchronized in session 3, out of a total of 93 time-windows (Table 2 and Figure 2). This corresponds to at least two participants' physiological arousal being synchronized in 31,2% of the total session time. More specifically, Demetra's ASV values led Costa's ASV The PDC analysis identified 29 time-windows in which the participants' ASV were synchronized in session 3, out of a total of 93 time-windows (Table 2 and Figure 2). This corresponds to at least two participants' physiological arousal being synchronized in 31,2% of the total session time. More specifically, Demetra's ASV values led Costa's ASV in one timewindow, and the therapist's ASV in four. In contrast, Costa's ASV led Demetra's ASV in eight time-windows, and the therapist's ASV in nine. Lastly, the therapist's ASV led Demetra's ASV in six time-windows, and Costas' in eight. Overall, in session 3, Costas' autonomic arousal was found to lead IPS to a greater degree than Demetra's; moreover, the therapist had a leading role in several parts of the session, while Demetra primarily had a pacing role.

**Table 2.** Number of time-windows showing significant PDC synchronization between clients and therapist in session 3.


**Note**: Number of time-windows in session = 93. Time-windows in which at least two participants show significant PDC = 29.

periods in the session that are significant for the process of therapy.

in one time-window, and the therapist's ASV in four. In contrast, Costa's ASV led Demetra's ASV in eight time-windows, and the therapist's ASV in nine. Lastly, the therapist's ASV led Demetra's ASV in six time-windows, and Costas' in eight. Overall, in session 3, Costas' autonomic arousal was found to lead IPS to a greater degree than Demetra's; moreover, the therapist had a leading role in several parts of the session, while

**Table 2.** Number of time-windows showing significant PDC synchronization between clients and

Therapist 4 9 **Note**: Number of time-windows in session = 93. Time-windows in which at least two participants

term to describe time-windows in which more than one of the six possible directed synchronizations were observed. We consider these time-windows as particularly significant. Specifically, four time-windows showed increased physiological synchrony. Notably, in one time-window, all three participants were physiologically synchronized, with a mutual IPS between the clients and both clients' arousal also leading the therapist's ASV. Moreover, as can be seen in Figure 2, the time-windows with IPS tended to cluster

In addition, in session 3, several time-windows showed increased IPS; we use this

Demetra 8 6 Costas 1 8

**Leading Role Demetra Costas Therapist**

Demetra primarily had a pacing role.

therapist in session 3.

Pacing role

show significant PDC = 29.

**Figure 2.** Time-windows showing significant PDC in session 3. The height of each line indicates the number of pairs that show significant PDC within each time-window. The vertical lines denote the boundaries of the Topical Episodes. **Figure 2.** Time-windows showing significant PDC in session 3. The height of each line indicates the number of pairs that show significant PDC within each time-window. The vertical lines denote the boundaries of the Topical Episodes.

3.1.2. IPS in Session 14 The PDC analysis of the penultimate session identified 10, out of a total of 58, time-windows in which the participants' ANS arousal was synchronized (Table 3 and Figure 3). This corresponds to at least two participants' physiological arousal being synchronized in 17.2% of the total session time. More specifically, Demetra's arousal led Costas' ASV in four time-windows, and the therapist's ASV in one. Costas' arousal led Demetra's ASV in one time-window, and the therapist's in two. Lastly, the therapist's ASV led Demetra's ASV in four time-windows, and Costas' in one. Overall, IPS in the In addition, in session 3, several time-windows showed increased IPS; we use this term to describe time-windows in which more than one of the six possible directed synchronizations were observed. We consider these time-windows as particularly significant. Specifically, four time-windows showed increased physiological synchrony. Notably, in one time-window, all three participants were physiologically synchronized, with a mutual IPS between the clients and both clients' arousal also leading the therapist's ASV. Moreover, as can be seen in Figure 2, the time-windows with IPS tended to cluster around specific points in the session. We consider this clustering of IPS as reflecting time periods in the session that are significant for the process of therapy.

#### 3.1.2. IPS in Session 14

The PDC analysis of the penultimate session identified 10, out of a total of 58, timewindows in which the participants' ANS arousal was synchronized (Table 3 and Figure 3). This corresponds to at least two participants' physiological arousal being synchronized in 17.2% of the total session time. More specifically, Demetra's arousal led Costas' ASV in four time-windows, and the therapist's ASV in one. Costas' arousal led Demetra's ASV in one time-window, and the therapist's in two. Lastly, the therapist's ASV led Demetra's ASV in four time-windows, and Costas' in one. Overall, IPS in the penultimate session was equally led by Demetra and the therapist, and both clients had similar pacing roles, with Demetra pacing the therapist's ASV and Costas pacing Demetra's. Three time-windows showed increased physiological synchrony in this session, and again, time-windows with PDC tended to cluster together.

**Table 3.** Number of time-windows showing significant PDC synchronization between clients and therapist in the penultimate session.


**Note:** Number of time-windows in session = 58. Time-windows in which at least two participants show significant PDC = 10.

and again, time-windows with PDC tended to cluster together.

**Figure 3.** Time-windows showing significant PDC in the penultimate session (session 14). The height of each line indicates the number of pairs with significant PDC within each time-window. The vertical lines denote the boundaries of the Topical Episode. **Figure 3.** Time-windows showing significant PDC in the penultimate session (session 14). The height of each line indicates the number of pairs with significant PDC within each time-window. The vertical lines denote the boundaries of the Topical Episode.

penultimate session was equally led by Demetra and the therapist, and both clients had similar pacing roles, with Demetra pacing the therapist's ASV and Costas pacing Demetra's. Three time-windows showed increased physiological synchrony in this session,

**Table 3.** Number of time-windows showing significant PDC synchronization between clients and

Therapist 1 2 **Note:** Number of time-windows in session = 58. Time-windows in which at least two participants

Demetra 1 4 Costas 4 1

**Leading Role Demetra Costas Therapist**

#### 3.1.3. IPS between Sessions 3.1.3. IPS between Sessions

therapist in the penultimate session.

Pacing role

show significant PDC = 10.

A comparison of the PDC values between the two measurement sessions shows that (i) the total time spent in interpersonal physiological synchrony was significantly lower in session 14 as compared to session 3, and (ii) the global architecture of the interpersonal physiological synchrony network was reorganized to become more balanced as therapy progressed (Figure 4). As mentioned above, the percentage of the total session time with IPS decreased from 31.2% in session 3 to 17.2% in session 14. The IPS between the therapist and each of the clients showed the most marked decrease, from twenty-seven time-windows (28.1% of the session time) in session 3 to eight time-windows (13.7% of total session time) in session 14. This reduction in IPS over the course of therapy can be seen to reflect the clients' reduced affective arousal and their gradual disengagement from the process. As therapy progressed, the clients' difficult feelings and conflicts were expressed, elaborated upon and gradually reconstructed, and the physiological synchrony between the clients and the therapist decreased. This is in line with the characteristics of the closing stages of therapy, which entail less affectively charged and more A comparison of the PDC values between the two measurement sessions shows that (i) the total time spent in interpersonal physiological synchrony was significantly lower in session 14 as compared to session 3, and (ii) the global architecture of the interpersonal physiological synchrony network was reorganized to become more balanced as therapy progressed (Figure 4). As mentioned above, the percentage of the total session time with IPS decreased from 31.2% in session 3 to 17.2% in session 14. The IPS between the therapist and each of the clients showed the most marked decrease, from twenty-seven time-windows (28.1% of the session time) in session 3 to eight time-windows (13.7% of total session time) in session 14. This reduction in IPS over the course of therapy can be seen to reflect the clients' reduced affective arousal and their gradual disengagement from the process. As therapy progressed, the clients' difficult feelings and conflicts were expressed, elaborated upon and gradually reconstructed, and the physiological synchrony between the clients and the therapist decreased. This is in line with the characteristics of the closing stages of therapy, which entail less affectively charged and more reflective conversations, as well as a process of gradual disengagement from the therapeutic relationship and the work of therapy. *Entropy* **2022**, *24*, x FOR PEER REVIEW 11 of 18 eight. In contrast, in session 14, the overall synchrony was more equally driven by all participants, producing a more balanced or 'democratic' network structure (Figure 4). This finding points to the co-creation of a more equally distributed and balanced embodied relatedness between participants as therapy was reaching termination.

In order to deepen our understanding of the relational meaning of IPS as it fluctuated through a session, the clinical process in session 3 was qualitatively analyzed

this thematic coding allows researchers to identify key themes in a session and track the process of meaning co-construction, thus obtaining a relatively fine-grained description of meaning making through the session. Next, a qualitative analysis was performed to identify the significant moments in the session, which were defined as those topical episodes where: (a) important issues in the couple's life were introduced and narratively elaborated; (b) associated emotions were recognized, explored and expressed; and (c) the meaning of these key issues began to be reconstructed. The central theme in this session concerned Demetra's low mood and her strong ambivalent feelings regarding her role as a mother. Two topical episodes were identified through the qualitative analysis as entailing the elaboration of this theme, accompanied by intense emotional expression; these

The theme of Demetra's conflicts in her role as a mother was first introduced in TE4, approximately ten minutes into the session (duration 10′40″). This episode started with the therapist asking how the couple would choose to spend their time together if they had the finances and caretaking support. Costas made several suggestions that Demetra firmly rejected as she felt 'bored' with everything. Through the therapist's gentle curiosity and empathic questioning, Demetra's boredom was gradually reconstructed as entailing intense sadness; Demetra cried as she described her low mood, exhaustion, and sense of suffocation in her role as a mother. Towards the end of the episode, Costas gently talked about Demetra's lack of interest in sex, and this led to the expression of more sadness by Demetra. This episode contained the elaboration of the key theme of the session along with nonverbal displays of affect, as well as several markers of a moderately strong therapeutic alliance; this was particularly evident in the relationship between Demetra and the therapist, as well as within the couple (Table 4). The PDC analysis for this episode shows a cluster of five time-windows with IPS, accounting for 39% of the episode time; two of these time-windows show increased IPS, whereby more than two participants are in-sync (Figure 3). In other words, the findings from the PDC analysis concur with those of the narrative analysis and with the coding of the therapeutic alli-

**Figure 4.** Time-windows showing significant PDC per leading participant. **Figure 4.** Time-windows showing significant PDC per leading participant.

*3.2. Physiological Synchrony and Clinical Process*

are briefly described below.

Furthermore, IPS in a multi-actor setting such as couple therapy is more complex, as there are six possible pairs of participants. A shift was observed in how the IPS was distributed between participants; in session 3, the IPS was mainly driven by Costas, who led Demetra's and the therapist's ASV in seventeen time-windows and paced the therapist in eight. In contrast, in session 14, the overall synchrony was more equally driven by all participants, producing a more balanced or 'democratic' network structure (Figure 4). This finding points to the co-creation of a more equally distributed and balanced embodied relatedness between participants as therapy was reaching termination.

#### *3.2. Physiological Synchrony and Clinical Process*

In order to deepen our understanding of the relational meaning of IPS as it fluctuated through a session, the clinical process in session 3 was qualitatively analyzed drawing upon narrative principles, and the findings were subsequently examined in relation to the PDC analysis. In brief, the session was segmented into topical episodes [83]; this thematic coding allows researchers to identify key themes in a session and track the process of meaning co-construction, thus obtaining a relatively fine-grained description of meaning making through the session. Next, a qualitative analysis was performed to identify the significant moments in the session, which were defined as those topical episodes where: (a) important issues in the couple's life were introduced and narratively elaborated; (b) associated emotions were recognized, explored and expressed; and (c) the meaning of these key issues began to be reconstructed. The central theme in this session concerned Demetra's low mood and her strong ambivalent feelings regarding her role as a mother. Two topical episodes were identified through the qualitative analysis as entailing the elaboration of this theme, accompanied by intense emotional expression; these are briefly described below.

The theme of Demetra's conflicts in her role as a mother was first introduced in TE4, approximately ten min into the session (duration 1004000). This episode started with the therapist asking how the couple would choose to spend their time together if they had the finances and caretaking support. Costas made several suggestions that Demetra firmly rejected as she felt 'bored' with everything. Through the therapist's gentle curiosity and empathic questioning, Demetra's boredom was gradually reconstructed as entailing intense sadness; Demetra cried as she described her low mood, exhaustion, and sense of suffocation in her role as a mother. Towards the end of the episode, Costas gently talked about Demetra's lack of interest in sex, and this led to the expression of more sadness by Demetra. This episode contained the elaboration of the key theme of the session along with nonverbal displays of affect, as well as several markers of a moderately strong therapeutic alliance; this was particularly evident in the relationship between Demetra and the therapist, as well as within the couple (Table 4). The PDC analysis for this episode shows a cluster of five time-windows with IPS, accounting for 39% of the episode time; two of these timewindows show increased IPS, whereby more than two participants are in-sync (Figure 3). In other words, the findings from the PDC analysis concur with those of the narrative analysis and with the coding of the therapeutic alliance, suggesting that this topical episode was important for the therapy process on both semantic and embodied levels.

The same theme was further elaborated with increased emotional expression in TE7. This long episode (duration 14040") took place in the middle of the session. It started with Demetra crying as she described feeling trapped and suffocating in her relationship with their baby; she vividly described her frustration and rage towards their baby, the wish to hurt him that she sometimes experienced, her angry outbursts towards him, and the intense guilt that she felt after such outbursts. As she listened to Demetra's emotional narrative, the therapist displayed many nonverbal signs of affiliation and empathy. She also introduced the hypothesis that Demetra's difficulties and sense of failure result from comparing herself to an idealized version of motherhood; this was followed by a productive conversation that challenged Demetra's idealization of her own mother, as well as the dominant discourse of ideal motherhood [60]. In terms of the alliance, this episode contained markers of a

moderately strong alliance between Demetra and the therapist, whereas Costas displayed some markers of difficulty in the alliance. Based on the PDC analysis, there were five time-windows with IPS in this episode (including two windows with increased IPS), and these account for 28% of the episode time. Notably, this topical episode contained one time-window during which all three participants were physiologically synchronized. This corresponds to the point in the session where Demetra cried as she talked about the rage and guilt she felt towards their baby. A cluster of time-windows with IPS can be seen at the end of the episode, as Demetra's sadness and sense of suffocation in her maternal role were expressed. Once again, in this topical episode, the PDC analysis identified a period in the session that entailed intense affective expression by Demetra, empathic responsiveness by the therapist and a strong therapeutic alliance.

**Table 4.** Comparison of SOFTA scores and number of PDC identified significant time-windows for each Topical Episode.


Next, we examined a topical episode identified as significant through the PDC, but not through the qualitative analysis: TE12 contained a cluster of time-windows with IPS that account for 66% of the episode time. The episode took place towards the end of the session (duration 5020") and focused on Demetra's lack of sexual desire, a delicate and affectively charged issue in the couple's relationship. In this instance, the PDC analysis identified a part of the session that was not identified as important from a qualitative perspective, but which proved to be significant later in the treatment.

In sum, the two topical episodes that were identified as important for the process of therapy through the qualitative analysis also entailed increased IPS and increased ratings of the therapeutic alliance, as compared to the rest of the session. These findings are in line with the literature that points to the role of IPS in empathy, affiliation, rapport and the therapeutic alliance [7,11,43]. At the same time, PDC analysis proved useful in identifying significant moments in the session.

In order to explore the relational significance of IPS, we examine the findings from the PDC analysis in relation to the interactions between participants in session 3. With regards to the therapist–client(s) interaction, periods of physiological synchrony between the therapist and Demetra account for 10.1% of the session time (corresponding to 27.8% of the total time with IPS); IPS between the therapist and Costas is significantly higher, accounting for 18.3% of the total session time (corresponding to 47% of the total time with IPS). This points to the presence of a more intense affective connection between Costas and the therapist, and this is in line with the clients' respective SRS scores (Table 1). This finding reflects the complex interplay between explicit and implicit aspects of interaction in psychotherapy. More specifically, the therapist was very responsive to Demetra's distress on an explicit level, as she openly expressed empathy and affiliation with Demetra's painful

conflicts regarding motherhood. At the same time, she was significantly more in-sync with Costas on an embodied level. This is in line with findings that behavioral and physiological synchrony seem to be independent processes that do not always co-occur [84]. The therapist was able to maintain a balanced therapeutic alliance with both members of the couple, and this was achieved through different modalities. In other words, implicit and explicit modalities of interaction were used to manage different interactional goals [30]. This is important, as split alliances, i.e., situations where the therapist takes sides by colluding with one partner, have a negative impact on the outcome of therapy [79]. In addition, there is some evidence that the therapeutic alliance with the male partner is critical for therapy outcome in heterosexual couples [70].

With regards to the couple's relationship, Costas' physiological arousal led Demetra's ASV significantly more than the opposite. Thus, although Demetra's difficulties were central on the level of talk, on an implicit embodied level, Costas had a more powerful influence on the interaction. Again, this is an observation that illustrates the complex interplay between the verbal and embodied aspects of psychotherapeutic interaction. A possible interpretation of this observation is that of affect co-regulation, where Costas' presence could be seen to regulate Demetra's affective arousal, as often happens in complementary, i.e., homeostatic, couple relationships [36]. A closer examination of the level of both partners' arousal, the valence of their affective displays, and an examination of the talk in these episodes would help contextualize and interpret these observations more fully.

#### **4. Conclusions**

The findings of this study point to the complex interplay between explicit and implicit levels of interaction and the potential added value of including physiological synchrony in the study of interactional processes in couple therapy [36,55,72,85,86]. In line with contemporary theories of therapeutic change, a key assumption of this work is that psychotherapy entails processes of intersubjective meaning making that take place through different modalities and, presumably, with different degrees of conscious awareness [23]. From this perspective, including measures of physiological activation in the study of psychotherapy sessions can help examine implicit, embodied interactional processes that contribute significantly to the formation of the therapeutic alliance, the co-creation of new meanings and, ultimately, therapeutic change. Although several research methods have been developed to study the talk in interaction [48,58], these methods generally fail to grasp the implicit, procedural level of interaction. Our attempt to include measures of autonomic arousal in studies of the therapy process and to operationalize implicit interactional processes of embodied responsiveness are in the spirit of exploring ways to include the implicit realm when studying the psychotherapy process [23].

Research to date suggests that IPS reflects different interactional processes, and these need to be disentangled for the field to progress [11,34,54]. In our study, we propose the use of a windowed Partial Directed Coherence-based approach as a metric to calculate physiological synchrony, as this allows a more nuanced examination of the dynamic nature of IPS in psychotherapy sessions. PDC analysis allows us to examine the therapy process in specific interactional events in the session, and this *micro-focus* provides a more finegrained description of interactional dynamics as they develop, thus allowing a more nuanced interpretation of the role of IPS in the therapy process. Importantly, PDC analysis allows us to examine the *directionality* of synchronous interactions, which again adds another layer of complexity to our understanding of the role of physiological synchrony in the therapy process. Therefore, the proposed approach models the couple's interactions within the setting of a therapy as a self-organizing system, a system that is both open and complex, exchanging energy and information between its component parts and with its surroundings [87]. This exchange may be synchronic and diachronic, in spatial distribution and time transitions, therefore demanding multidimensional theoretical models to represent its hybrid nature [88].

A key aim of this study was to explore the links between a quantitative approach to the study of IPS and the characteristics of interactional dynamics and the clinical process, and this mixed-method analysis produced promising initial findings. More specifically, it examined shifts in IPS between the start and end of therapy in a successful couple therapy and identified a reduction in IPS as therapy progressed. This decrease primarily concerned the therapist–client(s) interaction and was interpreted as a reflection of progress, in the sense of a decrease in the intensity of negative affects expressed by the clients and the need for therapist empathy, as well as the couple's gradual disengagement from the process of therapy in line with the termination phase. In addition, the network of IPS between the three participants became more balanced. Both these findings are in line with a good therapy outcome, and as such, they provide support for the clinical validity of PDC analysis.

The main limitation, inherent in this approach, is that only one couple is included; hence, the descriptive outcome of the study cannot be generalized. Nonetheless, we propose that this detailed analysis provides a necessary step for evaluating the usefulness of employing PDC analysis to examine IPS in therapy sessions, which can now be further elaborated. Studying IPS via a windowed PDC approach may lead to an even more detailed identification of the underlying processes if the characteristics of the ANS signals during significant time-windows are further investigated. In addition, calculating positive vs. negative correlations of ANS activity or specific patterns of ANS reactivity within the significant time-windows may be used in future studies to examine their associations with different intersubjective processes, such as empathy, alliance or affect contagion. We hope that future work in this field will exploit the strengths of the PDC analysis and further our understanding of the embodied, relational aspects of the therapy process.

**Author Contributions:** Conceptualization, E.A. and E.P.; literature review: E.A., E.P., C.L. and F.P.; methodology, E.A., E.P. and P.K.; formal analysis, E.A., E.P., P.K. and C.L.; original draft preparation, E.A.; writing—review and editing, E.A., E.P., C.L. and F.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of the Scientific Board of the Psychiatric Hospital of Thessaloniki in June 2013.

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

**Data Availability Statement:** The ANS data presented in this study are available on request from the corresponding author, and the data from the video material are not publicly available due to ethical restrictions.

**Acknowledgments:** We would like to thank Vasileia Lerou, Katerina Lazaridou and Kalliopi Papadopoulou for participating in the coding of the sessions with the SOFTA and the segmentation into Topical Episodes. We would also like to thank V. Escudero for useful comments on the use of the SOFTA analysis of the topical episodes.

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

#### **References**


## *Article* **Recurrence-Based Synchronization Analysis of Weakly Coupled Bursting Neurons under External ELF Fields**

**Aissatou Mboussi Nkomidio 1,†, Eulalie Ketchamen Ngamga <sup>2</sup> , Blaise Romeo Nana Nbendjo <sup>1</sup> , Jürgen Kurths 2,3 and Norbert Marwan 2,\***


**Abstract:** We 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. We find, by numerical simulations, that neurons can exhibit different spiking patterns, which are well observed in the structure of the recurrence plot (RP). We 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 by 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, phase synchronization of the coupled neurons occurs and, by adding an ELF electric field, this synchronization increases depending on the amplitude of the externally applied ELF electric field. We further suggest a novel measure for RP-based phase synchronization analysis, which better takes into account the probabilities of recurrences.

**Keywords:** neuron; electric field; weak coupling; gap junction; synchronization; recurrence plot

## **1. Introduction**

Action potentials, or spikes, are responsible for the transmission of information through the nervous system [1]. A neuron can generate various temporal patterns of spike signals when it is driven by stimuli or noise from both internal or external environments. Therefore, analyzing spiking patterns of neurons under different stimulations plays an important role in the exploration of the encoding and decoding mechanism of information for neurons. External environmental stimuli in the brain can be of various origins, such as a wide utilization of power lines or electrical equipment. Electromagnetic exposure in the environment today is nearly one hundred times stronger than in previous centuries and many neuronal diseases are probably caused by electromagnetic exposure, as reported by Huang et al. [2]. Experiments with transcranial electrical stimulation have shown that electric field magnitudes in the cortex can be as high as 0.4 mV/mm for a 1 mA stimulation current. For typical electrode configurations used in clinical trials, maximal field intensities of up to 0.8 mV/mm were found when applying 2 mA. More extended areas can reach values of 0.28 mV/mm (95th percentile) under 2 mA stimulation [3–5].

An electromagnetic field can affect the neuron sensibility [6–9]. It also exhibits the excitability of many nerve cells, such as hippocampal cells, or cortical neurons [10,11]. Neurons exposed to an electromagnetic field can change the normal firing properties, which may lead to many neural diseases such as amyotrophic lateral sclerosis, senile dementia, Parkinson's disease, and Alzheimer's disease [7,12–14].

**Citation:** Nkomidio, A.M.; Ngamga, E.K.; Nbendjo, B.R.N.; Kurths, J.; Marwan, N. Recurrence-Based Synchronization Analysis of Weakly Coupled Bursting Neurons under External ELF Fields. *Entropy* **2022**, *24*, 235. https://doi.org/10.3390/ e24020235

Academic Editors: Franco Orsucci and Wolfgang Tschacher

Received: 30 December 2021 Accepted: 1 February 2022 Published: 3 February 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/).

On the other hand, neurons are strongly coupled in the brain, and they need to synchronize information to encode and decode. Synchronization is a universal concept of nonlinear dynamics [15]. In the brain system, synchronization is a typical form of group motion rhythm, which means the neurons discharge at the same time or their discharge rhythms have at least some kind of relationship [16,17]. Neuronal synchrony activities can be found not only among coupled neuron groups in the same brain region but also among uncoupled neuron groups in the same brain region or among different cortical areas; moreover, synchronization can cross over two hemispheres of the brain [18]. Synchronization processes are crucially important for the neuronal system, and wellcoordinated synchrony within and between neuronal populations appears to play an important role in neuronal signaling and information processing.

To study synchronization between neurons, different models of neuron dynamics have been developed, such as the Hodgkin–Huxley (HH) model and all the models derived from it. One of the derived models is the Morris–Lecar (ML) model [19,20]. It has the advantage of exhibiting class I and II neurons. Most studies on neuron synchronization use the Morris–Lecar model under an external electric field. For example, Kitajima and Kurths [21] investigated forced synchronization of electrically coupled class I and class II neurons under different coupling strengths. It was found that class II neurons have a wide parameter region of forced synchronization. However, in general, such studies did not consider the effect of small variations of the coupling strength between neurons.

The assumption of weak neuronal connection is based on the observation that the typical size of a postsynaptic potential is less than 1 mV, which is small in comparison to the mean potential necessary to discharge a cell or the average value of the action potential [22]. In a study of synaptic organization and dynamical properties of weakly connected neuronal oscillators, Hoppensteadt and Izhikevich [23] showed some phase synchronization between neurons in this range of coupling. Moreover, Izhikevich [24] studied the synchronization of elliptic bursters in a range of weak connectivity and found that such weakly connected bursters need few bursts to synchronize and synchronization is possible for bursters having quite different quantitative features. These phenomena were found in different neuron models, such as the FitzHugh–Rinzel, ML, and HH models.

The important question is if, even in the range of small coupling strength, a pair of neurons weakly coupled with gap junction are able to synchronize under the effect of an electric field (EF). Because electromagnetic stimulation can cause many disorders in the neural system, the theoretical investigation of the impact of an external EF on the synchronization of weakly coupled neurons is an important step to understand what happens in the brain during this exposure. Thus, in this work, we study the synchronization of a pair of ML neurons weakly connected with gap junction under an externally applied extremely low frequency (ELF) EF. Here, extremely low frequency means a frequency range between 0 and 10 Hz. Mammalian neurons show intrinsic resonance with frequency selectivity for inputs within the range from 4 to 10 Hz [25–30]. Gap junctions (channels that physically connect adjacent cells) provide an efficient and extremely fast way to propagate those signals between neurons [31,32]. In contrast, signal transmissions via chemical synapses have a significant delay (in the order of milliseconds) [33] and are not fast enough to respond to the EF. Therefore, we consider here coupling via gap junction, allowing direct response to the ELF EF. Using recurrence plot-based time series analysis, we investigate how the applied EF affects the condition of synchronization of the coupled neurons. This specific method has the advantage of being able to compare the phases of chaotic weakly coupled systems, even within noncoherent regimes or for spiky signals [34].

Recurrence plots (RPs) represent manifold recurrence features of a dynamical system in phase space [35] and are widely applied in the field of neuroscience [36–42]. For example, RPs can differentiate the stochastic and deterministic dynamics of irregularly firing cortical neurons [43] and show the average dynamics within a network of synchronized neurons [44] or spontaneous activity in neuronal in vitro cultures [45]. They are also powerful tools to study inter-relationships, coupling directions, phase synchronization, and

generalized synchronization [34,46–48] and have been applied in different fields, such as chemistry, engineering, physiology, financial markets, and climatology [41,49–53]. Based on EEG measurements, the joint recurrence and the correlation of probability of recurrence were used to reconstruct brain networks [54,55]. A similar approach was used to study the synchronization between neurons based on the Hindmarsh–Rose model [39].

The correlation of probability of recurrence is a commonly used measure for recurrencebased phase synchronization analysis [34,50]. However, this measure is based on Pearson correlation and, thus, has a methodological concern because of the spiky nature of the probability of recurrence.

In this study, we first formulate the mathematical modeling of a single ML neuron and present typical neuron bursting patterns under varying ELF EFs and their corresponding recurrence features as obtained by RPs. We then study the synchronization of two weakly coupled chaotic bursting neurons with and without the influence of an ELF EF. We present the effect of EFs on the mismatch of the mean frequencies of both neurons, even when they are weakly connected. For this purpose, we suggest a slight modification of the recurrence-based phase synchronization measure.

#### **2. Model**

#### *2.1. Morris–Lecar Neuron Model under an Extremely Low Frequency Electric Field*

The ML neuron model is a model for electrical activity in the barnacle muscle fiber [19]. It is a simplified version of the HH neuron model for describing the discharge and the refractory properties of real neurons. It can explain the dynamical and biophysical mechanisms of the action potential initiation. This model is chosen as a compromise between a realistic representation of neuronal dynamics and an analytically tractable system. Furthermore, it has an advantage in that the excitability of types I and II can be obtained with a single parameter change. It can also exhibit a variety of bursting types involving regular bursting or irregular bursting and complex bifurcation structures [20,56–58].

The ML model has a fast activation variable *v* (membrane voltage) and a slower recovery variable *w*. *v* represents voltage (expressed in mV) and controls the instantaneous activation of fast currents (*i*fast); *w* is a function of *v* and controls the activation of slower currents (*i*slow). *c dv dt* is the current flowing through the capacitor related to the variation of ionic density between external and internal faces of the membrane. *i*fast, *i*slow, and *i*leak are ionic currents characterizing the movement of charged particles through the ion channels. This movement of charged particles is due to the opening and closing of each ion channel. *i*stim and *c* are the external input current and the membrane capacity, respectively. Finally, this model is given by the following equations:

$$c\frac{dv}{dt}\_{\cdot} = \left.i\_{\rm stim} - i\_{\rm fast} - i\_{\rm leak} - i\_{\rm slow} \right. \tag{1}$$

$$\frac{dw}{dt} \quad = \quad \varphi \frac{m\_2(v) - w}{b(v)}\tag{2}$$

with the currents

$$\begin{array}{rcl} i\_{\text{fast}} & = & g\_{\text{fast}} m\_1(v) (v - e\_{\text{Na}}), \\ i\_{\text{slow}} & = & g\_{\text{slow}} w (v - e\_{\text{K}}) \\ i\_{\text{leak}} & = & g\_{\text{leak}} (v - e\_{\text{leak}}). \end{array}$$

The parameters *e*Na, *e*K, and *e*leak represent the equilibrium potentials of Na+, K+, and leak ions, respectively, and *g*fast, *g*slow, and *g*leak are the maximal conductances of the corresponding ion currents. They reflect the ion channels' densities distributed over the

membranes. Control parameter *ϕ* is used to control the rate of change of *w*. The steady states *m*<sup>1</sup> and *m*<sup>2</sup> are nonlinear functions of *v*, given by

$$m\_1(v) \quad = \quad 0.5 \Big( 1 + \tanh \frac{v - u\_1}{u\_2} \Big) \tag{3}$$

$$m\_2(v) \quad = \quad 0.5 \Big( 1 + \tanh \frac{v - \mu\_3}{\mu\_4} \Big) \tag{4}$$

$$b(v) \quad = \frac{1}{\cosh \frac{v - \mu\_3}{2\mu\_4}}.\tag{5}$$

*u*<sup>1</sup> and *u*<sup>3</sup> are the activation midpoint potentials at which the corresponding currents are half activated. *u*<sup>2</sup> and *u*<sup>4</sup> denote the slope factors of the activation. The time constant of the potassium activation is *b*. When a time-varying ELF EF is applied to the brain, it can induce a charge movement in the brain tissue; in which case, the current flow occurs mostly in the extracellular medium [59]. Therefore, an external EF will induce a membrane depolarization ∆*v* which will modulate neuronal bursting behavior. For the sake of simplicity, we consider a steady external sinusoidal electrical field

$$v\_{\text{e}} = \frac{A}{\omega} \sin \omega t + V\_{\text{E}} \tag{6}$$

where *V*<sup>E</sup> is the direct voltage, *A* the amplitude, and *ω* the frequency of the ELF EF. The field-induced membrane depolarization ∆*v* can be expressed by [60]

$$
\Delta v = \frac{A}{\omega} \frac{\sin \omega t - \cos \omega t}{1 + \left(\omega t\_1\right)^2} + V\_{\text{E}} \tag{7}
$$

with *t*<sup>1</sup> significantly small and the frequency in the extremely low frequency area *ωt*<sup>1</sup> 1. Thereby, Equation (7) can be simplified to

$$
\Delta v = \frac{A}{\omega} \sin \omega t + V\_{\rm E}.\tag{8}
$$

According to Equation (8), the sinusoidal EF *v*<sup>e</sup> equals its field-induced membrane depolarization ∆*v*. Considering that ∆*v* acts as an additive perturbation to the membrane potential, the dynamics of a neuron during exposure can be described by [61]

$$c\frac{dv}{dt}\_{\!\!\!\/} = \left.i\_{\text{stim}} - \frac{d\Delta v}{dt}\_{\!\!\/} - i\_{\text{fast}} - i\_{\text{leak}} - i\_{\text{slow}}\right.\tag{9}$$

$$\frac{dw}{dt} \quad = \quad \varphi \frac{m\_2(v) - w}{b(v)}\tag{10}$$

with

$$\begin{array}{rcl} i\_{\text{fast}} & = & g\_{\text{fast}} m\_1(v)(v + \Delta v - e\_{\text{Na}}), \\ i\_{\text{slow}} & = & g\_{\text{slow}} w(v + \Delta v - e\_{\text{K}}) \\ i\_{\text{leak}} & = & g\_{\text{leak}} (v + \Delta v - e\_{\text{leak}}). \end{array}$$

We assume that the synaptic input current *i*stim = 0 in order to study the response of a cortical neuron model exposed to an external sinusoidal field. Throughout this paper, we use the same parameter values for the ML model as explained in Table 1 [62].


**Table 1.** Parameters used for the ML model.

#### *2.2. Bursting Patterns of a Neuron*

To explore how the neuron model responds to the externally applied ELF EF, we study the dynamics described by Equations (9) and (10) under the sinusoidal stimulus *v*e, Equation (6). The simulations are implemented using the 4th-order Runge–Kutta method with a time step of 0.01 ms. Initial conditions are chosen as the resting values of membrane voltage in the absence of stimuli, that is, *v*(0) = −65 mV and *w*(0) = 0. The length of the time series is 2000 ms. The response of a neuron induced by an EF depends on the EF's frequency *ω* (Figure 1 for *ω* in the range 0 ≤ *ω* ≤ 0.5 rad/ms). The amplitude of the external EF is set very small, *A* = 0.1.

The firing pattern of a neuron stimulated by an external EF varies when changing the frequency *ω* (Figures 1 and 2). For an ELF EF with very low frequency, the neuron fires periodically. We find *n* spike bursting states, and the number *n* can be large. *n* spike bursting means that we have *n* action potentials in every stimulus period (Figures 1A–D and 2) for *ω* = [0.001, 0.120] rad/ms. After this range *ω* of *n*-periodic bursting, the neuron bursts synchronously to the stimulus *ω* = [0.121, 0.280] rad/ms. After this range of *ω*, the neuron exhibits a chaotic response (Figure 1E) with *ω* = [0.281, 0.320] rad/ms where the membrane potential responses are aperiodic and irregular. As *ω* is further increased, a mode locking pattern of bursting appears (Figure 1F), finally followed by synchronized firing with only one action potential in every stimulus, which can maintain this state for a long-term frequency band (Figure 1G). Neuron dynamics are obviously very sensitive to the frequency of the stimulus by the ELF EF.

**Figure 1.** Spiking patterns of ML neuron membrane voltage under an external EF for different frequencies: (**A**) *ω* = 0.02 rad/ms, (**B**) *ω* = 0.05 rad/ms, (**C**) *ω* = 0.06 rad/ms, (**D**) *ω* = 0.10 rad/ms, (**E**) *ω* = 0.286 rad/ms, (**F**) *ω* = 0.32 rad/ms, and (**G**) *ω* = 0.5 rad/ms.

**Figure 2.** Bifurcation diagram of ML neuron dynamics under an external EF for varying frequencies *ω* (based on interspike intervals of the membrane voltage). For better visibility of the dynamics for larger *ω*, the *y*-axis was bounded to 250 ms.

#### **3. Recurrence Quantification Analysis (RQA)**

In the following, neuron dynamics will be studied using recurrence quantification analysis (RQA). This method quantifies certain recurrence features of the dynamical system in its corresponding phase space [35,63]. We define a recurrence of a trajectory <sup>~</sup>*x*(*t*) <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* (with *m* the dimension of the system) of a dynamical system by saying that the trajectory has returned at time *t* = *j* to the former point in phase space visited at *t* = *i* (with *i* ∈ [1, *N*] and *N* the length of time series) if

$$R\_{\vec{i},\vec{j}} = \Theta\left(\varepsilon - \left||\vec{x}(i) - \vec{x}(j)\right||\right) \tag{11}$$

where *ε* is a pre-defined threshold and Θ(·) is the Heaviside function. We have a matrix of (0, 1), where 1 at (*i*, *j*) means that ~*x*(*i*) and ~*x*(*j*) are neighbors and 0 means that they are not. The resulting black and white representation of this binary matrix is called a recurrence plot (RP). For the selection of the recurrence threshold *ε*, different strategies are available, depending on the research question [64–69]. Here, we use an approach to select *ε* in a way that ensures a certain recurrence point density. This allows a better comparability between RPs of different systems [68].

The RP method has been intensively studied and applied in the last years. Different measures of complexity have been proposed that can classify different dynamics, identify dynamical transitions, or detect couplings, causality, or synchronization [35].

If not all state variables of the state vector ~*x* are available, a phase space reconstruction has to be applied. Here, we use the recently proposed PECUZAL method to reconstruct the phase space trajectories [70]. This method allows us to use multiple embedding delays *τ*. The embedding parameters are listed in Table 2.


**Table 2.** Embedding parameters indicated by the PECUZAL algorithm.

RPs of the different bursting neurons represent a typical pattern (Figure 3, using an *ε* that ensures a recurrence point density of 0.15). Each "dashed-dotted" diagonal line in the RPs corresponds to a spike. For the alternating spiking behavior, we have a set of dashed lines followed by an extended black region (Figure 3A,B). The set of *n* spikes is well distinguished by the number of dashed lines (see some orange boxes marked in the figure). The block-like black region represents the silent state between each stimulus, which is a period for which the neuron cannot respond to a stimulus. On the small scale, the diagonal lines show some additional patterns, i.e., small structures sitting perpendicularly

at these lines or thickenings, similar to bumps or knobs. This is a typical feature of slow–fast systems [71].

**Figure 3.** RPs of the membrane voltage *v* of selected bursting neurons: (**A**) 4 spike burst (*ω* = 0.05 rad/ms), (**B**) 2 spike burst (*ω* = 0.10 rad/ms), (**C**) chaos (*ω* = 0.286 rad/ms), and (**D**) 1 spike (*ω* = 0.5 rad/ms). Diagonal lines represent the spikes, the larger extended structures represent the "silent" epochs, and the structures perpendicular to the diagonal lines and small thickenings represent the slow–fast dynamics (blue circle in (**D**)). The orange boxes in (**A**,**B**) mark a sequence of diagonal lines. The number of diagonal lines counted from the main diagonal of such a box towards the corner of this box represents the number of spikes within this period. Recurrence threshold *ε* is selected to ensure a recurrence point density of 0.15.

In order to go beyond the visual impression of the RP, we use recurrence quantification analysis (RQA) [35,72]. The RQA measures are based on the recurrence point density and the diagonal and vertical line structures of the RP. For example, the recurrence point density <sup>1</sup> *<sup>N</sup>*<sup>2</sup> ∑ *Ri*,*<sup>j</sup>* corresponds to the probability that a state will recur. The calculation of this measure can also be restricted to a diagonal-wise calculation, i.e., the recurrence point density along a diagonal with distance *τ* from the main diagonal *Ri*,*<sup>i</sup>* = 1 [35]. This gives us an estimator of the probability that the system returns to a previous state after time *τ* and is called the *τ*-recurrence rate,

$$RR\_{\tau} = \frac{1}{N - \tau} \sum\_{i=0}^{N-\tau} R\_{i, i+\tau} \tag{12}$$

where *τ* is the set time and *N* the total number of points in the phase space. The distance between the peaks in an *RR<sup>τ</sup>* plot corresponds to the period length of oscillations or the interspike intervals of spike trains similar to the neuron's spiking/bursting patterns.

The spike trains of 4 spikes, 2 spikes, chaos, and 1 spike have their specific probability distributions for recurrence after lag *τ* (Figure 4). Where the 1-periodic spike occurrence is clearly visible for 1 spike (Figure 4D), the 2 and 4 spikes produce more subtle probability distributions, revealing different periodicities and large blocks between the bursting periods (Figure 4A,B). The *RR<sup>τ</sup>* of the chaotic bursting exhibits a more complicated distribution of peaks corresponding to the unpredictable occurrence of spikes (Figure 4C).

**Figure 4.** Probability of recurrence after time *τ* (*τ*-recurrence rate) for the bursting neurons as shown in Figure 3: (**A**) 4 spike burst, (**B**) 2 spike burst, (**C**) chaos, and (**D**) 1 spike. The *n* bursts are visible as the rather thin side peaks of the main peaks (in addition to the main peak).

#### **4. Coupling of Two Bursting Neurons**

#### *4.1. Model and Numerical Simulation*

*c*

A coupling between ML neurons is realized by a gap junction. We suppose that the two neurons are slightly different by considering different values of *u*<sup>2</sup> and *u*<sup>3</sup> in Equations (4) and (5), i.e., *u*2,1 = −18.0 mV and *u*2,2 = −18.1 mV, and *u*3,1 = −12.8 mV and *u*3,2 = −10 mV for neurons 1 and 2, respectively. Moreover, both neurons start using different initial conditions *v*1(0) = −65.6 mV and *v*2(0) = −60 mV. We integrate the model for 50,000 time steps (with *dt* = 0.05) and remove the first 10,000 points as transients. Using the ELF EF frequency that leads to chaotic bursting *ω* = 0.286 rad/ms, the coupled chaotic bursting ML neurons under ELF EF exposure can be expressed as

$$c\frac{dv\_1}{dt}\quad =\quad i\_{\rm stim} - \frac{d\Delta v\_1}{dt} - i\_{1,\rm fast} - i\_{1,\rm slow} - i\_{1,\rm leak} - g(v\_1 - v\_2) \tag{13}$$

$$\frac{dw\_1}{dt} = -\varrho \frac{m\_2(v\_1) - w\_1}{b(v\_1)}\tag{14}$$

$$\frac{d\upsilon\_2}{dt} = \dot{\imath}\_{\text{stim}} - \frac{d\Delta\upsilon\_2}{dt} - \dot{\imath}\_{2,\text{fast}} - \dot{\imath}\_{2,\text{slow}} - \dot{\imath}\_{2,\text{least}} - g(\upsilon\_2 - \upsilon\_1) \tag{15}$$

$$\frac{dw\_2}{dt} = -\varrho \frac{m\_2(v\_2) - w\_2}{b(v\_2)},\tag{16}$$

where *g* is the gap which represents the electrical junction between the neurons. With these two different chaotic neurons, we will now study the phase synchronization between them and focus on the range of weak coupling, i.e., with 500 values of *g* within the interval *g* = [0, 0.15] (Figure 5).

**Figure 5.** Spiking pattern in the membrane voltage of two weakly coupled neurons in an ELF EF in chaotic regime (*ω* = 0.286) with (**A**) no synchronization with *g* = 0.01 and (**B**) phase synchronization with *g* = 0.04.

When bursting begins at the same time in the coupled neurons, we have bursting synchronization irrespective of the neurons' spiking behaviors within a given burst event. From a dynamical point of view, since we assign a phase that increases by 2*π* at each burst event, we regard bursting synchronization as a kind of phase synchronization [73]. Thus, first we will determine the phase of each chaotic bursting neuron. A frequently used approach to calculate the phase of a signal is using the Hilbert transformation [15]

$$\phi(t) = \arctan2\left(v\_H(t), v(t)\right) \tag{17}$$

where *v<sup>H</sup>* is the complex part of the Hilbert transform of the membrane voltage *v*(*t*) and *φ*(*t*) increases continually with time. Since chaotic neurons have chaotic spikes, the phase of chaotic neurons changes also chaotically. Unfortunately, this approach does not work well for spiky signals and can cause slipping of the instantaneous phases. Nevertheless, for long-term averages, it provides useful results.

To detect phase synchronization of chaotic coupled neurons and to evaluate the effect of an ELF EF on this synchronization, we first consider the absolute phase difference between the membrane voltage of both neurons without and with applied EF. Phase synchronization occurs if the difference *φ*1(*t*) − *φ*2(*t*) between the phases of the two neurons does not grow with time [74]. This means that the two neurons, on average, generate spikes almost simultaneously. With the knowledge of the phase *φ*(*t*), the frequency *ω*¯(*t*) = *<sup>d</sup>φ*(*t*) *dt* and the mean frequency Ω = D *dφ*(*t*) *dt* <sup>E</sup> can be defined. A weaker form of synchronization is frequency locking. Frequency locking between coupled systems can be measured by the mismatch between the average frequencies ∆Ω = Ω<sup>1</sup> − Ω2, with ∆Ω → 0 for phase locking.

The weakly coupled neurons show frequency locking without ELF EF when the coupling exceeds a critical value (Figure 6). The frequency mismatch ∆Ω between both neurons is constant between *g* = 0 (no coupling) and *g* = 0.025. After this value, ∆Ω is decreasing and vanishes around *g* = 0.066, indicating the onset of synchronization frequency locking between the neurons.

With ELF EF applied, the frequency difference is smaller, even for *g* = 0, and decreases much faster than without EF; the neurons become frequency-locked for *g* = 0.037 (Figure 6). Thus, in a range of weakly connected neurons, applying an external ELF EF on the chaotic coupled ML neurons enhances frequency-locked synchronization. This confirms earlier findings of synchronized neurons using a different model of weakly connected bursters [24].

**Figure 6.** Frequency mismatch between chaotic coupled neurons for increasing coupling *g* without external EF (red) and with external EF where *A* = 0.1 and *ω* = 0.286 (blue).

Since the firing pattern strongly depends on the amplitude of the ELF EF, we expect that the occurrence of frequency locking also depends on this external stimulus amplitude. In fact, we find that an amplitude value of *A* = 0.15 is strong enough to cause a complete synchronization of two neurons even without coupling (Figure 7). Therefore, we select a lower amplitude value of *A* = 0.1, where we still have a significant frequency mismatch between the uncoupled neurons. A weak coupling between the neurons leads, finally, to frequency-locked synchronization for lower ELF EF amplitudes.

**Figure 7.** Frequency mismatch between chaotic coupled neurons for increasing amplitude *A* and for different coupling strengths *g*.

To test for phase synchronization, i.e., whether the difference *φ*1(*t*) − *φ*2(*t*) remains constant, we will use an alternative method which can derive the phases of spiky signals in a more reliable way.

#### *4.2. Phase Synchronization Analysis Using Recurrence Features*

Phase synchronization is related to recurrence of states. Therefore, RPs are a natural tool to study phase synchronization [35]. The spiking pattern causes regular and almost periodic occurrence of diagonal line structures in the RPs (Figure 8). Here, we use a recurrence threshold *ε* to ensure a recurrence point density of 0.1. Although we notice a certain amount of similarity between the RPs of neuron 1 and neuron 2 in the nonsynchronized regime, we still see deviations in the line patterns of the RP of neuron 2 (Figure 8A,B). In contrast, the RPs of neuron 1 and neuron 2 for the in-phase synchronized regime show a striking similarity (Figure 8C,D).

**Figure 8.** Recurrence plots of the membrane voltage for weakly coupled neurons as shown in Figure 5 for (**A**,**B**) no synchronization, *g* = 0.01, and (**C**,**D**) phase synchronization, *g* = 0.04. Embedding parameters were estimated using the PECUZAL method [70]; the recurrence threshold is selected to ensure a recurrence rate of *RR* = 0.1.

The vertical distance between these diagonal line structures is related to the phase. Therefore, we can use the density of recurrence points along diagonals parallel to the main diagonal, the *τ*-recurrence rate *RRτ*, Equation (12), as an estimator of the phase distribution and compare it between different systems. For two nonsynchronized systems, the recurrence probabilities should differ significantly (Figure 9A). During phase synchronization, *RR<sup>τ</sup>* should have high probabilities at the same *τ* values; thus, the shape of *RR<sup>τ</sup>* should be very similar (Figure 9B). Therefore, *RR<sup>τ</sup>* has been used to construct a measure for phase synchronization between two signals *x*<sup>1</sup> and *x*<sup>2</sup> by calculating the Pearson correlation of probability of recurrence (CPR) between *RRx*<sup>1</sup> *<sup>τ</sup>* and *RRx*<sup>2</sup> *<sup>τ</sup>* [34]

$$\mathbb{C}PR^{\rm P} = \frac{\text{cov}\left(RR\_{\rm \tau}^{\rm x\_1}, RR\_{\rm \tau}^{\rm x\_2}\right)}{\sigma\_{\rm RR\_{\tau}^{\rm x\_1}}\sigma\_{\rm RR\_{\tau}^{\rm x\_2}}},\tag{18}$$

with *σRR<sup>τ</sup>* the standard deviation of the corresponding *RR<sup>τ</sup>* series. CPR values of 1 would then correspond to phase synchronization and 0 to no synchronization. Here, it is important to remove the first peak in *RR<sup>τ</sup>* close to *τ* = 0 because these values correspond to the main diagonal in the RP present in all systems [50]. Therefore, this first peak would indicate some kind of similarity between *RRτ*(*x*1) and *RRτ*(*x*2) even for completely desynchronized systems. Such exclusion of the first part of the *RR<sup>τ</sup>* series corresponds to applying a Theiler window [75]. Here we used a Theiler window of 25 mS.

**Figure 9.** *τ*-recurrence rate for weakly coupled neurons as shown in Figures 5 and 8 for (**A**) no synchronization, *g* = 0.01, and (**B**) phase synchronization, *g* = 0.04. For phase synchronization, the *τ*-recurrence rate series for both neurons are almost completely overlapping.

Another concern on calculating the CPR measure using Equation (18) is the spiky shape of the *RR<sup>τ</sup>* series, biasing the Pearson correlation estimation. As an alternative, we could use the Spearman rank correlation instead of the Pearson correlation,

$$\text{CPR}^{\text{S}} = \frac{\text{cov}\left(\text{R}(\text{RR}^{\text{x}\_{1}}\_{\text{\text{\text{\textdegree}}}}), \text{R}(\text{RR}^{\text{x}\_{2}}\_{\text{\text{\textdegree}}})\right)}{\sigma\_{\text{\textdegree R}(\text{RR}^{\text{x}\_{1}}\_{\text{\text{\textdegree}}})} \sigma\_{\text{\textdegree R}(\text{RR}^{\text{x}\_{2}}\_{\text{\text{\textdegree}}})} \,\text{},\tag{19}$$

with the *RR<sup>τ</sup>* series converted to the ranks R(*RRτ*). This correlation measure is expected to work better for non-normal distributed data, as the *RR<sup>τ</sup>* series would be.

Both CPR measures clearly show the onset of phase synchronization at *g* = 0.066 for neurons without ELF EF and at *g* = 0.037 for neurons with ELF EF (Figure 10). There are some differences between *CPR*<sup>P</sup> and *CPR*<sup>S</sup> . During phase synchronization, *CPR*<sup>P</sup> is almost 1, but *CPR*<sup>S</sup> is slightly below 1, even more obvious for the coupled neurons without ELF EF. However, for phase synchronization, we would expect to have a CPR value of 1. Moreover, the transition from a nonsynchronized regime to a synchronized regime is not as abrupt as indicated by ∆Ω (Figure 6), but *CPR*<sup>S</sup> changes almost abruptly from very low values to very large values, whereby *CPR*<sup>P</sup> shows a more gradual increase (and even step-wise increase for *A* = 0). This finding indicates that the Spearman-based CPR is obviously not a better choice than the Pearson-based CPR measure.

**Figure 10.** Correlation of probability of recurrence CPR based on Pearson (dotted) and Spearman (line) correlations indicating the onset of phase synchronization between chaotic coupled neurons without external EF (red) and with external EF where *A* = 0.1 and *ω* = 0.286 (blue).

The *RR<sup>τ</sup>* series represents probabilities of recurrence. Therefore, it seems more natural to use a measure that can directly quantify the difference between probability populations, such as Kullback–Leibler distance [76] or Hellinger distance [77]. Here we test the use of the Hellinger distance

$$H(RR\_{\tau}^{\mathbf{x}\_1}, RR\_{\tau}^{\mathbf{x}\_2}) = \frac{1}{\sqrt{2}} \left\| \sqrt{RR\_{\tau}^{\mathbf{x}\_1}} - \sqrt{RR\_{\tau}^{\mathbf{x}\_2}} \right\|,\tag{20}$$

which corresponds to the Euclidean norm of the square root distances between the *RR<sup>τ</sup>* series of the two signals. Values of *H* close to 0 indicate phase transition, whereas values close to 1 indicate nonsynchronized regimes.

To assess whether the variation of *H* indeed reveals phase synchronization, we use a simple block shuffling approach to test the null hypothesis that the signals are not synchronized. Block shuffling splits a time series into a number of blocks (here, we used five blocks) of equal width at random indices and randomly concatenates these blocks to create a new surrogate time series. Such surrogates preserve short-term temporal properties but destroy long-term dynamical information and, thus, correlations when compared with another signal. The distribution *p*(*H*) derived from the ensemble of surrogates is then used to define the confidence limit of 95% (simply by using the 95% quantile of this test distribution *p*(*H*)). Considering *A* = 0, we find the confidence limit by *H*0.95 as 0.17. Values of *H* below this value can be considered to represent phase synchronization.

The measure *H* indicates the transition from the nonsynchronous to the phase synchronization regime of the two weakly coupled neurons (Figure 11). The change in *H* is significant. Moreover, the variation in *H* over increasing coupling strength *g* reveals the more gradual change to the phase synchronization as well as the step-like transition to phase synchronization for the situation without ELF EF caused by phase jumps.

**Figure 11.** Hellinger distance of the *τ*-recurrence rate indicating the onset of phase synchronization between chaotic coupled neurons without external EF (red) and with external EF where *A* = 0.1 and *ω* = 0.286 (blue). A drop of *H* below the confidence limit of 95% (dotted line) represents the significance of this finding.

#### **5. Conclusions**

The synchronization of weakly coupled Morris–Lecar neurons under common external forcing has been studied previously. For example, Kitajima and Kurths [21] used interspike intervals (to study frequency locking) and Yi et al. [62] considered the average firing rate. In general, the numerical calculation of the phases of spiky signals using the Hilbert transform is problematic. An alternative way to identify the phases in dynamical systems is to use recurrence plots [34]. This method can find phases for noncoherent and spiky signals. We, therefore, used a recurrence-based approach, which decodes the phase in terms of specific recurrence patterns in the recurrence plot, and demonstrated its potential for the study of spiking patterns of neurons.

In this work, the spiking patterns of Morris–Lecar neurons under ELF sinusoidal EF and the synchronization of two neurons weakly coupled with gap junction under ELF EF were investigated using this recurrence-based approach. The representation of the dynamics of the neurons' membrane voltages by recurrence plots provided a convenient approach to compare the recurrence features of their spiking patterns. Various spiking patterns, such as periodic and chaotic bursting and periodic spikes, were observed. The spiking patterns were found to be very sensitive to changes of the stimulus frequency.

Moreover, the recurrence approach allows us to consider phase differences between the spiking patterns in a more robust way than the frequently used Hilbert transform. We have introduced an alternative measure for testing phase synchronization using recurrences. Instead of comparing the probabilities of recurrences (as represented by the *τ*-recurrence rate) by correlation coefficients, we suggest to use the Hellinger distance as a more natural measure because it quantifies the differences between probabilities. The typically used Pearson correlation is biased because the *τ*-recurrence rate does not follow a normal distribution. The Spearman rank correlation could be an alternative, but we found additional bias due the large number of zeros in the *τ*-recurrence rate series.

By using recurrence-based synchronization measures, we found that even without external EF, phase synchronization of two ML neurons can occur for a range of values of coupling strength. Moreover, phase synchronization can be enhanced by an additional external EF. This physiological behavior might be of importance for the functioning of the brain when exposed to electromagnetic fields, such as by power lines, electrical equipment, or cellular radio towers.

**Author Contributions:** Conceptualization, A.M.N., J.K. and N.M.; software, A.M.N., E.K.N. and N.M.; methodology, investigation, and writing—original draft preparation, A.M.N. and N.M.; evaluation, writing—review and editing, J.K., B.R.N.N. and E.K.N.; visualization, N.M.; supervision, N.M. and J.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was made possible by the financial support from the partnership between Deutsche Forschungsgemeinschaft (DFG) and the Academy of Science for the Developing World (TWAS) under the grant number 31405812. E.J.N. and J.K. acknowledge the Volkswagen Foundation (Grant No. 85391).

**Data Availability Statement:** Code used to perform the analysis of this study is available via Zenodo, doi:10.5281/zenodo.5910812 (accessed on 29 December 2022).

**Acknowledgments:** We thank Paul Woafo for support and supervision. We further acknowledge the help of K. Hauke Kraemer and Frank Hollmann for support in the implementation of the model in Julia.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

