New Framework for Understanding Cross-Brain Coherence in Functional Near-Infrared Spectroscopy (fNIRS) Hyperscanning Studies

Summary

Wavelet transform coherence (WTC) is a common methodology for assessing the coupling between signals that is used in functional near-infrared spectroscopy (fNIRS) hyperscanning studies. A toolbox for assessing the directionality of the signal interaction is presented in this work.

Abstract

Despite the growing body of functional near-infrared spectroscopy (fNIRS) hyperscanning studies, the assessment of coupling between two neural signals using wavelet transform coherence (WTC) seems to ignore the directionality of the interaction. The field is currently lacking a framework that allows researchers to determine whether a high coherence value obtained using a WTC function reflects in-phase synchronization (i.e., neural activation is seen in both members of the dyad at the same time), lagged synchronization (i.e., neural activation is seen in one member of the dyad prior to the other member), or anti-phase synchronization (i.e., neural activation is increased in one member of the dyad and decreased in the other). To address this need, a complementary and more sensitive approach for analyzing the phase coherence of two neural signals is proposed in this work. The toolbox allows investigators to estimate the coupling directionality by classifying the phase angle values obtained using traditional WTC into in-phase synchronization, lagged synchronization, and anti-phase synchronization. The toolbox also allows researchers to assess how the dynamics of interactions develop and change throughout the task. Using this novel WTC approach and the toolbox will advance our understanding of complex social interactions through their uses in fNIRS hyperscanning studies.

Introduction

In recent years, there has been a shift in the types of studies conducted to understand the neural bases of social behavior^1,2. Traditionally, studies in social neuroscience have focused on neural activation in one isolated brain during a socially relevant task. However, advances in neuroimaging technology now allow for the examination of neural activation in the brains of one or more individuals during social interaction as it occurs in "real-life" settings³. In "real-life" settings, individuals are able to move freely, and patterns of brain activation are likely to change as information is exchanged and social partners receive feedback from one another⁴.

Hyperscanning is a method that assesses this bidirectional information exchange by measuring the brain activity from two or more individuals simultaneously⁵. An emerging body of research has utilized functional near-infrared spectroscopy (fNIRS), a non-invasive neuroimaging technique that, in comparison to other neuroimaging techniques, is less susceptible to motion artifacts⁶. Hyperscanning via fNIRS allows for the assessment of inter-brain synchronization (IBS) in real-life settings while the interactive partners move freely and naturally. This is particularly relevant for work with infants and young children, who tend to be quite active. IBS has been reported to reflect mutual understanding between interactive partners, which serves as the foundation for effective social interaction and communication and mediates shared intentionality^1,7,8.

Several methods are used to evaluate the IBS of two brains. Such methods include time series correlations, such as cross-correlation and the Pearson correlation coefficient^9,10 (see a review by Scholkmann et al.¹⁰). Other methods involve evaluating the strength of coupling in the frequency domain. Such methods include phase-locking value (PLV) and phase coherence (see a review by Czeszumski et al.¹¹). One of the most common methods in fNIRS studies uses wavelet transform coherence (WTC)-a measure of the cross-correlation of two time series as a function of frequency and time¹⁰.

WTC uses correlational analyses to calculate the coherence and phase lag between two time series in the time-frequency domain. FNIRS hyperscanning studies have used WTC to estimate IBS in many domains of functioning, including action monitoring¹², cooperative and competitive behaviour^5,13,14,15, imitation¹⁶, mother-infant problem solving¹⁷, and teaching-learning behavior^18,19,20,21. Typically, in hyperscanning studies, cross-brain coherence, as measured by WTC, during an experimental task is compared to cross-brain coherence during a control task. These findings are usually presented with a WTC "hot plot", which shows the coherence across the two brains at each time point and frequency (see Figure 1).

