Hybrid Harmony is an open-source software package written in Python (https://github.com/RhythmsOfRelating/HybridHarmony), accompanied by a visualization module and a sonification module (Figure 2). The software consists of a backend that handles data acquisition and performs analyses, and a Graphical User Interface (GUI) made with PyQt5 (https://www.riverbankcomputing.com/software/pyqt/), where users can control parameters for the analyses (Figure 2). We introduce the software and discuss compatible hardware (EEG systems) in section Hardware, the processing pipeline including preprocessing and connectivity analysis in section Data Preprocessing, Connectivity Analysis, Normalization, and the visualization and sonification modules in section Visualization and Sonification. Data transfer protocols, i.e., LabStreamingLayer (LSL, https://github.com/sccn/labstreaminglayer) and Open Sound Control (OSC, Wright and Momeni, 2011), are described in Supplementary Material, Data Transfer Protocol. Detailed instructions can be found on the GitHub page.

The tool is compatible with any EEG device that interfaces with LabStreamingLayer, and has been tested with MUSE (https://choosemuse.com/; Bhayee et al., 2016), emotiv (EPOC and EPOC+; https://www.emotiv.com/; Williams et al., 2020), the SMARTING system from mBrainTrain (https://mbraintrain.com/; Grennan et al., 2021), and Brain Vision LiveAmp systems (https://www.brainproducts.com/; Fang et al., 2019), and can be expanded to accommodate other systems that export data to LabStreamingLayer.

The first stage of data preprocessing is a buffering procedure that holds and segments incoming data time-series. Incoming streams from LSL are stored in a 30 s buffer updating at the EEG data's sampling rate (e.g., 250 Hz for Brain Vision LiveAmp system). Hybrid Harmony then selects the most recent time window to perform the signal processing procedure. The time window is determined by “window size” and is 3 s by default. The rate of the analysis depends on the computation bandwidth of the system Hybrid Harmony is running on. For example, on many systems we tested data are analyzed roughly 3.5 times per second.

The time window is filtered with the infinite impulse response (IIR) filter (Oppenheim, 1999) into frequency bands of interest (e.g., 8–12 Hz for the alpha frequency band), and then Hilbert transformed to generate the instantaneous analytic signal (Oppenheim, 1999). Users can choose to output spectral power concurrently by selecting the “sending power values” checkbox: power spectral density will be computed and sent to LSL along with connectivity values.

The GUI allows users to change processing parameters via “Frequency bands for analysis,” “Input data streams,” and “window size” (Figure 3). “Frequency bands for analysis” is an editable table specifying the frequency bands of interest, where “Freq. Band” denotes the frequency band name, “Min. Freq” and “Max. Freq” denote the lower and upper bounds of the frequency band, and “weight” determines how the different bands are weighted relative to one another. Figures 3, 4 show the default setup for this table. “Input data streams” displays the name, channel count, and sampling rate of the incoming streams, and its editable cells (e.g., “theta channels,” “alpha channels,” etc.) determine the specific channels to use for each frequency band. Lastly, the “window size” text field determines the length of the data segment for the analysis.

Connectivity Analysis takes the analytic signal from Data Preprocessing as input, and computes one connectivity value for each electrode pair per frequency band in every participant pair. Then, connectivity values are averaged across electrode pairs, so the output is one connectivity value per frequency band for every participant pair. The exact computations are adapted from the python-based inter-brain analysis pipeline HyPyP (Ayrolles et al., 2021), and listed in Table 1. The user may choose from a list of connectivity measures by changing the “Connectivity metric” parameter (Figure 5). Currently implemented metrics include coherence, imaginary coherence, envelope correlation, power correlation, phase-locking value (PLV), and Circular Correlation Coefficient (CCorr). The mathematical equations and references of these metrics are provided in Table 1. In the equations, X(i, t) denotes the analytic signal for subject x at channel i for time point t, and Y(j, t) is that for subject y at channel j for time point t, and the star sign denotes the complex conjugate. The result of each equation is one synchrony value per electrode pair, written as r_{i, j}, and the computation is carried out for all electrode pairs (i, j) where i belongs to subject x 's channels and j belongs to subject y 's channels. Note that we are computing all metrics using the analytic signal from the previous step to streamline the computation. We compute synchrony values from the analytic signals using Hilbert Transform (HT), an alternative analysis to the windowed Fast Fourier Transform (FFT). The analytic signal from HT and the spectra from the windowed FFT both represent the amplitude and phase of the signal in their real and the imaginary parts, respectively, except that the analytic signal is “instantaneous,” while the windowed FFT is an average value over a period (Kovach, 2017). Therefore, while the cross-spectra in coherence is usually the expected value of the product between the two signals' spectra, we used a computation (Table 1) adapted from Equation (1) (Kovach, 2017), i.e., the cross-spectra are expressed as the expectation of the dot product between X(i, t) and the complex conjugate of Y(i, t). Using this formulation is appropriate when investigating the synchronization of signals as it allows us to measure their similarity on a sample-by-sample basis, not just as an average over a relatively long time-window.

Currently, only one metric at a time can be employed, but a user can run multiple instances of the software and thus output multiple metrics simultaneously. The user can then record these streams using LabRecorder (Supplementary Material, Saving Data Through LSL). For visualization, it is possible to differentiate the streams based on the unique source_id in the metadata of the LSL stream, and choose only one to display. However, this feature is not yet developed in our visualization module.

The metrics, which, as mentioned above, are a subset of those implemented in the hyperscanning analysis pipeline HyPyP (Ayrolles et al., 2021), include two variations of correlation (envelope and power correlation) and two of coherence (coherence and imaginary coherence) measures, which are traditional linear methods to estimate brain connectivity, as well as two measures of phase synchrony (PLV and CCorr). While correlation methods are predominantly employed in hyperscanning fMRI studies to characterize joint action and shared attention (Koike et al., 2016), as shown in the case study below, we have validated power correlation as a neurofeedback synchrony signal. Coherence is more commonly used in fNIRS and EEG hyperscanning studies (Liu et al., 2016; Dikker et al., 2017; Miller et al., 2019). Imaginary coherence, i.e., the imaginary part of the coherence, was developed in response to coherence's susceptibility to zero-lagged spurious synchrony, and has been found to reflect personality traits and social closeness in hyperscanning studies (Nolte et al., 2004; Dikker et al., 2021). For phase synchrony, we included PLV, which has been widely used in hyperscanning studies to capture joint action (Dumas et al., 2010), verbal interaction (Perez Repetto et al., 2017), decision-making (Tang et al., 2016) and other tasks. CCorr measures the covariance of phase variance between two data streams and is more robust to coincidental synchrony (Burgess, 2013) compared to PLV, and has been used to investigate touch (Goldstein et al., 2018), learning (Bevilacqua et al., 2019) and language (Perez Repetto et al., 2017) in hyperscanning studies.

In addition to the connectivity metric, users are also able to choose “connectivity type” (Figure 3), which determines the electrode pair combination for connectivity. For each frequency band, the computation is carried out for every possible electrode pair between the participants, and then averaged based on the “connectivity type” parameter. If “connectivity type” is “one-to-one,” only electrode pairs in the matching position are considered (e.g., Fp1 channel of participant A is only paired with Fp1 of participant B and C, etc.); alternatively, if it is set to “all-to-all,” all electrode pairs are considered in v the averaging.

Connectivity Analysis outputs data chunks to LSL as a “Marker” stream under the name “Rvalues.” The size of this data chunk depends on the number of subjects and the number of frequency bands chosen for analysis. For example, if there are 4 subjects and 4 frequency bands, the data chunk will be a vector of length 24 (6 combinations of pairs times 4 frequency bands). Additionally, if the checkbox “sending through OSC” is selected, the same data chunks will simultaneously be transmitted through the OSC protocol with parameters in “OSC IP address” and “OSC port.”

As part of Connectivity Analysis, normalization of connectivity values is implemented with two options: manual normalization (labeled as “Manual” in Figure 3) and baselining with a pre-recorded file (labeled as “from file” in Figure 3). With a Min-Max normalization method, the user can use either of the options, or a mixture of both with a weighting factor adjusted by the slider “Weight.” The minimum and maximum limits are then weighted between the “Manual” and “from file” options.

Beyond inter-brain coupling, interpersonal synchrony has been examined in more depth in other aspects of behavior, including movement (Oullier et al., 2008; van Ulzen et al., 2008; Varlet et al., 2011), language (Pickering and Garrod, 2004) and physiological rhythms such as heart rate and respiration (Müller and Lindenberger, 2011; Noy et al., 2015). As described above, movement and physiological synchrony may be both cause and effect of inter-brain coupling. In future iterations of Hybrid Harmony, we hope to extend the software to incorporate multiple data streams including physiological and movement data. This would allow users to compare the social relevance of various forms of synchrony, and possibly to tease apart interrelationships between (neuro)physiological and behavioral coupling (e.g., Dumas et al., 2010; Mayo and Gordon, 2020; Pan et al., 2020). Similarly, while we prioritized EEG research in Hybrid Harmony, given the increasingly rich fNIRS hyperscanning literature (Liu et al., 2016; Nozawa et al., 2016; Miller et al., 2019) and recent successful work in fNIRS neurofeedback (Gvirts and Perlmutter, 2020; Kohl et al., 2020) we believe an extension of Hybrid Harmony to include metrics suitable for fNIRS data would be a very welcome and fruitful future direction for the software.

EEG systems' susceptibility to movements and environmental noises can greatly compromise data quality and introduce spurious synchrony in our measure. In controlled lab studies, motion artifacts are often carefully removed manually and through data decomposition (e.g., principal component analysis). In the neurofeedback setting, however, such procedures haven't been widely implemented. In our practice, we tried to address the issue empirically with several solutions. For example, we piloted a version of the Mutual Wave Machine where sudden motion-related fluctuations in the data were removed, but this dramatically influenced the experience: participants often react enthusiastically to a sudden increase in light, only to get “punished” for facial expression, which would discourage them from naturally engaging with each other. We considered patching the data with correlations from non-contaminated stretches of data, but this would lead to arbitrary choices. We therefore instead opted for an alternative solution where participants were told explicitly that because extensive head and facial movements can dramatically affect the EEG signal, what they were seeing could also be caused by synchronous noise or synchronous movement. While this option sufficed for the experiential side of the neurofeedback, it is suboptimal with regard to data fidelity. In future releases of the software, we will incorporate support for online data cleaning procedures such as toolboxes (Mullen et al., 2015) and EEG systems that provide built-in data cleaning options in their software, e.g., the SMARTING system by mBrainTrain (Lee et al., 2020).

In addition to challenges related to data cleaning, the real-time nature of the analysis procedure poses challenges in terms of its interpretation. For instance, data fidelity is much higher when applying filtering and correlation analysis on larger stretches of data, but in the type of analysis employed here, this would compromise the immediacy of the neurofeedback. Therefore, it is important to note that our real-time approach might not be able to characterize those types of synchrony that may not be temporally aligned over short intervals. For example, fMRI studies have suggested there may be delays of up to 8 seconds in inter-brain synchrony between speakers and listeners (Stephens et al., 2010a; Dikker et al., 2014; Misaki et al., 2021). Indeed, while successful joint action is typically associated with the coupling of motor movements (Dumas et al., 2010), being “in sync” or “on the same wavelength” is often taken to imply interactive alignment at the level of mental representations (Garrod and Pickering, 2009; Pickering and Garrod, 2013), usually involving more “abstract” constructs such as sharing viewpoints (Van Berkum et al., 2009). These mental representations may or may not be linked to convergence at the temporal level. In line with this dissociation, in two instances of the Mutual Wave Machine we asked participants to reflect on their “connection strategies.” Pairs who used either eye contact or joint action as a connection strategy (mimicry, laughter, motion coordination) exhibited an increase in inter-brain synchrony over time as measured by Imaginary Coherence and Projected Power Correlation (Dikker et al., 2021). Such an increase in synchrony was not observed for pairs who tried “thinking about the same thing.” While these results validated our approach in capturing synchrony in joint action, they do not exclude the possibility of synchrony on the abstract mental representations, given that such synchrony may entail more complex temporal dynamics.

In future iterations, we hope to incorporate alternatives to the time-aligned approach, such as introducing a temporal delay between the data streams and real-time adaptive normalization. These additions will hopefully increase sensitivity to endogenous synchronizers and facilitate more controlled experimental designs (Stephens et al., 2010b).

Finally, it is worth reiterating that very little is known about the correspondence between different synchrony analysis metrics and socio-psychologically relevant factors. Although some metrics are better than others in theory (e.g., CCorr is more robust to spurious synchrony than PLV), and some are more common in hyperscanning studies than single brain studies (e.g., phase synchrony is more common than coherence/correlation), the exact pros and cons of each metric require further investigation. As more becomes known about the mapping between psychological processes as inter-brain synchrony metrics, we will add guidance for users with respect to the choice of metrics in different situations.