A new approach for the quantification of synchrony of multivariate non-stationary psychophysiological variables during emotion eliciting stimuli

Emotion eliciting situations are accompanied by changes of multiple variables associated with subjective, physiological and behavioral responses. The quantification of the overall simultaneous synchrony of psychophysiological reactions plays a major role in emotion theories and has received increased attention in recent years. From a psychometric perspective, the reactions represent multivariate non-stationary intra-individual time series. In this paper, a new time-frequency based latent variable approach for the quantification of the synchrony of the responses is presented. The approach is applied to empirical data, collected during an emotion eliciting situation. The results are compared with a complementary inter-individual approach of Hsieh et al. (2011). Finally, the proposed approach is discussed in the context of emotion theories, and possible future applications and limitations are provided.


INTRODUCTION
Researchers agree that emotion eliciting situations are accompanied by changes of multiple variables associated with subjective, physiological and behavioral responses. The reasons for a possible coupling of the response variables is a topic of ongoing discussion in various emotion theories. There is no uniform terminology to describe the simultaneous changes of the response variables. The term "coherence" is frequently used (e.g., Rosenberg and Ekman, 1994;Reisenzein, 2000;Mauss et al., 2005;Sze et al., 2010;Herring et al., 2011;Hsieh et al., 2011;Dan-Glauser and Gross, 2013) to describe the simultaneity of changes in the response variables. Further terms that are used in the research community to describe the interrelation of the responses are "synchronization," "organization of response systems," or "concordance," Throughout this paper, the terms synchronization and synchrony are interchangeably used to describe the simultaneous changes of response variables. This article introduces a new approach for the quantification of the synchrony of the response variables that is able to account for non-stationarity. A signal is nonstationary if its mean and covariance function are time-varying (Brillinger, 2001). First, different assumptions regarding the functionality of the synchrony concept in different emotion theories are explained. Subsequently, the new approach, which basically consists of two steps, is introduced. In a first step, time-frequency based bivariate coherence measures are derived (Muma et al., 2010). In a second step, these measures are used in an state space modeling approach to obtain an overall synchrony measure of the simultaneous activation of psychophysiological responses. The approach is then applied to empirical data collected during an emotion eliciting situation and compared with a complementary approach of Hsieh et al. (2011). Finally, the proposed approach is discussed in the context of emotion theories, and possible future applications and limitations are provided.

