Hyper-Brain Networks Support Romantic Kissing in Humans

Coordinated social interaction is associated with, and presumably dependent on, oscillatory couplings within and between brains, which, in turn, consist of an interplay across different frequencies. Here, we introduce a method of network construction based on the cross-frequency coupling (CFC) and examine whether coordinated social interaction is associated with CFC within and between brains. Specifically, we compare the electroencephalograms (EEG) of 15 heterosexual couples during romantic kissing to kissing one’s own hand, and to kissing one another while performing silent arithmetic. Using graph-theory methods, we identify theta–alpha hyper-brain networks, with alpha serving a cleaving or pacemaker function. Network strengths were higher and characteristic path lengths shorter when individuals were kissing each other than when they were kissing their own hand. In both partner-oriented kissing conditions, greater strength and shorter path length for 5-Hz oscillation nodes correlated reliably with greater partner-oriented kissing satisfaction. This correlation was especially strong for inter-brain connections in both partner-oriented kissing conditions but not during kissing one’s own hand. Kissing quality assessed after the kissing with silent arithmetic correlated reliably with intra-brain strength of 10-Hz oscillation nodes during both romantic kissing and kissing with silent arithmetic. We conclude that hyper-brain networks based on CFC may capture neural mechanisms that support interpersonally coordinated voluntary action and bonding behavior.

There is neurophysiological evidence that coordinated behavior is accompanied by synchronized brain activity [7][8][9][10]18,19] and oscillatory coupling of other biological functions, such as respiration and cardiac activity [11]. To investigate these phenomena, various synchronization or coupling measures have been proposed and used [7,8,11,20]. Usually, when using timefrequency decompositions, brain networks are constructed and considered for specific frequencies [8][9][10]18,20]. However, based on conceptual considerations [21,22], it seems important to also consider coupling across frequencies. For instance, large-scale theories of neural organization (e.g., [23]) and prior research [22,24,25] strongly suggest that cross-frequency coupling (CFC) plays a critical rule in neural information exchange. To overcome the limitations of single frequency representations, we extended the phase synchronization measures used in our earlier work [8,11] to CFC. We then used these coupling indexes to construct hyper-brain networks that represent intra-and inter-brain synchronization within and across frequencies. Our novel approach allows researchers to represent the complex interplay among different frequencies in the context of a hyper-brain network (i.e., a network consisting of interacting nodes across two or more brains of interacting people). In comparison to all earlier approaches, where different brain sites (different electrodes in the case of the EEG) are defined as nodes, we considered single brain sites oscillating at different frequencies as different nodes by using a CFC measure. In other words, each node is a combination of an electrode location and of an oscillation frequency. Thus, the hyper-brain networks considered in this article consist of electrodefrequency nodes of two brains.
Measures derived from graph theory are increasingly being used in the neurosciences [26][27][28][29][30][31][32][33]. In modularity analyses, when using novel CFC approach, the same electrode can participate in different modules of the hyper-brain network dependent on the oscillation frequency. It is well known that single neurons and also brain areas can be involved in multiple overlapping cell assemblies [34][35][36][37]. Such assemblies oscillating synchronously at different frequencies provide an efficient basis for integrative processes in the brain [38]. CFC, allowing accurate timing between different oscillatory rhythms, may indicate integration or communication between different cell assemblies (cf. [39][40][41]). We assume that hyper-brain networks constructed by using CFC both within and

Psychological assessment
We assessed the partners' feelings during the test session, partnership and kissing satisfaction as well as the quality of the kissing. Most items were based on a 5-point rating scale ranging from 1 (not at all) to 5 (very much). The items are summarized in Table 1 and Table S1. The assessment was carried out during and after the EEG session. Kissing quality was assessed immediately after the RK and K-SA sessions, correspondingly. The other items were assessed after the entire EEG session.
EEG data acquisition and preprocessing EEG measurement took place in an acoustically and electromagnetically shielded cabin. Separate amplifiers with separate grounds were used for each individual, optically coupled to a computer. Under all task conditions, the EEG was simultaneously recorded from both participants using two electrode caps with 64 Ag/AgCl electrodes placed according to the international 10-10 system, with the reference electrode at the right mastoid. Vertical and horizontal electrooculograms (EOGs) were recorded to control for eye blinks and eye movements. Additionally, a bipolar lip EMG was obtained from two EMG electrodes placed over the orbicularis oris muscle. The sampling rate was 5000 Hz. Recorded frequency bands ranged from 0.01 to 1000 Hz. EEG recordings were re-referenced to an average of the left and right mastoid separately for each subject, resampled at 1000 Hz, and filtered with a band pass ranging from 1 to 70 Hz. Eye-movement correction was accomplished by independent component analysis [42]. Using the lip EMG channels, which were filtered with a highpass filter of 4 Hz, we manually pasted markers on the kiss onset. Spontaneous EEG activity was then divided into epochs of 3 sec, starting 500 ms before the kiss onset and ending 2500 ms after it. Artifacts from head and body movements were rejected by visual inspection, after an artifact rejection based on a gradient (a maximum admissible voltage step of 50 mV), and a difference criterion (a maximum admissible absolute difference of 200 mV between two values in a segment) had not rendered satisfactory results. Epochs that were artifact-free in both participants were selected for further analysis. To reduce the amount of data and to overcome the problem of volume conduction between neighboring electrodes, we selected 21 electrodes based on the 10-20 system (Fp1, Fpz, Fp2, F7, F3, Fz, F4, F8, T7, C3, Cz, C4, T8, P7, P3, Pz, P4, P8, O1, Oz, and O2). These electrodes are distributed across the entire cortex so that the information of the remaining electrodes would be rather redundant.
In the EEG and lip EMG data, power spectra were calculated using the Fast Fourier Transform (FFT) with the Hanning window. In the EEG spectra, we then determined the individual alpha peak frequency and the spectral power in the low (6.8-9.7 Hz) and high (9.7-12.6 Hz) alpha frequency bands separately for each kissing condition. In the lip EMG spectra, we determined average power at the individual maximum (610 Hz) and at 60 Hz (610 Hz). Power values were normalized using a logarithmic transformation.

