Phase Locking

Calculating Phase locking with bmtool¶

By Gregory Glickert

In [1]:

RunningInCOLAB = 'google.colab' in str(get_ipython())
if RunningInCOLAB:
    %pip install bmtool &> /dev/null

RunningInCOLAB = 'google.colab' in str(get_ipython()) if RunningInCOLAB: %pip install bmtool &> /dev/null

In [2]:

import numpy as np
from scipy import signal as ss
import matplotlib.pyplot as plt
from bmtool.analysis.lfp import fit_fooof
from bmtool.analysis.entrainment import calculate_spike_lfp_plv, calculate_ppc, calculate_ppc2
import warnings
warnings.filterwarnings("ignore")

import numpy as np from scipy import signal as ss import matplotlib.pyplot as plt from bmtool.analysis.lfp import fit_fooof from bmtool.analysis.entrainment import calculate_spike_lfp_plv, calculate_ppc, calculate_ppc2 import warnings warnings.filterwarnings("ignore")

Warning: no DISPLAY environment variable.
--No graphics will be displayed.

In [3]:

np.random.seed(9) # lucky number 9
fs = 1000  # 1 kHz sampling rate
duration = 200  # 200 seconds
t = np.arange(0, duration, 1/fs)

# Define oscillation frequencies
beta_freq = 15  # Hz
gamma_freq = 40  # Hz

# Create a simulated LFP with multiple frequency components
lfp = (0.5 * np.sin(2 * np.pi * gamma_freq * t) + 0.2 * np.sin(2 * np.pi * beta_freq * t) +  0.3 * np.random.randn(len(t)))

# Generate phase information for each frequency
beta_phase = np.angle(ss.hilbert(np.sin(2 * np.pi * beta_freq * t)))
gamma_phase = np.angle(ss.hilbert(np.sin(2 * np.pi * gamma_freq * t)))

# Generate beta-synchronized spikes
# Strong preference for specific beta phase, weak response to gamma
beta_spike_probability = (0.8 * (1 + np.cos(beta_phase - np.pi/4)) + 
                            0.1 * (1 + np.cos(gamma_phase - np.pi/3)))

beta_spike_train = np.random.rand(len(t)) < beta_spike_probability * 0.02
beta_spike_times = t[beta_spike_train]

# Generate gamma-synchronized spikes
# Strong preference for specific gamma phase, weak response to beta
gamma_spike_probability = (0.1 * (1 + np.cos(beta_phase - np.pi/4)) + 
                            0.8 * (1 + np.cos(gamma_phase - np.pi/3)))

gamma_spike_train = np.random.rand(len(t)) < gamma_spike_probability * 0.02
gamma_spike_times = t[gamma_spike_train]

print(f"Generated {len(beta_spike_times)} beta-synchronized spikes")
print(f"Generated {len(gamma_spike_times)} gamma-synchronized spikes")

np.random.seed(9) # lucky number 9 fs = 1000 # 1 kHz sampling rate duration = 200 # 200 seconds t = np.arange(0, duration, 1/fs) # Define oscillation frequencies beta_freq = 15 # Hz gamma_freq = 40 # Hz # Create a simulated LFP with multiple frequency components lfp = (0.5 * np.sin(2 * np.pi * gamma_freq * t) + 0.2 * np.sin(2 * np.pi * beta_freq * t) + 0.3 * np.random.randn(len(t))) # Generate phase information for each frequency beta_phase = np.angle(ss.hilbert(np.sin(2 * np.pi * beta_freq * t))) gamma_phase = np.angle(ss.hilbert(np.sin(2 * np.pi * gamma_freq * t))) # Generate beta-synchronized spikes # Strong preference for specific beta phase, weak response to gamma beta_spike_probability = (0.8 * (1 + np.cos(beta_phase - np.pi/4)) + 0.1 * (1 + np.cos(gamma_phase - np.pi/3))) beta_spike_train = np.random.rand(len(t)) < beta_spike_probability * 0.02 beta_spike_times = t[beta_spike_train] # Generate gamma-synchronized spikes # Strong preference for specific gamma phase, weak response to beta gamma_spike_probability = (0.1 * (1 + np.cos(beta_phase - np.pi/4)) + 0.8 * (1 + np.cos(gamma_phase - np.pi/3))) gamma_spike_train = np.random.rand(len(t)) < gamma_spike_probability * 0.02 gamma_spike_times = t[gamma_spike_train] print(f"Generated {len(beta_spike_times)} beta-synchronized spikes") print(f"Generated {len(gamma_spike_times)} gamma-synchronized spikes")

Generated 3596 beta-synchronized spikes
Generated 3661 gamma-synchronized spikes

We can see that we generated an LFP with two different frequencies¶

In [4]:

hz,pxx = ss.welch(lfp,fs)
_,_ = fit_fooof(hz,pxx,plot=True,report=True,plt_range=(1,100),freq_range=[1,100])

hz,pxx = ss.welch(lfp,fs) _,_ = fit_fooof(hz,pxx,plot=True,report=True,plt_range=(1,100),freq_range=[1,100])

==================================================================================================
                                                                                                  
                                   FOOOF - POWER SPECTRUM MODEL                                   
                                                                                                  
                        The model was run on the frequency range 3 - 98 Hz                        
                                 Frequency Resolution is 3.91 Hz                                  
                                                                                                  
                            Aperiodic Parameters (offset, exponent):                              
                                         -3.8986, -0.0807                                         
                                                                                                  
                                       2 peaks were found:                                        
                                CF:  15.00, PW:  1.239, BW:  7.81                                 
                                CF:  39.85, PW:  2.120, BW:  7.81                                 
                                                                                                  
                                     Goodness of fit metrics:                                     
                                    R^2 of model fit is 0.9808                                    
                                    Error of the fit is 0.0539                                    
                                                                                                  
