Quantifying team cooperation through intrinsic multi-scale measures: respiratory and cardiac synchronization in choir singers and surgical teams

A highly localized data-association measure, termed intrinsic synchrosqueezing transform (ISC), is proposed for the analysis of coupled nonlinear and non-stationary multivariate signals. This is achieved based on a combination of noise-assisted multivariate empirical mode decomposition and short-time Fourier transform-based univariate and multivariate synchrosqueezing transforms. It is shown that the ISC outperforms six other combinations of algorithms in estimating degrees of synchrony in synthetic linear and nonlinear bivariate signals. Its advantage is further illustrated in the precise identification of the synchronized respiratory and heart rate variability frequencies among a subset of bass singers of a professional choir, where it distinctly exhibits better performance than the continuous wavelet transform-based ISC. We also introduce an extension to the intrinsic phase synchrony (IPS) measure, referred to as nested intrinsic phase synchrony (N-IPS), for the empirical quantification of physically meaningful and straightforward-to-interpret trends in phase synchrony. The N-IPS is employed to reveal physically meaningful variations in the levels of cooperation in choir singing and performing a surgical procedure. Both the proposed techniques successfully reveal degrees of synchronization of the physiological signals in two different aspects: (i) precise localization of synchrony in time and frequency (ISC), and (ii) large-scale analysis for the empirical quantification of physically meaningful trends in synchrony (N-IPS).

