Multi-Phase Locking Value: A Generalized Method for Determining Instantaneous Multi-frequency Phase Coupling
MMulti-Phase Locking Value: A Generalized Method for Determining InstantaneousMulti-frequency Phase Coupling
Yuan Yang,
1, 2, ∗ Bhavya Vasudeva, Hazem H. Refai, and Fei He † Stephenson School of Biomedical Engineering, The University of Oklahoma, Tulsa, Oklahoma-74135, USA Department of Physical Therapy and Human Movement Sciences,Feinberg School of Medicine, Northwestern University, Chicago, Illinois-60611, USA Indian Statistical Institute, Kolkata, West Bengal 700108, India Department of Electrical and Computer Engineering,The University of Oklahoma, Tulsa, Oklahoma-74135, USA Centre for Data Science, Coventry University, Coventry CV1 2JH, UK
Many physical, biological and neural systems behave as coupled oscillators, with characteristicphase coupling across different frequencies. Methods such as n : m phase locking value and bi-phase locking value have previously been proposed to quantify phase coupling between two resonantfrequencies (e.g. f , 2 f/
3) and across three frequencies (e.g. f , f , f + f ), respectively. However,the existing phase coupling metrics have their limitations and limited applications. They cannot beused to detect or quantify phase coupling across multiple frequencies (e.g. f , f , f , f , f + f + f − f ), or coupling that involves non-integer multiples of the frequencies (e.g. f , f , 2 f / f / I. INTRODUCTION
Complex systems such as the human brain behave asa series of oscillators with their instantaneous phasesdynamically coupled over multiple frequency bands [1–5]. Methods such as n : m phase locking value (PLV)[6], bi-phase locking value (bPLV) [7] and their vari-ants [8–10] have previously been proposed to detectand quantify different types of phase coupling. The n : m PLV measures phase coupling between two res-onant frequencies when n cycles of one oscillatory sig-nal are phase locked to m cycles of another oscillatorysignal, i.e. mφ ( f n , t ) − nφ ( f m , t ) ≤ const [6, 11]. ThebPLV quantifies quadratic phase coupling among threefrequencies, where a pair of frequencies, f and f arecoupled to a third frequency f = f + f or f − f , i.e. φ ( f , t ) ± φ ( f , t ) − φ ( f , t ) ≤ const [7].However, phase coupling could be shown in more com-plicated patterns involving more than three frequencies(e.g. f , f , f , f , f + f + f − f ) as well as their non-integer multiples (e.g. f , f , 2 f / f / ∗ [email protected] † [email protected] approach for quantifying integer multi-frequency phasecoupling [12]. This method has been applied to the hu-man nervous system to advance our understanding ofnonlinear neuronal processes and their functions in move-ment control [13, 14] and sensory perception [15]. TheMSPC is a straightforward extension of bPLV based onhigh-order spectra [16]; however, it does not cover eitherthe non-integer multi-frequency phase coupling (e.g. f , f , 2 f / f /
3) or non-integer resonant coupling (e.g.2:3 coupling [17] revealed by n : m PLV) problems.Thus, this paper aims to introduce a more generalizedapproach, namely multi-phase locking value (M-PLV),that integrates the concepts of MSPC and n : m PLV toallow the detection and quantification of various types ofphase coupling, including integer and non-integer, multi-frequency and resonant phase coupling. The proposedM-PLV provides us with a tool to explore the unre-ported non-integer multi-frequency phase coupling thathas never been captured by existing phase coupling meth-ods. Furthermore, different from commonly used instan-taneous phase coupling metrics, the proposed methodalso allows the detection of delayed phase coupling andthe associated time lag between coupled oscillators. Wetested M-PLV on two scenarios where synthetic coupledsignals are generated using white Gaussian signals, and asystem comprised of multiple coupled R¨ossler oscillators.The rest of this paper is organized as follows: SectionII describes M-PLV, Section III summarizes the experi-ments used to validate the method, Section IV presents a r X i v : . [ q - b i o . N C ] F e b the results and discussion, and Section V concludes thepaper. II. MULTI-PHASE LOCKING VALUE (M-PLV):THEORY AND CALCULATION
The proposed M-PLV is generalized approach that in-tegrates the concept of MPSC [12] and n : m PLV [6]. Itnot only provides us with a formulated mathematical de-scription for the phase coupling problems separately de-scribed by MSPC and n : m PLV, but also permits the de-tection and quantification of non-integer multi-frequencyphase coupling that cannot be assessed by using existingphase coupling methods.
A. M-PLV
The MSPC considers the case where multiple inputfrequencies f , f , ..., f L are coupled to an output fre-quency f Σ based on an integer combination, such that f Σ = (cid:80) Ll =1 m l f l , m l ∈ N : (cid:32) L (cid:88) l =1 m l φ ( f l , t ) (cid:33) − φ ( f Σ , t ) ≤ const (1)The MSPC does not cover the case where non-integermultiples of input frequencies are coupled to the outputfrequency. To address this gap, the proposed M-PLVgeneralizes the relation between frequencies as nf Σ = (cid:80) Ll =1 m l f l or f Σ = (cid:80) Ll =1 m l n f l . It can be seen that al-though m l , n are integers, their ratio can give real num-bers. This idea is in line with the concept of n : m PLV [6], but allows assessment of phase coupling betweenmultiple input frequencies and one targeted output fre-quency.Moreover, there may exist a delay τ in the systembetween the input and the output, such that the cou-pling can be detected only after this delay has been com-pensated by aligning the indices of all the instantaneousphases. Incorporating these factors, the proposed M-PLVaims to detect and quantify a more generalized phasecoupling phenomenon that can be described as: (cid:32) L (cid:88) l =1 m l φ ( f l , t − τ ) (cid:33) − nφ ( f Σ , t ) ≤ const (2)Based on this theoretical definition, the formula of M-PLV (Ψ) is given as follows for the calculation:Ψ( f , f , ..., f L ; m , m , ..., m L , n ; t, τ ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) K K (cid:88) k =1 exp (cid:32) j (cid:32)(cid:32) L (cid:88) l =1 m l φ k ( f l , t − τ ) (cid:33) − nφ k ( f Σ , t ) (cid:33)(cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (3) where K is the number of observations, φ k ( f l , t ) is aninstantaneous input phase at the k th observation, whichcan be obtained from the Hilbert transform of narrow-band filtered time series with the spectrum centered atfrequency f l [18].When τ = 0, the proposed method can be used for de-tecting and quantifying the simultaneous multi-frequencyphase coupling. Additionally, if n = 1, M-PLV is furtherreduced to MSPC, for measuring simultaneous integermulti-frequency phase coupling: M SP C ( f , f , ..., f L ; m , m , ..., m L ; t ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) K K (cid:88) k =1 exp (cid:32) j (cid:32)(cid:32) L (cid:88) l =1 m l φ k ( f l , t ) (cid:33) − φ k ( f Σ , t ) (cid:33)(cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (4)Noteworthy, bPLV [7] is basically a special form ofMSPC or M-PLV when the interest is in determiningquadratic phase coupling: bP LV ( f , f ; t ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) K K (cid:88) k =1 exp ( j (( φ k ( f , t ) + φ k ( f , t )) − φ k ( f + f , t ))) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (5)When L = 1, M-PLV can also be reduced to n : m PLV:
P LV n : m ( f n , f m ; m, n ; t )= (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) K K (cid:88) k =1 exp ( j ( mφ k ( f n , t ) − nφ k ( f m , t ))) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (6)As such, M-PLV not only allows the detection andquantification of non-integer multi-frequency coupling,but also provides a generic mathematical framework thatcan accommodate all common forms of phase coupling inthe existing literature. B. Detecting significant M-PLV
In order to detect the time window and frequency atwhich phase coupling is significant, a reference thresh-old value of M-PLV is required. For this purpose, the95% significance threshold is obtained by a Monte Carlosimulation [12]. The null hypothesis is that the phase dif-ference ∆ φ ( t ; k ) is completely random so that the cyclicphase difference ∆ φ ( t ; k ) = ∆ φ ( t ; k ) mod 2 π will be uni-formly and randomly distributed in the interval [0 , π ].The M-PLV corresponding to other frequency combina-tions for all instants t as well as those corresponding tothe combination of interest for the instants t (cid:48) = t − t c ( t c is the estimated coupling window) are taken as surro-gate data of uniformly and randomly distributed phasevalues of ∆ φ ( t ; k ). This procedure is repeated N (e.g.1000) times to obtain the statistical distribution of M-PLV values for the different number of observations, K ,(e.g. K = 50 , , . . ., C. Delay Estimation
In order to estimate the delay τ , the M-PLV for differ-ent values τ i within a given range calculated. The valueof τ i corresponding to the maximum value of M-PLV isthe estimated delay ˆ τ of the system. III. EXPERIMENTS
We tested M-PLV on two scenarios where syntheticcoupled signals are generated satisfying (1) using whiteGaussian signals alone, as well as (2) from a system com-prised of multiple coupled R¨ossler oscillators. In the fol-lowing simulations, the sampling frequency is 1 kHz.
A. Coupled white Gaussian signals
In this case, x ( t ) and y ( t ) are two independent whiteGaussian signals (zero mean and unit variance). Thesynthetic signal, y c ( t ) is generated as follows: y c ( t ) = y ( t ) − y ( f Σ , t c )+ x | m | ( f , t c ) x | m | ( f , t c ) A | m | x ( f , t c ) A | m | x ( f , t c ) A y ( f Σ , t c ) (7)where t is in the range of [0.001,10] s, x ( t c ) represents x ( t ) in the phase coupling time window t c = [2 . , . x ( f , t c ) is a narrowband signal with the spectrumcentered at frequency f , which is obtained after x ( t c )is passed through a Butterworth band-pass filter [19]centered at frequency f (bandwidth: 2 Hz, 6 th order). A x ( f , t c ) is the envelope of the Hilbert transform of x ( f , t c ). In order to eliminate the effect of filter on thesignal phase, zero-phase shift filter (Matlab function: filt-filt.m) is used in this study. The normalization of thesignal x ( f , t c ) by its envelope A x ( f , t c ) prevents abruptchanges in its amplitude.In these designed signals, there is phase coupling be-tween y c ( t ) and x ( t ) in the time interval t c , following therule f Σ = m f + m f , serving as the ground truth inthis “white” box problem for testing the M-PLV for in-teger ( n = 1) multi-frequency phase coupling with zerodelay ( τ = 0).In order to check for the phase coupling between x ( t )and y c ( t ), M-PLV is calculated based on Eq.(3), and theset of input frequencies includes f and f . B. Coupled R¨ossler oscillators
In this case, y ( t ) is white Gaussian signal, while x i ( t )are obtained from a system comprised of coupled R¨ossleroscillators, which consists of N − N th oscillator. The system is character-ized by the following equations:˙ x i = − N − (cid:88) j =1 m j ω j n x i − z i + ε i N − (cid:88) j =1 m j x j n − x i (8)˙ u i = ω i x i + au i (9)˙ z i = c + z i ( x i − b ) (10)where ε i = 0 for i < N and ω j = 2 πf j . These coupledoscillators are designed to mimic a multiple-input-single-output (MISO) system. In this case, Eq.(7) can be gen-eralized to include a larger number of signals coupled atdifferent frequencies, so that nf Σ = (cid:80) Nl =1 m l f l and thecoupled signal can be obtained as follows: y c ( t ) = y ( t ) − y ( nf Σ , t c )+ A y ( nf Σ , t c ) (cid:32) x | m | ( f , t c ) A | m | x ( f , t c ) ... x | m N | N ( f N , t c ) A | m N | x N ( f N , t c ) (cid:33) (11)where x i ( f j , t c ) is obtained after x i ( t c ) is passedthrough a Butterworth band-pass filter [19] centered atfrequency f j (bandwidth: 2 Hz, 6 th order). In order tointroduce a delay τ in the system, t c can be replacedby t c − τ in the above equation. The coupling is evalu-ated between x N ( t ) and y c ( t ) by calculating the M-PLVaccording to Eq.(3).The 95% significance threshold and delay τ can be es-timated through the procedure described in Section II.Band II.C. IV. RESULTS AND DISCUSSIONA. Coupled white Gaussian signals: integer ( n = 1 )multi-frequency phase coupling with zero delay( τ = 0 ) The results are shown for f = 29 Hz, f = 13 Hz, m = 2, and m = −
1, so that f Σ = 2 × − ×
13 = 45Hz. Fig.1 shows M-PLV plotted as a function of timeand frequency for varying numbers of epoches K ( K =500 , , f = 29 and f = 13 Hz toexamine whether the significant M-PLV is only detectedon the target frequency 45 Hz rather than other frequen-cies. It is observed that M-PLV shows significant valuesfor f Σ = 45 Hz in the time window t (cid:48) c ∼ t c = [2 . , . t (cid:48) c = [2 . , . . , . . , .
6] s for K = 500 , , and 900, respectively. Theerror of time window estimation can be defined as thedifference between t (cid:48) c and t c and divided by the windowsize. The errors are below 5 % for all tested K values. Tofurther demonstrate the performance of M-PLV, Fig. 2shows a few of example plots of M-PLV for K = 600for some possible combination frequencies of f = 29 and f = 13 Hz (e.g. 29 − ×
13 = 3 Hz, 0 × ×
13 = 39 Hz,etc). Significant M-PLV is only detected at the targetedfrequency f Σ = 45 Hz within the coupled time window. B. Coupled R¨ossler oscillators: multi-frequencyphase coupling with a delay
In these simulations, we set K = 400, N = 3, and theparameters of the coupled R¨ossler oscillators (Eq.(9) -(10)) as a = 0 . c = 0 . b = 10, and ε N = 0 . n = 1) multi-frequency phase coupling with zerodelay ( τ = 0) in a MISO system, the oscillators are sim-ulated for 80 seconds and two sets of 30 000 samples areobtained from the simulated signals, with t = [10 . , t c = [17 . , .
5] s for the first set and t = t + 40 s, t c = t c + 40 s = [57 . , .
5] s for the second set. Inthis case, f = 3 Hz, f = 5 Hz, m = − m = 2, sothat f Σ = − × × f = 3 Hz, f = 5 Hz to examine whether the significant M-PLV isonly detected on the target frequency 7 Hz rather thanothers (e.g. 2 × − × t (cid:48) c = [17 . , . t (cid:48) c = [57 . , . f = 7 Hz, f = 13 Hz, m = 1, m = 1, and n = 5 so that f Σ = (1 × × = 4 Hz. Also, t = [10 . ,
40] s and t c = [17 . , .
5] s. Fig. 5 shows the results obtainedfor f Σ = 4 Hz. Using the 95% significance threshold, t (cid:48) c = [17 . , . τ is set as 1 s. In this case, t = [10 . ,
40] s, t c = [17 . , .
5] s, f = 29 Hz, f = 13 Hz, m = 2,and m = −
1, so that f Σ = 45. Fig. 6 shows the aver-age M-PLV obtained for varying τ i . The estimated localmaxima over 10 such simulations is ˆ τ = 0 . ± . V. CONCLUSION
In this paper, a new method for quantifying multi-frequency phase coupling has been proposed. Thismethod addresses the limitation of existing approachesthat only allow the detection of coupling between tworesonant frequencies (i.e. n : m PLV) or quadratic cou-pling between three frequencies (i.e. bPLV). The M-PLVallows us to quantify various types of phase coupling,including both integer and non-integer phase couplingacross multiple frequencies, so as to permit the explo-ration of more complicated, even unreported phase cou-pling phenomena in the real world. Simulation studieshave been performed on synthetic coupled signals gener-ated using white Gaussian signals, and a complex systemcomprised of multiple coupled R¨ossler oscillators. Ourresults suggest that the proposed method can achieve areliable estimate of the frequency combination as well asthe time window during which phase coupling is present.Furthermore, this method can be used for a precise es-timation of the delay between the input and the outputwhen delayed phase coupling is present between the os-cillators.
ACKNOWLEDGEMENT
This work was supported by NIH/NICHD1R21HD099710 and the Dixon Translational Re-search Grants Initiative from Northwestern MemorialFoundation and Northwestern University Clinical andTranslational Sciences Institute (UL1TR001422), North-western University. B. Vasudeva received stipends fromS. N. Bose Scholars Program 2019 and Departmentof Physical Therapy and Human Movement Sciences,Northwestern University. [1] F. Varela, J.-P. Lachaux, E. Rodriguez, and J. Mar-tinerie, The brainweb: phase synchronization and large-scale integration, Nature reviews neuroscience , 229(2001).[2] M. Breakspear, Dynamic models of large-scale brain ac-tivity, Nature neuroscience , 340 (2017).[3] R. T. Canolty and R. T. Knight, The functional role ofcross-frequency coupling, Trends in cognitive sciences ,506 (2010). [4] O. Jensen and L. L. Colgin, Cross-frequency coupling be-tween neuronal oscillations, Trends in cognitive sciences , 267 (2007).[5] F. He and Y. Yang, Nonlinear system identification ofneural systems from neurophysiological signals, Neuro-science, in press (2021).[6] G. B. Ermentrout, n: m phase-locking of weakly cou-pled oscillators, Journal of Mathematical Biology , 327(1981). FIG. 1: M-PLV as a function of time and frequency for K = (a) 500, (b) 750 and (c) 900.FIG. 2: M-PLV for K = 600 as a function of time (unit: ms)for the set of frequencies (a) 3, (b) 39, (c) 45, (d) 55, (e) 71,and (f) 87 Hz.FIG. 3: M-PLV for the first set of values of the coupledR¨ossler oscillators, as a function of time, for the set offrequencies (a) 1, (b) 7, (c) 9, (d) 11, (e) 13, and (f) 15 Hz.[7] F. Darvas, J. Ojemann, and L. Sorensen, Bi-phase locking— a tool for probing non-linear interaction in the humanbrain, NeuroImage , 123 (2009).[8] B. Schelter, M. Winterhalder, R. Dahlhaus, J. Kurths,and J. Timmer, Partial phase synchronization for mul-tivariate synchronizing systems, Physical review letters FIG. 4: M-PLV for the second set of values of the coupledR¨ossler oscillators, as a function of time, for the set offrequencies (a) 1, (b) 7, (c) 9, (d) 11, (e) 13, and (f) 15 Hz.FIG. 5: M-PLV as a function of time for n = 5 (non-integermultiples) at frequency 4 Hz. , 208103 (2006).[9] M. Vinck, R. Oostenveld, M. Van Wingerden,F. Battaglia, and C. M. Pennartz, An improved indexof phase-synchronization for electrophysiological data inthe presence of volume-conduction, noise and sample-sizebias, Neuroimage , 1548 (2011).[10] Z. Vahabi, R. Amirfattahi, F. Shayegh, and F. Ghassemi,Online epileptic seizure prediction using wavelet-based FIG. 6: Average M-PLV as a function of delay τ . The local maxima occurs at ˆ τ = 1 .
02 s in this case.bi-phase correlation of electrical signals tomography, In-ternational journal of neural systems , 1550028 (2015).[11] M. Rosenblum, A. Pikovsky, J. Kurths, C. Sch¨afer, andP. Tass, Chapter 9 phase synchronization: From theoryto data analysis, in Neuro-Informatics and Neural Mod-elling , Handbook of Biological Physics, Vol. 4, edited byF. Moss and S. Gielen (North-Holland, 2001) pp. 279 –321.[12] Y. Yang, T. Solis-Escalante, J. Yao, A. Daffertshofer,A. C. Schouten, and F. C. T. van der Helm, A generalapproach for quantifying nonlinear connectivity in thenervous system based on phase coupling, InternationalJournal of Neural Systems , 1550031 (2016), pMID:26404514.[13] Y. Yang, T. Solis-Escalante, J. Yao, F. C. Van Der Helm,J. P. Dewald, and A. C. Schouten, Nonlinear connectivityin the human stretch reflex assessed by cross-frequencyphase coupling, International journal of neural systems , 1650043 (2016). [14] Y. Yang, J. P. Dewald, F. C. van der Helm, and A. C.Schouten, Unveiling neural coupling within the sensori-motor system: directionality and nonlinearity, Europeanjournal of neuroscience , 2407 (2018).[15] N. Gordon, N. Tsuchiya, R. Koenig-Robert, and J. Ho-hwy, Expectation and attention increase the integrationof top-down and bottom-up signals in perception throughdifferent pathways, PLoS biology , e3000233 (2019).[16] C. L. Nikias and J. M. Mendel, Signal processing withhigher-order spectra, IEEE Signal processing magazine , 10 (1993).[17] A. J. Langdon, T. W. Boonstra, and M. Breakspear,Multi-frequency phase locking in human somatosensorycortex, Progress in biophysics and molecular biology ,58 (2011).[18] B. Boashash, Estimating and interpreting the instanta-neous frequency of a signal. i. fundamentals, Proceedingsof the IEEE80