ON THE FUNCTIONALITY OF SYNCHRONY OF RESPONSES IN EMOTION THEORIES
Emotion theories make different assumptions regarding the functionality of a synchrony of response variables: Basic emotion theories state that different emotions have distinct and coordinated patterns of physiological responses. According to basic emotion theories, the specific psychophysiological response variables are activated simultaneously during an emotional experience but are less associated with each other during rest (Tomkins, 1962;Izard, 1977;Ekman, 1992;Levenson, 1994). A synchronized response is desirable as it prepares the body for an adequate reaction to a stimulus and leads to an appropriate reaction to environmental demands (Tomkins, 1962;Izard, 1977;Ekman, 1992;Levenson, 1994). Different emotional responses are organized by central mechanisms in the brain, such as the amygdala (e.g., Whalen et al., 1998;LeDoux, 2000;Murphy et al., 2003), orbitofrontal cortex (e.g., Hornak, 2003;Murphy et al., 2003;Goodkind et al., 2012), the insular (e.g., MacLean, 1990;Damasio et al., 2000;Murphy et al., 2003), and other brain regions (e.g., Ekman, 1992;Panksepp, 2008;Izard et al., 2011). The synchronized specific responses in time and intensity have been interpreted as evidence for the existence of a causal mechanism (e.g., Kettunen et al., 1998;Levenson, 2003). Further, emotions can also activate so called 'affect programs' that include behavioral and physiological changes, that might be similar for different individuals (Tomkins, 1962;Ekman and Cordaro, 2011). For example, fear, anger, or amusement may be accompanied by a specific reaction pattern in terms of subjective experience, physiology, and behavior (Ekman, 1992;Murphy et al., 2003;Levenson, 2011;Panksepp and Watt, 2011).
A different viewpoint is taken by so-called dimensional approaches that do not classify the emotional experience in distinctive categories, such as anger or fear. Instead, these approaches discriminate between different emotional states by introducing a two-dimensional space that is spanned by valence (pleasure/displeasure) and arousal (activated/deactivated) as sufficient means to discriminate between different emotional states (Russell, 1980(Russell, , 2003Barrett and Russell, 1998;Barrett, 2006). Higher dimensional spaces have been proposed by e.g., Bradley and Lang (1994) or Fontaine et al. (2007). The assumption that an event causes an emotion and the emotion causes a synchronized, specific change in cognition, behavior, and physiological reaction is criticized by Russell (2003). Instead, in the conceptual framework of the core affect, the hypothesis is made that that emotions do not have a common cause. From this perspective, a synchrony between the responses is not required (see also Barrett, 2006;Russell, 2009). On grounds of genetic differences, personal experiences (e.g., Lykken and Tellegen, 1996), responsiveness to stimuli, attributions, and other factors, individuals can respond differently to the same emotion eliciting situation (Russell, 2003). Hormonal changes, endocrine dysfunction, illness, satiety and diurnal rhythms can internally influence the emotional response (Russell, 2003), which results in very specific patterns of psychophysiological variables on the intra-and inter-individual level.
A third body of research relies on the concept of appraisal. Appraisal indicates an evaluation of the situation regarding its personal significance (Lazarus, 1991) and therefore, in the same situation there can be a great variation between the emotional state of individuals (Scherer et al., 2006;Kuppens et al., 2009). Still, a synchronous response during an emotional episode is expected (Scherer, 2009). Some of the appraisal theories view the appraisal of a situation, and not the event itself (as projected by the basic emotion theory; Russell, 2003), as a causal mechanism responsible for the elicitation of an emotion (e.g., Schachter and Singer, 1962;LeDoux, 1989;Roseman et al., 1990;Scherer, 1993). Lewis (2005), on the other hand, assumes a recursive, complex relationship between subsystems of the nervous system and hence, not a linear relationship of one causal mechanism followed by a cascade of responses. Grandjean et al. (2008) consider various feedback loops between the synchronized response of the peripheral-, motivational-, monitor-, cognitive-, and motorsystem to be responsible for the conscious awareness of an emotion (e.g., Scherer, 1984Scherer, , 1987Scherer, , 2009Fries, 2005). The synchrony of multiple response variables, however, is considered to be necessary for a conscious emotional experience (Grandjean et al., 2008). Also, several appraisal approaches support discrete categories of emotion (e.g., Oatley and Johnson-laird, 1987;Roseman et al., 1990) while others approve the dimensional perspective (e.g., Scherer et al., 2006;Kuppens et al., 2012;Wu et al., 2012). The above approaches are theoretical and the question remains how to quantify the synchrony of the response variables. The specific response variables constitute multivariate time series with time-varying distributions. Additionally, (Dan-Glauser and Gross, 2013) recently stated that defining a measure of synchrony is a challenging and timely topic in emotion theory. In the following section, the results concerning the synchrony of peripheral physiological response variables during emotion eliciting situations are reviewed.

THE SYNCHRONY OF PERIPHERAL PHYSIOLOGICAL RESPONSE VARIABLES DURING EMOTION ELICITING SITUATIONS: A BRIEF REVIEW OF RESULTS
In early psychological research of emotion and stress, the synchrony of peripheral physiological measures was a major focus of the analyses (Wenger, 1942;Lacey and Lacey, 1958;Lazarus et al., 1963;Nesse et al., 1985). For example, the autonomic response system (ANS) plays an important role during stress (e.g., Carroll et al., 2009;Bibbey et al., 2013), posttraumatic stress disorder (e.g., Zucker et al., 2009;Ehlers et al., 2010) or anxietydisorders like panic-attacks (e.g., Roth et al., 1988;Meuret et al., 2011). Information on the temporal interdependences of peripheral physiological measures can enhance the understanding of the underlying functioning of the ANS (Kettunen et al., 1998;McAssey et al., 2013) and can provide important information about the psychophysiological processes during an emotional episode (McAssey et al., 2013). However, it is unclear under which conditions, and to what exact quantitative extent, a synchronous physiological reaction occurs (Hsieh et al., 2011;McAssey et al., 2013). Applying intra-individual time series models, (Kettunen et al., 1998) reported that the synchronization between electrodermal activity and heart rate within an individual is associated with a higher level of arousal and behavioral activity (see also Lazarus et al., 1963). In their stochastic network configuration approach, (Hsieh et al., 2011) found three clusters within an overall system that is formed by 15 psychophysiological signals: one behavioral cluster and two physiological clusters (blood pressure and cardiovascular parameters). They reported a higher synchronization between the different clusters as the intensity of an emotional stimuli was measured from a neutral condition, suggesting that there is a higher association during an emotional experience episode. A higher synchrony was also reported within each cluster. Based on brain activity analysis, (Costa et al., 2006) found by using a synchronization index, a higher synchrony of various EEG channels during emotional film stimuli than during neutral film clips. Their results indicate a higher information exchange during emotional responses (for similar results see also e.g., Miskovic and Schmidt, 2010).
Numerous studies on psychophysiological correlates of emotional stimuli have been undertaken. However, reactions from emotion eliciting stimuli are not universal on the inter-and intraindividual levels. Individuals vary in the intensity and duration of an emotional episode in terms of subjective experience, physiological and behavioral reactions (e.g., Grandjean et al., 2008;Kuppens et al., 2009). Thus, physiological response patterns during emotion tend to differ between individuals (e.g., Marwitz and Stemmler, 1998;Kristjansson et al., 2007). Kristjansson et al. (2007) applied a two level growth curve model, in which the first level explains the variance within participants and the second level the variance between participants. Marwitz and Stemmler (1998) used correlations and ANOVA to analyze individual response specificity. The physiological responses were found to be influenced by the individual appraisal (see Scherer, 2009, for an overview) , emotion regulation (e.g., Levenson, 1993, 1997;Dan-Glauser and Gross, 2011), and context specific attributes (see Cacioppo et al., 1992, for an overview). Dan-Glauser and , Gross and Levenson (1993), and Gross and Levenson (1997) used ANOVAs to detect the effect of emotion regulation on physiological data. Dan-Glauser and Gross (2013) showed in their study, by applying cross-correlations, that synchrony within the physiological channel decreased, if participants were instructed to suppress their emotions. Lacey and Lacey (1958), on the other hand, showed that the physiological response can also vary within an individual (see also Cacioppo et al., 1992). As described in the previous section, basic emotion theory, as well as some of the appraisal approaches, suggest a higher intra-individual synchrony during emotion. Nevertheless, some responses, e.g., respiratory and cardiovascular measures (respiratory sinus arrhythmia; RSA) are also synchronized in order to assist biological functions within our body (Yasuma, 2004;Ben-Tal et al., 2012;Garcia et al., 2013).
In contrast to the large number of studies that examine the correlates of emotional response systems to emotional stimuli (Rosenberg and Ekman, 1994;Calvo and Miguel-Tobal, 1998;Reisenzein, 2000;Bonanno and Keltner, 2004;Mauss et al., 2005;Sze et al., 2010;Hsieh et al., 2011;Dan-Glauser and Gross, 2013), there are only a few empirical studies that provide a quantitative measure of synchrony (Wenger, 1942;Lacey and Lacey, 1958;Lazarus et al., 1963;Nesse et al., 1985;Kettunen et al., 1998). There may be several reasons for this: Firstly, the analysis of physiological response variables is difficult, since they require multivariate, nonlinear, and non-stationary analysis methods (Zong and Chetouani, 2009). Non-stationarity arises when the joint probability density function (pdf) of the response variables changes over time (see next section for more details). Most of the approaches applied so far e.g., cross-correlation, implicitly rely on the stationarity of the physiological signals, and such an assumption is not fulfilled in practice (Muma et al., 2010).
Secondly, specifying a model that quantifies a time-varying synchrony of multiple response variables is not trivial (Dan-Glauser and Gross, 2013).
Thirdly, not only the psychophysiological responses, but also the emotions, impose additional challenges and the demand for sophisticated analytical methods. Emotions have often been treated as static phenomena, similar for different individuals, and have been analyzed by using nomothetic approaches, neglecting the intra-individual variability and the dynamic process of an emotional experience (Kuppens et al., 2009). Inter-individual analysis of the mean values of physiological responses overlooks the possible synchronous response within an individual over time. Therefore, an intra-individual analysis is more appropriate for the analysis of synchrony (Lazarus et al., 1963;Mauss et al., 2005;Hsieh et al., 2011). In the following section, a new approach, which takes the non-stationarity of the data into account, is introduced for analyzing multivariate synchrony of peripheral physiological measures in an individual during an emotional event.

A NEW APPROACH FOR THE QUANTIFICATION OF SYNCHRONY OF MULTIVARIATE PSYCHOPHYSIOLOGICAL SIGNALS
In this section, we introduce a new time-frequency based latent variable approach for the quantification of the synchrony of peripheral physiological responses, such as the activity of the heart, respiration, and the electrodermal activity level. The concepts of spectral bivariate coherence and time-frequency bivariate coherence from a signal processing perspective are first discussed. These concepts provide the basis for the quantification of synchrony of non-stationary psychophysiological response variables. Bivariate coherences are used as indicators of a latent state-space model, which quantifies one latent synchronized variable 1 .

SPECTRAL BIVARIATE COHERENCE
A frequently used function to examine the linear relation at frequency f between two signals x(t) and y(t), which are a function of time t, is the spectral coherence C XY (f ) (Marple, 1987;Brillinger, 2001), which is defined according to (1) Here, S XY (f ), S XX (f ), and S YY (f ) denote the cross-spectrum and the auto-spectra of x(t) and y(t), respectively. The signals x(t) and y(t) can, e.g., be two different psychophysiological time series measurements. In practice, spectra can be estimated, e.g., based on the periodogram, which is computationally efficient, since it uses the fast Fourier transform (FFT; Brigham, 2002). The raw periodogram is not a consistent spectral estimate, since its meansquared error does not decrease to zero as the number of samples used in the computation increases to infinity (Marple, 1987). Consistent spectral estimates can be obtained, e.g., by Welch's and Bartlett's methods, by approaches that: (i) split the original measurement into K segments, for which the periodogram is computed and (ii) obtain the overall spectral estimate by averaging over the K periodograms of the segments. Alternatively, consistency can be achieved by the Blackman-Tukey estimator which performs a smoothing of the periodogram. This is achieved by a convolution of the periodogram with a spectral window to reduce the variance. In the time-domain, this operation corresponds to a multiplication of the sample covariance sequence with a lag window of smaller size than the data size (Stoica and Moses, 2005).
C XY (f ) displays the consistency-of-phase-relationship between x(t) and y(t) and provides a frequency selective measure of the phase coupling between the two signals (Brillinger, 2001).
The values of the magnitude squared coherence |C XY (f )| 2 will always satisfy the relationship 0 ≤ |C XY (f )| 2 ≤ 1. Since C XY (f ) is normalized by the product of the auto-spectra (see Equation 1), 1 In this paper a distinction between spectral coherence and synchrony is made. Spectral coherence is a well defined bivariate measure in the signal processing literature (Marple, 1987;Brillinger, 2001). The term "synchrony" defines an underlying latent state of an overall system-wise synchrony. it is independent, e.g., of different amplitudes of x(t) and y(t). If x(t) and y(t) are completely uncorrelated, their coherence will be zero. If, on the other hand, |C XY (f )| 2 = 1 then y(t) can be fully predicted from x(t) by a linear and time-invariant system (Marple, 1987). The case |C XY (f )| 2 < 1, as described in Bendat and Piersol (1990), may be either due to: (i) noise entering the measurements, (ii) a non-linear functional relationship between the signals, or (iii) further inputs in addition to x(t) contributing to the output. It should be emphasized, however, that the above mentioned statements about C XY (f ) are only valid for second order stationary signals (Marple, 1987;Bendat and Piersol, 1990;Brillinger, 2001). Second order stationarity implies a constant mean and variance, as well as an autocovariance function that does not depend on t. To overcome the limitations that arise from the assumption of stationary signals, coherence measures for non-stationary signals have been developed. These are based on the wavelet transform (Grinsted et al., 2004) and time-frequency distribution based methods (White and Boashash, 1990;Matz and Hlawatch, 2000;Muma et al., 2010;Orini et al., 2011). The following example illustrates some limitations of classical spectral coherence analysis and provides motivation for the use of timefrequency representations. Figure 1 displays the time, frequency, and joint time-frequency representations of two non-stationary signals x(t) and y(t).
• x(t) consists of a superposition of a linear frequency modulated signal (a signal that consists of a single frequency component whose instantaneous frequency changes linearly over time; starting at 1 Hz for t = 0 and reaching 2 Hz at t = 20 s) and two sinusoidal signals with frequencies of 3 Hz and 4 Hz. • y(t) consists of a superposition of a linear frequency modulated signal (an instantaneous frequency starting at 1 Hz for t = 0 and reaching 2 Hz at t = 20 s) and a sinusoidal signal with a frequency of 4 Hz and of 10 s duration.
The time-frequency representation, like a sheet of music, displays the frequency content of a signal evolving over time. This wealth of information is lost when only the frequency domain is considered, since in the spectrum estimation process, averaging over time is performed. The loss of information is passed on to the spectral coherence measure C XY (f ). Figure 2 (first panel) displays the spectral coherence estimate for the above example.
C XY (f ) only indicates that there is a high synchrony between x(t) and y(t) in the frequency region between 1 and 2 Hz, but the evolution/progression over time is averaged out. Furthermore, C XY (f ) takes a value of about 0.5 at a frequency of 4 Hz, suggesting a moderate coupling between the signals. This average value neglects the fact that the coherence is nearly equal to one for half of the time and zero for half of the time. These examples illustrate, if the signals are non-stationary, i.e., their spectra evolve over time, that it is important to take the time variation information into account when establishing measures of coherence.

TIME-FREQUENCY BIVARIATE COHERENCE
The notion of time-frequency (TF) coherence was first defined by White and Boashash (1990) and then extended by Matz and Hlawatch (2000). The definition is where S XY (t, f ), S XX (t, f ) and S YY (t, f ) are the cross and auto time-frequency distributions of x(t) and y(t), respectively. There exist some conditions, which are stated in Matz and Hlawatch (2000), that are necessary for C XY (t, f ) to be well defined and produces meaningful results. In particular, these conditions guaran- and y(t) are related via a linear time-invariant filter of sufficiently short length compared to the stationarity width (for details see White and Boashash, 1990). The time-frequency distribution used in this paper, is the spectrogram which satisfies the conditions for C XY (t, f ) to be well defined and has been suggested for use in TF coherence estimation by White and Boashash (1990). The spectrogram is based on the short-term Fourier transform (STFT), and is simply a sequence of FFTs of windowed data segments, where the windows overlap in time. The spectrogram is defined as the magnitude square of the elements of the STFT. The spectrogram yields a time-frequency plot which contains columns of spectral estimates for a specific moment in time. In addition to the choice of the nature of the time-frequency distributions that underpin S XX (t, f ), S YY (t, f ) and S XY (t, f ), it is also necessary to perform a smoothing operation on the distributions, (Matz and Hlawatch, 2000). In this paper, as in Muma et al. (2010), the smoothing is performed by a 2 dimensional filtering with a Gaussian kernel as shown in Figure 3. The parameters of the Gaussian kernel define the time and frequency resolution capability of C XY (t, f ) and have been chosen such that both sharp changes in time (time resolution) and individual frequency regions of interest (frequency resolution) can be resolved. Figure 2 (second panel) illustrates the resolution capabilities based on the signals displayed in Figure 1. It can be seen that C XY (t, f ) is able to display the time-varying coherence of the two non-stationary signals. As discussed earlier, this is not possible with the spectral coherence C XY (f ), alone. Similarly, time domain methods, that rely on stationarity, e.g., correlation coefficients, cannot adequately describe the timevarying relationship between these two exemplary non-stationary signals. In particular, computing the correlation coefficient relies on a non time-varying correlation function (i.e., stationarity).
To illustrate the usefulness of time-frequency based methods for analyzing psychophysiological data, three synchronously measured signals have been considered: An electrocardiogram (ECG) signal, a respiratory signal, and a signal measuring the electrodermal activity level (EDA). The results are shown in Figure 4. Visual inspection quickly reveals that the spectra change over time, indicating the non-stationarity of the signals. A more formal testing for stationarity of some physiological signals, among them ECG signals and respiratory signals, has been performed in Muma et al. (2010) using a frequency domain test, suggested by Brcich and Iskander (2006). The test is able to determine if stationarity exists at a given frequency region of interest. The test, in general, rejected the hypothesis of stationarity for ECG and respiratory signals. As can be seen in Figure 4, the power of the signals is not evenly spread out in the time-frequency plane, but instead is concentrated in delineated regions. If only one such region exists, the signal is referred to as "mono-component," while if multiple delineated regions exist, the signal is referred to as "multicomponent." For example, a linear frequency modulated signal is a mono-component signal. An ECG signal is a multi-component signal, which consists of the pulse frequency region (varying between about 1.2 and 1.4 Hz for the example shown in Figure 4, first panel) and its harmonics, which are located in regions given by integer multiples of the pulse region. The power of the EDA (Figure 4, second panel) and respiration (Figure 4, third panel) signals is primarily concentrated in one region, with only a small amount of power in other regions. Unlike the ECG and respiration signals, the EDA, in general, does not have a cyclic nature, but contains trends and abrupt changes. Most of its signal power is concentrated in the lowest frequencies (see Figure 4, second panel) and forms the delineated region of the EDA in the timefrequency plane. However, analysis of the EDA data showed that the EDA signal also contains power in other regions during some time intervals. For example, frequency domain analysis shows that both the EDA data of Figure 5  coupling depends on an emotional state. The existence of a small amount of power in higher frequency regions of the EDA signal, is evident in Figure 6 where the element-wise logarithm of the timefrequency distribution of the EDA signal depicted in Figure 4, second panel, is shown. For this example, the pulse frequency region and its harmonics become clearly visible in the EDA signal. For a thorough analysis of EDA signals, the interested reader is referred to Lim et al. (1997) . Figure 7 plots the pairwise coherences (as given by Equation 2) of the three signals. For this example, the coherence between the ECG and EDA signals (first panel) is "high" in the pulse frequency region and its harmonics, while it is "low" in the frequency region that contains most of the EDA signal power. The coherence between the ECG and respiration signals (second panel), on the other hand, takes maximal values at the respiration frequency region, is "moderate" at the pulse frequency region and its harmonics, and "low" elsewhere. The coherence between the EDA and the respiration signals (third panel) is low for the respiration and EDA frequency regions and "moderate" elsewhere. Instead of a narrative description of the pairwise relationships, we briefly describe a scalar measure which quantifies this information. For this, the concept of coherences within delineated regions in the time-frequency plane, as detailed in Muma et al. (2010), is followed. Figure 8 shows the pulse and respiration frequency regions. Muma et al. (2010) proposed an algorithm to detect these regions within the time-frequency plane. The idea of the algorithm, in brief, is: First, the frequency of the auto-spectrum that contains maximal energy f max,0 is determined for a data segment, in our case, of length 10 s. The frequency of maximal energy f max (t) of the non-stationary signal varies around f max,0 and can be found by searching near the maximal value of the periodogram at each time instant. The delineated regions include all neighboring frequencies for which the power drops less than 3 decibel compared to f max (t). By applying a mask onto C XY (t, f ), which is equal to one within the detected regions and zero elsewhere (see Figure 8), analysis is restricted to the regions of interest. Note that, unlike classical spectrum based heart rate variability (HRV) analysis, (e.g., Niskanen et al., 2004;Acharya et al., 2006) and references therein, fixed frequency bands (VLF, LF, HF), which are set a-priori and remain constant over time, are not assumed. Our algorithm automatically adapts the frequency region of interest over time to the data at hand. Based on this, scalar time-varying coherence measures between two signals, such as, for example, the average coherence at a frequency region of interest for each time instant, can be determined. Figure 9 plots the pairwise (bivariate) coherence measure for the three signals in the respective frequency regions. It is clearly visible from these examples that a coupling between the signals exists and that the amount of synchronization varies over time. In the next section, a systemwise coherence measure based on the pairwise coherences, is defined.

A LATENT STATE SPACE MODELING APPROACH TO A SYSTEM-WISE SYNCHRONY MEASURE
In this subsection, the well-known state space modeling approach (e.g., Durbin and Koopman, 2001) is used to specify an overall synchrony measure. The basic idea is to use the multivariate time series of pair-wise coherence measures within delineated frequency regions (see previous subsection and Figure 9) and to specify a measurement model that operationalizes a latent variable which represents an underlying state of an overall systemwise synchrony. In addition to the measurement model, a structural model describes the regression of a state relative to its previous states.
For a given individual i, the specification of the measurement and structural model is given in the so-called non-innovation form (e.g., Gilbert, 1993): In the measurement model (Equation 3), we assume that a given p-dimensional observed variable z it , at time point t, can be regressed on a q-dimensional (latent) state vector ξ it , where e it is a p-dimensional residual vector (white noise). The variable z it hereby includes the pair-wise coherence measures in selected frequency regions (delineated regions), whereas ξ it represents the overall synchrony for an individual i at time point t, which is the variable of interest. The estimation of the pair-wise coherence measures is described in the previous section. H t and R t are timevarying coefficient matrices, respectively of dimensions (p × q) and (p × p) that reflect the time-varying relationship between the latent variable vector ξ it and the observed variable vector z it , and the time-varying relationship between the residual vector e it and the observed variable vector z it , respectively. In the structural model (Equation 4), it is assumed that the state vector ξ it can be regressed on a previous state vector ξ i(t−1) and on a m-dimensional covariate vector u it (e.g., an intervention), where η it is a q-dimensional latent residual (white noise) vector. Again, F t (q × q), G t (q × m), and Q t (q × q) are time-varying coefficient matrices. The coefficient matrix H t indicates the timedependent strength of the relationship of the state variable ξ it and its indicators/measures. Predictions of latent states (ξ i(t+1) ) can be conducted, for example, by applying the Kalman filter (e.g., Grewal and Andrews, 2001;Shumway and Stoffer, 2006). By using a smoothing algorithm (e.g., Durbin and Koopman, 2001), scores for the latent variable ξ it can be derived for the time series. The model parameters can be, for example, determined via a maximum likelihood estimator and imposing additional distributional assumptions (for example, normality of the latent variable ξ it ). It is not possible to freely estimate all the parameters in Equations (3) and (4). In order to have an identified model, the typical additional restrictions of latent variable modeling (e.g., scaling variables) are needed (Bollen, 1989). For a given individual i the timedependent pattern of the parameters reflects the changing structural relationship of the individual reactions in the observed variables over the course of time; specifically, here, the psychophysiological responses (e.g., Marwitz and Stemmler, 1998). The time series of the resulting latent states ξ it represents the course of the overall system-wise synchrony. Figure 10 gives an example of a time series of the system-wise synchrony measure (latent state variable ξ it ) during an emotional episode. The first 10 s (100 points) show the synchrony while a participant was watching a affective neutral picture (barstool: #7025) of the International Affective Picture System (IAPS; Lang et al., 2008). During the last 10 s (100 points) the participant was confronted with a disgust eliciting picture (half ripped-off finger: #3150). As a descriptive result, it can be seen the overall synchrony measure was, on average, higher during the disgust eliciting picture than during the neutral picture.

APPLICATION OF TWO APPROACHES FOR THE QUANTIFICATION OF SYNCHRONY USING DATA FROM AN EMOTION REGULATION STUDY
In this section, we illustrate the application of the proposed measure of synchrony, and the application of a measure proposed by Hsieh et al. (2011), to data from a larger emotion regulation study where psychophysiological signals were collected while participants were watching funny film clips. First, a short description of the larger study, is given. Second, assuming that the induction of emotions leads to person-specific changes in synchrony of the psychophysiological signals, the time series of the overall synchrony measure is analyzed by applying the proposed approach. Third, the results are compared with those obtained by using the approach for the quantification of the (aggregated) synchrony of Hsieh et al. (2011).

Participants
The sample of the emotion regulation study consisted of 58 German male undergraduate students of engineering. The mean age of the overall sample was 23.11 years (SD = 3.72). The procedure was fully explained to participants before written informed consent was obtained. Participants took part in a lottery where they could win a portable media player.

Procedure
Participants watched a random sequence of two funny clips (each 10-min long), which were cuts of the German version of the slapstick comedy film "You Don't Mess with the Zohan," starring Adam Sandler. In a small pilot-study, ratings of the clips showed that one film clip was slightly funnier than the other. Blood pressure was measured 5 and 10 min after each clip. Physiological measures were obtained continuously, and self-reported ratings of experienced amusement were collected for each film clip.

Measures
The data consisted of physiological and subjective responses during the film clips. The subjective responses are not considered here. Physiological measures were continuously sampled with a frequency of 256 Hz. The following signals were obtained using a BIOPAC MP150 System (Biopac Systems Inc., 2011): ECG, EDA, and respiration signals. In addition to these continuous signals, discrete measures of systolic and diastolic blood pressure were taken every 5 min. In this article, analysis was restricted to the ECG, the EDA, and the respiration signalw. Motion artifacts in the ECG signals were removed using a method proposed by Strasser et al. (2012).

APPLICATION OF THE PROPOSED APPROACH FOR THE QUANTIFICATION OF SYNCHRONY
In order to obtain an overall quantity of the synchrony, the following two step procedure was applied. First, the time-frequency based pairwise coherence measures were obtained. Using three signals (here: ECG, EDA, and respiration) results in six pairwise measures (two measures for each pair of physiological signals; see for example Figure 9). Figure 11 gives an example of the six resulting time series for one participant over a period of 120 s. During the first 60 s, the participant was watching the film clip, which was rated as moderately funny. During the last 60 s, the participant was shown the funnier film clip. As can be seen from Figure 11, some of the pairwise measures are very similar (for example Series 5 and 6). Second, the obtained pairwise measures z it = (z 1it , z 2it , . . . , z 6it ) T were indicators of a single overall latent variable ξ it time series of the system-wise synchrony using the state-space modeling approach as described above (cp. Equation 3). For a given person i, the measurement variable vector is given as: where a 1i , . . . , a 6i represent additional intercepts, h 1i , h 2i , . . . , h 6i are loadings, and e 1it , . . . , e 6it ∼ N(0, σ e ). N(0, σ e ) is the zero Gaussian distribution with a standard deviation equal to σ e . The latent variable ξ it was assumed to be uni-dimensional (see Equation 4), reflecting the overall synchrony. For a given person i, the structural model is given as: where g 1i is an intercept, f 1i is a regression coefficient, and η it ∼ N(0, σ η ). Figure 12 provides examples of overall synchrony measures for two participants. Again, during the first 60 s, participants were watching a film clip, which was rated as moderately funny. During the last 60 s the participants were shown a funnier film clip. Table 1 presents parameter estimates, standard error estimates, and confidence intervals for two participants (see Equations 5, 6). The results between the two participants vary substantially. As can be seen from Table 1, for participant #1 the loadings h 1 and h 3 are not significant and thus the measures z 1 and z 3 are not reliable indicators of synchrony. The other indicators (z 2 , z 4 , z 5 , z 6 ) indicate the degree of overall synchrony. In contrast, for participant #2 a different set of indicators (z 1 , z 4 , z 5 , z 6 ) imply an overall level of synchrony. Furthermore, the sign of the loading of indicator z 4 is different, which means that for participant #1 an increase in the overall synchrony leads to an increase in the relationship of the two variables associated with respiration and EDA (h 4 = 0.97), while for participant #2 an increase in the overall synchrony leads to an decrease in the relationship of the two variables associated FIGURE 7 | Example of coherences between ECG, EDA, and respiration signals in the time-frequency domain. The first panel shows the coherence between the ECG and EDA signals, the second panel shows the coherence between the ECG and respiration signals, and the third panel shows the coherence between the EDA and respiration signals.
with respiration and EDA (h 4 = −0.39). These simple examples show that the patterns of synchrony of the psychophysiological signals are very person-specific.

APPLICATION OF A COMPLEMENTARY APPROACH FOR THE QUANTIFICATION OF THE SYNCHRONY PROPOSED BY HSIEH ET AL. (2011)
To provide a comparison to existing procedures, the recently proposed method by Hsieh et al. (2011) is applied to the dataset. The method estimates an aggregated system-wise synchrony by creating a stochastic network, where high connectivity corresponds to high synchrony of the J = 3 signals (here: ECG, EDA, and respiration signals). First, receiver operating characteristics (ROC) were obtained. For a random variable X, the ROC displays the probability of correct detection, i.e., the probability of correctly accepting the null hypothesis, as a function of x, vs. the probability of false alarm, i.e., the probability of falsely accepting the null hypothesis, as a function of x. In our experiment the null hypothesis (H 0 ) is associated to the time interval with the moderate funny film clip, i.e., the baseline, while the alternative (H 1 ) is associated with the funnier (emotionally more activating) film clip. The ROC area A is defined as the area between the ROC curve and a diagonal line given by an identity function: Here, c(x) is a linear function, that maps the interval of x onto the interval [0, 1], which is equivalent with the diagonal line in this context. P(H 0 |H 0 , x) and P(H 0 |H 1 , x) are the probabilities of correct detection and false alarm given a threshold at x, respectively (Hsieh et al., 2011). For the signals acquired while the participants were watching both film clips, in accordance with Hsieh et al. (2011), are divided into K 2 = 12 non overlapping segments. The ROC areas were then calculated by defining P(H 0 |H 0 ) as the probability of deciding for a particular film clip given that the participants were actually watching this film clip. In this case, we used the cumulative distribution function (CDF) of the single segments, and P(H 0 |H 1 ) as the probability of deciding for a particular film clip given that the other film clip was watched. P(H 0 |H 0 ) and P(H 0 |H 1 ) are used to form the CDFs of the baseline and of the activation phase, respectively.
In the second step, the Spearman rank coefficients of the ROC areas of 24 non-overlapping blocks were calculated. These were formed by 12 blocks of 50 s each from the emotionally neutral phase and 12 blocks of 50 s each from the emotionally eliciting phases. It should be noted that the acquired signal blocks exhibited a non-stationary character. By calculating the Spearman rank coefficients of the ROC area sequences of the non-stationary signal blocks for each individual and signal-pair, a single scalar value is obtained as an averaged measure of coherence. This scalar does not take into account any stochastic changes of the signals, and, hence, it disregards significant information. Applying this procedure to every subject and signal pair, yields a J × J rank coefficient matrix for each subject and phase. Applyinĝ Here the 1st, 2nd, and 3rd rows and columns belong to the ECG, EDA, and respiratory signals, respectively. In the third step, a stochastic network was defined, where the nodes correspond to the ECG, EDA, and respiratory signals, and the edges were defined by values from the partition matrixp(h). To introduce randomness, a matrix u was formed whose entries were generated by the uniform distribution on the interval [0, 1]. In general, an edge between two nodes j, j was established, if b con (j, j ) = p(j, j ) > u(j, j ), i.e., if the Boolean matrix b con had a "1" entry at (j, j ). If p(j, j ) ≤ u(j, j ), there was no connection between j and j nd the correspondent entry in the b con matrix was 0. Each node of the stochastic network is either activated or inactive. To test for system-wise synchrony, in the beginning, a node must be selected for activation. After activation, the node sends a signal, via its edges, to other directly connected nodes, which then activate their neighboring nodes. After sending out the signal, the node that was activated in the first place becomes inactive. This activation, and inactivation, of nodes was iterated until (i) all nodes became activated, in which case system-wise synchrony was achieved, or (ii) a system state was obtained for which system-wise synchrony was impossible to achieve. For J nodes, this was done J times, each time activating a different node at the beginning. In our case of J = 3 nodes, system-wise synchrony could only be achieved by a fully connected network. Based on a Monte Carlo simulation with 10,000 repetitions, the stochastic network achieved system-wise synchrony 511 times. Thus, the resulting stochastic network had a probability of about 5.1% to achieve system-wise synchrony.
Due to the previous averaging over the subjects, the calculated probability gives us an expected value for achieving system-wise coherence for any of the given subjects. It is important to realize that the result of the simulation is not dependent on an individual subject and that the results depend on the assumption of stationarity.

DISCUSSION
A new approach for the quantification of synchrony of multivariate non-stationary psychophysiological signals has been proposed. After calculating bivariate time-frequency based coherence measures, a state space modeling procedure was applied to obtain an overall measure of synchrony. The approach gives information on the intra-individual level about the course of synchrony of psychophysiological reactions.

METHODOLOGICAL AND SUBSTANTIVE CONSIDERATIONS
The use of multivariate time series of physiological reactions for the quantification of an overall synchrony has several methodological and substantive implications. Firstly, the FIGURE 10 | Example of an overall synchrony measure while a participant was watching a neutral picture and a disgust eliciting picture for 10 s each.
approach provides time-sensitive information about the intraindividual level of synchrony. Thus, it is complementary to alternative approaches, which give information on cross-subjects based (aggregated) synchrony. See, for example the approach proposed by Hsieh et al. (2011) 2 .
Secondly, the proposed approach also provides information about the inter-individual difference of the structure of the synchrony measure. Specifically, the person-specific reaction patterns of how the physiological signals are related with the (latent variable) synchrony have been quantified. The factor loadings within the state space approach provide person-specific information on how the simultaneous activation of two given signals, e.g., ECG and EDA signals, is determined by an overall synchrony and whether the same signals are reliable indicators of synchrony across individuals. As we have seen from the above discussed empirical example, this is not necessarily the case.
Thirdly, by applying the time-frequency based procedure from Muma et al. (2010), the proposed approach directly addresses problems of non-stationarity. In general, non-stationarity is a challenge in the examination of time series Scharf (1991);Vaseghi (2008). When a relationship between longitudinal data has to be quantified, non-stationarity leads to a underestimation of the time-dependent relationship (see examples from above).
Fourthly, a practical implication of the proposed approach is that researchers are given a simple two step tool, which is able to answer foundational research questions from emotion theory. The (real time) temporal resolution of synchrony of (nonstationary) measures allows for a person-orientated identification FIGURE 11 | Example of six time-frequency based pairwise coherence measures (series) for one individual watching funny film clips. Series 1: ECG-EDA coherence (EDA region), Series 2: ECG-EDA coherence (ECG region), Series 3: RE-ECG coherence (ECG region), Series 4: RE-ECG coherence (RE region), Series 5: RE-EDA coherence (EDA region), Series 6: RE-EDA coherence (RE region). During the first 60 s, the participant was watching a film clip, which was rated as moderately funny. During the last 60 s, the participant was shown a funnier film clip.

www.frontiersin.org
January 2015 | Volume 5 | Article 1507 | 13 of response systems that have changing distributional characteristics over time (for example respiration rate). By addressing substantive questions on the intra-individual level the construct validity of concepts based on psychophysiological data can be strengthened (Borsboom et al., 2003). Lastly, there are implications for emotion theory. The proposed approach facilitates assessing whether emotion eliciting stimuli lead to a stronger synchrony of psychophysiological reactions during positive and negative emotional states. In principle, within the time series representation of the state space model, treatment effects can be specified. With our empirical data collected from a study on humor we only found weak FIGURE 12 | Example for two overall synchrony measures for two participants. During the first 60 s, the participants were watching a film clip, which was rated as moderately funny. During the last 60 s, the participants were shown a funnier film clip.
(descriptive) effects between two funny film clips. However, recent results from other studies on disgust indicate that differences between phases can be found. Based on our research, it is possible to conclude that testing for treatment effects in a time series is a precise and powerful tool which is in some cases superior to examining treatment effects with aggregated data.

LIMITATIONS AND FUTURE DIRECTIONS
The proposed approach provides novel possibilities for analyzing psychophysiological signals, nevertheless it also has some limitations. Both are briefly discussed: Although being able to incorporate non-stationary signals, the proposed approach is not capable of analyzing signals of constant amplitude. For example, when subjective data are measured using a potentiometer, participants can continuously adjust their subjective experience during emotional stimuli. Unfortunately, participants do not change their ratings continuously. Between changes in the position of the potentiometer, the signals are constant. For such time intervals, the proposed approach can not be used for the quantification of synchrony of physiological and subjective data. Nevertheless, participants have an enduring/varying emotional experience in this period. This general problem, which also occurs for (retrospective) paper-and-pencil ratings, remains unsolved. A potential solution could be an extension of the proposed approach, which operates on a feature space of the signals. Specifically, instead of looking directly at the bivariate coherence of measured signals, a feature estimation step could be added that extracts psychophysiological features from the signals (on a higher data level). The bivariate coherence from the first step may be based on the synchrony of these features. Enhanced robustness and a higher flexibility would possibly be the result. Also utilizing the feature space may save computational resources: Features change more slowly, which allows for reducing the sampling frequency in the coherence computation. An additional limitation of the proposed approach is the assumed linearity of the bivariate coherences. In future work, the linear coupling between the signals could be relaxed, which would reflect a more realistic specification of underlying processes.
A possible future research direction consists in a detailed examination of coherence, and overall synchrony patterns, for specific individuals during specific emotion eliciting situations. On the one hand, individuals and their emotions could be recognized in applied settings (e.g., man-machine interaction). On the other hand, basic research questions of emotion theory could be addressed on a more detailed level. For this, further evaluation of more psychophysiological signal types in different emotion eliciting situations will be necessary. For example, the proposed approach works well in situations, in which emotions are elicited during a relatively long period of time (here several seconds). However, it is unclear how the proposed approach performs in the case of, for example, emotional priming, in which emotions are elicited for a very short period of time. The applicability of the proposed approach strongly depends on the responsiveness of physiological signals that were chosen for the analysis. Respiration, ECG, and EDA, for instance, are very slow or inadequate for the representation of affective responses in the context of priming. Therefore, adequate signals such as EEG waves should be chosen for the analysis.
From a practical perspective, it would be interesting to extend the use of the proposed approach to applications in the context of the strongly developing research field of affective computing (e.g., Picard, 1995Picard, , 1997Scherer, 2010). Computers, which recognize, interpret, and process emotions could be used not only in health services, but also in education science, human-computerinteraction and other applications. Regarding the performance of the affect detection, the comparison of different algorithms (in terms of classifiers), would be an important field of research (Hudlicka, 2003;Kolodyazhniy et al., 2011), which also depends on the integration of different measured response signals, for example speech, physiological, and behavioral/mimic measures. Clearly, response signals differ in the sense that some are more informative than others, depending on which emotions need to be separated (e.g., Larsen et al., 2003).
Furthermore, as an important application, biofeedback might be considered (Thompson and Thompson, 2003). Since biofeedback addresses numerous signal types, such as behavior, muscle tone, brainwaves, breath, skin conductance, heart rate, temperature and pain perception, and many more, an adaptation of the proposed approach needs to be conducted in future research in order to be applicable in different settings (e.g., LaVaque, 2003;Dawson et al., 2007;Tassinary et al., 2007).
Finally, although beyond the scope of this article, an interesting field of research has emerged that examines the synchrony of psychophysiological signals and empathy (Marci and Orr, 2006;Hulsman et al., 2011;Oliveira-Silva and Gonçalves, 2011;Reed et al., 2013). In this research, couples of subjects have been examined with respect to their synchrony of psychophysiological responses in the context of empathy (Levenson and Ruef, 1992). Although, the neurobiological mechanisms underlying the synchronous empathic responses are mostly unclear, an adaption of the proposed approach might be helpful for further research in the field. The extension is not straightforward, because the interindividual level is added to the current setting. Therefore, it might be interesting from a psychometric perspective to examine the synchrony of responses that stem from related but not identical subjects.