As suggested by Czesumaski et al.¹¹, WTC has become the standard analytical approach for analyzing fNIRS hyperscanning. WTC analysis is a flexible, "tool-agnostic" method for data visualization and interpretation²². The coherence coefficient heatmap, which provides a narrative form of analysis allowing for the easy identification of periods of synchronous or asynchronous behavior as well as the intensity of brain activity during the completion of a task, is the main advantage of WTC and makes it a strong tool for applied research²². WTC has an advantage over correlation techniques. Correlations are sensitive to the hemodynamic response function's (HRF) shape, which is thought to differ between individuals (particularly in terms of age) and among different brain areas. In contrast, WTC is unaffected by interregional changes in the (HRF)²³. Researchers have used the wavelet approach to study fMRI time series. Zhang et al.²⁴ compared the commonly used functional connectivity metrics, including the Pearson correlation, partial correlation, mutual information, and wavelet coherence transformation (WTC). They carried out classification experiments using large-scale functional connectivity patterns derived from resting-state fMRI data and natural-stimulus fMRI data of video viewing. Their findings indicated that WTC performed best in classification (specificity, sensitivity, and accuracy), implying that WTC is a preferable functional connectivity metric for studying functional brain networks, at least in classification applications²⁴.

Figure 1: Wavelet transform coherence (WTC). WTC shows the coherence and phase angle between two time series as a function of both time (x-axis) and frequency (y-axis). The coherence increase is depicted by the red color in the graph, and the small arrows in the graph show the phase angle of the two time series. The rightward-pointing arrow represents in-phase synchronization; the downward-pointing and upward-pointing arrows represent lagged synchronization; and the leftward-pointing arrow represents anti-phase synchronization³⁰. This figure was adapted from Pan et al.¹⁹. Please click here to view a larger version of this figure.

Recently, Hamilton²⁵ articulated several limitations to the interpretation of cross-brain coherence data in fNIRS hyperscanning studies. One of Hamilton's primary concerns was that coherence measures (e.g., WTC) only report effects as symmetrical (i.e., two brains are correlated, showing the same pattern of change). However, many social interactions are asymmetric (e.g., information flow between a speaker and a listener) in that two participants may play different roles, and it is not clear that WTC can capture this information. Here, this concern is addressed by a new framework that allows for a straightforward interpretation of the cross-wavelet power by using the cross-wavelet phase to detect directionality. This framework will also allow examination of how the dynamics of interactions develop and change throughout a task.

While WTC and correlation methods assess functional connectivity, other methods assess effective connectivity, attempting to extract the causal influences of one neural element over another. Transfer entropy is a measure from the field of information theory that describes the transfer between jointly dependent processes²⁶. Another related method is Granger causality analysis (GCA), which has been described as equivalent to transfer entropy²⁶.

In the existing literature of fNIRS hyperscanning studies, Granger causality analysis (GCA) has been widely used to estimate the coupling directionality between fNIRS time series data obtained during a variety of different tasks, such as cooperation⁵, teaching¹⁹, and imitation¹⁶. GCA employs vector autoregressive models to assess the directionality of coupling between time series in brain data. Granger causality is based on prediction and precedence: "a variable X is said to 'G-cause' variable Y if the past of X contains the information that helps to predict the future of Y over and above information already in the past of Y"²⁷. Accordingly, the G-causality is analyzed in two directions: 1) from subject A to subject B and 2) from subject B to subject A.

While GCA analysis serves as a complementary analysis aimed at determining whether a high coherence value obtained using a WTC function reflects IBS or lagged synchronization (one signal leading the other), it does not allow for the determination of whether anti-phase synchronization has occurred. In traditional neuroimaging studies, in which only one participant is scanned (i.e., the "single-brain" approach), an anti-phase pattern means that the activity in one brain region is increased while the activity in the other brain region is decreased²⁸. In the hyperscanning literature, the presence of anti-phase synchronization may suggest that neural activation is increased in one subject, and at the same time, neural activation is decreased for the other subject. Therefore, there is a need to provide a comprehensive model that can detect the directionality. More specifically, this model will be able to detect anti-phase synchronization (in which the direction of activity in one individual is opposite to that of their partner) in addition to in-phase synchronization and lagged synchronization.

In an attempt to address the concern that WTC shows only symmetrical effects, where both brains show the same pattern of change²⁵, a new approach to identify the type of interaction by examining the phase of synchronization (i.e., in-phase, lagged, or anti-phase) is presented (see Figure 2). To this end, a toolbox using the WTC method to classify the different types of interactions was developed. The types of interactions are classified by using relative phase data from cross-wavelet transform analysis.

