Phase Locking Value revisited: teaching new tricks to an old dog

1 Phase Lock ing Value revisited : teaching ne w tricks to an old dog Ricardo Bruña 1,2 , Fernando Maestú 1,2 , and E rnesto Pereda 1 ,3 1 Laboratory for Cognitive an d Computational Neuroscienc e, Canter for Biomedical T echnology, Technical University of Madrid, P ozuelo de Alarcón, Madrid , Spain 2 Departamento de Psicología Experimental, Procesos Psicológicos y Logopedia, Faculty of Psycholog y, Complutense Univer sity of Madr id, Pozuelo de Alarcón, Ma drid, Spain 3 Electrical Engineeri ng and Bio engineering Group, D epartment of I ndustrial Engineering & IUNE, University of La Lagu na, San Cristobal de La Laguna, T enerife, Spain Keywords : Brain connectiv ity, Computationa l efficiency, Neuro science, Synchroniza tion Abstract: Despite the increase in calculation power in the l ast decades, the estimation of brain connectivity is still a tedious task. T he high computational cost of the algorithms esca lates with the square of the number of signal s evaluated, usua lly in the rang e of thousands. I n this work w e propose a re - formulation of a w idely used algorithm that allows the estimation of whole brain connec tivity in much smaller t imes. We start from the original implementation of Phase Locking Value (PLV) and re - formulated it in a highly com putational eff icient way. Besi des, this formulation stresses its strong sim ilarity with coherence, which we used to introduce t wo new m etrics insensitive to zero lag synchronization, the imaginary part of PLV (iPLV) and its corrected coun terpart (ciPL V). The new implementation of PLV avoids some highly CPU-expensive operati ons, and achieved a 100-fold speedup over the original algorithm . T he ne w deri ved metrics were h ighly robust in t he presence of volum e conduction. ciPLV, in particular, proved capable of ignoring zero -lag connectivity , while correctly estimating nonzero-lag connectivity. Our implem entation of PL V mak es it possible to calculate whole-brain connectivity in much shorter times. The results of the simulations using ciPLV suggest that this metr ic is ideal to measure synchronization in the presence o f volume conduc tion or source le akage effects. 2 I. Introduction The study of the brain caught the interest of research ers before neuros cience was even a discipline. For a long time, the brain was considered an organ, or rather a set of org ans, with distinct parts taking care of distinct duties, following the so -called phrenologica l point of view. However, in recent years some have begun to departu re from this traditional approach in which distinct parts of the brain are believed to have different functions. Instead, they argue than cognition, thought and action are supported by the collective action of sets of brain areas (network s) (Friston 1994). Within this parad igm, t he brain areas form a network (th e c onnectome ), entwined by white matter tracts (anatomical connectom e), in which different sub -netw orks communicate dynam ically with each other to perform different functions (functional conne ctome). Whereas the anatomical connectome is assessed from diffusion tensor imaging (DTI) ( Le Bihan et al. 2001), the functional conn ectome is constructed f rom neuroimag ing techniques such as functional magnetic reson ance imaging ( fMRI) and electro- and magneto-phy siology (Friston et al. 2003; Brookes et al. 2011 ). In the latter case, howev er, the de scription of the connectom e is not straightforward, as the underlying mechanics of the brain are unknown. Instead, we have to trust on measures of brain a ctivity at differen t recording sites . In this framework, functional connectivity (FC) i s defined as the existence of statistical dependence bet ween the activities in t wo or more sites above chance lev el, a dependence that c an be evaluated in d ifferent ways. In t he case of signals such as M/EEG, one of the m ost studied connectivity hypothesis is that of phase synchronization (PS) (Lachaux et al. 1999; Fries 2015; Garcés et al. 2016; López et al. 2014) . Under this sc enario, two brain areas show the activity of different os cillators a nd, in the ca se of connectivity between different regions, their oscillatio n propertie s (i.e. phase or frequency) should be related. I f this rela tion can be m athematically evaluated, w e get an e stimator of phase connectivity . Of the many PS measures available in t he l iterature, one of the most used, mainly due to its simplicity, is the Phase Lock ing V alue (PLV, (Lachaux et al. 1999)), sometimes called Mean Phase Coherence (MPC) (Mormann et al . 2000). T his measure evaluates t he instantaneous phase difference of the signals under the hypothesis that connected areas generate signals whose instantaneous phases evolve together. In this case, t he phases of the signals are said to be “locked” , and the ir difference is therefore constant. However, real-world signals are inherently noisy, and it is not always possib le to be sure that th e evaluated signa l only comes from one oscillator. The problem is solv ed by allow ing som e deviation from the condition of a constant phase difference. Thu s, PLV evaluates the spread o f the distribu tion of phase differences, and the connectivity estimation is l inked t o t his spread. The narrower the distributi on of the phas e 3 difference, the higher the PLV value, which ranges between zero ( no phase dependence) and one (complete phase dependence). We will elaborate on the mathematics behind the estimation of the PLV in the next s ection. Here we b riefly describe the m ain steps involved. The most common method to calcula te PL V is based on t he instantaneous phase of the signa ls obtained using th e Hilbert or the wavel et transform s. In both cases, the calculations a re fast and reliable, which has undoubted ly contributed to its use. The calculation of the phase difference and its spread is also mathem atically simple, and, in principle, not very time-consum ing. However, even when the c omputation o f PLV is f ast and easy for a fair am ount of data, the c omputational co st grows wit h the square of t he number of signals. This is especially important in the case of M/EEG, distributed source- space analysis, where such num ber ranges from t he thousands to the hundreds of thousands. It is also noteworthy that, des pite the estim ation of the PLV and other related indices from data has a long history, from time to time new studies appear in the litera ture (Ewald et al. 2012; Kovach 2017) in which different as pects of the calculation procedures and t he i ndices themselves are analyzed with a f resh eye to refine the existing metho dologi es. T his paper intends to be one of such studies t o shed new li ght on PS estimation, interpretat ion and use. Thus, we start fr om the ori gina l PLV formulation of Lachaux (Lachaux et al. 1999) and rewrite it to obtain an equivalent expression that is far easier to compute, reducing t he required tim e for the PLV calculation up to a f actor of 100. Besides, we also show that this new formulation is closely related to that of coherency, thereby allowing an interpretation of this well-studied functi on in terms of the PLV (see also (Kovach 2017) ). As an addi tional advantage, such interpretation allows the formulation of two PLV -based zero-lag- insensitive measures: t he imaginary PLV (iPLV) and the corrected iPLV (ciPLV), which can be us ed as alternatives to the imaginary part of t he coherency (Nolte et al. 2004) and its corrected version (Ewald et al. 2012), res pectively , i n the assessm ent of direct PS from M/EEG data. II. Methods A. Computational op timization In the ir original paper, Lachaux (Lachaux et al. 1999) defined the PLV as a time-dependen t connectivity measured tailored t o st udy evoked activity . The i dea behi nd their definition i s that the stimulus resets t he phase of the neural oscillators so that signals connected i n a given t ime should have a stable phase- difference along trials. Its mathematica l formulation reads : 4    󰇛  󰇜        󰇡  󰇛   󰇜   󰇛   󰇜 󰇢    (1) where  is the number of trials and   󰇛    󰇜 is t he instantan eous phase for signal  in trial  at time  . This definition can be extended to resting- state data, by assessing phase locking as a stable phase - difference over t ime, thereby obtain the so-called MPC 1 (Mormann et a l. 2000):         󰇡  󰇛  󰇜   󰇛  󰇜 󰇢     (2) where  is the data length. In either case, one has to extract t he i nstantaneous phase 󰇛 󰇜 of each signal. Besi des, for the phase to be physi cally m eaningful, it is necessary that only one oscillator is present in each si gnal. This is achieved, e.g., by m eans of a narrow -band pass fil tering or, equivalen tly, the convolutio n with a narrow band com plex wavelet such as that of M orlet (Bruns 2004). After the filte ring process, w e obtain a band-pass versi on of the Hilbert ana lytical signal:  󰇝  󰇛  󰇜 󰇞     󰇛  󰇜    󰇛  󰇜      󰇛  󰇜   󰇛  󰇜      󰇛  󰇜 (3) where   represent the Hilbert transfo rm of  , and BP stands for band pass . The instantaneo us phase is the angle between the real and the imaginary parts of the Hilbert analytical signal, or the angle between the original ( band-pass) signal and its H ilbert transform . The instantaneous phase is usua lly extracted from t his analytical signal, the phase differenc e estimated and, finally, the exponentiation calculated to get t he unit phase difference vector. However, these t wo operations (phase extraction and exponentiation) are computationally expensive, but, as we w ill show, they can be easily circumvented by using the pro perties of exponentials. First, let us obta in the oscillatory pa rt of the analyt ical signal by norm alizing ( 3):  󰇗    󰇛  󰇜      󰇛  󰇜     󰇛  󰇜     󰇛  󰇜     󰇛  󰇜   󰇛  󰇜      󰇛  󰇜 (4) From this, we eas ily derive the e xponential of the p hase difference: 1 Henceforth, for the sake of simplicity, we will use PLV i,j to refer to both the time varying phase locking value as derived by Lachaux et al., and the mean phase coherence . In any ca se, both indices c an be easily distinguished since the former one ( PLV i,j (t) ) explicitly depends on time , where as the latter one, aver aged over the w hole data segme nt, does not. 5  󰇗   󰇛  󰇜  󰇡 󰇗    󰇛  󰇜 󰇢      󰇛  󰇜      󰇛  󰇜       󰇛  󰇜     󰇛  󰇜   󰇡  󰇛  󰇜   󰇛  󰇜 󰇢 (5) where 󰇛  󰇜  represents com plex conjugate. Thus, expression (2) can b e rewritten as:         󰇗    󰇛  󰇜  󰇡 󰇗    󰇛  󰇜 󰇢      (6) or, using v ector algebra:        󰇗    󰇗    (7) where  󰇗  is a vector version of  󰇗   󰇛  󰇜 and  󰇗  , its transpose conjugate of  󰇗 . T his calculation is computationally very efficient, and allows a considerable speed up of the estim at ion with low memory penal ization. The calculat ion efficiency is discusse d in the Results sect ion. The efficient f ormulation can be extended to Lachaux’s PLV i,j (t) by constructing t he vectors with the  -th sample of each trial :    󰇛󰇜     󰇗     󰇗     (8) B. Relation to cohe rence If we expand (1) us ing the oscillatory part of  , we get:   󰇛  󰇜        󰇛  󰇜     󰇛  󰇜      (9) Phase synchronization only mak es sense for sig nals composed of a si ngle oscillatory component . Yet, i t is mathem atically possible to calculate the PLV from both a br oadband signal and an arbitrarily narrow ba nd one, ev en when these calcu lations do not hav e, in principl e, physical sense. If we take the extreme case of a sing le oscillator whose spectrum is non-zero only at a given frequency,   , the signal can be written as an out -of- phase cosine, where the phase at the initial time is equal to the Fourier phase, and the phase at any ot her time is det ermined by the delay and the frequency. 6   󰇛  󰇜     󰇛       󰇜      󰇛      󰇜    󰇛   󰇜        ( 10 ) where  󰇛  󰇜 is the Fourier transform of  󰇛  󰇜 at   and  󰇝  󰇞 stands for the real part of  . In this case, Lachaux’s definition of PLV is no longer tim e depe ndent, and its formulation can be simplified as:             󰇛  󰇜      󰇛  󰇜         󰈏   󰇛    󰇜  󰇡  󰇛    󰇜 󰇢     󰇛    󰇜    󰇛    󰇜    󰈏 ( 11 ) This formula closely r esembles that of coherence (Nunez et al. 1997) . In f act, coherence can be rewritten as:    󰇛  󰇜  󰇻    󰇛    󰇜  󰇡  󰇛    󰇜 󰇢    󰇻      󰇛    󰇜      󰇛    󰇜        󰈏     󰇛    󰇜    󰇛    󰇜       󰇛    󰇜  󰇡  󰇛    󰇜 󰇢     󰇛    󰇜    󰇛    󰇜  󰈏      󰇛    󰇜      󰇛    󰇜       ( 12 ) With this form ulation, it is clear that coherenc e is a weighted averag e of the unit phase vectors, i.e., a version of PLV weig hted by t he joint ampl itude of the signals at a giv en frequency. Alternatively, coherenc e can be a regarded as a version of PLV weighted by the signal- to - noise (SNR) ratio of each trial, because if the environmental conditions (noise) ar e stable, the amplitude of a signal is pr oportional to its SNR. In any case, coherence and PLV are ti ghtly related and both are different fo rmulations of the same pr inciple. It can be argued that PLV and coher ence are inherently different, as the PLV must be calculated over a whole oscillator (i.e. a frequency band) and coherence is ca lculated independently for each frequency. However, follow ing the pr evious logic, it is triv ial to understand that PLV in a given band is related to the average coherence in the frequency range encompassed by the band weighted by their relativ e amplitude. We wi ll further elabora te on this relationsh ip in the Results section. C. Zero-lag-insensitiv ity versions Despite its popularity , the PLV presents an important l imitation when used to assess functional brain connectivity: its sensitivity to v olume conduction (Stam et al. 2007) and source-leakage 7 effects. I n fact, in the case of M/EEG data at the sensor level, several sensors can simultaneously pick up the activity from the same source (volume conduction ). In turn, at the source level, due to the low spatial resolution of the data, different neighboring sources may share som e act ivity (source leakage). Luckily, due to the low capacitance of the tissues of the head for the physiological frequencie s and the sm all distance that t he currents have to travel, t he propagation of the si gnals of interest ca n be considered instantaneo us (Nolte et al. 2004; Stam et al. 2007). Under this assumption, volum e conduction/source leak age occurs with zero -lag propag ation. In other words, the phase difference of the part of the sig nals related to such s purious connecti vity must be zero. Following this logic, Cornelis Stam (Stam et al. 2007) and Guido No lte (Nolte et al . 2004) ca me with two different PS metrics that discard zero -lag connectivity and are, thus, insensitive t o volume conduction. They are the phase lag index (PLI) and the imaginary part of coherency , r espectively . Due to the t ight relation of PLV with coherency, it is possible to extend t he imaginary part of coherency to PLV to obtain a PLV -based measure insensitive to vol ume conduc tion effects, which, for symm etry, we will term im aginary PLV (iPLV ):           󰇗    󰇗     ( 13 ) where  󰇝  󰇞 stands for the im agi nary part of  . This measure is insen sitive to zero-lag effects, as i t remov es the contribu tion of the zero phase differences that, due to the complex exponentiation, gives real PLV values. However, (13) is not normalized, as its upper bound, corresponding to two signals with a phase difference  >0, is  󰇛  󰇜 . This can be corrected analogously to what (Ewald et al. 2012) did for the imaginary part of coherency , to define a correc ted imag inary PLV (ciPLV):          󰇗     󰇗        󰇡    󰇗     󰇗    󰇢  ( 14 ) This definition of c iPLV is sim ilar to the definition of lagged coherenc e as introduced by Pascual- Marqui and cowo rkers (Pascual-Marqui e t al. 2011). III. Results A. Speedup achieved us ing the proposed a lgorithm In order to study the behavior of the proposed formulati on of the PLV algorithm , we calculated PLV using three different im plementations f or differe nt data lengths and number of signals. Even 8 when th e com put ation i s c ompletely deterministic, and any data set with the sam e characteristics would produce th e same results, fo r the sake of fide lity, we used real s ource space data. Test data consisted of five minutes of eyes -closed resting-state magnetoencephal ographic (MEG ) activity acquired using a 306-sensors Elekta Vectorview system (Elekta, AB, Stockholm, Sweden), located inside a magnetically shielded room (VacuumSchmelze GmbH, Hanau, Germany) at t he Laboratory for Cognitive and Computational Neuro science (Madrid, Spain). Data was acquired with a sampling rate of 1000 Hz (anti -alias band -pass filter of 0.1 -330 Hz) and filtered using a spatiotemporal signal space separa tion method (Taulu & Simola 2006). Data were segmented in 20 or 40 four-second artifact-free segments and band-pass filtered in t he classical al pha band (8 to 12 Hz) usi ng a 2000t h order FIR fi lter in two pas ses with 2 seconds of real data as padding on both si des. The instantaneous phase of t he signal was determined using Hilbert’s analyt ical signal. In order to avoid edge effects, the analytica l signal was calculated prior to the r emov al of the padding. In order to get usable timescales for the original implementations , data was downsam pl ed by a factor of ten. We then placed 2459 dipoles inside the head of the subject, in a 1 cm ho mogeneous three - dimensional grid. Source space data was calculated using a realistic single shell as a f orward model (Nolte 2003) and a beamfor mer as inverse method (van Veen et al. 1997). The obtained spat ial filter was applied t o the sensor space data to obtain up to 2459 time series in source space. Last, all - to - al l PLV connectivity matrix was calculated for differe nt numbers of sources, ranging from 500 to 2459. PLV calculation was performed using Matlab (The Mathwo rks I nc., Natick, Massachusetts) Figure 1. Execution time ( left) and RAM use ( right) of the three evaluated PLV algorithms for a different number of signals and 40 trials. The plotted values for each data point are mean and standard deviation of the execution tim e and RAM use over 5 executions. The propo sed algorithm achieves a 100 -fold speed-up with a sma ll in crease in memory use. Please note the logarithmic scale of the y-axis. 9 and the three functions wh ose codes are in the Appendix . Time and memory perform ances were measured using Matlab’s tic/ toc and FieldTrip’s (Oostenv eld et al. 2011) memtic/m emtoc functions, respectively. Figure 1 shows the r esults of th e performance test using 40 trials and different numbers of signals. Every confi gurat ion was evaluated 5 times, and the figure shows the mean and the standard deviation. The original implementation is the slowest, but the m o st m e mory-efficient of the thre e evaluated algorithms. The optimized implementation achieves almost a 3 -fold speedup with a slight increase in m emory consumption. Howev er, the pr oposed algorithm is able to get a 100 -fold speedup over the optimized im plementation, with only a marginal (especially at a high nu mber of signals) increase in mem ory use. Results using 20 trials of data are not shown here, but the behavior was similar, with a pproximate ly half the tim e and memory consum ption. These results show that the proposed algorithm can achieve a dramatic speed -up in the calculation of t he PLV values. The algorithm even allows the ca lculation of PL V from non- downsampled data, calculating the full connectivity matrix of 2459 time -series in an average of 113 seco nds, using an average of around 8 GB of RAM memory , over 100 executions. Extrapo lating the results shown in the previous p aragraph, this execution, even w ith the op timized implementation , would take 3 hours. B. Comparison of PLV and c oherence As noted i n the Methods section, both PLV and coherence are phase-synchronization estimation algorithms. Both metrics measure similar pr operties of the data, but while coherence gives one result per frequency, PLV gives one result per frequency band. This difference aris es from the fact that coherence estimates phase synchronization from Fourier’s phase, whereas PLV estim ates it from Hilbert’s phase. In the extreme case of an infinitely narrow band, both phase definition s converge, but in a norm al case where t he band is finite, the results of coherence and PLV w ould diverge. In order t o evaluate the behavior of the coherence and the PLV in signals with a known degree of coupling, we used a pair of coupled chaotic systems: a Rössler system and a Lorenz one (Quian- Quiroga et al. 2000). In this setup, the Rössler system acts as driver, and an oscilla tion frequency can be defined from it. The slave Lorenz system is driven by t he Rössler to an extent determined by the coupling parameter C , ranging from zero (com pletely i ndependen t systems) to one. The mathem atical definition of the coup led systems is: 10 󰇱   󰇗    󰇛      󰇜   󰇗    󰇛        󰇜   󰇗    󰇛             󰇜 󰇱   󰇗    󰇛      󰇜   󰇗    󰇛                  󰇜   󰇗              ( 15 ) where subscripts 1 and 2 refer to the Rössler and t he Lorentz system, res pectively ,  is a scaling parameter det ermining the fundamental frequency of the Rössler oscillator and  is the coupling parameter. The scaling parameter was s et to 10, establishing the oscillatory f requency at 0.585· π radians per sample, and the oscillatory band wa s set to 0.570·π to 0.600·π radians per sample. With this setup, we generated 50 pairs of signals of 20,00 0 samples for different values of C ranging between zero and one. PS between the systems was finally estimated using the first variable of each chaotic system  . In order to get a band-wise coherence value, we adopted tw o different approaches using eith er the maximum or t he mean v alue of coherence over the e valuat ed band. T he values of PLV and both Figure 2. Valu es of the synchronization indices analyzed for different coup lings for a pair of Rössler and Lorenz s ystems. The dark line represents th e mean value over 50 executions, and the shadow line represents the standard deviation . (a) PLV, maxim al and ave rage coherence for the b and betwe en 0.570•π and 0 .600•π radian s per sample. Coherence was calculated u sing a Hamming window of length 40 0 with 200 samples of overlapping. (b) The same that in (a), using a window length of 1 00 sam ples with 50 samples of overlappin g. Note th at only coherence -based metrics are affected by this change. 11 coherence approache s are shown in Figure 2a. Clearly, both PLV and coherence e volve closely with increasing coupling, with coherence values slightl y overestimated for low couplings. This overestimation is due to the l ower number of phases used for the coherence calcul ation as compared to the PLV. Indeed, for the calculation of coher ence data was split in to 400 samples segments with 200 samples of overlapping , giving a total of 99 segments, and 99 independent phases. For the calculation of PL V, with a bandwidth of 0.030·π r adians per sam ple, there is thre e effective sample for every 100 sam ples, giving a total o f 600 effective s amples, and thus 600 independent phases. To address this problem, Figure 2b shows the results usi ng a smaller window for the coher ence calculation. In this case, with 100 sam ples windows and 50 sam ples overlap, coherence is calculated using 399 phases. The overestimation for low coupling i s clearl y lower than in t he previous case, but the m aximal value of coherence drops to 0.85. In this case , the reduction in the Figure 3. V alues of the synchronization indices analyzed for different couplings for on e Rössler and one Lorenz system. The dark line represents th e mean value over 50 executions, and the shadow line represents the standard deviation. (a) Synchronization estimated usin g PLV, iPLV and ciPLV. (b) The same that in (a), a fter adding a 1 0% of li near mixing to the signals. (c) Idem after adding a 20% of linear mix ing. 12 calculated synchron ization can be expl ained from the frequency sm oothing associated with a shorter time window. Due t o t he narrowness of the frequency peak of the R össler oscillator, a l arger amount of frequency sm oothing forces the coherence to be determined with a phase that include s both the peak it self and the surround ing frequencie s. The effect is similar when evaluating the mean coh erence over a band, as coupled frequencies are averaged together with non-coupled ones. Howev er, PLV does not show this effect, as the instantaneous ph ase of the si gnal is c alculated ov er the whole band, and not from an average of frequencies. In this sense, PLV is superior to coherence, as the position and narrowness of the oscillator over the ob served band are not relev ant. PLV could ev en handle an oscill ator with varying frequency, given that the oscillatory frequency never goes outside of t he observ ed b and. On the other hand, maxim al coherence would only return the coupling at the most comm on frequency, ignoring the coupling at other frequencies, and mean coherence would dramatically underestim ate the connectiv ity. C. Effects of volume con duction As commented above, one of the main criticisms to PLV and coherence as applied to brain signals is their sensitiv ity to volume conduction. It can be modeled as an instantaneous pro jection of one signal onto the other, giving rise t o zero-lag synchronization. As shown in the Methods section, Guido Nolte’s proposed imag inary part of coherency (Ewald et al. 2012; Nolte et al. 2004) can be extended to the PLV, giving a set of zero-lag insensitive measures. In order to evaluate the behavio r of these PL V derived metrics, we u sed the sam e pair of chaotic system s defined in th e previous section. T he systems show nonzero -lag phase synchroniz ation, and the volum e conduction can be introduced using ins tantaneous linear m ixing (Porz et al. 2014 ; Haufe et al. 2013) :                 ( 16 ) where  is a p arameter determ ining the amount of mixing. Figure 3 shows the synchronization estimated between t he signals described in the previous section using the PL V, iPLV and ciPLV, and values of  of 0, 0.1 and 0.2. Figure 3a shows t he estimated synchronization for V=0 (without zero lag mixing), measured with PLV, iPLV and ciPLV. Values of PLV and ciPLV are almost i dentical, indicating that the synchronization of the pair of the chaotic oscillator is nonzero lag. However, iPLV values ar e slightly lower, especially for high couplings, indicating that the lag bet ween both si gnals is almost, but not exactly , one-quarter of a cy cle, so the PLV com plex vector has on ly a small re al component. 13 Figure 3b, corresponding to V=0.1, shows the effect of 10% volume conduc tion between the signals, introducing low spurious zero-lag synchronization. PLV values increase especially at lower couplings, whereas iPLV and ciPLV values remain unaffected. iPLV shows lower synchronizat ion values that ciPLV, w hich remains alm ost identical to the V=0 cas e. This is bec ause the complex PLV vector now has an important real component that is completely ignored by i PLV. This shows that ciPLV is superior, as it uses this real component to norm alize the im aginary one. Figure 3 c shows the effect of 20% volume conduction, with sim ilar results, but exacerbating t he estimation errors in PLV an d iPLV. ciPLV, however, continues to ext ract the correc t synchronization value . IV. Discussi on The formulation of PLV described in this paper allows for a speedup of two orders of mag nitude in the calculation of this index. This result is espec ially relevant in the study of EEG/MEG source- space connectiv ity. Source space models based in volum es, as the one used in the Results section above, include around 2500 sources for t he whole brain, and around 1500 for gray matter. In the case of sources based i n surfaces, which are norm ally used in m inimum norm est imates (Pascual- Marqui et al. 1994; Gram fort et al. 2013), t he num ber of sources ranges bet ween 1000 and 10000 per hemisphere. In both cases, a whole brain connectivity study, even with subsam pling, would take up to several hou rs per subject (and band). The usual approach in those ca ses is the parcellation of the brain and the ext raction of a representative tim e cour se per parcel (Farahibozorg et al. 2017; Korhonen et al. 2014) . This reduces the num ber of time courses to the order of the hundreds, thereby allowing the estimation of whole- brain connectivity in a reasonable am ount of time. However, this method forces the defini tion of homogeneous parcels, defined by one unique time ser ies, thus discarding som e information associated to inter -area variability. With the proposed algorithm the calculation of source - to -sourc e PLV would take only a few m inutes per subject (and band), allowing the estimation of whole -brain connectivity for a whole study in a matter of hour s. In additi on, our reform ulation of the original PLV definition stressed the similarities between PLV and coherence. T he com parison of both metrics as applied to a pair of chaotic systems with different coupling levels showed t hat PLV and coherence values are closely rel ated, with PLV showing a better behavior in the (very realistic) situation where the couplin g tak es place over a whole (possibly narrow) band instead of at a single frequency. This is, howev er, not a flaw of coherence, but t he direct effect of its definition. PLV seems best to evaluat e synchronization over a whole band, wherea s coherence does to ev aluate it at fi xed frequencies. 14 Finally, we evaluated the behavior of PLV in the relevant case of volum e conduction/sou rce leakage as simulated by instantaneous linear mixing of the signals corresponding to the t wo systems. In analogy to co herency, we introduc ed the i mag inary part of PLV and i ts corrected version. Both algorithms proved insensitive to volume conduction, but iPLV was flawed by the effect of a real component of the complex PLV vector. As for the imaginary part of coherency , iPLV only reaches t he maximal value when the phas e difference of the two si gnals is exactly π/2 radians, whereas it dr ops to almost zero when the phase difference is small but consistently nonzero. ciPLV corrects this behavior in a similar way to the corrected imaginary part of coherency (Ewald et a l. 2012) , being both un biased and insensitive to volume c onduction/source leakage effects. This is not the fi rst time that a metric based on the im aginary part of PLV is used (Dimitriadis et al. 2017; Palva et al. 2017; Wang et al. 2018) , but, to our knowledge, the expected behavior had not be en thoroughly tested yet. Both sensor and source space EEG/MEG connectivity have to deal with the burden of zero - lag connectivity . Zero lag connectivity between t wo sensors is usually interpreted as the effect of a common source i n MEG and an effect of volum e conduction in EEG. At the source level, it is usually interpret ed as an effect of sour ce leakage. B oth interpretations a re based in the idea th at brain signals cannot travel instantaneously between different parts of the brain, and any real connectivity should imply a delay (Stam et al. 2007). Even when some scenarios allow zero-lag connectivity in the brain (Kovach 2017), it is not possib le t o distinguish real zero-lag connectivity from volum e c onduction or source leak age. Several measures have been developed t o try to overcome t his problem. T he most notable of these measures is the imag inary part of coherence, and its corrected counterpart. Howev er, even when coherence is effectiv ely a phase-synchronizatio n measure, t he results do not completely fit in the phase synch ronization hy pothesis in the brain. P LI is another m easure developed t o estim ate phase connectivity ignoring the contribution of zero lag (Stam et al. 2007), but the metric ha s showed a low tes t-retest reliability (Gar cés et al. 2016 ; Colclough et al. 20 16). The variations of PLV proposed i n this paper mimic their coherence counterparts, and have proved effective t o remove zero-lag connectivity while keeping intact nonzero-lag connectivity. This allows the s tudy of sensor or source space E EG/MEG connec tivity e valuating only real synchronizations, and completely ignoring volume conduction or source leakag e ghost synchronizations. This method, as di scussed, is not perfect (Kovach 2017), but in the same fashion that imaginary coherency (Nolte et al. 2004), allows the evaluation of real connectiv ity without the influence of z ero-lag interference. 15 Note, however, that i n any case, and according to the latest results in t he literature (Palva et al., 2017; Wang et al., 2017) no bivariate FC index (whether zero lag insensitive or not) is free from the effect of spurious detection of connectivity due t o source leakage. Nevertheless, for the detection of tru e connectivity, we believ e that both iPLV and especia lly ciPLV are excellent choices that can be very efficiently es timated u sing the algorithm s we introduced here. V. Conclusion In this paper, we p ropose a new form ulation of the PLV a lgorithm that allow s its fast computa tion by bypassing CPU-demanding calculations. The direct use of t his new algorithm in Matlab allows for a speedup by a f actor of 100 from a vector-optimized implementation. Ev en though this implem entation can be improv ed, it is un likely that a simila r speedup can be achieved using the original form ulation. Noteworthy, the optimization is only based in a mathem atical reform ul ation of the original definition of PLV, in which the phase extraction and th e exponentiation necessary for its calculatio n are replaced by a simple matrix multiplication. Besides, the algorithm can be implemented in any language, and the code prov ided in the Appendix , created in Matlab, is com pletely machine- independent. Thus, differently to a toolbox we recently released (García- Prieto et al. 2017), t here is no need to compile the source code in C language, which mak es the present formulation much easier to use. Note also that this paper only deals wi th the mathematical defin ition, and there m ay be room for com putational improv ement. In conclusion, we have show n that, despite its longevity, it is still possible to refine further th e estimation of the (allegedly) m ost popular method of PS, the PLV index. T he new formulation no t only al lows for a much faster estimation of the index but also, by highlighting the similarities between PLV and coheren ce, prompts the definition of two new measures insensitiv e to zero lag. These measures, derived from the vector product of the real and imag inary part of PLV, ar e analog to the (co rrected) im aginary part of coherence and are therefore robust again volum e conductio n and source leakage effects. We hope these new developm ents will encourage furt her the use PLV- related m easures in the analy sis of neurophysiolog ical data at the sensor and the s ource level. Appendix In the R esults section, we t ested three different implem entations for the calculation of PLV. In this Appendix, we show the three different codes developed in Matlab to evaluate the behavior of each algorithm. In t he codes pr ovided, all three implem entations are fed with the Hilbert analytic 16 signal, with shap e signals per sam ples per trials. The first implementation includes two nested l oops, so every pair of si gnals is evaluated independently . This option is the most memory conservative but also the slowest one, as only two signals are ev aluated in each inter action. The code is: [ nc, ns, nt ] = size( data ); phs = angle( data ); plv = zeros( nc, nc, nt ); for t = 1: nt for c1 = 1: n c for c2 = c1 + 1: n c dphs = phs( c1, :, t ) - phs( c2, :, t ); plv( c1, c2, t ) = abs( mean( exp( 1i * dphs ) ) ); end end end The se cond implementation is a vectorized and memory effi cient version of the algorithm. In this implem entation every signal is compared against all the other signals, avoiding t wo loops. This option uses slightly more memory but achieves a speedup of factor 2.5 over the ori ginal one . The code is: [ nc, ns, nt ] = size( data ); phs = angle( data ); plv = zeros( nc, nc, nt ); for t = 1: nt tplv = complex( zeros( nc ) ); for s = 1: ns dphs = bsxfun( @minus, phs( :, s, t ), phs( :, s, t )' ); tplv = tplv + exp( 1i * dphs ); end plv( :, :, t ) = abs( tplv / ns ); end The last implementation uses the proposed alg orithm. This algorithm allows the direct calculation of all source pairs at once, with negligible memory in crease. In addition to avoiding the angle an d exp functions, th e implementation removes the need to use a loop ov er the signals. Th e code is: [ n c , ns, nt ] = size( data ); ndat = data ./ abs( data ); plv = zeros( n c , n c , nt ); for t = 1: nt 17 plv( :, :, t ) = abs( ndat( :, :, t ) * ndat( :, :, t )' ) / ns; end Acknowledgment The authors would like to thank X. Bello for her review of the style of this manuscript. This study was partially supported by one project from t he Spanish Ministry of Economy and Competitiveness (TEC2016- 80063- C3 -2-R) to EP and one pre-doctoral fellowship from the Spanish Ministry of Education (FPU13 /06009) to RB. 18 References Le Bihan, D. et al. , 2001 . Diffusion tensor imagin g: concepts and app lications. Jou rnal of mag netic resonan ce imaging : JMRI , 13(4), pp .534 – 546. Brookes, M.J. et al. , 2011. Measuring functional connectivit y using ME G: Me thodology a nd co mparison with fcMRI. NeuroIma ge , 56, pp. 1082 – 11 04. Bruns, A., 2004. Fourier -, Hilbert- and wavelet-based sign al analysis: Are they rea lly different app roaches? Journal o f Neuroscience Methods , 137(2), pp.321 – 332. Colclough, G.L. et al., 20 16. How reliable are MEG resti ng -state con nectivity metrics? Neu roImage , 138, pp.284 – 293. Dimitriadis, S.I . et al., 2017. T opological Filtering o f Dynamic Functional Brain Networks Unfolds Informative Chron nectomics: A Novel Data -Driven Thresholding Scheme Based on Orthogonal Minimal Spanning T rees (OMST s). Frontiers in Neuroin formatics , 11 . Ewald, A. et al., 2012 . Estimating true brain connectivity from EEG/MEG data invariant to linear and static transformations in sensor space. NeuroImage , 60(1), pp.476 – 488. Farahibozorg, S.-R., Henso n, R.N. & Hau k, O., 2017. Adaptive cor tical parcellations for source reconstr ucted EEG/MEG connecto mes. Neuro Image . Fries, P., 2015 . Rhythms for Cognition: Com munication through Co herence. Neu ron , 8 8(1), pp.220 – 235. Friston, K.J., 1994. Functional and effective connectivity in neuroimaging : A s ynthesis. Human Brain Mapping , 2(1 – 2), pp.56 – 78. Friston, K.J., Harr ison, L. & Penny, W., 2 003. Dynamic causal modelling. NeuroImage , 19(4 ), pp.1273 – 1302. Garcés, P., Martín-Buro, M.C. & Maestú, F., 2016. Quantifying the T est -Retest Reliability of Mag netoencephalograph y Resting-State F unctional Connecti vity. Brain Connectivity , 6(6), pp. 448 – 460. García-Prieto, J ., Bajo, R. & Pereda, E., 2017. Efficient co mputation of f unctional brain net works: to wards real-time functional co nnectivity. Fron tiers in Neu roinformatics , 11 , p.8. Gramfort, A. et al., 2013. MEG and EEG data analysis with MNE -Python. Fron tiers in Neuroscience , (7 DEC). Haufe, S. et al. , 2013. A critical assessment of connectivit y measure s for EEG data: A si mulation stud y. NeuroImage , 64(1), pp .120 – 133. Korhonen, O., Palva, S. & Palva, J.M., 2014. Spar se weightings for collap sing inverse solutions to cortical parcellations opti mize M/EEG source reco nstruction accurac y. Jo urnal of Neurosc ience Meth ods , 226, pp.147 – 160. Kovach, C.K., 201 7. A B iased Look at Phase Locki ng: Brief Critical Review and P roposed Remedy. IE EE Transactions on Signal Processing , 65(17), pp.4468 – 4 480. Lachaux, J.P . et al., 1999. Measuring p hase synchrony in brain signals. Human brain mapping , 8, pp. 194 – 208. López, M.E. et al., 2014. Alpha-band hypersynchronization in p rogressive mild cognitive impair ment: a 19 magnetoencephalograp hy study. The Jo urnal of neuroscience : the official jou rnal of the Society for Neuroscience , 34 (44), pp.14551 – 9. Mormann, F. et al., 200 0. Mean pha se coherence as a measur e for pha se synchronization a nd its app lication to the EEG of epilepsy patient s. Ph ysica D: Nonlin ear Phenomena , 144(3), pp.3 58 – 369 . Nolte, G. et al., 2 004. I dentifying true b rain interac tion fro m EEG d ata usi ng the imaginary part of coherenc y. Clinical Neurophysiolo gy , 115(10), pp.2292 – 2307. Nolte, G., 2003. T he m agn etic lead field theore m in the quasi -static app roximation and it s use for magnetoencephalograp hy forward calculation in realistic volume conductor s. Physics in medicine a nd biology , 48(22), pp.3637 – 3652. Nunez, P.L. et al., 1997 . EEG coherency. E lectroencephalog raphy and Clinical Neurophysiology , 103(5), pp.499 – 515. Oostenveld, R. et al., 201 1. FieldT rip: Open source soft ware for advanced analysis o f MEG, EEG, a nd invasive electrophysiological data. Computa tional intelligence and neuroscience , 2011, p.15 6869. Palva, J.M. et al., 2 017. Ghost interactions in MEG/EE G source spac e: A note of caution on i nter -areal coupling measures. bioRxiv . Pascual-Marqui, R.D. et a l., 2 011. Assessing interact ions in t he brain with exact low -resolution electromagnetic to mography. Philo sophical Tran sactions o f the R oyal Society A: Mathema tical, Physical and Engineering Sciences , 369(1952) , pp.3768 – 3784. Pascual-Marqui, R.D., Michel, C.M. & Lehmann, D., 1 994. Lo w r esolution electro magnetic tomography: a new method f or localizing electrical activity in the brain. Internationa l Jou rnal of Psychophysiology , 18(1), pp.49 – 65. Porz, S., Kiel, M. & Lehn ertz, K., 2014. Can spurious indicatio ns for phase s ynchronizatio n due to superimposed signals be avoid ed? Chao s , 24(3). Quian-Quiroga, R., Arnhold, J. & Grassberger, P ., 2000. Learning driver -response relationships fro m synchronization patterns. P hysical review. E, Sta tistical physi cs, plasmas, fluid s, and related interdisciplinary top ics , 61( 5 Pt A), pp.5142 – 5148. Stam, C.J., Nolte, G. & Da ffertshofer, A., 2 007. P hase lag i ndex: Assessme nt o f functio nal connectivit y from multi cha nnel EEG and MEG with di minished bias fro m co mmon sou rces. Hu man Brain Map ping , 28(11), pp.1178 – 1193. Taulu, S. & Simola, J., 200 6. Spatiote mporal signal space separ ation method for rejecting nearby inter ference in MEG measure ments. Physics in medicine an d biology , 51, pp.1759 – 1 768. van Veen, B.D. et al., 1997. Localization of br ain electrical ac tivity via linearly constrained minimum variance spatial filtering. IEEE tran sactions on bio-med ical engineering , 44, pp.867 – 8 80. Wang, S.H. et al., 2018 . Hyperedge bundling: A practica l solution to spurious i nteracti ons i n ME G/EEG source connectivity anal yses. Neu roImage .