Description of phase synchronization or coupling measures
To investigate phase coupling in a directed and frequencyresolved manner, we first applied an analytic or complex-valued Morlet wavelet transform to compute the instantaneous phase in the frequency range from 0 to 100 Hz in 0.5-Hz steps. The complex mother Morlet wavelet, also called Gabor wavelet, has a Gaussian shape around its central frequency f: in which s is the standard deviation of the Gaussian envelope of the mother wavelet. The wavelet coefficients were calculated with a time step of 5, leading to a time resolution of 5 ms and frequency resolution of 0.5 Hz. In order to identify the phase relations within and between any two channels or frequencies, the instantaneous phase difference was then computed from the wavelet coefficients for all possible electrode and frequency pairs within and between the brains. On the basis of instantaneous phases for two signals (X and Y) given as: and W Y (f n ,t) = arg[Q Y (f n ,t)], correspondingly, the n:m phase synchronization between two oscillations at the frequencies f m and f n were determined. The generalized phase difference (DW) according to n?f m = m?f n was calculated by: The n:m phase synchronization index (PSI) was then defined by: where ,N. denotes the averaging across time. In the case of within-frequency coupling (WFC) with f m = f n , PSI is calculated in the same way by setting m = n = 1. During calculation of the PSI, we not only determined the mean direction or the length of the vector but also the angle of this vector (h) in complex space: To determine the directed cross-frequency coupling, we used the adaptive Integrative Coupling Index (aICI) in this work. In contrast to our earlier studies [8,11], we used an adaptive algorithm, which allowed us to calculate coupling depending on the angle of phase differences determined in a given time window.
In other words, aICI no longer reflects in-phase synchronization, where the angle of phase differences is close to zero, but is suitable for the determination of phase coupling at any chosen or previously determined phase angles (e.g., h). Given the estimates of the phase difference between two signals, it is possible to ascertain how long the phase difference remains stable in defined phase angle boundaries by counting the number of points that are phase-locked at a defined time window. We slightly modified the procedure described in Müller and colleagues [8] in that we defined phase angle boundaries not related to phase zero but to the phase angle h. The further procedure was similar, as depicted in Figure 1. After the complex wavelet transform of the signals (Fig. 1A) and determination of PSI and h (Fig. 1B), we divided the range between h2p/4 and h+p/4 into two ranges and distinguished between positive and negative deviations from phase h (Fig. 1C). As shown in Figure 1D, we marked negative deviations in the range between h2p/4 and h in blue (coded with ''21'') and the positive deviations in the range between h and h+p/4 in red (coded with ''+1''). Phase difference values beyond these range were marked with green (coded with ''0'') and represent non-synchronization. In the case of two channels, X (e.g., Fz) and Y (e.g., Cz), a blue stripe in the diagram would mean that the phase of channel Y precedes the phase of channel X and a red stripe would mean that the phase of channel X precedes the phase of channel Y. We then counted the number of data points that are phase-locked separately in each of these two ranges. Before counting, successive points in the defined range (between h2p/4 and h+p/4) with a time interval shorter than a period of the corresponding oscillation at the given frequency (T i = 1/f i ) were discarded from the analysis. In the case of CFC, the lower frequency was considered. This cleaning procedure effectively eliminated instances of accidental synchronization. Figure 1E represents synchronization pattern of several electrode pairs after this cleaning procedure. On the basis of the counting described above, we obtained several synchronization indices: (1) the Positive Coupling Index, PCI, or the relative number of phase-locked points in the positive range (between h and h+p/4); (2) the Negative Coupling Index, NCI, or the relative number of phaselocked points in the negative range (between h2p/4 and h); (3) the Absolute Coupling Index, ACI, or the relative number of phase-  Coding of the phase difference of two signals, S1 (e.g., Fz) and S2 (e.g., Cz), at a given frequency (h2p/4,S1-S2, h: blue stripes; h,S1-S2,h+p/4: red stripes; S1-S2,2p/4 or S1-S2.+p/4: green stripes = nonsynchronization The aICI measure ranges from 0 to 1. It is an asymmetric measure (i.e., aICI AB ?aICI BA ) expressing both common (absolute) and ''positive'' contributions to phase synchronization. Phase differences and corresponding phase synchronization measures were determined for six different frequencies of interest (FOI): 5, 10, 20, 30, 40, and 60 Hz, and all corresponding combinations between them, reflecting different frequency relations, such as 1:2, 1:3, 1:4, 1:6, 1:8, 1:12, 2:3, and 3:4. In addition to the EEG-EEG coupling, we also determined lip-lip EMG and lip EMG-EEG coupling. In the case of lip EMG, only prominent 60-Hz oscillations were used. We restrict the description of our study results to the aICI measure, which is most informative due to its directionality.

Graph-theoretical approach (GTA)
Network construction. The coupling measures determined as described above were used to construct a connectivity matrix or a graph determining the network properties. In comparison to all earlier approaches, where different brain sites (different electrodes in the case of the EEG) were defined as nodes in such a graph, we were able to define single brain sites oscillating at different frequencies as different nodes by using a CFC measure. This means that the same electrode site at the six FOI (5, 10, 20, 30, 40, and 60 Hz) represents six different nodes that communicate with other nodes at the same or different frequencies. In other words, each node is a combination of spatial representation (electrode location) and of the oscillation frequency. The structure of such a graph is represented in Figure S1. The advantages of a network architecture allowing for CFC are: (1) not only connections between but also within the brain areas can be captured; (2) different brain areas can communicate with each other at multiple frequencies. In addition to the 21 EEG electrodes of each of the two kissing partners at six different frequencies, two EMG lip responses of female and male partners at the prominent frequency of 60 Hz were used for the construction of our hyper-brain networks. There were 254 nodes altogether (262166+2 = 254) in each network.
To investigate the network topology of the real networks, we also constructed regular (lattice) and random networks that have the same number of nodes and mean degree as our real networks. For these purposes, we randomized the edges in the real network to achieve a random network with the same number of nodes and edges. Lattice networks were configured like random networks, but in addition edges were redistributed after an initial random permutation such that they lay close to the main diagonal with increasing order of their weights. For these purposes, each column in the adjacency matrix was split into two parts around the diagonal. Further, all edges in these two parts were redistributed in increasing order, and then merged again into the column. Lattice networks reconstructed in such a way have the same number of nodes and edges as the initial real network but are characterized by ring or lattice topology incorporating nearest-neighbor connectivity [43]. These network reconstructions for random and regular networks were carried out 10 times for each individual network. Average network topology was then determined for these repeated reconstructions. Real and correspondingly reconstructed networks are displayed in Figure 2 for comparison. For this representation, grand averages of connectivity values were used.
Threshold determination. Before determining the network properties by means of the GTA in a next step, the threshold of the synchronization or coupling measures had to be calculated. For this purpose we generated surrogate data by (a) computing the amplitude and phase spectrum of a real signal using a Fourier transformation; (b) phase shuffling, whereby the phase values of the original spectrum are used in random order and the sorted values of the surrogate sequence are replaced by the corresponding sorted values of the reference sequence; and (c) inverse Fourier transformation back to the time domain. In this way, the real and the surrogate data retain the same power spectrum but a different time course. Surrogate data computed in such a way for all epochs at all considered channels were then used for the calculation of the corresponding synchronization measures. Thereafter, we applied a bootstrapping procedure with 1,000 resamples of the coupling measures gained from the surrogate data set and determined the threshold as the bootstrapping mean plus the confidence interval at a significance level of p#0.0001 separately for each frequency combination. Only coupling values greater than these threshold values were considered for network construction and further analyses.
In general, the choice of a threshold plays an important and nontrivial role in network construction, but is necessarily always arbitrary. In our case, the use of CFC makes this problem even more complex. To overcome this, at least partly, we (1) used different threshold values for different FOI and their combinations, as they were determined using the surrogate data procedure, and (2) applied a range of thresholds by the multiplication of initial threshold values by 10 different equidistant multiplication factors from 1.0 to 1.045 with a lag of 0.005. Thresholds determined in this way were used to construct different graphs with sparse connections covering a sparsity-range between 16% and 43% of the strongest edges preserved. Thus, we could compare the topological network properties at different sparsity or costs levels. As we were interested in coupling strengths, we used weighted networks at different sparsity levels.

Network metrics
Degrees and strengths. As aICI is a directed measure, we obtained the node in-and out-degrees in the network, in which the in-degree is the sum of all incoming connections of the node i k in in-and out-strength, respectively. Thus, the strength can be considered as the weighted degree [32]. For statistical evaluation, we determined out-strengths for each of the nodes of the whole hyper-brain network of each kissing couple. Additionally, we determined the within-brain out-strength for each of the kissing partners, and then calculated the between-brain out-strength by subtracting the within-brain out-strength from the hyper-brain out-strength.
Clustering coefficient (CC) and characteristic path length (CPL). If the nearest neighbors of a node are also directly connected to each other, they form a cluster. For an individual node, the clustering coefficient (CC) is defined as the proportion of the existing number of connections to the total number of possible connections. In the case of a weighted directed graph the mean CC is calculated by the formula [44]: being the number of weighted directed triangles around a node i. Another important measure is the characteristic path length (CPL). In the case of an unweighted graph, the shortest path length or distance d i,j between two nodes i and j is the minimal number of edges that have to be passed to go from i to j. This is also called the geodesic path between the nodes i and j. The CLP of a graph is the mean of the path lengths between all possible pairs of vertices [45]: where L i = CPL i is the average distance or average shortest path length between node i and all other nodes. In the case of a weighted and directed graph, the weight and direction of the links are considered.
Global (E glob ) and local (E loc ) efficiency. Global efficiency (E glob ) is defined as the average inverse shortest path length and is calculated by [46]: Like CPL, E glob is a measure of the integration of a network, but whereas CPL is primarily influenced by long paths, E glob is primarily influenced by short ones. Calculating E glob is advantageous over distance in disconnected networks: The efficiency between disconnected pairs of nodes is set to zero (the inverse of infinity).
Local efficiency (E loc ) is similar to the CC and is calculated as the harmonic mean of neighbor-neighbor distances [46]: Figure 2. Representation of real, lattice, and random networks. A: Coupling (aICI) matrices covering within-frequency coupling (WFC) and cross-frequency coupling (CFC) between the 254 nodes of the network. In the real hyper-brain network (left), the nodes are organized by electrode location (Fp1, Fpz, Fp2, F7, F3, …, O2), oscillation frequency (5, 10, 20, 30, 40, and 60 Hz), and brains (female, male); the last two nodes are lip EMG channels oscillating at 60 Hz for a female and a male, correspondingly. The lattice network (middle) was configured by the randomization of the edges in the real network and consecutive redistribution in such a way that the strongest edges lay close to the main diagonal. The random network (right) was configured by the randomization of the edges only. The lattice and random networks were reconstructed in such a way that they have the same number of nodes and edges as the initial real network, but are characterized by ring (lattice) or random network topology. B: The same networks as in A, represented in the form of a circle, where the nodes are in clockwise order beginning with 0u, marked with arrow. Note: more information about network construction can be found in Figure S2. doi:10.1371/journal.pone.0112080.g002 Like CC, E loc is a measure of the segregation of a network, indicating efficiency of information transfer in the immediate neighborhood of each node and showing how fault-tolerant the system is.
Small-worldness. To investigate the small-world (SW) properties of a network it has become common to compare its clustering coefficient and characteristic path length to those of regular lattices and random graphs. At least two specific properties of small-world network (SWN) related to control networks (random and lattice) are significant: (1) The CC of the SWN (CC SWN ) is much higher than that of random networks (CC SWN ..CC rand ), but the CPL of the SWN (CPL SWN ) is only slightly higher than that of the random network (CPL SWN $ CPL rand ), and (2) the CC of the SWN is lower than that of lattice networks (CC SWN #CC latt ), but the CPL of the SWN is much lower than that of the lattice network (CPL SWN ,,CPL latt ). Specific quantitative SW metrics were developed in addition to these main graph metrics. Foremost, the so-called SW coefficient s, is related to the main metrics of a random graph (CC rand and CPL rand ) and is determined on the basis of two ratios c~CC=CC rand and l~CPL=CPL rand [47]: The SW coefficient s has been used in numerous networks showing SW properties and has been found to be greater than 1 in the SWN (s.1).
The second SW metric was defined by comparing the CC of the network of interest to that of an equivalent lattice network and comparing the CPL of the network to that of an equivalent random network [48]: This metric normally ranges between 21 and +1 and is close to zero for SWN (CPL SWN <CPL rand and CC SWN <CPL latt ). In addition, positive values of v indicate a graph with more random characteristics (CPL SWN <CPL rand and CC SWN ,,CPL latt ), while negative values indicate a graph with more regular (lattice-like) characteristics (CPL SWN ..CPL rand and CC SWN <CPL latt ). The clear advantage of the v metric as compared to s is the possibility to define the extent to which the network of interest is like its lattice or random equivalents [48].
Community structures and definition of node roles within the brain networks. To further investigate the topological properties of the hyper-brain networks, community structures as well as indices of modularity (M), the within-module degree (Z i ) and the participation coefficient (P i ) were determined (cf. [32]). For this calculation, the modularity optimization method for directed networks [49] as implemented in the Brain Connectivity Toolbox [32] was used (for limitations regarding this method, see [50]). The optimal community structure is a subdivision of the network into non-overlapping groups of nodes in a way that maximizes the number of within-module edges, and minimizes the number of between-module edges. The modularity (M) is a statistic that quantifies the degree to which the network may be subdivided into such clearly delineated groups or modules. It is given for directed networks by the formula [49]: where l~P ij a ij is the number of edges in the graph, and a ij is defined to be 1 if there is an edge from j to i and zero otherwise, k in i and k out i are the in-and out-degrees of the node i, and d m i ,m j is the Kronecker delta. High modularity values indicate strong separation of the nodes into modules. M = 0 if nodes are placed at random into modules or if all nodes are in the same cluster [49]. To test the modularity of the empirically observed networks, we compared them to the modularity distribution (N = 100) of random networks, that is, to simulated networks with the same number of nodes and edges as the original network [51].
The within-module degree Z i indicates how well node i is connected to other nodes within the module m i . It is determined by [52]: where k i (m i ) is the within-module degree of node i (the number of links between i and all other nodes in m i ), k k(m i ) and s k(m i ) are the mean and standard deviation of the within-module degree distribution of m i . The within-module degree Z i is zero if all the nodes of the module have the same number of edges (e.g., if all the nodes within the module are fully interconnected with each other); otherwise it has negative or positive values depending on the number of links at the different nodes.
The participation coefficient P i describes how well the nodal connections are distributed across different modules [52]: where M is the set of modules, k i (m i ) is the number of links between node i and all other nodes in module m i , and k i is the total degree of node i in the network. Correspondingly, P i of a node i is close to 1 if its links are uniformly distributed among all the modules, and zero if all of its links lie within its own module. Z iand P i -values form a so-called Z-P parameter space and are characteristic for the different roles of the nodes in the network [52]. The community structures and corresponding measures (M, Z i , and P i ) were determined for the common networks of the kissing couples for each kissing condition separately.

Statistical analyses
Partner-oriented kissing satisfaction was analyzed using a twoway repeated measures ANOVA with a between-subject factor Sex (female vs. male) and a within-subject factor Items (3 Items, indicated in Table 1). Immediate kissing quality was assessed separately for RK and K-SA conditions and analyzed using a three-way repeated measures ANOVA, Sex6Items6Kissing (RK vs. K-SA). Alpha peak frequency was analyzed using a four-way repeated measures ANOVA with a between-subject factor Sex and three within-subject factors Kissing (RK, K-SA, and HK), Anterior-Posterior (frontal, central, parietal, and occipital), and Left-Right. Spectral power in the low and high alpha frequency band was subjected to a five-way repeated measures ANOVA with a between-subject factor Sex and four within-subject factors Kissing, Alpha Band, Anterior-Posterior, and Left-Right. Lip EMG at the individual maximum (610 Hz) and at 60 Hz (610 Hz) were analyzed using a two-way repeated measures ANOVA with a between-subject factor Sex and a within-subject factor Kissing (RK, K-SA, and HK). For statistical analyses of network properties, individual electrodes were collapsed into three regions along the anterior-posterior axis (frontal, central, and parieto-occipital) for each FOI separately. Strengths were statistically evaluated for hyper-brain networks and separately for within-and between-brain connections. Strengths, CC and CPL were analyzed using a four-way repeated measures ANOVA with a between-subject factor Sex (female vs. male) and three withinsubject factors Frequency (5,10,20,30,40, and 60 Hz), Site (frontal, central, and parietal), and Kissing (RK, K-SA, and HK). Greenhouse-Geisser epsilons were used in all ANOVAs for nonsphericity correction when necessary. The Scheffe test was employed for the post-hoc testing of kissing condition differences.
To assess correlations between kissing satisfaction and quality and electrophysiological data, Pearson product correlations were computed between the averaged scores of partner-oriented kissing satisfaction or kissing quality and the network or other EEG measures. The difference between the correlation coefficients was tested using Steiger's procedure [53]. Table 1 summarizes psychologically assessed data (means and standard deviations) of partner-oriented kissing satisfaction and kissing quality during the experiment in the female and male partners. Further results of the psychological assessment can be found in the Table S1. Immediate kissing quality was assessed separately for RK and K-SA conditions. A two-way repeated measures ANOVA Sex6Items for partner-oriented kissing satisfaction revealed only significant item differences (F(2,56) = 7.30, P = 0.003, g 2 = 0.21). A three-way repeated measures ANOVA Sex6Items6Kissing for immediate kissing quality revealed significant main effects of all three factors: Sex (F(1,28) = 5.83, P = 0.023, g 2 = 0.17), Items (F(2,56) = 8.62, P = 0.001, g 2 = 0.24), and Kissing (F(1,28) = 29.90, P,0.0001, g 2 = 0.52), but no significant interactions. As shown in Table 1, women generally reported higher immediate kissing quality during the experiment than men, and kissing quality was rated higher for RK than for K-SA.

EEG and lip EMG spectral power analyses
First, we determined alpha peak frequency and spectral power in low and high alpha frequency bands in each of the female and male partners during the three different kissing conditions. The alpha peak frequency differed by kissing condition (F(2,56) = 3.77, P = 0.029, g 2 = 0.12). A posthoc Scheffe test showed significant differences between K-SA and HK (K-SA.HK: Mean Diff. = 0.217, Crit. Diff. = 0.205, P,0.05). Spectral power in the low and high alpha frequency bands varied as a function of the kissing condition, as indicated by significant interactions, Kissin-g6Alpha Band (F(2,56) = 3.62, P = 0.041, g 2 = 0.12) and Kissin-g6Alpha Band6Anterior-Posterior (F(6,168) = 3.00, P = 0.025, g 2 = 0.10), and as a function of Sex (Sex6Alpha Band6Anterior-Posterior, F(3,84) = 8.33, P = 0.003, g 2 = 0.23 and 6Alpha Band6Left-Right, F(2,56) = 4.59, P = 0.025, g 2 = 0.14). In general, spectral power was highest during RK and lowest during K-SA in the low alpha frequency band, and lowest during HK in the high alpha frequency band, especially at parieto-occipital regions. The spectral power in the low alpha frequency band was higher in men than in women, especially at parieto-occipital sites, and higher in women than in men in the high alpha frequency band at occipital and frontal regions.
Second, we determined the average power of lip EMG at the individual maximum (610 Hz) and at 60 Hz (610 Hz). Lip EMG oscillations at 60 Hz were used together with EEG channels for the construction of hyper-brain networks (see Methods). A twoway repeated-measures ANOVA (Sex6Kissing) revealed only a significant main effect of Sex for average spectral power at the individual maximum, F(1,28) = 4.6, P,0.05, g 2 = 0.14, and a marginally significant main effect of Sex for average spectral power across 60 Hz, F(1,28) = 3.9, P = 0.058, g 2 = 0.12, indicating a higher spectral lip EMG power in women than in men. There were no significant differences in the frequency of the maximum amplitude.

Network analyses
As described in Methods and displayed in Figure S1, we constructed our hyper-brain networks using CFC measures between 21 EEG electrodes of each of the two kissing partners at six different frequencies and two EMG lip responses of female and male partners at the prominent frequency of 60 Hz. There were 254 nodes altogether (262166+2 = 254) in each network (see Figure S1 for details). To investigate the SW properties of the real networks, we constructed two different control networks: regular (lattice) and random networks containing the same number of nodes and edges. Figure 2 displays network structures for real, lattice, and random networks using grand average coupling values. As expected, the lattice network showed a regular structure with high connectivity between neighbors, while random networks had a random structure with equally distributed short-and long-range connections. Real networks showed SW structure as indicated by high connectivity between neighbors and also a portion of longrange connections resembling the random network. Oscillations at the alpha frequency appeared to fulfill a pace-setting function in the real network, showing CFCs with all the frequencies used here. This strong connectivity of the alpha oscillations to all other oscillation frequencies seems to be a characteristic property of the real network constructed on the basis of CFC. Whether this network topology is characteristic for kissing only or has a broader applicability remains to be seen. Interestingly, lattice and random networks showed a comparable number of modules (4 modules in this simulation that are indicated by the different colors in Figure 2) with a low modularity value for the random network (M = 0.05) and a high modularity value for the lattice (M = 0.40). The real network was divided into 9 modules with relatively high modularity value (M = 0.29).
The real and control networks for each couple were constructed for 10 different adaptive thresholds (see Methods for details), and the SW metrics were measured as a function of the threshold ( Figure 3). As expected, increasing thresholds resulted in lower costs, which indicate sparser networks (Fig. 3A). Sparser networks have lower global but higher local efficiency ( Fig. 3B and 3C), which are correspondingly related to a higher CPL and also a higher CC (Fig. 3D and 3E). Thus, sparsity in hyper-brain networks leads to higher segregation but lower integration of information flow. The small-worldness coefficient s was always greater than 1 and increased with lower costs (Fig. 3F) indicating SW properties for all networks independently of the threshold or sparseness. The other small-worldness coefficient v ranged between 20.3 and +0.3 for individual networks and decreased with lower costs (Fig. 3G), also indicating SW properties of the observed networks and a tendency to become more regular with higher sparseness. Further, we compared network characteristics at the eighth threshold level (f 8 = 1.035) with high sparsity and optimal SW parameters. Figure 4 displays the network structures of one of the kissing couples under the three different kissing conditions. As mentioned above, alpha frequency interacting with all the other frequencies plays an exceptional role in the hyper-brain networks (Fig. 4A). Alpha frequency oscillations showed also strong CFC to theta oscillations, which remain therefore in a specific connection to each other and constitute together the so-called theta-alpha subnetwork binding the brains of the kissing couple together (Fig. 4B). This subnetwork has been observed practically in all kissing couples participating in the study. Through strong interconnection of alpha-frequency nodes with other frequencies, this subnetwork is strongly incorporated in the whole hyper-brain network structure and has a particular importance.

Statistical Evaluation of GTA measures
Out-strength, Characteristic Path Length, and Clustering Coefficient. Statistical analyses carried out for the three graphtheoretical measures (out-strength, CPL i , and CC i , see Table 2) showed that: (i) strengths differed between kissing conditions, in that they were higher during RK and K-SA than during HK, (ii) CPL i also differed between kissing conditions, with the shortest path length during RK and K-SA than during HK, (iii) CC i showed only significant interaction Kissing6Site (F 4,112 = 4.39, P, 0.01) with higher clustering during RK at frontal sites and lower clustering at parietal sites. Additionally, all three measures showed significant main effects of Frequency and Site, with the strongest effect for alpha frequency (stronger coupling strength but shortest path length and lower clustering) and a predominance of coupling strength, CC i and CPL i (shortest path length) at parieto-occipital sites (s. Table 2). The kissing effect was also modulated by Frequency and Site as shown by the significant interaction Kissing6Frequency6Site for out-strength (F 20,560 = 3.12, P, 0.005, g 2 = 0.10), with stronger strengths at 10 and 40 Hz, and at parieto-occipital sites.
Intra-and inter-brain strengths. To test whether the strength differences between the kissing conditions were not only due to enhanced within-brain coupling, we determined strengths separately for intra-and inter-brain connections (s. Methods for details). We found a significant main effect of Kissing for both intra-and inter-brain strengths (s. Table 2). The strengths were higher during RK and K-SA than during HK for the intra-brain and higher during RK than during HK for the inter-brain coupling. The differences between kissing conditions for interbrain coupling were modulated by electrode site, F 4.112 = 3.62, P, 0.05, g 2 = 0.12, with stronger differences between kissing conditions at parietal sites. We also found significant main effects of Frequency and Site, with stronger coupling strength at alpha frequency and at parieto-occipital sites.
Additionally, we tested the strength of lip EMG channels using a two-way repeated-measures ANOVA (Sex6Kissing). There were no significant differences between sex groups or kissing conditions. Correlation between GTA measures and subjectively assessed kissing satisfaction. Strength and CPL i for thetaoscillation nodes (5 Hz) in the hyper-brain network and especially for inter-brain (but not for intra-brain) connections during RK and K-SA conditions correlated significantly with items of partneroriented kissing satisfaction (see Table 3 and Figure 5). Correlations during HK did not differ reliably from zero, and the correlation between partner-oriented kissing satisfaction and interbrain strength during RK was significantly different from the corresponding correlation during HK (Zd = 1,97, P,0.05). Intrabrain strength for alpha-oscillation nodes (10 Hz) during RK and K-SA correlated significant positively with immediate kissing quality assessed after K-SA. Furthermore, there were significant correlations between kissing quality assessed after K-SA and hyper-brain strength (positive) as well as CPL i (negative) during RK. Unexpectedly, kissing quality assessed after RK was negatively correlated to hyper-brain strength and CPL i during RK and to parieto-occipital hyper-brain strength during K-SA. . Hyper-brain network properties under the three kissing conditions. A: Changes in hyper-brain network costs dependent on the coupling threshold. B: Changes in global efficiency (E glob ) in the hyper-brain, regular (lattice), and random networks dependent on the coupling threshold. C: Changes in local efficiency (E loc ) in the hyper-brain, regular (lattice), and random networks dependent on the coupling threshold. D: Changes in characteristic path length (CPL) in the hyper-brain and random networks dependent on the coupling threshold. The CPL of regular networks was always equal infinity and is, therefore, not presented in the diagram. E: Changes in the clustering coefficient (CC) in the hyper-brain, regular (lattice), and random networks dependent on the coupling threshold. F: Changes in the small-worldness coefficient (s) in the hyper-brain network dependent on the coupling threshold. G: Changes in the small-worldness coefficient (v) in the hyper-brain network dependent on the coupling threshold. RK = romantic kissing, K-SA = kissing while performing silent arithmetic, and HK = hand kissing. Hyper-brain network: red line; regular network: blue line; and random network: green line. doi:10.1371/journal.pone.0112080.g003 Modularity structures during kissing Next, we investigated modularity structures of the 15 individual hyper-brain networks, which could be partitioned well into different modules or communities, with modularity values ranging between 0.24 and 0.47 (M = 0.32, SD = 0.06). Modularity values were reliably higher than values observed in random networks (p, 0.0001) indicating the presence of nonrandom community structures. The number of modules varied between 4 and 14 with an average of 8.4 (SD = 2.2) modules for all individual hyper-brain networks and kissing conditions. There were no significant differences in the number of modules and modularity between kissing conditions. To define how a node was positioned within its own module and with respect to other modules, we calculated the within-module degree (Z i ) and participation coefficient (P i ) of the node i for the hyper-brain networks of the kissing couples. The within-module degree measures how 'well-connected' node i is to other nodes in the module, whereas participation coefficient reflects how 'welldistributed' the links of the node i are among the other modules. Together, Z i and P i form the so-called Z-P parameter space, with different regions indicating specific universal roles of the nodes positioned in this space or these regions. Although some attempts to define these regions or roles were reported in the literature [54][55][56][57], in our opinion, the boundaries of the roles in the given Z-P parameter space can not always be identical when using different measures and having different network structures. Therefore, they should be determined individually based on some specific rules. Figure 6 displays the Z-P parameter space for the three kissing conditions including all kissing couples. It can be seen that the topology of the Z-P parameter space is symmetric to the horizontal axis at Z i = 0 and have a specific bulb-like form. Nodes with the within-module degree Z i $1.4 were defined as hubs and nodes with Z i ,1.4 as non-hubs. We adopt this separation Z-value from our previous study with guitarists [8], which corresponds to the definition of hubs as nodes containing many more edges than most of the nodes in the module (cf. [54]). The number of hubs in our hyper-brain networks calculated across all kissing couples in the three kissing conditions was 4.8%. The boundary for provincial and connector nodes was set to 0.65, separating 48.6% of connector nodes. This separation value for provincial and connector nodes is close to the value suggested by Guimerà and Amaral [54]. In addition, we distinguished ultra-peripheral nodes (P#0.05) and kinless nodes (0.9#P#1.0). Corresponding to this separation, the Z-P parameter space was divided into six different roles: R1 (P i #0.05) -ultra-peripheral nodes, R2 (Z i ,1.4; 0.05,P i ,0.65) -non-hub peripheral nodes, R3 (Z i ,1.4; 0.65,P i ,0.9) -non-hub connector nodes,  Figure 7 displays the structure of the Z-P parameter space of a kissing couple under the three kissing conditions, with the nodes belonging to different modules coded by colors. It can be seen that most of the connector and also peripheral hubs share the same two or three largest modules in the network (Fig. 7A). In most cases of the kissing couples investigated in the study, the largest module represents the aforementioned theta-alpha subnetwork. This module also has the strongest connections to the other modules (Fig. 7B). However, a more thorough view on the modularity structure of the kissing couple presented in Figure 7 indicates that the modular organization of hyper-brain networks have a more complex structure (see Fig. 7C). During RK, five out of nine (in total) modules (marked in blue, red, green, yellow, and aquamarine) are relatively large and contain alpha-oscillation nodes together with other frequency nodes distributed across the two brains. The blue module shares theta-and alpha-frequency nodes in the two brains and represents the above-mentioned theta-alpha subnetwork. Furthermore, fronto-temporal alpha-frequency nodes in male brain share the beta-frequency nodes distributed across the female brain in the common module marked in red. The yellow module shares the frontal (and one occipital) alphafrequency nodes in the male brain with 30-Hz-oscillation nodes distributed across the female brain. The green and aquamarine modules belong practically to the female brain and bind together alpha and gamma oscillations in the female brain. All this leads to a very intertwined hyper-brain structure during RK. During K-SA, there are two (marked in blue and red) out of seven (in total) largest modules, whereby all the nodes in the female brain belong to these two largest modules and share these modules with alphafrequency nodes in the male brain. Overall, these two modules represent two very strong hyper-brain subnetworks: (i) alphagamma subnetwork (blue) and (ii) theta-alpha-beta subnetwork (red). During HK, there are two large modules (marked in blue and red) that share alpha-frequency nodes in the male brain with those of the female brain, whereby all the frequencies, with exception of 60-Hz oscillations in the female brain, join in these two modules; moreover, the nodes belonging to the first largest module (marked in blue) lie in the fronto-central regions in both the female and male brains, whereas the nodes from the second largest module (marked in red) are localized in the temporal and parieto-occipital regions, also in both brains. It thus appears that hyper-brain networks in kissing couples have a complex modular organization, in which alpha-frequency oscillations and their subnetworks, especially the theta-alpha subnetwork, play a crucial role. Figures 7D and 7E show the strongest connections within and between a couple's brains, respectively. It should be noted here that intra-brain networks based on CFC have very strong large-scale connections binding distributed cell assemblies, especially between the frontal and parieto-occipital areas. These regions are also strongly interconnected or synchronized between the brains. Inter-brain coupling generally reach out from the frontal regions of one partner to the parieto-occipital or central regions of the other partner and vice versa, whereas frontal-tofrontal connections are mostly reduced or attenuated.
Further, we calculated the numbers of nodes lying in the different regions of the Z-P parameter space for women and men separately under the three kissing conditions. (see Table 4). A twoway repeated measures ANOVA (Kissing6Sex) showed no significant differences between kissing conditions or sex groups   Table 3. Correlation coefficients for strength (hyper-, Intra-, and inter-brain) and characteristic path length (CPL i ) with average scores of subjective partner-oriented kissing satisfaction and immediate kissing quality during the task. at all. There is a tendency for a higher number of hubs and nonhub connectors and a smaller number of non-hub peripheral nodes during RK, compared with control conditions, but these differences were not statistically reliable.

Discussion
In this article, we reported a network analysis based on graph theory to examine cross-frequency coupling within and between Figure 5. Correlation plots displaying the correlation between strengths (hyper-, intra-, and inter-brain) and subjectively assessed partner-oriented kissing satisfaction and immediate kissing quality. A: Correlation between the hyper-brain strength (5-Hz oscillation nodes) and partner-oriented kissing satisfaction under the two kissing conditions (RK and K-SA). B: Correlation between the inter-brain strength (5-Hz oscillation nodes) and partner-oriented kissing satisfaction under the two kissing conditions (RK and K-SA). C: Correlation between the intra-brain strength (10-Hz oscillation nodes) and kissing quality (assessed after K-SA) under the two kissing conditions (RK and K-SA). doi:10.1371/journal.pone.0112080.g005 the brains of kissing couples. Our key finding is that CFC hyperbrain networks consisting of the two brains of the interacting people have SW properties, and that alpha frequency oscillations play an essential role in these networks, as they set the pace for other frequencies in the network. Furthermore, we identified theta-alpha subnetworks that bind together the two brains of kissing partners, especially during partner-oriented kissing, and that play a crucial role during interpersonal action coordination. Our analyses introduce a new method for capturing CFC in hyper-brain networks.

Advantages of a network architecture based on cross-frequency coupling
The human brain is a complex system that shows temporally coherent activity at multiple scales of time and space. This activity constitutes functional networks of different topology intermediate between highly regular lattices and random graphs. This topology is called SW topology [45]. We have found that our complex hyper-brain networks based on CFC within and between the brains of interacting people also possess SW topology, with the high degree of clustering found in a lattice and the short path length found in a random graph. We investigated our real and control networks for each participating couple by using 10 different adaptive thresholds corresponding to different sparsity levels that increase by threshold. Sparser networks showed higher local efficiency but lower global efficiency, which are related to higher CC and also higher CPL. Despite the fact that sparsity in hyper-brain networks led to higher segregation, but also to lower integration of information flow, the small-worldness coefficient s increased with lower costs and was always greater than 1, indicating that real networks retained small-world network topology at all threshold levels. Moreover, the other smallworldness coefficient v ranged between 20.3 and +0.3, also indicating that the observed networks belong to SWNs. About half of the kissing couples showed more random properties (v.0), whereas the other half showed a more regular network properties (v,0). A decrease of the coefficient v with a higher threshold or lower costs indicates a tendency of networks to become more regular. For further analyses, we chose the threshold level that provided a high sparsity of networks with about equal local and global efficiency. Modularity analyses showed that hyper-brain networks at this threshold level had mean modularity values at about 0.3, which were statistically always higher than that in random networks. This indicates nonrandom community structure in hyper-brain networks [56]. Furthermore, community structures were organized by combining electrode location and oscillation frequency, whereby each electrode oscillating at different frequencies participated in multiple modules. This gives rise to overlapping community structures within and between the brains. Normally, low-frequency nodes (e.g., theta and alpha) were distributed across the two brains and composed the so-called hyper-brain modules sharing electrodes/nodes from two brains (c.f. [8,9]). The largest hyper-brain module was normally the theta-alpha subnetwork, which could be found in practically all kissing couples. Besides this largest hyper-brain module, there were also smaller ones composed from other also high-frequency nodes. Such hyper-brain modules were also found in our earlier studies based on within-frequency analyses [8,9]. Hyper-brain networks based on CFC, however, have the advantage that they show how different frequencies interact with each other and build up overlapping community structures. This allows a better understanding of network topology and its organizational principles [58,59]. Neural cell assemblies are normally organized on the principle of overlapping communities, with neural units being shared by different overlapping cell assemblies [34][35][36][37]. The networks observed here follow the same principles and are therefore suitable for the investigation of hyper-brain networks arising during interpersonal action coordination.

Kissing satisfaction, lip EMG and EEG alpha activity
In terms of subjective experience, we found that women and men did not differ regarding partner-oriented kissing satisfaction but that women generally reported higher immediate kissing quality during the experiment than men. Not surprisingly, kissing quality was rated higher during RK than during K-SA. Women Figure 7. Modular organization of hyper-brain networks under the three kissing conditions. A: Scatterplots of Z-P parameter space with corresponding role regions (cf. the. 5). Different modules are coded with color of the circles. B: Circle modularity structure. The size of the circle (module) represents the common connectivity strength of the module, and connectivity strength between the modules is coded by line thickness. In both cases, the out-strengths were used for calculation. C: Modular organization of the female and the male brains. Each electrode contains six nodes representing different oscillation frequencies (5,10,20,30,40,and 60 Hz) in clockwise order, beginning from the top. The size of the circle corresponds to the out-strength of the node, and modules are represented by color, which is the same as in A and B. D: Within-brain connections in the female and the male brains, correspondingly. Coupling strength range from blue (low coupling) to red (high coupling). E: Between-brain connections between the female and the male brains. Coupling strength range from blue (low coupling) to red (high coupling). doi:10.1371/journal.pone.0112080.g007 and men did not differ in peak alpha frequency. However, alpha peak frequency was higher during K-SA than during RK or HK, which may reflect the cognitive demand of executing both tasks (kissing and silent arithmetic) simultaneously. Additionally, spectral power in the low and high alpha frequency bands varied as a function of the kissing conditions and as a function of sex. The spectral power was generally highest during RK and lowest during K-SA in the low alpha frequency band, and lowest during HK in the high alpha frequency band, especially at parieto-occipital regions. This means that RK enhances spectral power compared with HK in both low and high frequency bands, whereas K-SA is associated with higher spectral power in the high frequency band only. Likewise, men showed enhanced spectral power in the low alpha frequency band, especially at parieto-occipital sites, whereas women showed enhanced power in the high alpha frequency band, specifically at occipital and frontal regions. Despite this variation of alpha spectral characteristics by the sex and kissing conditions, neither alpha peak frequency nor alpha spectral power correlated significantly with kissing satisfaction or quality.
We also determined the average power of lip EMG across the individual maximum (610 Hz) and at 60 Hz (610 Hz), and the latter was used together with EEG channels for the construction of hyper-brain networks. Women showed generally higher lip EMG power than men, indicating stronger lip activity during kissing, independently of kissing condition. Lip activity also did not show any reliable correlations with kissing satisfaction or immediate quality of kissing during the experiment. It indicates that neither alpha band activity nor lip muscle activity are responsible for kissing satisfaction, at least not responsible alone.

Kissing from the viewpoint of network properties
Network analyses showed that hyper-brain networks during partner-oriented kissing (especially during RK) as compared with hand kissing have elevated strength and shortest path length indicating higher connectivity between different brain regions within and between the brains, and also optimized information transfer within and between them. The fact that higher strength during partner-oriented kissing was found for both intra-and inter-brain connections indicates that neural cell assemblies synchronize stronger not only within the brains of kissing partners but also between their brains. Thus, this evolutionarily important activity is associated with enhanced synchronicity between the brains of kissing partners. Such elevated inter-and also intra-brain connectivity was previously found with guitar duets [7][8][9][10], during spontaneous imitation of hand movements [19,60], and during card playing [18,20]. In these studies, cross-frequency coupling was not considered, and there was a preference for low (delta, theta) communication frequencies (cf. [7][8][9]) to bind two interacting brains. In other studies, especially when synchronization was measured across time [8,9,[18][19][20], higher frequencies (e.g., alpha, beta, and gamma) were involved as well. The important role of alpha observed in this study, however, has not been observed previously. Our results suggest that the important role of alpha is due to CFC, in the sense that the alpha frequency sets the pace for other frequencies. This interpretation is consistent with the available evidence. Alpha oscillations serve an important function in top-down attentional processes [61][62][63], and they synchronize with oscillations at other frequencies in response to cognitive demands [59,62,64,65]. The theta-alpha subnetwork identified in this study is in good agreement with these earlier results. This subnetwork comprises nearly all electrode sites in the two brains oscillating at theta or alpha frequencies, whereby the electrodes oscillating at alpha frequency (mostly at parietooccipital sites) fulfill the role of a connector hub, with strong connections to the other nodes of the subnetwork as well as to nodes or electrodes belonging to other modules or subnetworks. Based on these observations, we tentatively conclude that parietooccipital sites in the kissing partners may have played a leading or integrating role during kissing.
Kissing differs from other types of social interaction, such as music or card playing, by stronger reciprocal sensory and motor connections, including afferent and efferent feedback loops. Additionally to EEG couplings, we observed significant coupling between the lip EMGs of the two partners (lip-lip), as well as couplings between lip EMG and EEG within and across interaction partners. Within each partner, the EEG oscillations with the same frequency as the lip EMG (i.e., 60 Hz) were mostly strongly connected with the lip EMG. On the other hand, alpha frequency EEG oscillations were also strongly coupled, both with the own lip EMG and the lip EMG of the kissing partner. Couplings between EMG and EEG have been reported previously [66][67][68], but, for the first time, we show that EMG activity is synchronized with EMG and EEG activity of the interacting partner.
A further finding of this study is the significant correlation between strength (and also CPL i ) and subjectively assessed kissing satisfaction and quality. The coupling strength and integrative capacity of the frontal nodes and also of other brain regions may underlie this association. We only tested 5-and 10-Hz oscillations, which play an important role in the between-brain connectivity binding two brains together mostly in terms of theta-alpha subnetworks distributed across the two brains of kissing partners. Frontal theta hyper-brain strength during RK and K-SA and CPL i during RK showed reliable relationship to partner-oriented Table 4. Number of nodes (mean and standard deviation in %) lying in the specific regions of the Z-P parameter space separately for women and men under the three kissing conditions. kissing satisfaction. Importantly, this relationship was strongest for inter-brain connections covering all brain regions during RK, and frontal and parietal regions during K-SA. Given the fact that intrabrain strengths did not show a reliable relationship, it can be concluded that partner-oriented kissing satisfaction is associated with distributed inter-brain connectivity and frontal nodes in hyper-brain networks, showing strong connectivity and integration processes (shortest path length) in common hyper-brain network. Given the relatively low spatial resolution of EEG, we can only speculate about the neural circuitry below the frontal electrodes that form the substrate for this hyper-brain partner-oriented kissing-satisfaction effect. Neural activity in the medial prefrontal cortex is selectively enhanced during theory-of-mind tasks and mentalization [69][70][71]. Strong involvement of frontal regions in the within-brain, and especially the between-brain, synchronization has also been found with guitar duets [7][8][9]. Thus, these coupling strengths may reflect a synchronization of cell assemblies representing the coordination of one's own behavior with the behavior of the interaction partner [2,3,17,72]. This line of reasoning is further supported by the observation that kissing satisfaction correlated reliably with inter-brain but not with intrabrain strength. Thus, the significant relationship of kissing satisfaction with strength and CPL i during partner-oriented kissing (RK and K-SA) but not during hand kissing suggests that orientation toward the partner during kissing not only enhances frontal hyper-brain circuitry and distributed inter-brain coupling strength, but also induces neural integration (reduced CPL i ) in the whole hyper-brain network. This, in turn, may enhance the coordination of kissing behavior and enhance kissing satisfaction. Surprisingly, the immediate effect of kissing quality during the experiment showed a different relationship with network indicators dependent on whether kissing quality was assessed after the RK or K-SA test session. Kissing quality assessed after K-SA showed a positive association with hyper-brain strength and a negative association with CPL i for 10-Hz oscillation nodes during RK, whereas kissing quality assessed after RK showed an inverse relationship with hyper-brain strength and CPL i during RK, and also with parieto-occipital hyper-brain strength during K-SA. At the same time, intra-brain strength for alpha-oscillation nodes during RK and K-SA correlated positively with immediate kissing quality assessed after K-SA. This relationship with kissing quality assessed after RK was not reliable and rather negative. It seems that the assessment of kissing quality after a K-SA session is more veridical than after RK. Most of the participants reported after this session (K-SA) that combining kissing with performing silent arithmetic has been circumstantial, and they probably paid more attention to kissing and assessed them more differentially.
Kissing quality correlated reliably with intra-brain strength, but not with inter-brain strength, whereas the relationship to partneroriented kissing satisfaction was inverse, that is, it correlated reliably with inter-brain, but not with intra-brain strength. It seems that kissing quality is related above all to within-brain dynamics, whereas partner-oriented kissing satisfaction addresses above all the dynamics between the brains. Moreover, the former is related to alpha oscillations, whereas the latter to theta oscillations, which are slower. In a previous study [8], we also found that the interbrain connectivity was operating at lower frequencies than intrabrain connectivity. It remains to be seen whether these observations form part of a more general phenomenon.

Limitations
The present experiment has limitations and leaves room for questions to be addressed in future research. First, we used only six frequencies of interest for network construction. Using more finegrained frequency components would lead to a more differentiated representation of hyper-brain networks. At the same time, our results clearly show the benefits of including CFC in network construction. Second, our analyses were limited to phase-to-phase CFC. Other types of CFC (e.g., power to power, phase to power, power to frequency, or phase to frequency) are likely to provide further information about functionally relevant network properties [22]. Third, the coupling measures used in this study were linear and bivariate. Nonlinear and multivariate couplings [73,74] may also contribute to inter-brain dynamics and should be investigated in the future. In conjunction with high-density behavioral assessments, such measures may shed further light on the behavioral and neuronal dynamics of interpersonal action coordination. Finally, the order of RK and K-SA was fixed instead of counterbalanced to safeguard the ecological validity of RK. Hence, an unknown portion of the effects attributed to condition may have been due to sequence effects.

Conclusion
Methodologically, the results of this study show that hyper-brain networks constructed on the basis of CFC have clear advantages in the investigation of neural mechanisms of interpersonal action coordination. Such networks consider the interactions of different frequencies and incorporate overlapping community structures, and thereby provide a more complete representation of network topology and organization. We showed that CFC-based hyperbrain network topology differ between partner-oriented and solitary kissing. Significant relations of network properties, such as strength and average shortest path length (CPL i ), to subjectively assessed partner-oriented kissing satisfaction, point to the functional significance of hyper-brain network properties for interpersonal action coordination. Finally, this study shows that data acquisition and analysis methods for simultaneous EEG and EMG recordings from multiple persons are important tools to examine inter-brain oscillatory couplings during interpersonal interactions. Figure S1 Construction of the hyper-brain network using a CFC approach. A: Coupling (aICI) matrix covering within frequency coupling (WFC) and cross frequency coupling (CFC) between the 254 nodes of the kissing couple's hyper-brain network. The nodes are organized by electrode location (Fp1, Fpz, Fp2, F7, F3, …, O2), oscillation frequency (5,10,20,30,40, and 60 Hz), and brain (female, male); the last two nodes are lip EMG channels oscillating at 60 Hz for a female and a male, correspondingly. B: The same network as in A, represented in the form of a circle, where the nodes are in clockwise order for the female and the male nodes, representing 21 electrodes for each of the six frequencies used for network construction. The last two nodes are the lip EMG channels. It can be seen that most of the long-range connections are between the theta (5 Hz) and alpha (10 Hz) frequency nodes of the female and the male brains, representing the so-called theta-alpha subnetwork. (TIF)

Supporting Information
Table S1 Psychological assessment of partner-and relationship-oriented satisfaction. (DOCX) during data acquisition. The authors thank Julia Delius for language assistance.