Figure 2: Illustration of the different phase relationships of simple sine waves. (A) When the two signals, Signal 1 (blue lines) and Signal 2 (orange lines), reach their respective maximum, minimum, and zero values at the same time point, they are said to be showing in-phase synchronization³². (B) When one signal reaches its maximum value and the other signal reaches zero value at the same time point, they are said to be showing lagged synchronization (one is leading by 90°)^32,33,34. (C) When two time series shift in opposite directions, meaning one signal reaches the maximum and the other reaches the minimum value at the same time point, this is referred to as anti-phase synchronization²⁸. (D-P) In all other phase relationships between two time series, one signal is leading the other. In all positive phases, Signal 2 is leading Signal 1 (e.g., panels E, F, M, and N), whereas in all negative phases, Signal 1 is leading Signal 2 (e.g., panels D, G, H, O, and P). Notably, when the absolute value of the phase is higher, it becomes more distinctive which time series is leading the other (e.g., the leadership is more distinctive in panel J than in panel I, and in panel K, the leadership is more distinctive than in panel L). Please click here to view a larger version of this figure.

Protocol

The study was conducted at Florida Atlantic University (FAU) and was approved by the FAU Institutional Review Board (IRB).

1. Using Homer3 software (Table of Materials) to perform the pre-processing of the fNIRS hyperscanning data