A highly localized data-association measure, termed intrinsic synchrosqueezing transform (ISC), is proposed for the analysis of coupled nonlinear and non-stationary multivariate signals. This is achieved based on a combination of noiseassisted multivariate empirical mode decomposition and short-time Fourier transform-based univariate and multivariate synchrosqueezing transforms. It is shown that the ISC outperforms six other combinations of algorithms in estimating degrees of synchrony in synthetic linear and nonlinear bivariate signals. Its advantage is further illustrated in the precise identification of the synchronized respiratory and heart rate variability frequencies among a subset of bass singers of a professional choir, where it distinctly exhibits better performance than the continuous wavelet transformbased ISC. We also introduce an extension to the intrinsic phase synchrony (IPS) measure, referred to as nested intrinsic phase synchrony (N-IPS), for the empirical quantification of physically meaningful and straightforward-to-interpret trends in phase synchrony. The N-IPS is employed to reveal physically meaningful variations in the levels of cooperation in choir singing and performing a surgical procedure. over time. Standard algorithms, on the other hand, employ rigid basis functions with fixed frequencies (cosine in STFT, wavelet function in DWT), and therefore linear superpositions of additional harmonic components are required to represent nonlinear signals, thus spreading energy over a wider, higher frequency range and not admitting physical interpretation. Applications of EMD range from biosignal analysis [13,14], through to mechanical systems [15] and seismology [16].
Owing to the empirical nature of EMD, its direct component-wise application to multivariate signals may result in: (i) IMFs with different intrinsic oscillatory components (modes) across multiple data channels for a given IMF index-a phenomenon known as mode mixing, and (ii) multiple IMFs containing similar oscillatory modes for a given data channel-a phenomenon referred to as mode splitting. To mitigate these problems in multivariate scenarios, several extensions to EMD have been proposed, which include the complex EMD [17], rotation-invariant complex EMD [18], bivariate EMD [19], trivariate EMD [20], multivariate EMD (MEMD) [10,21] and noise-assisted MEMD (NA-MEMD) [9]. The general multivariate MEMD has found applications in brain-computer interface [22,23], image processing [24,25], nuclear engineering [26] and system characterization [8].
In spite of the undoubted usefulness of the EMD algorithm and its variants, a rigorous theoretical description for the underlying algorithms is still lacking. To this end, the synchrosqueezing transform (SST or WSST) was proposed in [27]. It is a post-processing technique, originally based on the continuous wavelet transform (CWT), for the generation of highly localized TF representations of nonlinear and nonstationary signals. It reassigns the wavelet coefficients in scale or frequency by combining the coefficients that contain the same instantaneous frequencies, such that the resulting energy is concentrated around the instantaneous frequency curves of the modulated oscillations. The CWT, however, has a limited TF resolution, and because the CWT and STFT are both extensively used to analyse and process multicomponent signals, a natural extension of the univariate CWT-based SST (WSST) to the STFT setting was proposed in [28] and is referred to as STFT-based SST (FSST) (see appendix A.2 for more details). To identify oscillations common to multiple data channels, an extension of the univariate WSST to multivariate cases was proposed in [29], the so-called multivariate synchrosqueezing transform (MSST or W-MSST). It employs the WSST channel-wise to obtain the concentrated coefficients, and then estimates the multivariate instantaneous frequency by combining, for each frequency band, the instantaneous frequencies across all the channels using the joint instantaneous frequency (see [27,[29][30][31] for more details). However, the performance of the W-MSST degrades with noise power, because: (i) the joint instantaneous frequency estimator is sensitive to noise and (ii) the CWT produces mathematical artefacts (additional noise)-wavelet coefficients which correspond to undesired frequency components. The additional noise (artefacts) generated by the W-MSST can be reduced by STFT-based MSST (F-MSST) (see appendix A.3 for more details).
To mitigate the aforementioned problems posed by both the intrinsic data association measures-poor time and frequency localization in IPS and loss of time information in ICoh-we, therefore, propose a highly localized data association measure which is achieved based on the combination of NA-MEMD and STFT-based SST and MSST (F-M/-SST) algorithms. The NA-MEMD is first employed to obtain physically meaningful intrinsic oscillations of a given multivariate signal. Owing to the resolution and noise problems in the CWT-based SST and MSST, we employ the STFT-based SST and MSST algorithms to generate highly localized TF representations of signal dependence between the intrinsic oscillations produced by the NA-MEMD. This procedure is referred to as the intrinsic synchrosqueezing coherence (ISC).
However, for certain scenarios, in the presence of -collaborative tasks which take place over a long period of time, -multiple occurrences of long and complex events during the tasks, -physically meaningful and straightforward-to-interpret trends in the level of cooperation during the events being of interest, and -prior knowledge of periods (i.e. frequency ranges) of the trends for physically meaningful and straightforward interpretation being unavailable, it is imperative to have a data-driven data-association measure which empirically quantifies intrinsic trends, whereby only those components which contain physically meaningful interpretation can be combined. To this end, we propose an extension to the standard IPS, referred to as the nested intrinsic phase synchrony (N-IPS), which further decomposes time series of the degrees of synchrony between data channels obtained using the standard IPS into multiple scales of synchrony; then only certain scales which admit meaningful physical interpretation are empirically combined (i.e. summed), without any prior knowledge of the frequencies of such scales. The aim of this study is to build upon the enhanced discrimination capability of the ISC and N-IPS data association metrics, in order to precisely identify physiological synchrony in frequency and time, using the ISC, and to empirically quantify physically meaningful intrinsic trends in the level of cooperation associated with events during the course of long cooperative tasks, using the N-IPS. We employ the ISC to precisely reveal the synchronized respiratory and HRV frequencies over time among a subset of bass singers of the Eric Whitacre Choir during 4 min and 40 s rehearsal and performance of the same musical score. The N-IPS is employed to empirically obtain physically meaningful trends in the levels of cooperation of: (i) subsets of soprano and bass singers of the Imperial College Chamber Choir during a 1 h evensong performance, through trends of synchrony in their respiratory and HRV signals and (ii) three pairs of catheterization laboratory team members (cardiology consultant, cardiology registrar, and physiologist) during a 2 h invasive coronary procedure, through trends of synchrony in their HRV signals.

Related work
In addition to our recently proposed data association measures, IPS and ICoh, there also exist several other synchrony measures. Cross-correlation is a simple measure of linear synchronization between two signals, and hence it cannot effectively deal with the nonlinear coupling behaviour, thus resulting in an undesired low value of correlation coefficient. Phase synchronization index (PSI) proposed in [32] is obtained by considering time-averaged phase difference between two signals, instead of considering the distribution of phase differences as employed in IPS and the proposed N-IPS (see §3.2 for more details). This technique can underestimate synchrony if the distribution of phase differences between two signals has more than one peak, and by averaging over time, phase differences can be cancelled out, resulting in an undesired low value of PSI. Note that the estimation of PSI in IPS, N-IPS and [32] is achieved via the calculation of instantaneous phase of the analytic signal generated using the Hilbert transform. Waveletbased PSI was introduced in [33], whereby instantaneous phase is calculated by convolving each signal with a complex wavelet function and PSI is obtained in the same manner as in IPS and [32]. As a central frequency and a width of the wavelet function must be specified, this approach for estimating PSI is sensitive to phase synchrony only in a certain frequency band.
Synchrony can also be measured by means of information-theoric concepts [34], whereby the mutual information between two signals is defined as the indication of the amount of information of a signal which can be obtained by knowing the other signal and vice versa. The physical meaning or interpretation of synchrony quantified using this approach, however, does not exist.
General synchronization-the existence of a functional relationship between the systems generating the signals of interest-can be characterized by the conditional stability of the driven chaotic oscillator if the equations of the systems are known [35]. For real-world data, however, the model equations are typically unavailable. The non-parametric method of mutual false nearest neighbours [36], which is based on the technique of delay embedding and on conditional neighbours, therefore, has been proposed to characterize general synchronization, yet this technique might produce errors if the signals of interest have more than one predominant time scale [37]. Phase and general synchronization can also be quantified using the concept of recurrence quantification analysis, whereby two signals are deemed: (i) phase synchronized if the distances between the diagonal lines in their respective recurrence plots coincide and (ii) generally synchronized if their recurrence plots are very similar or approximately the same.
All of the described measures are limited to quantifying synchrony between the signals as a whole, and cannot yield TF representations of synchrony. Although such representations can be generated from the IPS algorithm, through the Hilbert transform, we have empirically found that for effective estimation of time-varying synchrony using IPS relatively long sliding windows should be used; hence its time localization is poor. Furthermore, a number of realizations of IPS must be performed for statistical relevance, thus inevitably blurring out TF representations of synchrony. On the other hand, the ISC proposed here generates highly localized TF representations of synchrony and is suitable for the analysis of synchrony in nonlinear and non-stationary multivariate signals.

Algorithm and background
The ISC and N-IPS algorithms proposed in this study are described below.

Intrinsic synchrosqueezing coherence
The proposed ISC is a data-association measure which exhibits precise time and frequency localization in the analysis of nonlinear and non-stationary multivariate signals. This is achieved through the combination of NA-MEMD, FSST and F-MSST algorithms. The NA-MEMD is first employed to decompose a given multivariate signal into a set of physically meaningful narrow-band intrinsic oscillations (IMFs). The FSST is next employed channel-wise to generate multiple univariate multicomponent TF planes with high localization in both time and frequency. The F-MSST is then used to construct multivariate highly localized TF representations. These univariate and multivariate TF representations are subsequently employed to generate TF representations of signal dependence (synchrony) in IMFs; see algorithm 1 for detail of the proposed ISC. The synchrosqueezing coherence index (SCI) ranges between 0 and 1, with 1 indicating perfect synchrony and 0 a non-synchronous state. (v) TF representation of the degree of signal dependence (synchrony) in frequency and time between channels i and j is then determined via the SCI, given by ; ∀f , ∀t (3.1) (vi) Perform 30 realizations of NA-MEMD, repeat steps (i)-(v) for each realization, and average the 30 TF representations of signal dependence between channels i and j, SCI i,j ( f , t), in order to obtain a highly localized TF representation of synchrony.

Nested intrinsic phase synchrony
The IPS was originally proposed in the so-called intrinsic multi-scale analysis framework in [8] and generalizes standard phase synchrony by equipping it with the ability to operate at the intrinsic scale level. It employs NA-MEMD to decompose a given multivariate signal into its narrow-band intrinsic oscillations (IMFs), which makes it possible to quantify the temporal locking of the phase information in IMFs using the standard PSI. We here introduce an extension to the IPS, referred to as N-IPS, for the empirical quantification of physically meaningful and straightforward-to-interpret trends in phase synchrony. The N-IPS first employs the conventional IPS to quantify the intrinsic phase relationship in IMFs, and further decomposes the resulting time series of IPS into multiple scales, whereby only certain scales which contain physically meaningful and straightforward-to-interpret information are then empirically combined (summed), as outlined in algorithm 2. Note that trends in synchrony obtained using the N-IPS algorithm can be negative, because IMFs of synchrony which contain a positive offset in the raw PSI values can be neglected. In such circumstances, the baseline for the trends in synchrony obtained using N-IPS is an imperative, because it is used to judge whether the trends in synchrony are significant or not-they are deemed significant if above the baseline.
(i) Obtain IMFs via NA-MEMD, c n,k , where n = 1, 2, . . . , N, k = 1, 2, . . . , K, and K is the number of IMFs. (ii) For each channel n, calculate the instantaneous phases for the IMFs through the analytic signals generated using the Hilbert transform. (iii) Calculate phase difference between the instantaneous phases for the IMFs of channels i and j, φ i,j (t), where i = 1, 2, . . . , N, j = 1, 2, . . . , N, and i = j . (iv) Phase synchrony between channels i and j is then defined in terms of the deviation from perfect synchrony via the PSI [38], given by where S = − M m=1 p m ln p m is the Shannon entropy of the distribution of phase differences φ i,j (t − W 2 : t + W 2 ) within a window of length W, M is the number of bins within the distribution of phase differences, and p m is the probability of φ i,j (t − W 2 : t + W 2 ) within the mth bin. The maximum entropy S max is given by (v) Obtain IMFs of the phase synchrony between channels i and j, ρ i,j (t), via NA-MEMD, c ρ i,j ,l , l = 1, 2, . . . , L, and L is the number of IMFs of the phase synchrony. (vi) Empirically combine (sum) certain IMFs of the phase synchrony, c ρ i,j ,l , which contain physically meaningful and straightforward interpretation.

Applications
We shall now demonstrate applications of the proposed ISC algorithm in: (i) estimating degrees of synchrony in a linear bivariate signal, (ii) estimating degrees of synchrony in a nonlinear narrowband bivariate signal and (iii) the precise identification of physiological synchrony among choir members. We also demonstrate applications of the N-IPS algorithm in large-scale analyses for the empirical quantification of the levels of cooperation in: (i) choir singing and (ii) performing a surgical procedure. magnitude A oscillating at 0.01 Hz, A 0.01 , given by

Estimating degrees of synchrony in a linear bivariate signal
where A 0.01 = 5, f 1 = 5 Hz, f 2 = 11 Hz, f 3 = 18 Hz, k 1 = 1, k 2 = 0.5, k 3 = 2, k 4 = 0.1 and the sampling frequency f s = 200 Hz. The SNR of the first channel corrupted by additive 1/f noise n 1 was 10 dB, while the SNR of the second channel governed by the 'physiological' 1/f noise n 2 was varied between −10 dB and 20 dB. Figure 1 shows the spectrograms of the clean linear bivariate signal without the low-frequency trend estimated using the standard STFT. In this bivariate signal, there was strong synchrony between the two channels at 5 Hz from 0 s to 1.667 s, and weak synchrony at 18 Hz from 3.334 s to 5 s. The bivariate signal was decomposed using NA-MEMD with 10 adjacent WGN channels. Three IMFs of each of the two channels which governed the frequency band of 4-19 Hz were combined and averaged over 30 realizations of NA-MEMD. Alternatively, the bivariate signal was band-pass filtered using an FIR BPF of order 50 with the passband set to 4-19 Hz, a bandwidth similar to the combined IMFs. A Butterworth IIR BPF of order 3 with the same passband was also applied to the bivariate signal. The different filtered bivariate signals obtained using the NA-MEMD, FIR BPF and IIR BPF were then fed into the W-M/-SST and F-M/-SST algorithms to generate TF representations of signal dependence (synchrony) calculated using equation (3.1) in algorithm 1. Phase synchrony estimation was also employed to estimate phase relationship in the filtered bivariate signal obtained using NA-MEMD-IPS. Figure 2 shows TF representations of degrees of synchrony in the linear bivariate signal estimated using the different combinations of algorithms, where the SNR of the second channel was set to 0 dB, 5 dB and 10 dB (results from the other SNRs are not shown). The W-M/-SST combined with either FIR BPF, IIR BPF or NA-MEMD, produced noticeable spurious synchrony (see the first, third and fifth rows), because CWT introduced mathematical artefacts (noise), and consequently they were localized by the joint frequency estimator in W-MSST. As the 0.01 Hz trend added to the signal was not filtered out The SNR of the second channel governed by 1/f noise was set to 0 dB, 5 dB and 10 dB (a, b and c, respectively). Observed mathematical artefacts produced by CWT and highlighted in the red boxes. Undesired synchrony at 11 Hz was produced by the IIR BPF and F-M/-SST as shown in the yellow box. The ISC algorithm, NA-MEMD and F-M/-SST, performed the best due to no spurious synchrony and the identified weak synchrony at 18 Hz (see the synchrony in the green box).
completely using the FIR BPF and the remaining trend was localized by F-M/-MSST, the SCI values at this frequency consequently dominated the TF representations of degrees of synchrony (second row). Undesired synchrony at f 2 was produced by the combination of IIR BPF and F-M/-SST (see fourth row). Using IPS (last row), the degrees of synchrony estimated at f 1 were not highly localized, spreading around f 1 , and the estimation of weak synchrony at f 3 was not achieved. As desired, the proposed ISC algorithm (NA-MEMD and F-M/-SST, sixth row) did not produce spurious synchrony, yielded desired synchrony at f 1 and achieved the estimation of weak synchrony at f 3 at a high SNR.

Estimating degrees of synchrony in a nonlinear narrow-band bivariate signal
We next evaluated the performance of the seven combinations of algorithms in the task of estimating synchrony in a nonlinear narrow-band bivariate signal extracted from a nonlinear wide-band bivariate AM/FM signal corrupted by additive WGN, given by and where The NA-MEMD was employed to extract a narrow-band bivariate signal of 7-15 Hz (f m ± 4 Hz) from the corrupted wide-band bivariate signal. The corrupted signal was also alternatively band-pass filtered to 7-15 Hz by FIR and IIR BPFs with the same orders as in the previous example. The degrees of synchrony in the filtered narrow-band bivariate signals were then estimated as previously described. Figure 3 shows TF representations of degrees of synchrony in the nonlinear narrow-band bivariate signal estimated using the different combinations of algorithms, where the SNRs of both the channels were set to 0 dB, 5 dB and 10 dB (results from the other SNRs are not shown). Using W-M/-SST combined with either FIR BPF, IIR BPF or NA-MEMD, there exist noticeable spurious synchronies outside the band of interest (see the first, third and fifth rows). Synchrony in the band of interest was not exhibited using the combination of FIR BPF and F-M/-SST (second row), and was poorly estimated using the combination of IIR BPF and F-M/-SST (see fourth row). The proposed ISC algorithm (NA-MEMD and F-M/-SST, sixth row) did not produce spurious synchrony outside the band of interest, and exhibited the desired synchrony at high SNRs, outperforming IPS (last row) at 5 dB and 10 dB SNRs. Figure 4 shows average SNRs of the TF representations of synchrony produced by the seven combinations. The degrees of synchrony estimated at the frequencies f 1 and f 3 in the linear bivariate signal and those in the band of interest in the nonlinear bivariate signal were deemed 'signals', and the degrees of synchrony outside these were deemed 'noise'. Note that the SNR of the FIR and F-M/-SST in the first task could not be computed, because no 'signal' was present in its TF representations. The Z-test at a significance level of 0.01 was performed, in order to reveal statistical differences in the SNRs between the proposed ISC algorithm and the second best combination in both the tasks. Average SNRs of the proposed ISC algorithm (3.593 dB and 13.66 dB) were 3.625 dB and 3.2 dB significantly higher (p-value = 0.0001 and 0.0004, respectively) than those (−0.032 dB and 10.46 dB) of the second best combination, which is the IIR BPF and F-M/-SST, in estimating degrees of synchrony in, respectively, the synthetic linear and nonlinear bivariate signals. The ISC algorithm outperformed the others in both of the tasks, because essentially it did not produce spurious synchrony, making it the most reliable among them. It should also be noted that: (i) both FIR and IIR filters generally introduce several problems such as group delay and nonlinear phase response; (ii) these filters are not data-driven and not designed specifically for nonlinear and non-stationary signals as opposed to NA-MEMD; and (iii) CWT produces mathematical artefacts which degrade the performance of W-MSST.

Highly localized heart rate variability-and respiration-based choral synchrony analysis 4.3.1. Motivation
Our recent work [1,39] has shown that during choral singing, HRV and respiration of a subset of five singers synchronized, with an increase from the rehearsal to the performance on the average, and the physiological synchrony in both situations, as a function of frequency, dominated the LF band. This exemplifies the better coordination among the singers under the situation that was more stressful. That study, however, examined the average level of cooperation by averaging the PSI and coherence of the respiratory and HRV signals of three different voices (one tenor, one soprano and three bass singers) over time. This did not allow for the quantification of the level of cooperation of specific voices and at different time instants where the tempo of the musical score varied. In addition, a highly localized  TF representation of synchrony in the respiration and HRV signals for thorough investigation of the frequential and temporal dynamics of the physiological synchrony was not available.
To this end, we here employ: (i) IPS to reveal the degrees of time-varying synchrony in HRV and respiration among only the three bass singers, since as the same voice they performed long or short inhalation or exhalation almost exactly at the same time over the course of the 4 min and 40 s rehearsal and performance; and (ii) ISC to precisely quantify synchrony in the physiological signals in both time and frequency. Furthermore, the time-varying synchrony estimated using IPS was used to verify patterns of temporal dynamics and time localization of the TF representations of the physiological synchrony obtained using the proposed ISC algorithm.   performance at Union Chapel, London, UK. Respiration of each participant was recorded using a custommade respiration belt placed around the chest. For all participants, the ECG was recorded with three electrodes placed on the skin, just below the collar bone. The respiration belt and the electrodes were connected to an eight-channel portable biosignal data logger powered by a rechargeable coin cell battery. The data logger sampled the signals at 1 kHz and saved the respiratory and ECG data onto a micro-SD card. The data logger also recorded timestamps onto the micro-SD card in order to guarantee the synchronization between the devices. The respiratory signals were downsampled to 4 Hz, and the HRV was estimated from the ECG data by band-pass filtering between 8 Hz and 30 Hz, and the subsequent R-peak detection to obtain the RR-interval (i.e. HRV) time series with a sampling frequency of 4 Hz [40]. The respiratory (or HRV) signals of the three bass singers during both the rehearsal and the performance were used to form three-channel data which was decomposed using NA-MEMD with 10 adjacent WGN channels. The IMFs with indices 3-6 produced by the NA-MEMD of the threechannel multivariate HRV (or respiratory) signal of the bass singers contained the physically meaningful frequency range 0.04 Hz to 0.4 Hz, that is, exactly the LF/HF frequency band of HRV. The full band of interest in both the HRV and respiration data was produced by summing up the corresponding IMFs, in order to obtain the desired scale in data. The time-varying PSI values among the singers were obtained by averaging PSI values calculated from three combined-IMF pairs of the three data channels (pair 1: first and second bass singers, pair 2: first and third bass singers, pair 3: second and third bass singers), whereby the PSI values for each pair were computed using 20 s sliding windows with 18 s overlap (2 s increment) for the respiratory signals, and 40 s sliding windows with 36 s overlap (4 s increment) for the HRV signals. The PSI indices between the combined IMFs of the noise channels were also estimated in order to provide the PSI of random signals as a baseline.

Data acquisition and synchrony analysis
The ISC was performed on each of the same three combined-IMF pairs of the three data channels of the respiratory (or HRV) signals to obtain three TF representations of synchrony in each of the signals for all the pairs. The three TF representations were next averaged in order to obtain a TF representation of synchrony in the respiratory (or HRV) signals among the bass singers. Figures 5a and 6a show the timings of long or short exhalation of the bass singers when they performed in unison, and when they remained silent or inhaled. Figures 5b and 6b show degrees of synchrony in, respectively, the respiratory and HRV signals during the rehearsal and the performance, estimated using IPS with 30 realizations of NA-MEMD. Observe that the respiratory synchrony in both situations (figure 5b) reached their peaks approximately during a series of long exhalation. Figures 5e-f and 6e-f show TF representations of, respectively, respiratory and HRV synchrony during the rehearsal and the performance estimated using the proposed STFT-based ISC. In both situations, the respiratory synchrony (figure 5e-f ) was highly localized at 0.1 Hz due to a series of long exhalations they had to perform in unison. This exemplifies that the singers were constrained by the musical score to breathe in unison at  a very slow rate of six breaths per minute (breaths per minute = 60 × breathing frequency). This rate is remarkably slower than the normal breathing rate in adults, which varies between 12 and 18 breaths per minute [41]. The HRV synchrony (figure 6e-f ) was dominant in the LF band of HRV, 0.04-0.15 Hz. This conforms with our recent finding in [1], where ICoh revealed a large proportion of coherence in the HRV signals in the LF band.

Results of the analysis
The sums of SCI values across frequencies as a function of time of the respiratory and HRV signals during both the rehearsal and the performance are, respectively, shown in figures 5c-d and 6c-d (see thin dotted lines). Observe that the sums of SCI values of both the signals varied dramatically during both situations, indicating that the proposed ISC algorithm effectively yielded highly localized timevarying respiratory and HRV synchrony. To perform straightforward verification of patterns of temporal dynamics and time localization of the relatively 'faster' temporal dynamics of the sums of SCI values against the 'slow' temporal dynamics of the corresponding PSIs shown in figures 5b and 6b, the sums of SCI values of both the signals during both situations were used to form four-channel data which was decomposed using MEMD. We empirically found that the multivariate IMFs with indices 1-7 of the four-channel data could be deemed 'fast oscillatory components', while the multivariate IMFs with indices 8-11 'slow oscillatory components', and that the combinations of the multivariate IMFs in the latter set yielded the reconstructions of the sums of SCI values which exhibited: (i) 'slow' temporal dynamics similar to the 'slow' temporal dynamics of the corresponding PSIs for both the signals and  We also examined the degrees of synchrony in the bass singers' respiratory and HRV signals during both the situations using CWT-based ISC, which combines the conventional CWT-based SST and MSST algorithms, as shown in figure 7a-d. As the joint frequency estimator in the CWT-based MSST is sensitive to noise [29] and CWT produced mathematical artefacts (additional noise), the TF representations produced by the CWT-based ISC exhibited notable spurious synchrony at several frequencies. the dynamics of the level of cooperation in response to events during the course of long cooperative tasks. Therefore, quantifying trends in synchrony, i.e. the level of cooperation, in choir singers during an evensong performance, which typically takes place over a long period of time and comprises several events (singing and pausing), is the subject of this study.
Owing to the nature of an evensong performance, where all the events have different lengths, time series of synchrony in singers' respiratory and HRV signals typically exhibit nonlinear behaviour (i.e. deformed wave profiles). Such time series of synchrony can be empirically decomposed using a second NA-MEMD into a set of physically meaningful nonlinear components (IMFs) of the synchrony without prior knowledge of the frequency ranges (periods) of trends. In this way, we obtain physically meaningful and straightforward interpretation in response to the events, as opposed to low-or bandpass filtering, where cut-off frequencies are rigid, must be determined and specified beforehand, and are person-specific (unlike the proposed approach), which causes physically meaningful synchrony events to be missed. To this end, the proposed N-IPS algorithm was employed to empirically obtain physically meaningful and straightforward-to-interpret trends in the degrees of synchronization in: (i) respiratory signals among a subset of three soprano singers and a subset of three bass singers; and (ii) HRV signals among a subset of five soprano singers and a subset of four bass singers, of the 20-member Imperial College Chamber Choir during a 1 h evensong performance at St Martin-in-the-Fields church at Charing Cross, London, UK.

Data acquisition and synchrony analysis
The data were acquired in the same manner as in §4.3, except for the group and number of singers and the duration of the recording; for this analysis the respiratory and HRV signals were recorded from 10 min before the performance, 1 h of the performance and 10 min after the performance.
The respiratory (or HRV) signals of the soprano and bass singers were used to form multichannel data which was decomposed using NA-MEMD with 10 adjacent WGN channels. The IMFs which were produced by the NA-MEMD with indices 3-6 of the multivariate HRV (or respiratory) signal of the singers and contained the physically meaningful frequency range 0.04 Hz to 0.4 Hz (the LF/HF frequency band of HRV) were then summed up to obtain the desired scale in data. The PSI estimation was next performed in the same manner as in §4.3. The time-varying PSI indices between the combined IMFs of the signal and noise channels (baseline) were next used to form multichannel data which were further decomposed using NA-MEMD with 10 adjacent WGN channels. Certain IMFs of the multivariate synchrony (synchrony-IMFs) produced by the NA-MEMD (IMFs 5-7 for the respiratory synchrony, 5-6 for the HRV synchrony) were then combined, because the combinations of the synchrony-IMFs admit capturing of physically meaningful and interpretable trends in the respiratory and HRV synchrony; that is, temporal dynamics of the combinations of the synchrony-IMFs closely resemble the timings of the songs and the pauses, where increases (during songs) and decreases (during pauses) in trends of the physiological synchrony, i.e. the level of cooperation, were physiologically expected to be observed. The other IMFs, on the other hand, were neglected due to no clear physical meaning or straightforward interpretation; their temporal dynamics did not resemble the timings of the songs and the pauses, and distinct increases and decreases in trends of the synchrony during the songs and the pauses, respectively, cannot be observed. Observe that all the 'raw' synchrony obtained using the standard IPS algorithm (solid thin lines) varied dramatically during the course of events (singing and pausing), thus being less amenable to physical interpretation. Trends in the synchrony obtained using the proposed N-IPS algorithm, on the other hand, exhibited smooth variations of synchrony during the events, and smooth transitions between them. This remarkably offers more physically meaningful and interpretable results; that is, the levels of cooperation markedly increased during most of the songs (depicted by the upward pointing arrows), and decreased during the long pauses (the downward pointing arrows). Variations in the degrees of synchrony (levels of cooperation) were because the breathing rhythms, and by virtue of RSA the cardiac activities too, of the singers were modulated by the pieces of music performed, thus exhibiting stronger dynamic coupling reflected in increases in the synchrony of their physiological responses. During the pauses, however, their breathing rhythms were not dictated by any piece of music, thus resulting in decreases in the synchrony. It should be noted that the HRV and respiratory signals were recorded from an amateur choir. Their breathing rhythms, therefore, may not have been perfectly synchronized. This is reflected by the transitional periods of the trends in synchrony between the songs and pauses occasionally not being perfectly aligned with the exact transitional timings between the songs and pauses. It also must be emphasized that we here focused on the empirical quantification of physically meaningful and straightforward-to-interpret trends in synchrony, which is a first step towards prediction.

4.5.
Large-scale heart rate variability-based surgical synchrony analysis 4.5

.1. Motivation
Invasive cardiology procedures are performed on a daily basis in dedicated centres. Every patient with chest pain and diagnosed myocardial infraction goes for primary coronary angioplasty, which is an established management giving the best outcomes in terms of recovery and survival. Real-life HRV recording and simultaneous acquisition of its changes between team members involved in the procedure will give an objective tool for assessment of the team performance and efficiency. Empirical evidence supports the view that the quality of surgical performance depends not only the technical skills of the surgical teams, but also on good collaboration and effective teamwork [42]. Failures in non-technical skills in the operating room have been frequently implicated in more frequent errors, adverse events in surgical patients and longer operative times [43,44]. Non-technical skills also have a direct impact on the technical performance of surgical teams [45].
Multidisciplinary interventional cardiology catheterization laboratory (Cath Lab) teams are composed of an interventional cardiology consultant, a cardiology registrar, a nurse, a physiologist and a radiographer, each with a different background in terms of training and competencies. These team members are expected to work and function in full synchrony to achieve optimal outcomes. Close cooperation between physicians performing the procedure (cardiology consultant and cardiology registrar) and physiologist outside the suite who record and monitor pressures during coronary intervention is essential.    Degrees of raw synchrony (thin solid lines) and the corresponding baselines (thin broken lines) were estimated using IPS. Trends of the raw synchrony (thick solid lines) and the corresponding baselines (thick broken lines) were estimated using N-IPS. Time instants when the singers sang during the nine songs (S1-S9) are depicted by the shaded areas. Observe increases in the synchrony of singers' physiological responses during most of the songs depicted by the upward pointing arrows. These are due to their breathing rhythms and HRV being modulated by the demands of pieces of music.
The case presented in the manuscript is part of a pilot study, which aims to show that a team, whose members work together on a daily basis and have developed mutual understanding and trust, can respond quicker and more efficiently to unexpected situations. Common experience of team members has been found to play a critical role in good collaboration in the operating theatre and is increasingly recognized as an important mechanism for enhancing the safety of delivered healthcare services. It has been shown that team familiarity has a threefold greater impact on the duration of the surgical procedure than the experience of the main surgeon [46].

Data acquisition and synchrony analysis
ECG signals were recorded from a Cath Lab team starting 10 min before, during and after a percutaneous coronary intervention procedure, which was 1 h and 56 min long, and was carried out in the Cath Lab at Hammersmith Hospital, Imperial College Healthcare NHS Trust, London. During the procedure, team members were closely observed by an independent observer, who noted times and types of any glitches that occurred. Based on the published literature, we have divided possible glitches into 13 categories [47]. The electrode placement for ECG recording, data acquisition, HRV estimation and PSI estimation were carried out in the same manner as in §4. 3  trends in synchrony in the surgeons' HRV signals in order to cater for: (i) the nonlinear behaviour of the synchrony measure and (ii) the lack of prior knowledge of the frequency ranges (periods) of trends, thus providing physically meaningful and straightforward interpretation in response to glitches. The timevarying PSIs between IMFs of the HRVs and noise channels (baseline) were used to form multichannel data, which were further decomposed using NA-MEMD with 10 adjacent WGN channels. The IMFs produced by the NA-MEMD with indices 10-12 of the multivariate PSI values were combined to produce trends in the synchrony, because these combinations best represented changes in trends of the synchrony in response to gliches. Figure 9 shows trends in synchrony in HRV data during the recorded procedure, which exhibited cooperation in responses to glitches between the following pairs: (i) cardiology consultant-cardiology registrar, (ii) cardiology consultant-physiologist and (iii) cardiology registrar-physiologist. Glitches, which were recorded during the procedure included the following groups: (i) group 3-distractionsanything causing distraction from task (phone calls/bleeps, loud music, alarms), and (ii) group 5-equipment design-issues arising from equipment design that would not otherwise be corrected with training or maintenance (e.g. compatibility problems with different implants or wires, or equipment failure). Figure 9 shows that physiological responses of pairs of professionals responded to the recorded glitches with increased synchrony in HRV. Observe that the first glitch at the 27th min caused immediate negative deflections (increases) in trends of HRV synchrony between (i) the cardiology consultant and the physiologist (thin solid line) and (ii) the cardiology consultant and the cardiology registrar (thick broken line). Also, the third glitch at the 72nd min caused another immediate negative deflection (increase) in a trend of HRV synchrony between the cardiology registrar and the physiologist (thick solid line). These findings suggest these glitches were impeding the procedure and may affect procedural duration as well as outcome.

Results of the analysis
Simultaneous acquisition of HRV of team members involved allows us to objectively measure the involvement of each specialist during the case and observe their response to different interruptions and unexpected problems during the invasive cardiology procedures.

Conclusion
We have introduced a data-association measure which exhibits high frequency and time localization, termed ISC, for the analysis of nonlinear and non-stationary multivariate signals. This has been achieved through the combination of multivariate signal analysis using NA-MEMD and the generation of highly localized TF univariate and multivariate representations of the multivariate intrinsic oscillations using the STFT-based SST and MSST algorithms. Such a measure has enabled precise identification of physiological dependence in both frequency and time. The performance of the proposed ISC algorithm has been evaluated against other combinations of algorithms in the tasks of estimating synchrony in linear and nonlinear bivariate signals. The proposed ISC algorithm has been shown to outperform the other combinations of algorithms and to exhibit significantly higher average SNRs of TF representations of signal dependence. We have demonstrated an application of the ISC algorithm to the quantification of inter-channel dependence in respiratory and HRV signals. The proposed algorithm has exhibited highly localized synchrony in the respiratory and HRV signals among a subset of three bass singers of a choir during both a rehearsal and a performance. The proposed STFT-based ISC algorithm has been shown to outperform CWT-based ISC, which was built upon CWT-based SST and MSST, in the localization of synchronized frequencies in respiratory and HRV signals.
In addition to the ISC algorithm, we have proposed an extension of IPS, referred to as N-IPS, as a meaningful and straightforward-to-interpret data-association metric for trends in the level of cooperation. This is achieved by first employing the standard IPS to quantify intrinsic phase relationship between data channels, and then further decomposing time series of the multivariate degrees of phase synchrony into multiple scales, whereby certain intrinsic scales which contain physically meaningful and straightforward interpretation are then combined. This algorithm allows for empirical quantification of physically meaningful and straightforward-to-interpret trends in phase synchrony. Two applications of the N-IPS algorithm to the empirical quantification of physically meaningful trends in the level of cooperation through the empirical estimation of trends in synchrony in respiratory and HRV signals have been demonstrated. The N-IPS algorithm has effectively quantified physically meaningful increases and decreases in trends in levels of cooperation among subsets of soprano and bass singers of a choir during a 1 h evensong performance. It has also revealed significant increases in trends in levels of cooperation between pairs of cardiologists and physiologists when certain types of glitches occurred during a surgery. This is a first attempt to empirical quantification of intrinsic, physically meaningful and interpretable trends in the level of cooperation of long collaborative tasks.
Future work will focus on the development of Panorama-based SST and MSST algorithms in order to enhance the generation of TF representations [48]. The result of the surgical synchrony analysis is preliminary, and the N-IPS algorithm will be employed in the assessment of HRV changes in response to glitches within different surgical teams, to show how team familiarity has an impact on patient outcomes, i.e. complications, reintervention rate, length of surgical or endovascular procedures and length of post-operative hospitalization. A growing body of the literature suggests that shared knowledge and understanding of each other's roles and objectives during the invasive and high-risk procedures in turn facilitates team cooperation and coordination. For a team to function efficiently, its members should share a 'mental model' of the team's tasks, objectives, means and environment [49][50][51]. We are hoping that our prospective research will help us to facilitate a better training programme across UK hospitals in order to deliver better care to the patients.
Ethics. Ethics approval was obtained from the Joint Research Office at Imperial College London, reference ICREC_12_1_1. Written consents of the subjects who participated in these studies were obtained, and they were not paid for the participation.
Data accessibility. Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.80cv0 [52]. It can also be accessed from http://www.commsp.ee.ic.ac.uk/∼mandic/research/emd.htm. aspect of an MEMD approach, via either standard MEMD or NA-MEMD, is that even if mode mixing occurs it is likely to happen simultaneously in every channel. In this way, multicomponent comparisons at the IMF level are matched and physically meaningful. It was also suggested in [8] that the number of WGN channels should be as large as possible-computation time permitting-in order to reduce the degree of overlap between IMF spectra, that is, impose a more pronounced dyadic filterbank structure within the IMFs. See algorithm 3 for details of the NA-MEMD algorithm.

A.2. Short-time Fourier transform-based synchrosqueezing transform
The synchrosqueezing transform (SST) [27] is a post-processing technique originally applied to the CWT in order to generate highly localized TF representations of nonlinear and non-stationary signals. It has been extended to the STFT setting in [28], and is referred to as STFT-based SST (FSST). It performs STFT, instead of CWT, on a given univariate signal x(t) to obtain STFT coefficients U(η, t), given by where g is a sliding window and τ is the time lag. The STFT coefficients, U(η, t), are next relocated to the instantaneous frequency, F, according to the map (η, t) → (F(η, t)), where F is defined by .
The STFT-based SST coefficient is then given by A summary of the FSST is outlined in algorithm 4. (i) Given a univariate signal x(t), perform STFT on the signal in order to obtain the STFT coefficients U(η, t), given by equation A 2. (ii) Estimate the instantaneous frequency F(η, t) of the STFT coefficients U(η, t), given by equation A 3. (iii) Combine the STFT coefficients in order to obtain the STFT-based SST coefficients T( f , t), given by equation A 4.

A.3. Short-time Fourier transform-based multivariate synchrosqueezing transform
The multivariate synchrosqueezing transform (MSST) [29] is an extension of the CWT-based SST to multivariate cases. It was proposed in order to identify oscillations common to multiple data channels. It first partitions the TF domain into M frequency bands {f m } m=1,2,...,M . This makes it possible to identify a set of matched monocomponent signals from a given multivariate signal. The instantaneous amplitudes and frequencies present within those frequency bands can then be determined. The implementation of F-MSST is straightforward; FSST is employed instead of SST to obtain SST coefficients of each channel. For an N-variate signal x(t) with the corresponding SST coefficients of each channel n, T n ( f , t) obtained using FSST, and a set of frequency bands {f m } m=1,2,...,M , the instantaneous frequency F( f m , t) for each channel n and frequency band m is given by To estimate the multivariate instantaneous frequency F( f m , t) for a given frequency band m, the instantaneous frequencies across the N channels are then combined using

A.4. Areal density
Trends in synchrony obtained using the proposed N-IPS algorithm offer a meaningful and straightforward interpretation in signal dependence. This can also be further simplified by calculating the areal density (AD)-an area between a trend in synchrony and a corresponding baseline per unit of time. Figure 10a shows areal densities of trends in synchrony in the respiratory and HRV signals of the soprano and bass singers considered in §4.4. Observe that the areal densities increased during the songs (S x , where x represents the song number) and decreased between two consecutive songs (S x -S y , where x and y represent the two consecutive song numbers), except for only two circumstances: (i) during Songs 1 and 2 (S1-S2), where areal densities of the trend in synchrony in the HRV signals of the soprano and bass singers remained relatively high; and (ii) during Songs 6 and 7 (S6-S7), where areal densities of the trends in synchrony in the sopranos' HRV signals and the basses' respiratory signals increased from the previous song (Song 6, S6). Note that 'S7 (pause)' represents the later period of Song 7 where several pauses occurred, and consequently the trends and the areal densities of the trends in synchrony decreased. For rigour, we also calculated average areal densities of trends in synchrony in both signals of the sopranos and basses, as shown in figure 10b. The average areal densities exhibited distinct increases during the songs and decreases between the songs.