A Guide to Parent-Child fNIRS Hyperscanning Data Processing and Analysis

The use of functional near-infrared spectroscopy (fNIRS) hyperscanning during naturalistic interactions in parent–child dyads has substantially advanced our understanding of the neurobiological underpinnings of human social interaction. However, despite the rise of developmental hyperscanning studies over the last years, analysis procedures have not yet been standardized and are often individually developed by each research team. This article offers a guide on parent–child fNIRS hyperscanning data analysis in MATLAB and R. We provide an example dataset of 20 dyads assessed during a cooperative versus individual problem-solving task, with brain signal acquired using 16 channels located over bilateral frontal and temporo-parietal areas. We use MATLAB toolboxes Homer2 and SPM for fNIRS to preprocess the acquired brain signal data and suggest a standardized procedure. Next, we calculate interpersonal neural synchrony between dyads using Wavelet Transform Coherence (WTC) and illustrate how to run a random pair analysis to control for spurious correlations in the signal. We then use RStudio to estimate Generalized Linear Mixed Models (GLMM) to account for the bounded distribution of coherence values for interpersonal neural synchrony analyses. With this guide, we hope to offer advice for future parent–child fNIRS hyperscanning investigations and to enhance replicability within the field.


Introduction
Children begin to understand themselves and their social surroundings by engaging in embodied social interactions with others. As this process is interpersonal by definition, it should ideally be studied by simultaneously obtaining data from all interaction partners. However, such a second-person social neuroscience approach has only gained momentum in developmental research in recent years [1]. Methodological challenges made it difficult to obtain neural measures from two (or more) individuals at the same time during naturalistic, live interactions [2]. The emergence of functional near-infrared spectroscopy (fNIRS) and particularly the synchronized measurement of brain activity from parents and their children while they interact with each other-so-called hyperscanning-have recently allowed us to take this step. We are now able to investigate the neurobiological underpinnings of socio-cognitive and affective processes underlying these early social interactions and, thereby, deepening our understanding of child development from a second-person social neuroscience perspective.
Children's social brain development and the associated question of how children come to understand others have long been studied in experiments utilizing pre-recorded and live social stimuli. The obtained results provided important insights into children's abilities to decipher social information but were devoid of the information's social context and especially the dynamic and reciprocal nature of real social exchanges. Traditional developmental neuroscience has, thus, investigated the neural correlates of social perception robust and cross-correlations) [21][22][23]. Here, we will primarily focus on WTC to assess the relation between the two time-series of interacting partners' brain activity. WTC is suggested to be more suitable in comparison to correlational approaches, as it is invariant to interregional differences in the hemodynamic response function (HRF) [24]. Correlations, on the other hand, are sensitive to the shape of the HRF, which is assumed to be different between individuals (especially regarding age) as well as different within distinct brain areas. Moreover, a high correlation may be observed among regions that have no blood flow fluctuations.
In the remaining parts of this paper, we outline the data analysis procedure we have developed to study fNIRS hyperscanning data from naturalistic parent-child interactions. We provide insight into our approach by providing an exemplary dataset of 20 motherchild dyads during a naturalistic, cooperative versus individual problem-solving task to showcase all pre-processing and analysis steps. Statistical analyses will allow preliminary inferences on whether parent-child dyads show significant interpersonal neural synchrony, and if so within which region of interest (ROI) and during which experimental condition. We will be using MATLAB and RStudio to pre-process and analyze the data. Interested readers can follow the steps of the described analyses by examining the associated files on OSF (https://osf.io/wspz4/, accessed on 12 June 2021).

Sample Description
Our exemplary dataset consists of 20 randomly sampled mother-preschool child dyads from a real study data set and is only used for illustration purposes here. The children were 5-6 years of age (M = 5 years; 3 months; SD = 1.5 months). Mothers' age averaged at 37.20 years (SD = 3.51 years). Families were recruited from a database of volunteers based in and around a mid-sized city in eastern Germany. All dyads were of European white origin and came from middle to upper-class families based on parental education and family monthly income. Sixty percent of mothers had a university degree and 60% of families had a monthly income higher than 3000 €. Participants were remunerated for their participation.

Experimental Procedure
The primary interest of the study was whether interpersonal neural synchrony between mother and child increases during a cooperative problem-solving task in comparison to an individual problem-solving task. In the cooperative problem-solving task, mother and child sat face-to-face and were instructed to take turns forming seven geometric shapes into predetermined templates (e.g., rocket, bridge, lamp etc.) using wooden blocks (i.e., Tangrams). In the individual condition, mother and child were given four of the same templates and asked to reconstruct these by themselves. Because they were seated at the same table, a portable wooden barrier was put in between them to induce a separated task context. Mother and child were also instructed to refrain from talking to one another if possible. Each task lasted 120 s and was repeated twice, resulting in four task phases (see Figure 1A). In between those four task phases, three 80-s resting phases were included. fNIRS was simultaneously measured in both participants in bilateral temporo-parietal junction (left hemisphere: Channel 9-12; right hemisphere: Channel 13-16) and dorsolateral prefrontal cortex (left hemisphere: Channel 1-4; right hemisphere: Channel 5-8; see Figure 1B). For a full description of the study design, please refer to [25,26].

General Information on fNIRS Data Acquisition and Analysis
A visual summary of all processing steps is depicted in Figure 2. fNIRS data were obtained with a NIRScout 16-16 system (NIRx GmbH, Berlin, Germany) offering 16 sources and 16 detectors, which were divided to measure two interactants. In other settings, two or more devices can be synchronized to enable hyperscanning. The absorption of near-infrared light was measured at the wavelengths of 760 and 850 mm and the sampling frequency was 7.81 Hz. The start of each condition was indicated via triggers sent through an experimental program, for instance, OpenSesame [27]. Varying task conditions can be indicated by triggering a second pin for the end of the condition, while preset and standardized task duration can be manually added at a later time-point.
The NIRStar recording program that came with the NIRScout device used to acquire the present data saved it in a NIRx-specific format-i.e., each participant's data was stored in a separate folder and included the following files: *.avg, *.dat, *.evt, *.hdr, *.inf, *.set, *.tpl, *.wl1, *.wl2, *config.txt, and *probeInfo.mat. Other devices save fNIRS data in different data formats, but all data formats always include a data file comprising the raw wavelength data, as well as a data file in which triggers and/or the optode configuration are saved.

Optode Configuration and Raw Data Conversion
In the first pre-processing step, we must store the individual optode configuration as

General Information on fNIRS Data Acquisition and Analysis
A visual summary of all processing steps is depicted in Figure 2. fNIRS data were obtained with a NIRScout 16-16 system (NIRx GmbH, Berlin, Germany) offering 16 sources and 16 detectors, which were divided to measure two interactants. In other settings, two or more devices can be synchronized to enable hyperscanning. The absorption of nearinfrared light was measured at the wavelengths of 760 and 850 mm and the sampling frequency was 7.81 Hz. The start of each condition was indicated via triggers sent through an experimental program, for instance, OpenSesame [27]. Varying task conditions can be indicated by triggering a second pin for the end of the condition, while preset and standardized task duration can be manually added at a later time-point.
The NIRStar recording program that came with the NIRScout device used to acquire the present data saved it in a NIRx-specific format-i.e., each participant's data was stored in a separate folder and included the following files: *.avg, *.dat, *.evt, *.hdr, *.inf, *.set, *.tpl, *.wl1, *.wl2, *config.txt, and *probeInfo.mat. Other devices save fNIRS data in different data formats, but all data formats always include a data file comprising the raw wavelength data, as well as a data file in which triggers and/or the optode configuration are saved.

Optode Configuration and Raw Data Conversion
In the first pre-processing step, we must store the individual optode configuration as a Source-Detector (SD) file. The SD file can be prepared using the SDgui as part of the Homer2 toolbox, in which the optode positions can be manually entered according to each study's optode layout. In our case, having used a NIRx system, all necessary information to do so is available in the *probeInfo.mat and *config.txt files located in the corresponding data folder. The probeInfo.mat file is automatically created after each recording and stems from a manual configuration using NIRSite set up before the start of data collection. In NIRSite, sources and detectors can be manually placed on MNI head models or imported from digitized optode coordinates. Optode localization only needs to be saved once-if the same optode layout was used for all participants and individual digitizer information is unavailable-and the thereby created SD file can then be used generically for all participants. In this specific case, the function createSD file will automatically convert the optode configuration into an SD file.
Once the SD file is prepared, we can run the NIRxtoSPM function to convert the stored data into a MATLAB structure with the *.nirs format, which comprises all needed information for further processing with Homer2. This can be done using both Homer2 and SPM for fNIRS, and their data conversion functions are freely available online. For other systems than NIRx, these conversion functions should be available from the respective manufacturer (e.g., Hitachi: https://www.nitrc.org/projects/hitachi2nirs, accessed on 12 June 2021, Shimadzu: https://www.nitrc.org/projects/shimadzu2nirs/, accessed on 12 June 2021).

Pre-Processing and Visual Quality Check
We start the pre-processing by loading the *.nirs data of an individual into the MAT-LAB workspace. Next, we go through a rough automatic data quality check (initial pruning) using the function enprunechannels. We specify the expected range of the data (V)-dRange = [0.03 2.5], a signal-to-noise threshold-SNRthresh = 10 (the lower, the more conservative; with a range from 1-13), and the range of the inter-optode distance-SDrange = [2. 5 3.5]. Please be aware that these parameters are different for each fNIRS recording device. The obtained output includes the variable SD.MeasListAct, which labels bad channels (per row) with a 0. This information can be used to guide the decision on which channels to exclude later.
Subsequently, we convert the raw wavelength data (illustrated in Figure 3) into raw optical density (OD; illustrated in Figure 4, upper panel) using the function spm_fnirs_calc_od. Following the initial pruning and data conversion, we proceed to motion correction. The choice of motion correction has been discussed in various other methodological papers [30][31][32]. We have found the MARA algorithm [32] to provide the most efficient motion correction when applied to both our adult and preschool child fNIRS data. MARA implements a smoothing based on local regression using weighted linear least squares and a 2nd-degree polynomial model. We use the default setting for the length of the moving window to calculate the moving standard deviation (L = 1), the threshold for artifact detection (th = 3), and the parameter that defines how much high-frequency information should be preserved by the removal of the artifacts (alpha = 5). This parameter also corresponds to the length of the LOESS smoothing window. Motion corrected OD data is illustrated in Figure 4, lower panel.   . The x-axis shows the time in frames (sampling rate: 7.81 Hz). The y-axis shows optical density units. deoxygenated hemoglobin (HbR) using the function hmrOD2Conc, with the partial path length factor set to its default [6 6]. We then extract the HbO time-series and plot it for all channels using the wavelet transform function wt. Following, we check plots for each channel for their integrity, a clearly visible heart band, and visually observable motion artifacts. Examples for different types of artifacts detected during this manual, visual quality check are shown in Figure 5. As a final result, the number of bad quality channels is identified and saved in an array to be later excluded from further analyses. Once the quality check is finalized, we go back to the OD data after motion correction and apply a band-pass filter, with low-and high-pass parameters of 0.5 and 0.01, respectively, and using a second-order Butterworth filter with a slope of 12 dB per octave. Spatial filtering at this point makes sure that the data is also usable for individual analyses (GLM) and excludes physiological noise to a certain extent.
In the last step, before interpersonal neural synchrony analyses, we once more convert the filtered OD data to concentration changes in HbO and HbR (µmol; illustrated in Figure 6). This time, we use the spm_fnirs_calc_hb function from the SPMforfNIRS toolbox, Subsequently, we continue with an additional manual, visual quality check. To do so, we convert the OD data to concentration changes (µmol/L) in oxygenated (HbO) and deoxygenated hemoglobin (HbR) using the function hmrOD2Conc, with the partial path length factor set to its default [6 6]. We then extract the HbO time-series and plot it for all channels using the wavelet transform function wt. Following, we check plots for each channel for their integrity, a clearly visible heart band, and visually observable motion artifacts. Examples for different types of artifacts detected during this manual, visual quality check are shown in Figure 5. As a final result, the number of bad quality channels is identified and saved in an array to be later excluded from further analyses.
Once the quality check is finalized, we go back to the OD data after motion correction and apply a band-pass filter, with low-and high-pass parameters of 0.5 and 0.01, respectively, and using a second-order Butterworth filter with a slope of 12 dB per octave. Spatial filtering at this point makes sure that the data is also usable for individual analyses (GLM) and excludes physiological noise to a certain extent.
In the last step, before interpersonal neural synchrony analyses, we once more convert the filtered OD data to concentration changes in HbO and HbR (µmol; illustrated in Figure 6). This time, we use the spm_fnirs_calc_hb function from the SPMforfNIRS toolbox, as it allows for age-dependent modification of the modified Beer-Lambert Law. In the current sample, children's age averaged at the rounded number of 5, and mother's age averaged at the rounded number of 36. The age-dependent parameters are automatically  In good quality channels (left panel), the heart band (located around period 0.5 s) is clearly visible and continuous and there are no major motion artifacts. In bad quality channels (right panel), the heart band is not visible or obscured by strong motion artifacts (channels 4, 8, 15), or there is a very low signal-to-noise ratio (channel 11).

WTC
To calculate interpersonal neural synchrony, several methods have been used, such as Pearson correlation and Robust Correlations [19,33]. We rely on Wavelet Transform Coherence (WTC), the most commonly used method to estimate interpersonal neural synchrony in fNIRS hyperscanning paradigms [21]. WTC is used to assess the relation between the individual fNIRS time series in each dyad and each channel as a function of frequency and time. We estimate WTC using the cross wavelet and wavelet coherence toolbox [22,23].
Based on earlier literature [25,26,34], visual inspection, and spectral analyses, a frequency band of interest needs to be identified. In previous publications on interpersonal neural synchrony during parent-child problem solving, the frequency band of 0.02-0.10 Hz (corresponding to ~10-50 s) was identified as most task-relevant. Accordingly, the same frequency band is used for the present sample. This frequency band should avoid high-and low-frequency physiological noise such as respiration (~0.2-0.3 Hz) and cardiac pulsation (~1 Hz).
Furthermore, coherence values outside the cone of influence need to be excluded in the WTC analysis due to likely estimation bias. Average neural coherence (i.e., interpersonal neural synchrony) can then be calculated for each channel combination and epoch. Epoch estimations, and, therefore, epoch length, are defined by the physical constraints of time-frequency analyses, i.e., the minimal time needed to estimate an appropriate coherence value for the indicated frequency range. Three oscillations in a certain frequency are required for WTC estimations [22]. For example, epoch lengths for a frequency of 2 s must, therefore, be at least 6 s long. To ensure reliable coherence value estimation, most studies have opted towards averaging coherence values over entire conditions. However,  (1) The factor condition had the cooperation condition as reference level. Estimates for the single predictors indicate the change from the response when the predictor changes from the reference level to the level of the predictor (in parentheses).
(2) Confidence intervals were derived using the Wald method using the function confint. Included are estimates, standard errors (SE), confidence intervals (CI), the standard deviation of the random effect (RE.SD), and likelihood ratio test outputs for the model comparisons.
In a post-hoc analysis, we contrast the conditions in each pairing condition, i.e., true vs. random (see Figure 7). In the true pairings, the cooperation condition has higher esti-  In a real and sufficiently powered data set, only those channels/averaged regions of interest where there is a significant difference between real and random pairings should be taken forward for interpersonal neural synchrony analyses.

Interpersonal Neural Synchrony
Next, we analyze whether there are condition and regional differences in interpersonal neural synchrony in mother-child dyads. We, therefore, estimate a second GLMM with the true WTC values only. Again, WTC is entered as the response variable. Condi- . The x-axis shows the time in frames (sampling rate: 7.8125 Hz). The y-axis depicts the oxygenation change in µmol.

WTC
To calculate interpersonal neural synchrony, several methods have been used, such as Pearson correlation and Robust Correlations [19,33]. We rely on Wavelet Transform Coherence (WTC), the most commonly used method to estimate interpersonal neural synchrony in fNIRS hyperscanning paradigms [21]. WTC is used to assess the relation between the individual fNIRS time series in each dyad and each channel as a function of Sensors 2021, 21, 4075 9 of 17 frequency and time. We estimate WTC using the cross wavelet and wavelet coherence toolbox [22,23].
Based on earlier literature [25,26,34], visual inspection, and spectral analyses, a frequency band of interest needs to be identified. In previous publications on interpersonal neural synchrony during parent-child problem solving, the frequency band of 0.02-0.10 Hz (corresponding to~10-50 s) was identified as most task-relevant. Accordingly, the same frequency band is used for the present sample. This frequency band should avoid high-and low-frequency physiological noise such as respiration (~0.2-0.3 Hz) and cardiac pulsation (~1 Hz).
Furthermore, coherence values outside the cone of influence need to be excluded in the WTC analysis due to likely estimation bias. Average neural coherence (i.e., interpersonal neural synchrony) can then be calculated for each channel combination and epoch. Epoch estimations, and, therefore, epoch length, are defined by the physical constraints of timefrequency analyses, i.e., the minimal time needed to estimate an appropriate coherence value for the indicated frequency range. Three oscillations in a certain frequency are required for WTC estimations [22]. For example, epoch lengths for a frequency of 2 s must, therefore, be at least 6 s long. To ensure reliable coherence value estimation, most studies have opted towards averaging coherence values over entire conditions. However, few studies have also calculated neural synchrony in smaller epochs to study how it behaves over time [35,36]. Here, we average coherence over the entire condition due to the broad frequency band and relatively short condition duration. Condition onset time can be extracted from the variable s included in the initial *.nirs files. Condition durations were manually added.
To run the WTC function, we extract each individual's time-series for each channel and combine that time-series array with another column including the timestamps of each row-and, thus, the information on the data's sampling rate. The timestamp array, if not already available, can be created using the following formula: t = (0:(1/sampling_rate):((size(time_series, 1) − 1)/sampling_rate))' We, then, input the following variables into the function wtc ([t, hboSub1(:, i)], [t, hboSub2(:, i)], 'mcc',0), which translates into timestamp and fNIRS time-series of HbO of Subject 1, timestamps and fNIRS time-series of HbO of Subject 2, and the count for a Monte Carlo Simulation. Generally used to estimate the significance of coherence values, we set the simulation to zero to increase computational speed. The output variables include Rsq (coherence values), period (period/frequency band), and coi (cone of interest). As we need to calculate coherence for each channel combination separately, the rows comprising the frequency band of interest are extracted once for each dyad and used for all channel combinations. The information in coi is used to set values outside of the cone of interest to missing values (NaN). To sum up the estimation, interpersonal neural synchrony is calculated by averaging the variable Rsq over a certain frequency of interest (y-axis) and task duration (x-axis), which results in a coherence value for each channel combination in each dyad.

Control Analysis
There are various approaches for additional control analyses. Here, we describe an approach in which we pair a random mother time-series to children's time-series and bootstrap the pairing 100 times with replacements. Accordingly, 100 WTC values are computed per channel combination and dyad combination with a fixed child, resulting in 2000 values overall. The thereby obtained WTC values are then averaged over each channel combination and pairing and result in a coherence value for each channel combination and child-caregiver dyad. Each true WTC value, therefore, is compared against a "randomly generated" value.

Statistical Analysis
For statistical analysis, the data is exported in a long data format and imported into RStudio (RStudio Team, 2015). Coherence values are bound by 0 and 1 and, therefore, assume a beta distribution. While a Fisher's r-to-z transformation would be appropriate and indicated for subsequent analyses with a linear (mixed) model assuming a gaussian data distribution, we use the untransformed coherence values to calculate generalized linear mixed models (GLMM) using the package and function glmmTMB [37]. We chose to cluster individual channels into regions of interest to enhance the reliability of interpersonal neural synchrony measures in specific cortical regions, even across minor differences in cap fit and, thus, probe settings [20,26,38]. Post-hoc contrasts are conducted using the function and package emmeans and corrected for multiple comparison using Tukey's Honest Significant Difference [39]. The distributions of residuals are visually inspected for each model. Models are estimated using Maximum Likelihood. Model fit is compared using a Chi-Square Test (likelihood ratio test; [40]).
In the current data set, we chose to focus on the analysis of HbO only, for both brevity and its consistent association with neural synchronization (e.g., [25,36,38]). However, we recommend repeating all analyses for HbR (and ideally Hb total) to derive a comprehensive picture of all possible neural synchronization processes [41].

Control Analysis
First, we analyze the true pair data in comparison to the random pair data using a GLMM. This initial analysis is to ensure that real coherence values are higher than mere spurious and autocorrelations in the data. WTC is entered as the response variable. Condition (cooperation vs. individual), region of interest (left dlPFC vs. right dlPFC vs. left TPJ vs. right TPJ), and pairing (true vs. random) are added as fixed and interaction effects.
Due to the small sample size, only random slopes for conditions and random intercepts for dyads are inserted to allow for model convergence. For a larger sample, additional random slopes can be added to the random effects structure.
wtc~condition + pairing + roi + condition: pairing + condition: roi + pairing: roi + condition: pairing: roi + (1 + condition | ID) The GLMM reveals that including an interaction effect between condition and pairing, as well as a fixed effect for region of interest significantly improves the model's explained variance, χ 2 (1) = 5.396, p = 0.020. Individual fixed effects and further interaction effects show no significant increase in explained variance, p > 0.056. The full details on model outputs are described in Table 1.
In a post-hoc analysis, we contrast the conditions in each pairing condition, i.e., true vs. random (see Figure 7). In the true pairings, the cooperation condition has higher estimated means Next, we also contrast the pairings in each condition. In the cooperation condition, true pairings have higher estimated means than random pairings. The difference, however, is not significant, p = 0.225. In the individual condition, true pairings have lower estimated means than random pairings. The difference is significant, t = −1.98, p = 0.047. The results are averaged over all regions of interest.

wtc ~ condition + roi + condition : roi + (1 + condition | ID)
The GLMM reveals that none of the included fixed and interaction effects between condition and region of interest significantly improve the model's explained variance, p > 0.130 (see Figure 8 and Table 2 for model outputs).  In a real and sufficiently powered data set, only those channels/averaged regions of interest where there is a significant difference between real and random pairings should be taken forward for interpersonal neural synchrony analyses.

Interpersonal Neural Synchrony
Next, we analyze whether there are condition and regional differences in interpersonal neural synchrony in mother-child dyads. We, therefore, estimate a second GLMM with the true WTC values only. Again, WTC is entered as the response variable. Conditions and regions of interest are entered as fixed and interaction effects. Additionally, we assume random slopes for conditions and random intercepts for each dyad.
wtc~condition + roi + condition : roi + (1 + condition | ID) The GLMM reveals that none of the included fixed and interaction effects between condition and region of interest significantly improve the model's explained variance, p > 0.130 (see Figure 8 and Table 2 for model outputs).     Note. (1) The factor condition had the cooperation condition as reference level. Estimates for the single predictors indicate the change from the response when the predictor changes from the reference level to the level of the predictor (in parentheses). (2) Confidence intervals were derived using the Wald method using the function confint. Included are estimates, standard errors (SE), confidence intervals (CI), the standard deviation of the random effect (RE.SD), and likelihood ratio test outputs for the model comparisons.

Exemplary Dataset
In the current exemplary dataset, we compared interpersonal neural synchrony in frontal and temporal brain regions in mother-child dyads during a cooperative versus individual problem-solving task.
In an initial analysis step, we compared data across all regions and both conditions from true pairs to data from random pairs. In doing so, we found that only true pairs show a trending increase in interpersonal neural synchrony during the cooperation condition as compared to the individual condition. Interpersonal neural synchrony in true pairings, however, is generally not higher than in random pairings in the current data set. Thus, we are unable to conclusively determine that interpersonal neural synchrony is significantly increased beyond spurious correlations in the signal, common task responses, or arousal [40]. Such initial control analysis is crucial when investigating interpersonal neural synchrony, as the results from this exemplary data set show that even random pairings of mothers' and children's fNIRS time-series can show high coherence.
To illustrate how to perform the subsequent interpersonal neural synchrony analyses in data from true pairs, we continued to assess the exemplary data, despite the control analysis not revealing a significant difference between true and random pairs. Accordingly, we investigated the true pairs' interpersonal neural synchrony data for differences between cooperative versus individual problem-solving, as well as between neural regions of interest (i.e., right versus left dorso-lateral prefrontal and temporo-parietal cortex). In this exemplary data set, however, such analysis did not reveal any significant effects.

General Considerations
In our approach, we chose to assess interpersonal neural synchrony across regions of interest (by clustering data from up to four paired channels) rather than from single channel pairings, in accordance with other literature [20]. It is, however, feasible to analyze and report interpersonal neural synchrony from single channel pairings, if appropriate correction for multiple comparisons is employed [42]. At the same time, we would advise caution when establishing interpersonal neural synchrony measures across all possible channel pairings. This is because when computing coherence between channels with different anatomical locations in the assessed individuals, results can be difficult to interpret without clear a-priori hypotheses. The above said, such an approach can yield meaningful results, for instance, if activity in a speaker's brain regions associated with speech production is linked with a listener's activity in regions associated with speech perception [43].
A related and more general consideration pertains to the issue of interpersonal neural synchrony data interpretation and, thus, the understanding of its (causal) implication for social interaction-as recently pointed out elsewhere [12,14,15,44]. Within this context, a theory of embodied mutual prediction was proposed [45], which intrinsically links the emergence of coherent brain signals across two (or more) interacting individuals with behavioral monitoring and action prediction. Along these lines, it is crucial that interpersonal neural synchrony data is assessed and interpreted in tight association with behavioral as well as potentially other physiological data obtained during the social interaction simultaneously with fNIRS hyperscanning data. Such a multi-level approach allows for a more comprehensive account of bio-behavioral synchrony across several behavioral, physiological, and neural systems [6]. Furthermore, what emerges from such considerations is the need for appropriate experimental control conditions, beyond comparing an active and interactive experimental task to a simple resting condition.
Most of the so far available studies on interpersonal neural synchrony-especially pertaining to fNIRS hyperscanning research in parent-child dyads-have focused on establishing that interpersonal neural synchrony is increased during cooperation. Such research has crucially advanced our understanding of the importance of bio-behavioral synchrony for social interaction by extending it from behavior and physiology to brain activation. As pointed out in the introduction, however, parent-child interaction is characterized by a fluctuation between aligned and misaligned states [5]. In other words, it is likely that not only short times of high interpersonal neural synchrony, but more extended periods of deand re-synchronization are determining the quality and future impact of social interactions for child development (see also [7]). Interpersonal neural synchrony should, therefore, be assessed over longer tasks or compared between different task epochs and probed for changes [35,36].
We also need to better understand differences in interpersonal neural synchrony patterns, such as across distinct brain regions, interaction partners, and age. For example, there is emerging evidence that patterns of interpersonal neural synchrony in parent-child dyads seem to diverge for mothers and fathers, with fathers showing a more anatomically localized synchrony pattern [25]. Mothers, on the other hand, seem to show a more anatomically widespread synchrony pattern in frontal and temporal brain areas [26,36]. Furthermore, in an older sample with parents and their school-aged children, more anatomically localized synchrony was evidenced in a cooperative button-press task, in comparison to an anatomically more widespread synchrony pattern in a competitive condition [38]. Differences in interpersonal neural synchrony within parent-child dyads have also been described in relation to the child's gender, with boys and girls demonstrating differential interpersonal neural synchrony patterns during cooperation with their mothers [46]. Another avenue for future research is to follow changes in patterns of interpersonal synchrony between parent-child dyads longitudinally and in relation to long-term child development and dyadic outcomes.

Limitations
The present guide is not meant to represent an exhaustive review of previous developmental fNIRS hyperscanning research (please refer to [18] for this purpose), nor a comprehensive overview on the current best practices for fNIRS publications (please refer to [47] for this purpose). We furthermore note the following limitations.
The exemplary data set provided here does not comprise an optode registration process, which is typically acquired using fMRI, 3D digitizers, and/or photogrammetry [48][49][50]. Although probe alignment is controlled to a certain degree by using standardized EEG caps aligned to participants' nasion and inion, variance in optode location could still cause some difficulties in estimating a reliable synchrony value. To avoid misalignment of channel estimations and, thus, synchrony estimates, we opt for a region of interest approach, meaning that synchrony is statistically assessed over groups of four adjacent homologous channels making up four regions of interest. Therefore, synchrony is derived and interpreted for a broader area of the cortex rather than a single homologous channel pair, similar to EEG hyperscanning [21]. Future investigations using fNIRS hyperscanning would, however, benefit from co-registration of optode localization for analyses, including optical reconstruction.
Early naturalistic interactions are still difficult to investigate with fNIRS, as motion artifacts can remain in the data, despite the use of motion correction algorithms. Especially rapid changes in probe distance and separation from the skin during data acquisition can cause spikes in the recorded signal, which commonly available algorithms cannot fully correct. Therefore, future studies need to take both more infant/child optimized optode settings as well as more flexible motion correction approaches into consideration.
Developmental fNIRS hyperscanning can also be affected by systemic effects, such as confounds through task-related arousal [51]. A potential remedy for this issue could be the concurrent assessment of autonomic nervous system and/or physiological responses in the periphery, but also using additional short channel measurements. Although these additional measurements can make hyperscanning with children even more challenging, advancements in wearable diffuse optical tomography might provide a feasible future solution [52].
How we best estimate interpersonal neural synchrony in parent-child dyads remains debated. Although WTC is emerging as a dominant method in fNIRS hyperscanning (see the methods section and [21] for more information), correlational measures are often employed as well [19,33]. If, however, differences in the hemodynamic response increase with age-due to different rates of synaptogenesis and angiogenesis in the developing brain-we might require a novel analysis approach, for instance, cross-frequency coherence/correlation. Moreover, the non-directional analysis with WTC does not allow for interpretations of directionality (i.e., leader-follower relationships) within interpersonal neural synchrony.

Conclusions
The goal of this article was to showcase a specific implementation to an existing fNIRS hyperscanning pipeline for parent-child dyads. We provided pre-processing steps, interpersonal neural synchrony estimation approaches, as well as the statistical analyses approach for future studies. We used openly available MATLAB toolboxes to build the pipeline, as well as R packages to analyze the data. By making the approach openly available, we hope to enable improvements in standardizing data processing approaches in the field and to enhance replicability of studies.