rsos.royalsocietypublishing.org
Research
Cite this article: Hemakom A, Powezka K,
Goverdovsky V, Jaer U, Mandic DP. 2017
Quantifying team cooperation through
intrinsic multi-scale measures:respiratory and
cardiac synchronization in choir singers and
surgical teams. R. Soc. open sci. 4: 170853.
http://dx.doi.org/10.1098/rsos.170853
Received: 6 July 2017
Accepted: 7 November 2017
Subject Category:
Engineering
Subject Areas:
electrical engineering/biomedical engineering
Keywords:
multivariate empirical mode decomposition,
multivariate synchrosqueezing transform,
intrinsic multi-scale analysis, coherence,
respiration, heart rate variability
Author for correspondence:
Danilo P. Mandic
e-mail:d.mandic@imperial.ac.uk
Quantifying team
cooperation through
intrinsic multi-scale
measures: respiratory and
cardiac synchronization
in choir singers and
surgical teams
Apit Hemakom
1
, Katarzyna Powezka
2
, Valentin
Goverdovsky
1
, Usman Jaer
2
and Danilo P. Mandic
1
1
Department of Electrical and Electronic Engineering, and
2
Department of Vascular
Surgery, Imperial College London, London SW7 2AZ, UK
AH, 0000-0003-0621-5240
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
2017 The Authors. Published by the Royal Society under the terms of the Creative Commons
Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted
use, providedtheoriginalauthorand source are credited.
2
rsos.royalsocietypublishing.org R. Soc. open sci. 4: 170853
................................................
frequency (ISC), and (ii) large-scale analysis for the empirical quantification of physically meaningful
trends in synchrony (N-IPS).
1. Introduction
Cooperative human activities require high degrees of mental and physical synchronization among
multiple participants, to the extent that synchrony underpins performance level in activities, such as
choir singing, playing music in an ensemble, rowing, flying an aeroplane with a co-pilot or performing
surgical procedures. When it comes to quantifying the degree of synchronization among participants,
synchrony in physiological responses has been reported in respiration and heart rate variability (HRV)
among the choir members [1,2].
Synchrony in respiration among the choral singers is a result of their breathing rhythm being dictated
by the tempo and demands of a musical score; that is, they typically perform the short inhalation and
long exhalation in unison. In addition to the voluntarily controlled breathing, the respiration is also
involuntarily controlled by the autonomic nervous system, which comprises the sympathetic (SNS) and
parasympathetic (PNS) nervous subsystems. The SNS also accelerates other functions, such as the arterial
blood pressure and heart rate [3–5], by dilating bronchioles in the lungs, and by regulating neuronal
and hormonal responses to stimulate the body. The PNS, on the other hand, slows down physiological
functions when the body is at rest.
The interplay between the SNS and the PNS, among other factors, manifests itself in variations of the
timing of the cardiac cycle—HRV—in response to both external and internal factors. Changes in HRV
are commonly evaluated in two frequency bands: (i) the low-frequency (LF) band, 0.04–0.15 Hz, which is
linked to the interaction of the SNS and PNS, and (ii) the high-frequency (HF) band, 0.15–0.4 Hz, which
primarily reflects the activity of the PNS [6,7]. In addition, it is well understood that breathing modulates
HRV via a phenomenon referred to as the respiratory sinus arrhythmia (RSA), whereby the heart rate
accelerates during inspiration and decelerates during expiration. The RSA is usually attributed to the
activity of the PNS, so that the HF component of HRV is dominated by the changes in heart rate induced
by breathing.
In an attempt to quantify the degrees of synchronization in the singers’ physiological responses
(respiration and HRV), a quantitative measure of the level of cooperation has been recently proposed
in [1]. This is achieved via the assessment of phase relationship between multiple physiological
responses, based on the intrinsic phase synchrony (IPS) and intrinsic coherence (ICoh) measures
proposed in our recent work [1,8] under a framework referred to as intrinsic multi-scale analysis.The
algorithms under this framework are capable of quantifying intra- and inter-component dependence of
a complex system, such as multiple synchronies and causalities. The IPS and ICoh are implemented
through a combination of the novel data-driven multivariate signal decomposition algorithm called
noise-assisted multivariate empirical mode decomposition (NA-MEMD) (see [9,10] and appendix A.1 for more
details) and two standard data-association measures, phase synchrony and coherence. As desired, the
IPS accounts only for phase information between dependent signals, regardless of the differences in
the magnitude of intrinsic oscillations between data channels obtained using NA-MEMD. The time and
frequency aspects of synchrony, however, are not highly localized using this algorithm. Conversely, using
ICoh, the degree of signal dependence can be quantified as a function of frequency, at a loss of the time
dimension, because ICoh is computed over the whole dataset.
Despite several limitations, when used in conjunction with IPS and ICoh, NA-MEMD is an efficient
multichannel data processing method. It is an extension of the empirical mode decomposition (EMD)
[11,12] to multivariate cases, with the assistance of noise in order to enforce the dyadic filterbank
property. The EMD is essentially an adaptive, data-driven method for the analysis of nonlinear and
non-stationary univariate time series. It employs the so-called sifting process to decompose a given
signal into its multiple physically meaningful narrow-band amplitude/frequency modulated (AM/FM)
components, which are referred to as intrinsic mode functions (IMFs) and are used as bases for signal
representation.
Unlike conventional projection-based time–frequency (TF) algorithms, such as the short-time Fourier
transform (STFT) and the discrete wavelet transform (DWT), the IMFs—the adaptive basis functions
within EMD—enable physical interpretation for nonlinear signals, because IMFs are theoretically
designated as intrinsic oscillations with intrawave amplitude and frequency modulation, with a physical
meaning that the amplitudes and frequencies of these basis functions can intrinsically and nonlinearly vary
3
rsos.royalsocietypublishing.org R. Soc. open sci. 4: 170853
................................................
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 non-
stationary 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–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.
4
rsos.royalsocietypublishing.org R. Soc. open sci. 4: 170853
................................................
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.
2. 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. Wavelet-
based 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.
3. Algorithm and background
The ISC and N-IPS algorithms proposed in this study are described below.
5
rsos.royalsocietypublishing.org R. Soc. open sci. 4: 170853
................................................
3.1. 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.
Algorithm 1: Intrinsic synchrosqueezing coherence.
Input: A multivariate signal x(t), x(t) = [x
1
(t), x
2
(t), ..., x
N
(t)]
T
.
(i) Obtain IMFs via NA-MEMD, c
n,k
,wheren = 1, 2, ..., N, k = 1, 2, ..., K,andK is the number of
IMFs.
(ii) For each channel n, compute the Fourier spectra of IMFs and combine (sum) the IMFs governing
the frequency band of interest, IMFs
n
.
(iii) Apply FSST channel-wise to the combined IMFs, IMFs
n
, to generate N univariate TF
representations, T
n
( f , t), where n = 1, 2, ..., N.
(iv) Apply F-MSST to the combined IMFs between channels i and j, IMFs
i
and IMFs
j
,where
i = 1, 2, ..., N, j = 1, 2, ..., N,andi = j, to generate a multivariate TF representation of the two
channels, T
i,j
( f , t).
(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
SCI
i,j
( f , t) =
|T
i
( f, t)|·|T
j
( f, t)|
|T
i,j
( f, t)|
max
f ,t
|T
i
( f , t)|·|T
j
( f , t)|
|T
i,j
( f , t)|
; ∀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.
3.2. 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.