NOTE: Homer3 is a MATLAB application that analyzes fNIRS data to obtain estimates and maps of brain activation²⁹. Homer3 can be downloaded and installed from the following link (https://openfnirs.org/software/homer/).

1. Open MATLAB, and navigate to the folder where the raw .nirs files are saved. Select and open the folder.
2. Type Homer3 into the command window of MATLAB to launch the Homer3 GUI. Homer3. will detect the .nirs files and ask to convert to .snirf format (a universal file format for storing and sharing NIRS data independently of any specific application-specific file format such as MATLAB) to proceed with the pre-processing of the data.
3. After importing the .nirs files into .snirf format in Homer3, click on the Tools option in the Homer3 GUI, and select Edit Processing Stream.
4. In the ProcStreamEdit GUI, select the pre-processing steps from the Registry Function column to the Current Processing Stream column by clicking on Add. The included pre-processing steps are as follows:
   1. Use hmrR_intensity2OD to convert the intensity data to optical density.
   2. Use hmrR_MotionCorrectWavelet to correct motion artifacts using the appropriate filtering function.
   3. Use hmrR_OD2conc to convert the OD data to concentration.
   4. Use hmR_BlockAvg to calculate the Block Average on concentration Data.
      NOTE: The selection of the pre-processing steps may vary depending on the type of dataset.
5. To save the current processing stream, click on the Save option, and then exit the ProcStreamEdit GUI.
6. To run the pre-processing stream in the main Homer3 GUI, click on the RUN option. After Homer3 finishes running the selected processing stream, it will save the pre-processed time series for each participant in a .mat file format containing Hbo, Hbr, and Hbt for all the channels and events. A folder named homer output will be created by Homer3 in the currently selected folder to store these files.
7. A folder named derivatives will be created by Homer3 in the selected folder to store these files. Select the homer folder located in derivatives folder. Choose the .mat file for each brain and export Hbo, Hbr, Hbt.
   NOTE: The name of the output folder created by Homer3 is dependent on the Homer3 version.

2. Getting started with the LeaderFollowerByPhase toolbox

1. To analyze the type of interaction that occurs in a hyperscanning recording, use the LeaderFollowerByPhase toolbox, as described in the process shown in Figure 3. In MATLAB, select the .mat files for each brain, and load the Hbo (or Hbr) data of the specific channel and specific event into a one-dimension vector as signal1 and signal2.
2. In the MATLAB command line, define the parameters
   1. lowFreq, highFreq: Type lowFreq = [low FOI], and highFreq = [high FOI]. The default values are lowFreq = 0.01 Hz, highFreq = 1 Hz.
      NOTE: The parameters of the functions lowFreq and highFreq define the frequency of interest (FOI) range. WTC calculates the coherence across the two brains at each time point and frequency. The coherence values are typically averaged within a particular FOI.
   2. Define the parameter phaseRange; type phaseRange = [range in deg].
      NOTE: The default value is phaseRange = 90°. The phase ranges between 0° to 360° due to the circular modulo nature of the phase. The phase ranges are divided according to a range that surrounds four points. In the presented toolbox, a new approach for classifying asymmetric interactions is presented (Figure 4) by examining the coupling directionality using the phase angle values according to the ranges corresponding to lagged synchronization with Signal 1 leading (a range surrounding −90°) or Signal 2 leading (a range surrounding 90°), Signal 1, Signal 2 in-phase synchronization (a range surrounding 0), and Signal 1, Signal 2 anti-phase synchronization (a range surrounding +180° or −180°).
   3. Define the Threshold parameter. Type threshold = [threshold rsq val]. The default value is Threshold = 0.
      NOTE: The toolbox allows for the specification of a threshold coherence value by specifying the threshold parameter. This allows the researcher to select time points with a specified minimum coherence value. Consequently, only time points with coherence values that are higher than the specified threshold are taken into consideration.
3. Download the LeaderFollowerByPhase toolbox from the following link (https://www.ariel.ac.il/wp/sns/download/ or https://github.com/Minisharmaa/Leader-Follower-By-Phase).
4. Execute the MATLAB function LeaderFollowerByPhase by entering the command cohervalues = LeaderFollowerByPhase(signal1, signal2, lowFreq, highFreq, phaseRange, threshold) into the command line.
   NOTE: The coherence and phase calculations are performed using the WTC and XWT functions in MATLAB, respectively³⁰.
5. Examine the in-phase, Signal 1 leading, Signal 2 leading, and anti-phase synchronization values:
   1. Inspect the plots in MATLAB. The toolbox generates one figure with four plots.
      1. Coherence by type of interaction: Inspect the box chart plot on the top-left part of the figure, which shows the R-Squared (Rsq) according to each type of interaction (in-phase, Signal 1 leading, Signal 2 leading, anti-phase).
         NOTE: For a detailed description of the MATLAB box chart built-in function, see the following link (https://www.mathworks.com/help/matlab/ref/boxchart.html).
      2. Central indices by type of interaction: Inspect the bar graph on the top-right part of the output figure, which displays the max mean and median according to each type of interaction (in-phase, Signal 1 leading, Signal 2 leading, anti-phase).
         NOTE: For a detailed description of the MATLAB bar built-in function, see the following link (https://www.mathworks.com/help/matlab/ref/bar.html).
      3. Coherence over time: Inspect the scatterplot at the bottom-left part of the output figure, which displays the values of coherence and the types of interaction over time. The colored dots represent different types of interaction (black dots represent in-phase synchronization, dark grey dots represent Signal 1 leading, light grey dots represent Signal 2 leading, and purple dots represent anti-phase synchronization).
         NOTE: The figure shows the dynamics of the interaction: the exchange between the four types of interaction throughout the full-time series. For a detailed description of the MATLAB scatter function, see the following link (https://www.mathworks.com/help/matlab/ref/scatter.html).
      4. Time percentage: Inspect the pie chart at the bottom right part of the output figure, which displays the division of time according to the different types of interactions.
         NOTE: For a detailed description of the MATLAB pie function, see the following link (https://www.mathworks.com/help/matlab/ref/pie.html).
6. Inspect the output table with the statistic values (i.e., mean, max, median, and standard deviation) for each type of interaction (in-phase synchronization, Signal 1 leading, Signal 2 leading, anti-phase synchronization). The table also presents the percentage of time for which each type of interaction occurred. Each type of interaction appears in a different column.
7. Examine the output value in the extracted spreadsheet file (i.e., datatable.xlsx located in the current folder).

Figure 3: Overview of the workflow. (A) Mother-child dyads engaged in free play while hyperscanning fNIRSdata were collected. (B) Illustration of a mother-child time series.(C) Pre-processing of the time series using Homer3.(D,E) The use of a toolbox to examine different types of interactions, such as in-phase synchronization, anti-phase synchronization, and lagged synchronization. Please click here to view a larger version of this figure.

Figure 4: Classification of four different types of interactions based on phase. This represents the phase difference of the two neural time series in a 360° modulo. The phase difference can be thought of as a time lag between two values and is measured in degrees and radians or fractions of the wavelength. Here, the 360° modulo is divided into four different ranges depicting four different phases of interaction: (A) Signal 1 leading (a range around 90°, between 45° to 135°), (B) anti-phase synchronization between Signal 1 and Signal 2 (a range around 180° or −180°, between 135 to −135°), (C) Signal 2 leading (between −135° to −45°), (D) in-phase synchronization (a range around 0, between −45° to 45°). This division is the default approach (45° around each point); however, the toolbox allows one to configure a different division. While other configurations may not cover all 360°, it may yield a more exact definition of each type of interaction. Please click here to view a larger version of this figure.

Representative Results

This section demonstrates the types of analyses that can be carried out with the toolbox (which can be downloaded at https://www.ariel.ac.il/wp/sns/download/ or https://github.com/Minisharmaa/Leader-Follower-By-Phase.). For these analyses, fNIRS data collected with a small sample of infant-parent dyads were utilized. Six pairs of mother-infant dyads were tested using a validated behavioral task, the free-play task³¹, which is as close to a real-life infant-mother interaction as possible. Prior to the experiment, the infants and parents were fitted with a custom-made optode set for the collection of fNIRS data. The optode set used to collect the fNIRS data in this study comprised 8 sources (red dots) and 8 detectors (blue dots) that were configured to create 18 channels covering the prefrontal and temporoparietal regions bilaterally (see Figure 5). An NIRScout acquired the optical imaging data using two wavelengths: 760 nm, which is more sensitive to deoxyhemoglobin (HbR), and 850 nm, which is more sensitive to oxyhemoglobin (HbO). All the parents were females (age range = 26-36 years), and the infants were healthy and full-term (two females, four males, age range = 1-2 years) with no known developmental delays. Dyads were recruited through advertisements. Each parent gave informed consent prior to the experiment, and they were paid for their participation. For simplicity, the analysis focuses on data obtained at channel 18 of Dyad A.

Figure 5: Optode set used in the preliminary study. The optode set used to collect the fNIRS data in the preliminary study consisted of 8 sources (red dots) and 8 detectors (blue dots) that were configured to create 18 channels (yellow lines) covering the prefrontal and temporoparietal regions bilaterally. Please click here to view a larger version of this figure.

The toolbox was used to identify changes in the types of interaction that could be identified over time for that channel for a specific dyad. The parameters for the function were as follows: Signal 1 = channel 18 mother, Signal 2 = channel 18 mother, lowFreq = 0.0067, highFreq = 0.1142, phaseRange = 90, Threshold = 0.

Figure 6: Classification analysis at a threshold of 0. The threshold is set to 0 (Threshold = 0). (A) Box plots representing coherence value associated with interactions. There is one for each type of interaction, and the median and interquartile range (IQR) are indicated. Higher scores indicate a greater level of coherence. (B) The central indices of the coherence values of all the types of interactions. (C) The dynamics of the type of interaction change throughout the task. (D) The percentage of scores for each of the four types of interactions. Please click here to view a larger version of this figure.

First, each time point was classified into one of the four types of interaction (in-phase, anti-phase, mother-leading, or infant-leading) (Figure 6). Figure 6D shows the percentage of time points that were classified into each of the four types of interactions. It is important to note that while the toolbox can detect the percentages and the coherence value associated with interactions that are led by one of the participants (Figure 6A), its unique contribution is that it also presents the percentages (Figure 6D) and central indices of the coherence values of all the types of interactions, including anti-phase synchronization (Figure 6B). Finally, the toolbox allows one to examine how the dynamics of the type of interaction change throughout the task (Figure 6C). It is important to note that, similar to GCA analysis, the toolbox calculates these indices for each dyad separately. Group-level analysis using these indices should be conducted to determine the type of interaction.

Here, to explore the influence of changing these minimum threshold values on the classification of the types of interaction within a dyad, the classification analysis was repeated with a threshold of 0.5 on Dyad A (Figure 7).

Figure 7: Classification analysis at a threshold of 0.5. The threshold is set to 0.5 (Threshold = 0.5). (A) Box plots representing coherence value associated with interactions. There is one for each type of interaction, and the median and interquartile range (IQR) are indicated. Higher scores indicate a greater level of coherence. (B) The central indices of the coherence values of all the types of interactions. (C) The dynamics of the type of interaction change throughout the task. (D) The percentage of scores for each of the four types of interactions. Please click here to view a larger version of this figure.

As shown in Figure 7D, when using this threshold, the distribution of the different types of relative phase relationships changed. The percentage of anti-phase synchronization increased (from 35% to 59%), and the percentage of in-phase synchronization decreased (from 26% to 3%). This suggests that anti-phase synchronization may be the type of interaction that is more representative of this dyad. In other words, defining a threshold allows for the ability to conduct a more sensitive analysis in which only time points with a minimum level of coherence are averaged. It is important to note that determining the optimal coherence value threshold is a complicated process, as the optimal threshold may vary from one experiment to another and across different environments. Although the toolbox provides the possibility to set a threshold, more studies are needed to develop a protocol for identifying the optimal coherence value. Moreover, it is important to select threshold and frequency of interest values that still intersect with the Rsq values. For example, the function with the parameters lowFreq = 0.0067, highFreq = 0.1142, phaseRange = 90, and Threshold = 0.5 showed interactions with Rsq values only above 0.5, but the same function with a threshold of 0.7 resulted in an error, as there were no values above 0.7 within the frequency range.

Supplementary File 1: Wavelet transform coherence (WTC). The overview of wavelet transform and cross-wavelet transform, which are utilized to analyze the time-frequency characteristics and interdependence of two time series. The wavelet transform decomposes a time series into time-frequency space³⁵, whereas the cross-wavelet transform reveals the common power and phase between two time series^9,30. The text also introduces the coherence of the wavelet transform, which quantifies the degree of synchronization between two time series. The R-square value derived from the coherence of the wavelet transform reflects interdependence but does not distinguish between positive and negative correlations³⁶. Positive and negative correlations are assumed to indicate interrelationships^37,38. Please click here to download this File.

Discussion

One of the most common methods used in fNIRS studies is wavelet transform coherence (WTC), which is a measure of the cross-correlation of two time series as a function of frequency and time¹⁰. WTC calculates the coherence and phase lag between two time series using correlational analyses (Supplementary File 1). FNIRS hyperscanning studies have used WTC to estimate IBS in many domains of functioning, including action monitoring¹², cooperative and competitive behaviour^5,13,14,15, imitation¹⁶, mother-infant problem solving¹⁷, and teaching-learning behavior^18,19,20,21. Hyperscanning studies often compare cross-brain coherence measured using the wavelet coherence transform (WTC) during an experimental task to that of a control task. These comparisons are typically presented using a WTC "hot plot" that displays the coherence between the two brains at each time point and frequency. In addition, as seen in Figure 1, the phase lag information is shown by the direction of small arrows in the WTC "hot plot". However, previous studies have neglected to take into account the phase lag information represented by the direction of the small arrows in the WTC "hot plot" and have only estimated the inter-brain synchrony (IBS) by examining the coherence in the WTC plot. This oversight may result in inaccurate or incomplete findings.

The limitations discussed by Hamilton²⁵ regarding the interpretation of cross-brain coherence data in fNIRS hyperscanning studies are addressed in the new framework allowing a straightforward interpretation of the cross-wavelet power by using the cross-wavelet phase to detect directionality and also includes a coherence analysis module for calculating the coherence values by averaging them directly³⁹. This approach enables examining the development and change in interactions throughout a task and provides a reliable measure of coherence between signals.

Such an approach has been demonstrated in behavioral studies of interpersonal synchronization, which have used the relative phase data that can be extracted from the cross-wavelet analysis. Some studies have used these data to distinguish between in-phase and anti-phase coherence values. For example, this approach has been used to evaluate the hand movements of two improvising musicians⁴⁰ and to examine social postural coordination⁴¹. Some studies have examined the distribution of phase angles in movement data to understand the dynamics of interactions by using cross-wavelet coherence during structured⁴² and unstructured⁴³ conversations.

The relative phase between two time series allows for the detection of temporal shifts between signals of the same frequency. Indeed, in the field of EEG hyperscanning, most of the methods aiming to determine the degree of synchronization of neural time series assess the relative phase relationship between the two time series^13,44.

The critical steps of using the LeaderFollowerByPhase toolbox in fNIRS hyperscanning data are demonstrated in the protocol. Specifically, the protocol involves the pre-determination of Signal 1 and Signal 2 in MATLAB before running the toolbox. It is noteworthy that the parameters such as the frequency of interest (FOI), phase range, and threshold are optional and may use default values if left unset. The filtering and detrending of raw signals are recommended⁴⁵. Additionally, caution must be exercised when conducting band-pass filtering, as this can affect the selection of the FOI.

The FOI parameters (lowFreq, highFreq) require careful selection, specifically excluding high-frequency and low-frequency physiological noise, such as respiration (~0.2-0.3 Hz) and cardiac pulsation (0.6-1.2 Hz). It is recommended to take the low and high frequencies of interest as between 0.01 and 0.7 Hz, respectively⁴⁶, as this range effectively eliminates high-frequency noise such as heart beats (0.8-1 Hz) as well.

The phaseRange parameter defines a range around the phase angle values according to the ranges corresponding to lagged synchronization with Signal 1 leading (a range surrounding −90°) or Signal 2 leading (a range surrounding 90°), Signal 1, Signal 2 in-phase synchronization (a range surrounding 0°), and Signal 1, Signal 2 anti-phase synchronization (a range surrounding +180° or −180°). The width of the surrounding range around these four points is defined by phaseRange For example, if phaseRange is set to 90°, then the range for in-phase synchronization will be surrounding 0°, between −45° to 45°; the range for Signal 2 leading (lagged synchronization) will be surrounding 90°, between 45° to 135°, the range for anti-phase synchronization will be surrounding 180° or −180°, between 135° to −135°; and the range for Signal 1 leading (lagged synchronization) will be surrounding 180°, between −135° to −45°. The phaseRange parameter must be between 0° to 90° degrees, as otherwise, the following message will be displayed: "The value of phaseRange variable must be between 0 to 90". Although the range can be any number from 0° to 90°, the minimum recommended value is 30° (±15°). The Threshold value needs to be any value between 0 to 1, as otherwise, the following message will be displayed: "The value of the Threshold variable must be between 0 to 1". It is recommended to choose a threshold that is between 0.25 to 0.75.

Whilst the LeaderFollowerByPhase toolbox presents a promising approach, it is not without limitations. As mentioned above, determining the optimal coherence value threshold is a complicated process, as the optimal threshold may vary from one experiment to another and between different tasks. Testing this toolbox on more diverse datasets is necessary to obtain more accurate information about the optimal values for the threshold.

The ability to understand complex human interactions using fNIRS hyperscanning has been limited by the fact that the current approaches used to detect coupling between two neural signals ignore the directionality of the signals. Here, a more sensitive approach for analyzing the coherence of two neural signals using wavelet transform coherence (WTC) is proposed. The toolbox allows researchers to examine the coupling directionality by classifying the phase angle values as representing in-phase synchronization, lagged synchronization, and anti-phase synchronization.

This novel approach using the toolbox will provide more detailed information about the nature of dyadic interactions, which, to date, has been lacking. For example, whereas phase synchronization and anti-phase synchronization have been treated as identical (Supplementary File 1)³⁶, investigators will now be able to identify the extent to which the neural signals of the dyad members move in the same direction (both increase or both decrease) or opposite directions (one increases and the other decreases). This will have a transformative impact on the understanding of how the brain mediates social processes and behavior.

The proposed framework has promising potential for future applications in the field of interpersonal neural synchronization research, as it allows for the classification of distinct types of interactions, including in-phase synchronization, lagged synchronization, and anti-phase synchronization. By reanalyzing the previous findings with the new proposed framework, researchers can gain a more comprehensive understanding of the nature of the synchronization between the participants. Specifically, the ability to differentiate between in-phase and anti-phase interactions provides a new level of clarity that was previously unavailable, which could lead to more precise interpretations of the previous findings. This functionality of the framework could be applied to a wide range of scenarios, including exploring the role of interpersonal neural synchronization in social behavior, communication, and decision-making processes. Overall, the proposed framework represents a valuable contribution to the field and holds significant potential for future applications.

Materials

Name	Company	Catalog Number	Comments
NIRScout	NIRx Medical Technologies, LLC	n.a.	8 sources, 8 detectors
MATLAB	The Mathworks, Inc.	Matlab 2022a	In this protocol, several toolboxes and buit in MATLAB functions were used: HOMER3 toolbox was used to convert Intensity to OD, to remove motion artifacts through its function hmrMotionCorrectWavelet with default parameters and to convert OD to Conc. Wavelet Toolbox was used to compute WTC.