Reconciling contradictory models of subthalamic nucleus contributions to basal ganglia beta oscillations.

Code ↔ Paper

12 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 12 matches

[1] § Models and methods › Population firing rate model ↔ 01_rate_model/figures_in_paper.ipynb, lines 298–394 · score 0.87 · fast spiking interneurons, subthalamic nucleus, prototypical neurons, population firing rate, receptors, spiny
[2] § Results › Full neuronal network model ↔ 01_rate_model/figures_in_paper.ipynb, lines 298–394 · score 0.85 · fast spiking interneurons, subthalamic nucleus, prototypical neurons, basal ganglia, firing rate models, spiny
[3] § Results › Full neuronal network model › Results on the two-loop spiking network. ↔ 03_spiking_networks/figures_in_paper.ipynb, lines 178–257 · score 0.78 · Error bars, EIF het, QIF het, Proto STN, Beta power, Beta frequency
[4] § Models and methods › Multi-population basal ganglia model › Beta power quantification. ↔ 03_spiking_networks/utils.py, lines 157–187 · score 0.71 · 12.5–30 Hz, Beta power, beta frequency, 12.5 Hz, welch, nperseg
[5] § Results › Full neuronal network model › Results on the two-loop spiking network. ↔ 03_spiking_networks/figures_in_paper.ipynb, lines 121–159 · score 0.69 · EIF het, QIF het, EIF STN, Proto STN, parameter sweeps, network
[6] § Models and methods › Phase difference analysis ↔ 01_rate_model/figures_in_paper.ipynb, lines 396–458 · score 0.64 · dominant frequency, power spectrum, fft, wrapped, firing rate, signals
[7] § Models and methods › Phase difference analysis ↔ 03_spiking_networks/utils.py, lines 283–344 · score 0.63 · dominant frequency, power spectrum, fft, wrapped, firing rate, signals
[8] § Models and methods › Multi-population basal ganglia model › Exponential integrate-and-fire (EIF) model. ↔ 02_tts_and_phase/tts_utils.py, lines 215–241 · score 0.61 · slope factor, exponential term, threshold potential, reset, EIF, fires
[9] § Models and methods › Multi-population basal ganglia model › Power spectral density. ↔ 03_spiking_networks/utils.py, lines 119–155 · score 0.60 · Power spectra, bin width, Welch, segment
[10] § Models and methods › Multi-population basal ganglia model › Average beta frequency. ↔ 03_spiking_networks/utils.py, lines 157–187 · score 0.55 · 12.5–30 Hz, beta frequency, 12.5 Hz, power, population
[11] § Results › Full neuronal network model › Results on the two-loop spiking network. ↔ 03_spiking_networks/sweep_utils.py, lines 184–232 · score 0.55 · multiple seeds, standard error, Beta power, firing rate, network, spiking
[12] § Results › Full neuronal network model › Results on the two-loop spiking network. ↔ 03_spiking_networks/utils.py, lines 415–499 · score 0.52 · Proto firing rates, STN Proto, spiking networks, spectrogram, activation, connection

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 598 lines · 23 KB · no license · 5 matches

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import welch
from scipy import signal
from network_simulation import NetworkSimulation
import yaml
import os
from scipy.fft import fft, fftfreq
class NetworkVisualization:
def __init__(self, simulation):
"""
Initialize the visualization class with a reference to the simulation.
Args:
simulation (EINetworkSimulation): The simulation object containing the data to visualize.
"""
self.simulation = simulation
def calculate_firing_rate(self, pop, bin_width=10, normalize=False, drop_initial=0):
"""
Calculate the firing rate over time for a population.
Args:
pop (str): The name of the population to plot ('exc' or 'inh').
bin_width (int): The width of each bin in milliseconds.
normalize (bool): Whether to normalize the firing rate.
drop_initial (int): The initial portion of data to drop in milliseconds.
Returns:
tuple: A tuple containing the bin edges and the firing rate.
"""
# Get the spike train from the simulation results
spike_train = self.simulation.get_spike_trains()[pop]
# Concatenate all spike times for the population
all_spikes = np.concatenate(spike_train)
if drop_initial > 0:
# Drop the initial portion of data
all_spikes = all_spikes[all_spikes >= drop_initial]
# Define the bins such that each bin edge represents the end of the look-back period
bins = np.arange(drop_initial, self.simulation.T, bin_width)
else:
# Define the bins such that each bin edge represents the end of the look-back period
bins = np.arange(0, self.simulation.T, bin_width)
# Calculate the histogram
spike_counts, _ = np.histogram(all_spikes, bins=bins)
# Calculate the firing rate (spikes per second)
firing_rate = spike_counts / (bin_width / 1000.0 * len(spike_train))
# Normalize if required
if normalize:
firing_rate = firing_rate / firing_rate.max()
return bins[1:], firing_rate
def plot_firing_rate(self, pop, bin_width=10, ax=None, color='blue', normalize=False, drop_initial=0, label=None):
"""
Plot the firing rate over time for a population.
Args:
pop (str): The name of the population to plot ('exc' or 'inh').
bin_width (int): The width of each bin in milliseconds.
ax (matplotlib.axes.Axes): The axes to plot on.
color (str): The color of the plot line.
normalize (bool): Whether to normalize the firing rate.
drop_initial (int): The initial portion of data to drop in milliseconds.
"""
if ax is None:
fig, ax = plt.subplots()
# Calculate the firing rate
bin_right, firing_rate = self.calculate_firing_rate(pop, bin_width, normalize, drop_initial)
if label is None:
label = f'{pop} Firing Rate'
# Plot the firing rate
ax.plot(bin_right, firing_rate, label=label, color=color)
# Set the title and labels
ax.set_xlabel("Time (ms)")
ax.set_ylabel("Normalized" if normalize else "Firing rate (Hz)")
if ax is None:
plt.show()
def plot_raster(self, pop):
"""
Plot a raster plot for the specified population.
Args:
pop (str): The name of the population to plot ('exc' or 'inh').
"""
spike_trains = self.simulation.get_results()[f'{pop}_spikes']
# Set up the figure
fig, ax = plt.subplots()
# Loop over each neuron and plot its spike times
for i, spike_train in enumerate(spike_trains):
# i is the neuron index, spike_train is the list of firing times
ax.scatter(spike_train, [i] * len(spike_train), marker='|', color='black')
# Label the axes
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Neuron Index')
ax.set_title(f'Raster Plot of {pop.capitalize()} Neuron Spike Trains')
# Set the y-axis limits to ensure each neuron has its own space
ax.set_ylim(-0.5, len(spike_trains) - 0.5)
# Show the plot
plt.show()
def plot_power_spectrum(self, pop, bin_width=10, nperseg=256, ax=None, normalize=False, drop_initial=0, legend=None):
"""
Plot the power spectrum of the firing rate for a given population.
Args:
pop (str): The name of the population to analyze ('exc' or 'inh').
bin_width (int): The width of each bin in milliseconds.
nperseg (int): Length of each segment for the Welch method.
ax (matplotlib.axes.Axes): The axes to plot on.
normalize (bool): Whether to normalize the firing rate.
drop_initial (int): The initial portion of data to drop in milliseconds.
"""
if ax is None:
fig, ax = plt.subplots()
# Calculate the firing rate
bin_right, firing_rate = self.calculate_firing_rate(pop, bin_width, normalize, drop_initial)
# Calculate the sampling frequency from bin_right
fs = 1000.0 / (bin_right[1] - bin_right[0]) # Convert bin width from ms to Hz
# Calculate the power spectrum using Welch's method
freqs, power = welch(firing_rate, fs=fs, nperseg=nperseg, nfft=nperseg)
# Plot the power spectrum
if legend is not None:
ax.plot(freqs, power, label=legend)
else:
ax.plot(freqs, power)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power')
ax.set_title(f'Power Spectrum of {pop} Firing Rate')
ax.set_xlim(0, 100)
if ax is None:
plt.show()
def calculate_beta_power(self, pop, bin_width=10, nperseg=256, normalize=False, drop_initial=0):
"""
Calculate the beta power (average power between 12.5 and 30 Hz) of the firing rate for a given population.
Args:
pop (str): The name of the population to analyze ('exc' or 'inh').
bin_width (int): The width of each bin in milliseconds.
nperseg (int): Length of each segment for the Welch method.
normalize (bool): Whether to normalize the firing rate.
drop_initial (int): The initial portion of data to drop in milliseconds.
Returns:
float: The average power in the beta frequency range (12.5-30 Hz).
"""
# Calculate the firing rate
bin_right, firing_rate = self.calculate_firing_rate(pop, bin_width, normalize, drop_initial)
# Calculate the sampling frequency from bin_right
fs = 1000.0 / (bin_right[1] - bin_right[0]) # Convert bin width from ms to Hz
# Calculate the power spectrum using Welch's method
freqs, power = welch(firing_rate, fs=fs, nperseg=nperseg, nfft=nperseg)
# Find the indices of the frequencies in the beta range (12.5-30 Hz)
beta_indices = np.where((freqs >= 12.5) & (freqs <= 30))[0]
# Calculate the average power in the beta range
ave_beta_freq = np.sum(freqs[beta_indices]*power[beta_indices])/np.sum(power[beta_indices])
beta_power = np.mean(power[beta_indices])
return ave_beta_freq, beta_power
def compute_spectrogram(self, pop, bin_width=10, nperseg=256, noverlap=128,
demean=True, drop_initial=0):
"""
Compute the spectrogram of the firing rate for a given population.
Args:
pop (str): The name of the population to analyze.
bin_width (int): The width of each bin in milliseconds.
nperseg (int): Length of each segment for the spectrogram.
noverlap (int): Overlap between segments.
demean (bool): Whether to subtract the mean from firing rate.
drop_initial (int): The initial portion of data to drop in milliseconds.
Returns:
tuple: (frequencies, times, Sxx, fs, times_ms)
"""
# Calculate the firing rate
bin_right, firing_rate = self.calculate_firing_rate(pop, bin_width, False, drop_initial)
# Z-score the firing rate if requested
if demean:
firing_rate = firing_rate - np.mean(firing_rate)
# Calculate the sampling frequency from bin_right
fs = 1000.0 / (bin_right[1] - bin_right[0]) # Convert bin width from ms to Hz
# Calculate the spectrogram
frequencies, times, Sxx = signal.spectrogram(firing_rate, fs=fs,
nperseg=nperseg,
noverlap=noverlap,
scaling='spectrum')
# Calculate time offset for proper alignment
time_offset = drop_initial + (bin_width / 2) # Add half bin width for center alignment
times_ms = times * 1000 + time_offset # Convert to ms and apply offset
return frequencies, times, Sxx, fs, times_ms
def plot_spectrogram(self, pop, bin_width=10, nperseg=256, noverlap=128, ax=None,
demean=True, drop_initial=0, max_freq=100, cmap='viridis',
global_vmin=None, global_vmax=None, title=None):
"""
Plot the spectrogram of the firing rate for a given population.
Args:
pop (str): The name of the population to analyze.
bin_width (int): The width of each bin in milliseconds.
nperseg (int): Length of each segment for the spectrogram.
noverlap (int): Overlap between segments.
ax (matplotlib.axes.Axes): The axes to plot on.
demean (bool): Whether to z-score the firing rate.
drop_initial (int): The initial portion of data to drop in milliseconds.
max_freq (int): Maximum frequency to display in Hz.
cmap (str): Colormap for the spectrogram.
global_vmin (float, optional): Global minimum value for color scaling.
global_vmax (float, optional): Global maximum value for color scaling.
Returns:
matplotlib.image.AxesImage: The spectrogram image
"""
if ax is None:
fig, ax = plt.subplots()
# Compute the spectrogram
frequencies, times, Sxx, fs, times_ms = self.compute_spectrogram(
pop, bin_width, nperseg, noverlap, demean, drop_initial
)
# Use global vmin and vmax if provided, otherwise use local min and max
vmin = global_vmin if global_vmin is not None else Sxx.min()
vmax = global_vmax if global_vmax is not None else Sxx.max()
# Plot the spectrogram
im = ax.pcolormesh(times_ms, frequencies, Sxx,
shading='gouraud', cmap=cmap,
vmin=vmin, vmax=vmax)
# Add a colorbar
plt.colorbar(im, ax=ax, label='Power (dB)')
# Set labels and title
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Frequency (Hz)')
if title is not None:
ax.set_title(title)
else:
ax.set_title(f'Spectrogram of {pop} Firing Rate')
ax.set_ylim(0, max_freq) # Set maximum frequency
if ax is None:
plt.show()
return im
def calculate_phase_difference(self, pop1, pop2, bin_width=2, drop_initial=100, min_power_threshold=1e-6):
"""
Calculate the phase difference between two populations' firing rates.
Args:
pop1 (str): The name of the first population.
pop2 (str): The name of the second population.
bin_width (int): The width of each bin in milliseconds.
drop_initial (int): The initial portion of data to drop in milliseconds.
Returns:
--------
phase_diff : float
Phase difference in radians (wrapped to [-π, π])
dominant_freq : float
Dominant frequency in Hz
"""
# FFT of both signals
_, signal1 = self.calculate_firing_rate(pop1, bin_width=bin_width, drop_initial=drop_initial)
_, signal2 = self.calculate_firing_rate(pop2, bin_width=bin_width, drop_initial=drop_initial)
fft1 = fft(signal1 - signal1.mean())
fft2 = fft(signal2 - signal2.mean())
fs = 1000.0 / bin_width # Sampling frequency in Hz
# Frequency array (positive frequencies only)
freqs = fftfreq(len(signal1), 1 / fs)
pos_freqs = freqs[:len(freqs) // 2]
pos_fft1 = fft1[:len(freqs) // 2]
pos_fft2 = fft2[:len(freqs) // 2]
# Calculate power spectra
power1 = np.abs(pos_fft1) ** 2
power2 = np.abs(pos_fft2) ** 2
# Find dominant frequency from signal1 (our reference)
dominant_idx = np.argmax(np.abs(pos_fft1))
dominant_freq = pos_freqs[dominant_idx]
max_power1 = power1[dominant_idx]
max_power2 = power2[dominant_idx]
# Check if signals have sufficient power at the dominant frequency
if max_power1 < min_power_threshold:
print(f"Warning: Signal 1 ({pop1}) has no significant frequency components (likely constant or noise)")
return np.nan, np.nan
if max_power2 < min_power_threshold:
print("Warning: Signal 2 ({pop2}) has no significant frequency components (likely constant or noise)")
return np.nan, np.nan
# Extract phases at dominant frequency
phase1 = np.angle(pos_fft1[dominant_idx])
phase2 = np.angle(pos_fft2[dominant_idx])
# Calculate phase difference with signal1 as reference
phase_diff = phase2 - phase1
# Wrap to [-π, π]
phase_diff = np.degrees(np.arctan2(np.sin(phase_diff), np.cos(phase_diff)))
return round(phase_diff, 1), round(dominant_freq, 1)
from scipy.optimize import minimize
class StopOptimization(Exception):
pass
def optimize_firing_rates(params_pop, params_conn_all, target_rates,
T=200, dt=0.5, use_gpu=True, drop_time=100, error_threshold=5):
populations = list(target_rates.keys())
last_results = [{'error': float('inf'), 'rates': {pop: 0 for pop in populations}}]
def objective_function(currents):
params_pop_copy = {pop: params_pop[pop].copy() for pop in params_pop}
for i, pop in enumerate(populations):
params_pop_copy[pop]['I'] = currents[i]
sim = NetworkSimulation(params_pop_copy, params_conn_all, T, dt, use_gpu)
sim.simulate(drop_time=drop_time)
error = 0
rates = {}
for i, pop in enumerate(populations):
actual_rate = sim.calculate_average_firing_rate(sim.pop_dict[pop].spike_train, start_time=drop_time)
rates[pop] = actual_rate
if pop == 'D2':
error += (actual_rate - target_rates[pop])**2 * 25
if pop == 'Proto':
error += (actual_rate - target_rates[pop])**2 / 2
else:
error += (actual_rate - target_rates[pop])**2
print(f"{pop}: current={currents[i]:.2f}, rate={actual_rate:.2f} Hz (target: {target_rates[pop]} Hz)")
print(f"Total error: {error:.2f}")
print("-----------------------------------")
last_results[0] = {'error': error, 'rates': rates}
if error < error_threshold:
raise StopOptimization # 💥 Force exit
return error
initial_currents = [params_pop[pop]['I'] for pop in populations]
try:
result = minimize(
objective_function,
initial_currents,
method='Powell',
options={
'maxiter': 20,
'disp': True,
'ftol': 5e-2
}
)
except StopOptimization:
print("Optimization stopped early due to error threshold.")
result = None
optimized_params = {pop: params_pop[pop].copy() for pop in params_pop}
if result and result.success:
for i, pop in enumerate(populations):
optimized_params[pop]['I'] = result.x[i]
else:
# fallback: use last known good currents
for i, pop in enumerate(populations):
optimized_params[pop]['I'] = initial_currents[i]
return optimized_params, result, last_results[0]['rates']
def simulate_stn_intervention(T, params_pop, params_conn_all, dt=0.5, STN_start=400, STN_end=1200, seed = 422, colorbar_range=None):
drop_time = 100
print('STN inactive')
network_sim_1 = NetworkSimulation(params_pop, params_conn_all, T, dt, seed=seed)
activate_population_1 = {'pop_name': 'STN', 'start': 6000, 'end': 7000}
network_sim_1.simulate(drop_time=drop_time, activate_population=activate_population_1)
print('STN active')
network_sim_2 = NetworkSimulation(params_pop, params_conn_all, T, dt, seed=seed)
activate_population_2 = {'pop_name': 'STN', 'start': STN_start, 'end': STN_end}
network_sim_2.simulate(sin_amp=5, drop_time=drop_time, activate_population=activate_population_2)
network_vis_1 = NetworkVisualization(network_sim_1)
network_vis_2 = NetworkVisualization(network_sim_2)
bin_width = 2
# Compute spectrograms first to determine global color scaling
spec1_freq, spec1_times, spec1_Sxx, _, _ = network_vis_1.compute_spectrogram(
pop='Proto',
bin_width=bin_width,
drop_initial=drop_time,
nperseg=64,
noverlap=32
)
spec2_freq, spec2_times, spec2_Sxx, _, _ = network_vis_2.compute_spectrogram(
pop='Proto',
bin_width=bin_width,
drop_initial=drop_time,
nperseg=64,
noverlap=32
)
# # Compute difference of spectrograms (signed)
spec_diff = spec2_Sxx - spec1_Sxx
# print('Difference min',spec_diff.min(), 'max', spec_diff.max())
# Create figure with custom axes
fig = plt.figure(figsize=(20, 10))
ax_spec_diff = fig.add_axes((0.1, 0.7, 1, 0.2))
ax_rate = fig.add_axes((0.1, 0.4, 0.8, 0.2))
ax_stn_proto_rate = fig.add_axes((0.1, 0.1, 0.8, 0.2))
vmin, vmax = (colorbar_range if colorbar_range is not None else (None, None))
im_diff = ax_spec_diff.pcolormesh(
spec2_times * 1000 + drop_time + (bin_width / 2),
spec2_freq,
spec_diff,
shading='gouraud',
cmap='RdBu_r',
vmin=vmin,
vmax=vmax
)
plt.colorbar(im_diff, ax=ax_spec_diff, label='Power Difference (dB)')
ax_spec_diff.set_xlabel('Time (ms)')
ax_spec_diff.set_ylabel('Frequency (Hz)')
ax_spec_diff.set_title('Difference Spectrogram (STN active - STN inactive)')
ax_spec_diff.set_ylim(10, 30)
network_vis_2.plot_firing_rate('Proto', ax=ax_rate, normalize=False, drop_initial=drop_time, color='green',
bin_width=bin_width, label='STN active')
network_vis_1.plot_firing_rate('Proto', ax=ax_rate, normalize=False, drop_initial=drop_time, color='red', bin_width=bin_width, label='STN inactive')
ax_rate.set_xlim(ax_spec_diff.get_xlim())
ax_rate.set_title('Proto Firing Rate')
ax_rate.legend()
ax_rate.axvspan(STN_start, STN_end, color='green', alpha=0.2)
ax_stn_proto_rate_twin = ax_stn_proto_rate.twinx()
network_vis_2.plot_firing_rate('Proto', ax=ax_stn_proto_rate, drop_initial=drop_time, color='green', bin_width=bin_width)
network_vis_2.plot_firing_rate('STN', ax=ax_stn_proto_rate_twin, drop_initial=drop_time, color='blue', bin_width=bin_width)
ax_stn_proto_rate.axvspan(STN_start, STN_end, color='green', alpha=0.2)
ax_stn_proto_rate.set_ylabel('Proto Firing Rate (Hz)')
ax_stn_proto_rate_twin.set_ylabel('STN Firing Rate (Hz)')
ax_stn_proto_rate.set_title('Proto and STN Firing Rate (STN active)')
ax_stn_proto_rate.set_xlim(ax_spec_diff.get_xlim())
plt.show()
phase_diff, ave_freq = network_vis_2.calculate_phase_difference('Proto', 'STN', bin_width=bin_width, drop_initial=STN_start)
print(f"Phase difference between Proto and STN during intervention: {phase_diff} degrees at {ave_freq} Hz")
return network_sim_1, network_sim_2
def calculate_synchronization_index(membrane_potentials):
"""
Calculate the synchronization index (chi-squared) from membrane potentials.
Parameters:
-----------
membrane_potentials : numpy.ndarray
Array of shape (n_neurons, n_timepoints) containing membrane potential traces
for each neuron over time.
Returns:
--------
chi_squared : float
The synchronization index, bounded to [0, 1] interval.
"""
# Number of neurons
N = membrane_potentials.shape[0]
# Calculate the LFP-like signal (average membrane potential across neurons)
lfp = np.mean(membrane_potentials, axis=0)
# Calculate variance of the LFP signal
var_lfp = np.var(lfp)
# Calculate variance of each neuron's membrane potential
var_individual = np.var(membrane_potentials, axis=1)
# Calculate the sum of individual variances
sum_individual_var = np.sum(var_individual)
# Calculate chi-squared
chi_squared = (N * var_lfp) / sum_individual_var
return chi_squared
def load_params_yaml(file_path='params1.yaml', dt=0.5):
"""
Load parameters from a YAML file.
Args:
file_path (str): Path to the YAML file
Returns:
dict: Loaded parameters with dt placeholders replaced
"""
# Check if file exists
if not os.path.exists(file_path):
raise FileNotFoundError(f"Parameter file not found: {file_path}")
# Load the YAML file
with open(file_path, 'r') as file:
# Load YAML content
params = yaml.safe_load(file)
# Replace ${dt} placeholders with actual dt value
params_str = yaml.dump(params)
params_str = params_str.replace('${dt}', str(dt))
params = yaml.safe_load(params_str)
params_pop = params.get('params_pop', {})
params_conn_all = params.get('params_conn_all', [])
return params_pop, params_conn_all
def get_selected_params(params_pop, params_conn_all, selected_pops):
"""
Get selected population parameters and their connections.
Args:
params_pop (dict): Dictionary of population parameters
params_conn_all (list): List of connection dictionaries
selected_pops (list): List of population names to include.
If None, returns Proto, FSI/FSN, and D2.
Returns:
tuple: (selected_pop_params, selected_connections)
"""
# Create a deep copy to avoid modifying the original
selected_pop_params = {}
for pop in selected_pops:
if 'STN' in pop:
selected_pop_params['STN'] = params_pop[pop]
elif pop in params_pop:
selected_pop_params[pop] = params_pop[pop]
else:
raise ValueError(f"Population '{pop}' not found in params_pop.")
# Only include connections between selected populations
selected_connections = []
for conn in params_conn_all:
source = conn['source']
target = conn['target']
# Include connection if both source and target are in selected populations
if source in selected_pop_params and target in selected_pop_params:
selected_connections.append(conn.copy())
return selected_pop_params, selected_connections

utils.py at commit d2babf6, no license · at the source

Overview

Authors: Ka Nap Tse¹, G. Bard Ermentrout¹, Jonathan E. Rubin¹

ORCID iDs: Jonathan E. Rubin

Department of Mathematics, University of Pittsburgh, Pittsburgh, Pennsylvania, United States of America

Institutions: University of Pittsburgh (United States)

Journal: PLoS computational biology, volume 22, issue 8, article e1013942

Dates: received 23 January 2026; accepted 2 July 2026; published online 17 August 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1371/journal.pcbi.1013942 · PMID 42607099 · PMCID PMC13541254 · OpenAlex W7125773214

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.64898/2026.01.26.701663

Status: code verified

Categories: EEG (modality), human (organism), Parkinson's (population), computational (subfield)

Methods: Spectral & time-frequency, Connectivity, Single-unit activity, calcium imaging

MeSH: Basal Ganglia*, Beta Rhythm*, Models, Neurological*, Subthalamic Nucleus*, Action Potentials, Animals, Computational Biology, Computer Simulation, Humans, Nerve Net, Neurons, Parkinson Disease (* major topic)

Journal subjects: Biology and Life Sciences, Cell Biology, Cellular Types, Animal Cells, Neurons, Neuroscience, Cellular Neuroscience, Computational Biology, Computational Neuroscience, Single Neuron Function, Computer and Information Sciences, Neural Networks, Anatomy, Brain, Basal Ganglia, Medicine and Health Sciences, Physiology, Electrophysiology, Membrane Potential, Action Potentials, Neurophysiology, Population Biology, Population Dynamics, Research and Analysis Methods, Crystallographic Techniques, Phase Determination

Topic: Neurological disorders and treatments (Neurology, Medicine), according to OpenAlex

Funding: NIH (R01NS125814, R01DA059993)

Citations: not cited yet (Europe PMC); 50 references in the paper

Abstract

Recent computational studies of Parkinson’s disease have yielded contradictory findings regarding the role of the subthalamic nucleus (STN) in pathological beta oscillations, with some models implicating STN as essential for beta generation and others suggesting that STN suppresses oscillations. This work addresses these discrepancies by systematically investigating how the specific features of the integrate-and-fire neurons used in these models influence simulated basal ganglia network dynamics. Using both rate models and spiking network simulations incorporating coupled subthalamopallidal and pallidostriatal circuits, we demonstrate that the choice between leaky integrate-and-fire (LIF) and quadratic integrate-and-fire (QIF) models to represent STN neurons fundamentally impacts the phase relationship between STN and external globus pallidus prototypical (Proto) neuron populations. QIF STN neurons establish in-phase coupling with Proto neurons, which enhances beta oscillation amplitude, while LIF STN neurons develop anti-phase relationships that suppress beta power. Through intervention experiments and parameter sweeps across physiologically relevant firing rates, we show that these phase-related effects persist robustly across network conditions, and we mathematically establish conditions under which these results are guaranteed to hold. Our findings reveal that the fundamental mathematical structure underlying spike generation, rather than other biophysical details, determines whether the subthalamopallidal loop acts as a beta amplifier or suppressor. This mechanistic insight reconciles contradictory findings in the literature, demonstrates that seemingly minor modeling choices can have profound consequences for understanding disease mechanisms and therapeutic targets, and offers predictions for determining which model framework reflects the biological reality.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 12 matches between paragraphs and lines of code.

kanaptse/BG-Oscillation-Synergy-Suppression

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: d2babf669d92ecac936c92ae194347df5021460a, 24 May 2026

Languages: Python (9), Jupyter (8)

Size: 109 files, 17 scripts

Software Heritage: not archived

Found in: “Data Availability”

Holds: README, environment (requirements.txt), 6 notebooks

Not found: license file, CITATION.cff, tests, continuous integration, documentation

Tools: NumPy (17 files), Matplotlib (13 files), SciPy (4 files), CuPy (3 files), pandas (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

18 files

01_rate_model/figures_in_paper.ipynb, Jupyter, 613 lines, 3 matches
02_tts_and_phase/debug_heatmap_holes.ipynb, Jupyter, 179 lines
02_tts_and_phase/figures_in_paper.ipynb, Jupyter, 1,346 lines
02_tts_and_phase/heatmaps.ipynb, Jupyter, 846 lines
02_tts_and_phase/interactive_phase_plane_qif_lif.py, Python, 599 lines
02_tts_and_phase/interactive_tts.py, Python, 506 lines
02_tts_and_phase/interactive_tts_v0.py, Python, 462 lines
02_tts_and_phase/sine_waveform_simulation.ipynb, Jupyter, 462 lines
02_tts_and_phase/tts_utils.py, Python, 415 lines, 1 match
03_spiking_networks/.ipynb_checkpoints/figures_in_paper-checkpoint.ipynb, Jupyter, 393 lines
03_spiking_networks/.ipynb_checkpoints/network_simulation-checkpoint.py, Python, 170 lines
03_spiking_networks/.ipynb_checkpoints/parameter_exploration-checkpoint.ipynb, Jupyter, 163 lines
03_spiking_networks/figures_in_paper.ipynb, Jupyter, 263 lines, 2 matches
03_spiking_networks/network_simulation.py, Python, 170 lines
03_spiking_networks/neuron_population.py, Python, 246 lines
03_spiking_networks/sweep_utils.py, Python, 273 lines, 1 match
03_spiking_networks/utils.py, Python, 598 lines, 5 matches
README.md, Text, 124 lines

At the source: github.com/kanaptse/BG-Oscillation-Synergy-Suppression

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
17 scripts, each with its path and the digest of its content;
12 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

The codes underlying our findings are openly available at: https://github.com/kanaptse/BG-Oscillation-Synergy-Suppression We did not generate any data for this study.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 12 MeSH terms, 1 funder, 50 references.

Cite

This paper

Tse, K. N., Ermentrout, G. B., & Rubin, J. E. (2026). Reconciling contradictory models of subthalamic nucleus contributions to basal ganglia beta oscillations. PLoS computational biology, 22(8), e1013942. https://doi.org/10.1371/journal.pcbi.1013942

BibTeX

@article{tse2026reconciling, author = {Tse, Ka Nap and Ermentrout, G. Bard and Rubin, Jonathan E.}, title = {{Reconciling contradictory models of subthalamic nucleus contributions to basal ganglia beta oscillations}}, journal = {PLoS computational biology}, year = {2026}, month = aug, volume = {22}, number = {8}, pages = {e1013942}, publisher = {PLOS}, issn = {1553-734X}, doi = {10.1371/journal.pcbi.1013942}, url = {https://doi.org/10.1371/journal.pcbi.1013942}, pmid = {42607099}, pmcid = {PMC13541254} }

RIS

TY - JOUR AU - Tse, Ka Nap AU - Ermentrout, G. Bard AU - Rubin, Jonathan E. TI - Reconciling contradictory models of subthalamic nucleus contributions to basal ganglia beta oscillations T2 - PLoS computational biology J2 - PLoS Comput Biol PY - 2026 DA - 2026/08/17 VL - 22 IS - 8 SP - e1013942 SN - 1553-734X PB - PLOS DO - 10.1371/journal.pcbi.1013942 UR - https://doi.org/10.1371/journal.pcbi.1013942 LA - en ER -

CSL-JSON

{ "id": "10.1371/journal.pcbi.1013942", "type": "article-journal", "title": "Reconciling contradictory models of subthalamic nucleus contributions to basal ganglia beta oscillations", "container-title": "PLoS computational biology", "author": [ { "family": "Tse", "given": "Ka Nap" }, { "family": "Ermentrout", "given": "G. Bard" }, { "family": "Rubin", "given": "Jonathan E." } ], "container-title-short": "PLoS Comput Biol", "volume": "22", "issue": "8", "page": "e1013942", "DOI": "10.1371/journal.pcbi.1013942", "PMID": "42607099", "PMCID": "PMC13541254", "ISSN": "1553-734X", "publisher": "PLOS", "URL": "https://doi.org/10.1371/journal.pcbi.1013942", "language": "en", "issued": { "date-parts": [ [ 2026, 8, 17 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.