==================================================================================================

Now we can see how our spike times are entrained to the lfp we made. There are a few different metrics that we can use.¶

1. Phase-Locking Value (PLV) - calculate_spike_lfp_plv¶

Equation:

Unbiased Phase Locking Value¶

The unbiased Phase Locking Value is calculated as follows:

1. Standard PLV calculation:

   $$PLV = \left| \frac{1}{N} \sum_{j=1}^{N} e^{i\phi_j} \right|$$

   where $\phi_j$ is the phase at each spike time.

2. Pairwise Phase Consistency (PPC) calculation:

   $$PPC = \frac{PLV^2 \cdot N - 1}{N - 1}$$

3. Unbiased PLV calculation:

   $$PLV_{unbiased} = \sqrt{\max(PPC, 0)}$$

2. Pairwise Phase Consistency (PPC) - calculate_ppc¶

Equation: $$ \text{PPC} = \frac{2}{n(n-1)} \sum_{j=1}^{n-1} \sum_{k=j+1}^{n} \cos(\phi_j - \phi_k) $$

Where:

• $\phi_j,\phi_k$ = Phases at spike times $j$ and $k$
• $n$ = Number of spikes

This calculates phase consistency through pairwise comparisons of spike-phase differences.

---

3. Optimized Pairwise Phase Consistency (PPC2) - calculate_ppc2¶

Derivation of PPC2 from PPC¶

1. Original PPC¶

$$ \text{PPC} = \frac{2}{n(n-1)} \sum_{i < j} \cos(\phi_i - \phi_j) $$

2. Rewrite $\cos(\phi_i - \phi_j)$ Using Euler’s Formula¶

$$ \cos(\phi_i - \phi_j) = \text{Re}\left( e^{i(\phi_i - \phi_j)} \right) $$

3. Sum Over All Pairs¶

$$ \sum_{i < j} \cos(\phi_i - \phi_j) = \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \cos(\phi_i - \phi_j) - \frac{n}{2} $$

4. Express in Terms of Complex Exponentials¶

$$ \sum_{i < j} \cos(\phi_i - \phi_j) = \frac{1}{2} \text{Re} \left( \left| \sum_{j=1}^n e^{i\phi_j} \right|^2 \right) - \frac{n}{2} $$

5. Final PPC2 Formula¶

$$ \text{PPC2} = \frac{\left|\sum_{j=1}^{n} e^{i\phi_j}\right|^2 - n}{n(n-1)} $$

This algebraic equivalent of PPC provides computational optimization by:

1. Converting phases to complex unit vectors
2. Using vector magnitude properties to avoid explicit pairwise comparisons

All three metrics quantify phase locking between spike times and LFP oscillations

• all metrics should be the same with enough spike times

In [ ]:

# spikes were generated on second time scale so spike_fs=1
# unbias plv is the same as the square root of ppc
print(calculate_spike_lfp_plv(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40))
# 3 different methods for ppc all give the same answer
print(np.sqrt(calculate_ppc(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40,ppc_method='numpy')))
print(np.sqrt(calculate_ppc(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40,ppc_method='numba')))
# can be annoying to install correctly to work with GPU i know that it works install cudatoolkit in conda if you wanna try
try:
  print(np.sqrt(calculate_ppc(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40,ppc_method='gpu')))
except:
  print("GPU not working passing for now")
print(np.sqrt(calculate_ppc2(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40)))

# spikes were generated on second time scale so spike_fs=1 # unbias plv is the same as the square root of ppc print(calculate_spike_lfp_plv(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40)) # 3 different methods for ppc all give the same answer print(np.sqrt(calculate_ppc(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40,ppc_method='numpy'))) print(np.sqrt(calculate_ppc(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40,ppc_method='numba'))) # can be annoying to install correctly to work with GPU i know that it works install cudatoolkit in conda if you wanna try try: print(np.sqrt(calculate_ppc(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40,ppc_method='gpu'))) except: print("GPU not working passing for now") print(np.sqrt(calculate_ppc2(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=40)))

0.42532179470701004
0.4253217947070101

0.4253217947070101
0.4253217947070101
0.4253217947070101

We can also check the entrainment at different frequencies¶

In [6]:

freqs = [1,5,15,25,30,40,50,60,70,80]
ppc_gamma = []
ppc_beta = []

for i in range(len(freqs)):
    ppc_gamma.append(calculate_ppc2(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=freqs[i]))
    ppc_beta.append(calculate_ppc2(spike_times=beta_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=freqs[i]))
    
plt.plot(freqs,ppc_gamma,'o-', linewidth=2,label="FSI")
plt.plot(freqs,ppc_beta,'o-', linewidth=2,label='LTS')
plt.legend()
plt.show()

freqs = [1,5,15,25,30,40,50,60,70,80] ppc_gamma = [] ppc_beta = [] for i in range(len(freqs)): ppc_gamma.append(calculate_ppc2(spike_times=gamma_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=freqs[i])) ppc_beta.append(calculate_ppc2(spike_times=beta_spike_times, lfp_data=lfp,spike_fs=1, lfp_fs=fs, filter_method='wavelet',freq_of_interest=freqs[i])) plt.plot(freqs,ppc_gamma,'o-', linewidth=2,label="FSI") plt.plot(freqs,ppc_beta,'o-', linewidth=2,label='LTS') plt.legend() plt.show()