Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences.

Code ↔ Paper

8 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 8 matches

[1] § Materials and Methods › Clustering of Theta‐Skipping Cells ↔ zenodo-sciad-share/scripts/analysis-2.py, lines 495–541 · score 0.76 · silhouette scores, Spectral clustering, affinity matrix, theta skipping, Cells
[2] § Materials and Methods › Simulation of Left–Right Sweeps ↔ zenodo-sciad-share/shared_src/ratinabox_utils.py, lines 210–314 · score 0.67 · theta frequency, theta phase, fraction, theta cycle, RatInABox, speed
[3] § Materials and Methods › Firing Rate Map/Head‐Direction Tuning Curve ↔ zenodo-sciad-share/shared_src/kde.py, lines 17–120 · score 0.59 · Gaussian kernel, kernel density, resolution, firing rates, bandwidth, dimensions
[4] § Results › Simultaneous MEC–Hippocampal Recordings Reveal Consistent Coupling and Left–Right Alternation at the T‐Maze Decision Point ↔ zenodo-sciad-share/scripts/analysis-2.py, lines 1155–1212 · score 0.58 · map correlation, tuning curve, Theta skipping cells, matrix, Firing rate, windows
[5] § Materials and Methods › Firing Rate Map/Head‐Direction Tuning Curve ↔ zenodo-sciad-share/shared_src/kde_fast.py, lines 76–129 · score 0.56 · von Mises kernel, Circular, smoothed, density, bins
[6] § Results › MEC‐Hippocampal Coupling Occurs Under Less Cognitively Demanding Conditions ↔ zenodo-sciad-share/scripts/analysis-1.py, lines 556–629 · score 0.55 · egocentric internal direction, MEC internal direction, MEC cell, theta phase, Head direction, frames
[7] § Results › MEC‐Hippocampal Coupling Persists Even at Forced Turn Corners ↔ zenodo-sciad-share/scripts/analysis-1.py, lines 556–629 · score 0.55 · egocentric internal direction, MEC internal direction, MEC cell, theta phase, Head direction, frames
[8] § Results › Simultaneous MEC–Hippocampal Recordings Reveal Consistent Coupling and Left–Right Alternation at the T‐Maze Decision Point ↔ zenodo-sciad-share/shared_src/ratinabox_utils.py, lines 210–314 · score 0.50 · theta sequence, pre, theta phase, agent, sweeps, position

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 · 1,664 lines · 54 KB · CC-BY-4.0 · 2 matches

import sys
from functools import partial
from pathlib import Path
import pandas as pd
from scipy.io import loadmat
import scipy.signal
import numpy as np
from sklearn.cluster import SpectralClustering
from sklearn.metrics import silhouette_score
from scipy.stats import circmean
from scipy.ndimage import gaussian_filter1d
from scipy import stats
import matplotlib
import matplotlib.path as mpath
from matplotlib.patches import Rectangle
from matplotlib import colors
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from kalmax.kde import poisson_log_likelihood
import jax.numpy as jnp
# root_dir = Path("/Users/nakanomasahiro/PycharmProjects/research/hpc-analysis")
root_dir = Path("/Users/nakanomasahiro/PycharmProjects/research/hpc-analysis")
source_dir = Path(
"/Users/nakanomasahiro/Desktop/research/research-projects/left-right-sweeps-paper/zenodo-sciad"
)
# Define the dataset
dataset = "left-right-sweeps"
# Setup
# fmt: off
# sys.path.append(str(root_dir))
sys.path.append(str(source_dir))
from shared_src.utils import get_S, apply_mask_lambda, chunker, find_nearest, chunker_v2
from shared_src.kalmax_utils import kalmax_place_field_pos
from shared_src.kde_fast import get_circular_tuning_curve_all, kernel_density_estimation, get_circular_tuning_curve_all_v2
from shared_src.decoding_utils import make_bool_visited_map, get_error_stats
from shared_src.spike_autocorrelation import spike_correlation
from shared_src.linearlize import (
linearize_for_mmaze_v2,
)
from shared_src.fig_utils import (
erase_topright_axes,
erase_toprightleft_axes,
)
# sys.path.append(str(root_dir / "projects" / dataset))
from src.dataset import (
get_xarray_dataset_v1,
get_range_mask_from_ds,
)
from src.figure.ds_fig import (
plot_spike_correlation_from_ds_v2
)
# fmt: on
data_dir = root_dir / "data" / dataset
root_fig_dir = root_dir / "figs" / dataset
# fig_dir = root_fig_dir / "papers-conferences" / "paper" / "ScienceAdvances"
root_fig_dir = source_dir / "figs"
fig_dir = root_fig_dir
fig_dir.mkdir(parents=True, exist_ok=True)
#### Processing data
xs = [
0,
0,
-0.5,
-0.5,
0.475,
0.475,
]
ys = [-0.6, 0.55, 0.55, -0.6, 0.55, -0.6]
task_name = "mmaze"
session_name = "29502_1"
file_path = data_dir / "navigation" / task_name / f"{session_name}.mat"
area_name = "hc"
# edit here add some other functions
ds, df = get_xarray_dataset_v1(
file_path=file_path,
area_name=area_name,
fast_kde=True,
run_decoding=False,
include_unvisited_location=False,
# 2026-01-19 use proper binning mode
bin_mode="center",
)
ds_mec, df_mec = get_xarray_dataset_v1(
file_path=file_path,
area_name="mec",
fast_kde=True,
run_decoding=False,
include_unvisited_location=False,
# 2026-01-19 use proper binning mode
bin_mode="center",
)
X_1d, direction, X_proj, path_ids, borders = linearize_for_mmaze_v2(
X=ds.X.values, xs=xs, ys=ys
)
ds["X_1d"] = (("frame",), X_1d)
ds["direction"] = (("frame",), direction)
points = np.stack([xs, ys], axis=1)
points_1d, _, _, _, _ = linearize_for_mmaze_v2(X=points, xs=xs, ys=ys)
borders["top_left"] = points_1d[2]
borders["top_right"] = points_1d[4]
ds.attrs["borders"] = borders
ds["path_ids"] = (("frame",), path_ids)
area_name = "mec"
data_ = loadmat(file_path, squeeze_me=True, struct_as_record=False)
data = data_["Dsession"]
lmt_comp = getattr(data.lmt, area_name)
lmt_pos = lmt_comp.pos.XA # shape: [n_time_bins × 2]
ds["lmt"] = (["frame", "space"], lmt_pos)
lmt_id = lmt_comp.id.XA
ds["lmt_id"] = (["frame"], lmt_id)
ds_mec["lmt_id"] = (["frame"], lmt_id)
phi, tuning_curves = get_circular_tuning_curve_all_v2(
ds.head_direction.values, ds_mec.S.values
)
ds_mec["bin_angles"] = (["n_hd_bin"], phi)
ds_mec["hd_tuning_curve"] = (["cluster", "n_hd_bin"], tuning_curves)
ds_mec["lmt_id_tuning_curve"] = (
["cluster", "n_hd_bin"],
get_circular_tuning_curve_all_v2(ds.lmt_id.values, ds_mec.S.values)[1],
)
ds["cosine_head_direction"] = np.cos(ds.head_direction)
assert np.all(lmt_comp.hd.XA == ds.head_direction.values)
ds["ego_id"] = (ds.lmt_id - ds.head_direction + np.pi) % (2 * np.pi) - np.pi
# color_lis = ['red', 'blue']
color_lis = ["darkblue", "crimson"]
# colors = {1: "C0", -1: "C1"}
# 1D population decoding
X_1d_expanded = np.stack([X_1d, np.zeros_like(X_1d)], axis=1)
border_central_left = (borders["central_arm_end"] + borders["left_arm_start"]) / 2
border_left_right = (borders["left_arm_end"] + borders["right_arm_start"]) / 2
barriers = np.array(
[
[border_central_left, border_central_left, -0.02, 0.02],
[border_left_right, border_left_right, -0.02, 0.02],
]
)
firing_rate, firing_maps, spike_maps, position_d, bins, shape, env = (
kernel_density_estimation(X_1d, S=ds.S.values, barriers=barriers)
)
arm_labels = np.zeros(100).astype(int)
max_ = borders["right_arm_end"]
i1 = int(border_central_left / max_ * 100) + 2
i2 = int(border_left_right / max_ * 100) + 2
arm_labels[i1:i2] = 1
arm_labels[i2:] = 2
log_likelihoods = poisson_log_likelihood(
spikes=ds.S.values, mean_rate=firing_rate
) # shape (T, N_bins)
likelihoods = np.exp(log_likelihoods) # shape (T, N_bins)
normed_likelihoods = likelihoods / (likelihoods.sum(axis=1, keepdims=True) + 1e-12)
ML_modes = np.argmax(likelihoods, axis=1)
ml_pos = np.array(bins[ML_modes])
S_count = ds.S.values.sum(axis=1)
use_frames = np.where(S_count > 0)[0]
bin_labels = np.zeros_like(bins)
bin_labels = np.where(
(bins > border_central_left) & (bins <= border_left_right), 1, bin_labels
)
bin_labels = np.where((bins >= border_left_right), 2, bin_labels)
cluster_colors = ["lightgray", "darkblue", "crimson"]
# cluster_colors = ['lightgray', 'darkblue', 'crimson']
stripe_y = np.zeros((100, 1, 3))
for i, lbl in enumerate(arm_labels):
stripe_y[i, 0] = colors.to_rgb(cluster_colors[lbl])
chunk_dicts = {}
chunks, sizes = chunker(np.where((direction == 1) & (path_ids == 0))[0], 50)
chunks = [chunk for chunk, size in zip(chunks, sizes) if size > 200]
sizes = [size for size in sizes if size > 200]
chunk_dicts["outbound_central_arm"] = chunks
arm_ids = {"left": 1, "right": 2}
for arm_name, arm_id in arm_ids.items():
chunks, sizes = chunker(np.where((direction == -1) & (path_ids == arm_id))[0], 50)
chunks = [chunk for chunk, size in zip(chunks, sizes) if size > 200]
sizes = [size for size in sizes if size > 200]
chunk_dicts[f"inbound_from_{arm_name}_arm"] = chunks
arm_ids = {"left": 1, "right": 2}
for arm_name, arm_id in arm_ids.items():
chunks, sizes = chunker(np.where((direction == 1) & (path_ids == arm_id))[0], 50)
chunks = [chunk for chunk, size in zip(chunks, sizes) if size > 200]
sizes = [size for size in sizes if size > 200]
chunk_dicts[f"outbound_to_{arm_name}_arm"] = chunks
ds_dicts = {"hc": ds, "mec": ds_mec}
print(f"HPC cluster size: {ds.cluster.size}")
print(f"MEC cluster size: {ds_mec.cluster.size}")
range_dicts = {
"outbound_central_arm": {
"X_1d": (0.2, 1.1),
"direction": (0.9, 1.1),
"head_direction": (np.pi / 3, np.pi * 2 / 3),
},
"inbound_from_left_arm": {
"X_1d": (1.2, 1.7),
"direction": (-2, 0),
"head_direction": (-np.pi / 3, np.pi / 3),
},
"outbound_topleft_corner": {
"X_1d": (1.3, 1.7),
"direction": (0, 2),
"cosine_head_direction": (-2, -0.5),
},
}
range_dicts_supp = {
"inbound_from_right_arm": {
"X_1d": (3, 3.5),
"direction": (-2, 0),
"cosine_head_direction": (-2, -0.5),
},
"inbound_topright_corner": {
"X_1d": (3.6, 4.2),
"direction": (-2, 0),
"head_direction": (np.pi / 3, np.pi * 2 / 3),
},
"outbound_topright_corner": {
"X_1d": (3, 3.5),
"direction": (0, 2),
"head_direction": (-np.pi / 3, np.pi / 3),
},
"inbound_topleft_corner": {
"X_1d": (1.8, 2.4),
"direction": (-2, 0),
"head_direction": (np.pi / 3, np.pi * 2 / 3),
},
}
def get_theta_skipping_index_v2(
spike_times,
window_size=40,
spike_times_2=None,
method="numba_v2",
):
bins, counts = spike_correlation(
spike_times_1=spike_times,
spike_times_2=spike_times_2,
bin_size=1,
window=window_size,
method=method,
)
p1 = counts[window_size + 6 : window_size + 18].mean()
p2 = counts[window_size + 18 : window_size + 30].mean()
p_even = p2
p_odd = p1
skipping_index = (p_even - p_odd) / max(p_even, p_odd, 1e-12)
return skipping_index, p_even, p_odd
def show_correlogram(
spike_times,
window_size=40,
color="lightgray",
ax=None,
alpha=0.7,
label=None,
spike_times_2=None,
):
"""f'shuffled within mask, {skipping_index:.2f}'"""
bins, counts = spike_correlation(
spike_times_1=spike_times,
spike_times_2=spike_times_2,
bin_size=1,
window=window_size,
method="numba_v2",
)
if not ax:
fig, ax = plt.subplots(figsize=(1, 1))
ax.bar(
bins[window_size:],
counts[window_size:],
width=1,
alpha=alpha,
color=color,
label=label,
)
ax.set_xticks([0, 10, 40])
ax.set_xticklabels([0, 100, 400], rotation=45)
erase_topright_axes(ax)
return ax
def make_theta_cycle_chunks(theta_phase_in_mask):
"""
theta_phase_in_mask = ds.theta_phase.sel(frame=mask).values
"""
chunks, sizes = chunker_v2(np.where(np.diff(theta_phase_in_mask) > 0)[0], 1)
theta_cycle_chunks = []
last_idx = -1
for chunk in chunks:
if chunk[0] - last_idx > 1:
chunk.insert(0, chunk[0] - 1)
chunk.append(chunk[-1] + 1)
theta_cycle_chunks.append(chunk)
last_idx = chunk[-1]
sizes = [len(c) for c in theta_cycle_chunks]
# assert np.concatenate(theta_cycle_chunks).shape == mask.sum()
return theta_cycle_chunks, sizes
def run_shuffling_test_for_latent_position_cycling(
ds,
ML_modes,
bin_labels,
range_dict,
target_arm_1=1,
target_arm_2=2,
window_size=40,
n_samples=1000,
show_plots=False,
):
range_dict["S_count"] = (0, 100)
mask = get_range_mask_from_ds(ds, range_dict)
frames = np.where(mask)[0]
theta_cycle_chunks, _ = make_theta_cycle_chunks(
ds.theta_phase.sel(frame=mask).values
)
# Real data latent position
latent_pos_masked = bin_labels[ML_modes[mask]]
spike_times = frames[np.where(latent_pos_masked == target_arm_1)[0]]
spike_times_2 = frames[np.where(latent_pos_masked == target_arm_2)[0]]
real_skipping_index, _, _ = get_theta_skipping_index_v2(
spike_times, window_size=window_size, spike_times_2=spike_times_2
)
# Getting shuffling idx list
shuffled_idx_list = []
for i in range(n_samples):
shuffled_idx = np.concatenate(
[
theta_cycle_chunks[i]
for i in np.random.permutation(len(theta_cycle_chunks))
]
)
shuffled_idx_list.append(shuffled_idx)
latent_pos_masked = bin_labels[ML_modes[mask]]
n_samples = n_samples
# Perform shuffling
lis = []
for i in range(n_samples):
if shuffled_idx_list is not None:
shuffled_idx = shuffled_idx_list[i]
else:
shuffled_idx = np.concatenate(
[
theta_cycle_chunks[i]
for i in np.random.permutation(len(theta_cycle_chunks))
]
)
latent_pos_shuffled = bin_labels[ML_modes[shuffled_idx]]
spike_times_shuffled = frames[np.where(latent_pos_shuffled == target_arm_1)[0]]
spike_times_shuffled_2 = frames[
np.where(latent_pos_shuffled == target_arm_2)[0]
]
shuffled_skipping_index, _, _ = get_theta_skipping_index_v2(
spike_times_shuffled,
window_size=window_size,
spike_times_2=spike_times_shuffled_2,
)
lis.append(shuffled_skipping_index)
is_cycling = real_skipping_index < np.percentile(lis, 5)
print("Real skipping index:", real_skipping_index)
print("5th percentile of shuffled skipping index:", np.percentile(lis, 5))
print("Is cycling between latent positions:", is_cycling)
print("")
if show_plots:
titles = ["center arm", "left arm", "right arm"]
fig, axes = plt.subplots(1, 4, figsize=(5, 1))
for i, j in enumerate([0, 1, 2]):
ax = axes[i]
show_correlogram(
frames[np.where(latent_pos_masked == j)[0]],
ax=ax,
color=cluster_colors[j],
)
ax.set_title(titles[i])
ax.set_yticks([])
ax = axes[3]
show_correlogram(
frames[np.where(latent_pos_masked == 1)[0]],
spike_times_2=frames[np.where(latent_pos_masked == 2)[0]],
color="#262626",
ax=ax,
)
ax.set_title("left vs right")
ax.set_yticks([])
plt.figure(figsize=(2, 1.5))
plt.hist(lis, bins=30, color="gray", alpha=0.7)
plt.axvline(real_skipping_index, color="red")
erase_topright_axes()
################################
### Get theta skipping cells ###
################################
def compare_real_and_shuffled_correlogram(
ds,
cluster_idx,
mask,
theta_cycle_chunks,
ax=None,
c1="lightgray",
c2="#262626",
alpha1=0.7,
alpha2=0.3,
):
frames = np.where(mask)[0]
shuffled_idx = np.concatenate(
[theta_cycle_chunks[i] for i in np.random.permutation(len(theta_cycle_chunks))]
)
S = ds.S.sel(cluster=cluster_idx).sel(frame=mask).values
S_shuffled = S[shuffled_idx]
ax = show_correlogram(
frames[np.where(S)[0]], alpha=alpha1, color=c1, label="real", ax=ax
)
show_correlogram(
frames[np.where(S_shuffled)[0]],
ax=ax,
color=c2,
alpha=alpha2,
label="shuffled",
)
ax.legend(loc="upper right", bbox_to_anchor=(2.4, 1))
return ax
################################################
### Apply clustering to theta skipping cells ###
################################################
def get_affinity_matrix(ds, mask, cluster_lis, window_size=40):
affinity_mat = np.zeros((len(cluster_lis), len(cluster_lis)))
for i, clu1 in enumerate(cluster_lis):
spike_times_1 = np.where(ds.S.sel(cluster=clu1).sel(frame=mask))[0]
for j, clu2 in enumerate(cluster_lis):
spike_times_2 = np.where(ds.S.sel(cluster=clu2).sel(frame=mask))[0]
affinity_mat[i, j] = get_theta_skipping_index_v2(
spike_times=spike_times_1,
spike_times_2=spike_times_2,
window_size=window_size,
)[0]
return affinity_mat
def cluster_affinity_matrix(
affinity_mat, cluster_lis, random_state=42, cluster_range=range(2, 9)
):
assert affinity_mat.shape[0] == len(cluster_lis)
min_cluster = cluster_range.start
max_cluster = cluster_range.stop - 1
cluster_range = range(min_cluster, min(max_cluster, affinity_mat.shape[0] - 1) + 1)
M = affinity_mat.copy()
np.fill_diagonal(M, 0)
M = np.nan_to_num((M + M.T) / 2, nan=0.0)
M_min, M_max = M.min(), M.max()
if M_min < 0:
M = (M - M_min) / (M_max - M_min + 1e-12)
np.fill_diagonal(M, 0) # spectral doesn’t need self-similarity
D = 1 - M # simple distance from affinity in [0,1]
np.fill_diagonal(D, 0)
best = {"k": None, "score": -np.inf, "labels": None}
for k in cluster_range:
sc = SpectralClustering(
n_clusters=k,
affinity="precomputed",
assign_labels="kmeans",
random_state=random_state,
).fit(M)
score = silhouette_score(D, sc.labels_, metric="precomputed")
if score > best["score"]:
best = {"k": k, "score": score, "labels": sc.labels_}
labels = best[
"labels"
] # cluster label per cell (same order as real_theta_skipping_cells)
labels_dict = {}
for i in range(0, best["k"]):
labels_dict[i] = cluster_lis[labels == i]
order = np.argsort(labels)
M_ord = M[order][:, order]
cells_ord = [cluster_lis[i] for i in order]
return {
"best_k": best["k"],
"best_score": best["score"],
"labels": labels,
"labels_dict": labels_dict,
"M_ord": M_ord,
"cells_ord": cells_ord,
}
# affinity_mat = get_affinity_matrix(
# ds, mask, theta_skipping_cell_result["real_theta_skipping_cells"], window_size=40
# )
# clustering_result = cluster_affinity_matrix(
# affinity_mat, cluster_lis=theta_skipping_cell_result["real_theta_skipping_cells"]
# )
########################################################
### Overlap of place cell rate maps within the group ###
########################################################
def within_across_values(M, labels):
"""Extract within- and across-group metric values from symmetric matrix M."""
n = len(labels)
ut = np.triu(np.ones((n, n), bool), k=1) # upper triangle only (no diagonal)
same = labels[:, None] == labels[None, :]
diff = ~same
within_vals = M[same & ut]
across_vals = M[diff & ut]
return within_vals, across_vals
def permutation_test_for_within_across_group_representation_similarity(
M, labels, n_perm=5000, random_seed=None
):
if random_seed is None:
rng = np.random.default_rng()
else:
rng = np.random.default_rng(random_seed)
real_within, real_across = within_across_values(M, labels)
real_diff = np.nanmean(real_within) - np.nanmean(real_across)
n = len(labels)
null_diffs = np.empty(n_perm)
for i in range(n_perm):
perm = rng.permutation(labels)
w, a = within_across_values(M, perm)
null_diffs[i] = np.nanmean(w) - np.nanmean(a)
p = (1 + np.sum(null_diffs >= real_diff)) / (n_perm + 1)
z = (real_diff - np.mean(null_diffs)) / (np.std(null_diffs, ddof=1) + 1e-12)
return {"real_diff": real_diff, "p": p, "z": z, "null_diffs": null_diffs}
def show_within_across_distribution(M, labels):
real_within, real_across = within_across_values(M, labels)
fig, ax = plt.subplots(figsize=(1.5, 1.5))
ax.hist(real_within, bins=np.linspace(-1, 1, 25), alpha=0.5, label="within")
ax.hist(real_across, bins=np.linspace(-1, 1, 25), alpha=0.5, label="between")
ax.legend()
erase_topright_axes(ax)
ax.set_title(f"{np.nanmean(real_within) - np.nanmean(real_across):.3f}")
def show_matrix_with_labels(M, labels, vmin=-1, vmax=1):
order = np.argsort(labels)
M_ord = M[order][:, order]
fig, ax = plt.subplots(figsize=(2, 2))
im = ax.imshow(M_ord, cmap="viridis", vmin=vmin, vmax=vmax)
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
erase_topright_axes(ax)
def show_correlation_matrix(
M,
skipping_cells,
labels=None,
color_for=["blue", "red"],
ax=None,
plot_labels=True,
group_order=None,
vmin=None,
vmax=None,
):
if labels is None:
labels = np.arange(M.shape[0])
if color_for is None:
cmap = plt.get_cmap("hsv")
color_for = {i: cmap(i % cmap.N) for i in range(len(np.unique(labels)))}
if group_order is None:
order = np.argsort(labels)
else:
order = np.concatenate([np.where(labels == g)[0] for g in group_order])
M_ord = M[order][:, order]
cells_ord = [skipping_cells[i] for i in order]
labels_ord = np.asarray(labels)[order]
if ax is None:
fig, ax = plt.subplots(figsize=(1.2, 1.2))
im = ax.imshow(M_ord, aspect="equal", interpolation="nearest", vmin=vmin, vmax=vmax)
if plot_labels:
ax.set_xticks(np.arange(len(cells_ord)))
ax.set_xticklabels(cells_ord, rotation=90)
ax.set_yticks(np.arange(len(cells_ord)))
ax.set_yticklabels(cells_ord)
else:
ax.set_xticks([])
ax.set_yticks([])
# --- find contiguous cluster blocks along the diagonal ---
edges = np.flatnonzero(
np.r_[True, labels_ord[1:] != labels_ord[:-1], True]
) # block boundaries
erase_topright_axes(ax)
# ax.axis("off")
# --- draw rectangles ---
for i, (a, b) in enumerate(zip(edges[:-1], edges[1:])): # block = [a, b)
n = b - a
cid = labels_ord[a]
rect = Rectangle(
(a - 0.5, a - 0.5), n, n, fill=False, lw=2, ec=color_for[i], zorder=3
)
ax.add_patch(rect)
ax.set_xlim(-len(cells_ord) * 0.1, len(cells_ord) * 1.1 - 1)
ax.set_ylim(len(cells_ord) * 1.1 - 1, -len(cells_ord) * 0.1)
ax.spines.left.set(visible=False)
ax.spines.bottom.set(visible=False)
return ax
def get_percentile_shuffle_with_theta_preserved_v2(
frames,
S_in_mask,
theta_cycle_chunks,
percentile=95,
n_samples=10000,
window_size=40,
return_lis=False,
shuffled_idx_list=None,
):
lis = []
p_even_lis = []
p_odd_lis = []
for i in range(n_samples):
if shuffled_idx_list is not None:
shuffled_idx = shuffled_idx_list[i]
else:
shuffled_idx = np.concatenate(
[
theta_cycle_chunks[i]
for i in np.random.permutation(len(theta_cycle_chunks))
]
)
S_shuffled = S_in_mask[shuffled_idx]
spike_times_shuffled = frames[np.where(S_shuffled)[0]]
skipping_index, p_even, p_odd = get_theta_skipping_index_v2(
spike_times_shuffled, window_size=window_size
)
lis.append(skipping_index)
p_even_lis.append(p_even)
p_odd_lis.append(p_odd)
if return_lis:
return (
np.percentile(lis, percentile),
lis,
np.percentile(p_even_lis, percentile),
np.percentile(p_odd_lis, percentile),
)
else:
return (
np.percentile(lis, percentile),
np.percentile(p_even_lis, percentile),
np.percentile(p_odd_lis, percentile),
)
def get_theta_skipping_cells_with_and_wo_hard_threshold(
ds,
mask,
n_samples=1000,
firing_rate_threshold={"hc": 1, "mec": 5},
percentile=95,
hard_threshold_theta_skipping_index=0.5,
):
"""
Skip cells that have less than 1 hz for HPC and 5 hz for MEC firing rate within the mask
"""
theta_cycle_chunks, sizes = make_theta_cycle_chunks(
ds.theta_phase.sel(frame=mask).values
)
shuffled_idx_list = []
for i in range(n_samples):
shuffled_idx = np.concatenate(
[
theta_cycle_chunks[i]
for i in np.random.permutation(len(theta_cycle_chunks))
]
)
shuffled_idx_list.append(shuffled_idx)
shuffle_test_result = {}
loose_theta_skipping_cells = []
strict_theta_skipping_cells = []
fake_theta_skipping_cells = []
# bool_is_theta_skipping = []
frames = ds.frame.sel(frame=mask).values
for cluster_idx in ds.cluster.values:
S_in_mask = ds.S.sel(cluster=cluster_idx).sel(frame=mask).values
if S_in_mask.sum() < ((mask.sum() / 100) * firing_rate_threshold[ds.area_name]):
continue
spike_times_1 = frames[np.where(S_in_mask)[0]]
real_skipping_index, real_p_even, real_p_odd = get_theta_skipping_index_v2(
spike_times_1
)
(
shuffled_percentile_skipping_index,
shuffled_percentile_p_even,
shuffled_percentile_p_odd,
) = get_percentile_shuffle_with_theta_preserved_v2(
frames,
S_in_mask,
theta_cycle_chunks,
n_samples=n_samples,
return_lis=False,
percentile=percentile,
shuffled_idx_list=shuffled_idx_list,
)
# is_skipping = (real_skipping_index > shuffled_percentile_skipping_index) & (real_p_even > shuffled_percentile_p_even) & (real_p_odd < shuffled_percentile_p_odd)
is_skipping = real_skipping_index > shuffled_percentile_skipping_index
is_hard_threshold_skipping = (
real_skipping_index > hard_threshold_theta_skipping_index
)
shuffle_test_result[cluster_idx] = {
"real_skipping_index": real_skipping_index,
"shuffled_percentile_skipping_index": shuffled_percentile_skipping_index,
"real_p_even": real_p_even,
"shuffled_percentile_p_even": shuffled_percentile_p_even,
"real_p_odd": real_p_odd,
"shuffled_percentile_p_odd": shuffled_percentile_p_odd,
"is_skipping": is_skipping,
"is_hard_threshold_skipping": is_hard_threshold_skipping,
}
if is_skipping:
loose_theta_skipping_cells.append(cluster_idx)
if is_hard_threshold_skipping:
strict_theta_skipping_cells.append(cluster_idx)
else:
fake_theta_skipping_cells.append(cluster_idx)
# bool_is_theta_skipping.append(is_skipping)
result = {
"shuffle_test_result": shuffle_test_result,
"loose_theta_skipping_cells": np.array(loose_theta_skipping_cells),
"strict_theta_skipping_cells": np.array(strict_theta_skipping_cells),
"fake_theta_skipping_cells": np.array(fake_theta_skipping_cells),
# "bool_is_theta_skipping": np.array(bool_is_theta_skipping),
}
return result
def show_all_theta_skipping_cells(ds, theta_skipping_cell_result, mask, n_cols=10):
n_cells = len(theta_skipping_cell_result["real_theta_skipping_cells"])
n_rows = int(np.ceil(n_cells / n_cols))
fig = plt.figure(figsize=(n_cols * 1, n_rows * 1))
for i, cluster_idx in enumerate(
theta_skipping_cell_result["real_theta_skipping_cells"]
):
ax = fig.add_subplot(n_rows, n_cols, i + 1)
plot_spike_correlation_from_ds_v2(
ds, cluster_idx, mask=mask, ax=ax, window_size=40
)
ax.axvspan(0, 6, color="lightgreen", alpha=0.5)
ax.axvspan(6, 18, color="pink", alpha=0.5)
ax.axvspan(18, 30, color="lightgreen", alpha=0.5)
ax.set_xlim(0, 40)
ax.set_title(f"clu {cluster_idx}", fontsize=6)
erase_topright_axes(ax)
ax.set_xticks([])
ax.set_xlabel("")
ax.set_ylabel("")
# ax.spines['left'].set_fontsize(4)
plt.tight_layout()
return fig
def show_all_theta_skipping_cells_v2(
ds, theta_skipping_cells, mask, n_cols=10, yticks=False
):
n_cells = len(theta_skipping_cells)
n_rows = int(np.ceil(n_cells / n_cols))
fig = plt.figure(figsize=(n_cols * 1, n_rows * 1))
for i, cluster_idx in enumerate(theta_skipping_cells):
ax = fig.add_subplot(n_rows, n_cols, i + 1)
plot_spike_correlation_from_ds_v2(
ds, cluster_idx, mask=mask, ax=ax, window_size=40
)
ax.axvspan(0, 6, color="lightgreen", alpha=0.5)
ax.axvspan(6, 18, color="pink", alpha=0.5)
ax.axvspan(18, 30, color="lightgreen", alpha=0.5)
ax.set_xlim(0, 40)
ax.set_title(f"clu {cluster_idx}", fontsize=6)
erase_topright_axes(ax)
ax.set_xticks([])
ax.set_xlabel("")
ax.set_ylabel("")
if not yticks:
ax.set_yticks([])
# ax.spines['left'].set_fontsize(4)
plt.tight_layout()
return fig
def get_hpc_group_direction(
ds,
mask,
labels_dict,
):
assert list(labels_dict.keys()) == [0, 1]
res = []
for i in [0, 1]:
cum_place_field_center = ds.bin_positions.values[
np.nanargmax(
ds.firing_rate_map_mask.sel(cluster=labels_dict[i])
.mean(dim="cluster")
.values.flatten()
)
]
mean_position = ds.X.sel(frame=mask).mean(dim="frame").values
mean_hd = circmean(ds.head_direction.sel(frame=mask).values)
offset_direction = np.arctan2(
cum_place_field_center[1] - mean_position[1],
cum_place_field_center[0] - mean_position[0],
)
is_relative_left = np.sin(offset_direction - mean_hd) > 0
res.append(is_relative_left)
assert (res == [0, 1]) or (
res == [1, 0]
), "The two groups are not opposite in direction"
left_then_right = np.argsort(res)[::-1]
return left_then_right
def get_mec_group_direction(
ds_mec,
mask,
labels_dict,
):
assert list(labels_dict.keys()) == [0, 1]
res = []
for i in [0, 1]:
# tuning_curve_peak = ds_mec.bin_angles.values[np.argmax(ds_mec.hd_tuning_curve.sel(cluster=labels_dict[i]).values.mean(axis=0))]
y = ds_mec.hd_tuning_curve.sel(cluster=labels_dict[i]).values.mean(axis=0)
# plt.plot(ds_mec.bin_angles, y)
vector_sum = np.sum(y * np.exp(1j * ds_mec.bin_angles.values))
tuning_curve_peak = np.angle(vector_sum) - np.pi
mean_hd = circmean(ds.head_direction.sel(frame=mask).values)
is_relative_left = int(np.sin(tuning_curve_peak - mean_hd) > 0)
res.append(is_relative_left)
assert (res == [0, 1]) or (
res == [1, 0]
), "The two groups are not opposite in direction"
left_then_right = np.argsort(res)[::-1]
return left_then_right
def get_group_direction(ds, mask, labels_dict, area_name):
if area_name == "hc":
return get_hpc_group_direction(ds, mask, labels_dict)
elif area_name == "mec":
return get_mec_group_direction(ds, mask, labels_dict)
else:
raise ValueError("area_name must be either 'hc' or 'mec'")
def show_cumulative_firing_rate_map(ds, labels_dict, group_order=[0, 1]):
fig = plt.figure(figsize=(2, 1))
for i, j in enumerate(group_order):
clusters = labels_dict[j]
n_clusters = len(clusters)
ax = fig.add_subplot(1, 2, i + 1)
# cumulative_rate_map = np.zeros_like(
# ds.firing_rate_map_mask.isel(cluster=0)
# )
# for j, cluster_idx in enumerate(clusters):
# cumulative_rate_map += ds.firing_rate_map_mask.sel(
# cluster=cluster_idx
# )
cumulative_rate_map = ds.firing_rate_map_mask.sel(cluster=clusters).mean(
dim="cluster"
)
ax.imshow(cumulative_rate_map.T[::-1], interpolation="none")
ax.axis("off")
def show_overlaid_hd_tuning_curve(
ds_mec,
mask,
labels_dict,
plot_var="head_direction",
group_order=[0, 1],
plot_hd_and_id=False,
density=False,
):
angles = ds_mec.head_direction.sel(frame=mask)
mean_angle = circmean(angles, high=np.pi, low=-np.pi)
fig = plt.figure(figsize=(1.4, 0.6))
for i, j in enumerate(group_order):
ax = fig.add_subplot(1, 2, i + 1, polar=True)
ax.grid(True, linewidth=0.25) # default is around 1.0
ax.spines["polar"].set_linewidth(0.25) # default is around 1.0
clusters = labels_dict[j]
y = ds_mec.hd_tuning_curve.sel(cluster=labels_dict[j]).values
if density:
y = y / np.sum(y, axis=1, keepdims=True)
# for cluster_idx in clusters:
# plot_head_direction_selectivity_from_ds(
# ds_mec, cluster_idx, plot_hd_and_id=plot_hd_and_id,
# plot_var=plot_var,
# ax=ax, linewidth=0.5
# )
# ax.set_title("")
ax.plot(
np.append(ds_mec.bin_angles, ds_mec.bin_angles[0]),
# np.append(ds_mec.hd_tuning_curve.sel(cluster=labels_dict[i]).values, ds_mec.hd_tuning_curve.sel(cluster=labels_dict[i]).values[0]).T,
np.concatenate([y, y[:, 0][:, None]], axis=1).T,
color="#262626",
linewidth=0.5,
)
ax.set_theta_zero_location("W") # optional: 0° at top
ax.set_theta_direction(1) # optional: increase clockwise
ax.set_yticks([])
# ax.set_xticks([])
ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2])
ax.set_xticklabels([])
ax.arrow(
0,
0,
np.cos(mean_angle) * 0.5 * ax.get_rmax(),
np.sin(mean_angle) * 0.5 * ax.get_rmax(),
transform=ax.transData._b,
facecolor="#2BB04F",
edgecolor="#A7F2BC",
linewidth=0.4,
width=0.1 * ax.get_rmax(),
head_length=0.25 * ax.get_rmax(),
zorder=100,
)
def show_group_response(
ds, area_name, mask, labels_dict, group_order=[0, 1], density=True
):
if area_name == "hc":
show_cumulative_firing_rate_map(ds, labels_dict, group_order=group_order)
elif area_name == "mec":
show_overlaid_hd_tuning_curve(
ds, mask, labels_dict, group_order=group_order, density=density
)
else:
raise ValueError("area_name must be either 'hc' or 'mec'")
def get_mec_hc_group_average_lag(
ds_dicts, stats_all, range_dicts, condition, threshold, window_size=15
):
mask = get_range_mask_from_ds(ds_dicts["hc"], range_dicts[condition])
hc_labels_dict = stats_all[condition]["hc"][threshold]["clustering_result"][
"labels_dict"
]
mec_labels_dict = stats_all[condition]["mec"][threshold]["clustering_result"][
"labels_dict"
]
hc_group_order = get_hpc_group_direction(ds_dicts["hc"], mask, hc_labels_dict)
mec_group_order = get_mec_group_direction(ds_dicts["mec"], mask, mec_labels_dict)
group_lags = {}
for i, (hc_group, mec_group) in enumerate(zip(hc_group_order, mec_group_order)):
peak_frame_lis = []
for hc_clu in hc_labels_dict[hc_group]:
for mec_clu in mec_labels_dict[mec_group]:
spike_times_1 = np.where(
ds_dicts["hc"].S.sel(cluster=hc_clu).sel(frame=mask)
)[0]
spike_times_2 = np.where(
ds_dicts["mec"].S.sel(cluster=mec_clu).sel(frame=mask)
)[0]
counts = spike_correlation(
spike_times_1=spike_times_1,
spike_times_2=spike_times_2,
bin_size=1,
window=window_size,
method="numba_v2",
)[1]
peak_frame_lis.append(int(np.argmax(counts)))
lags = (np.array(peak_frame_lis) - window_size) * 10
group_lags[f"group_{i}_peak_frame_lis"] = peak_frame_lis
group_lags[f"group_{i}_mean_lag_ms"] = lags.mean()
t_stat, p_val = stats.ttest_1samp(lags, popmean=0, alternative="less")
group_lags[f"group_{i}_p_value"] = p_val
group_lags[f"group_{i}_t_stat"] = t_stat
group_lags[f"group_{i}_sample_size"] = len(lags)
group_lags["hc_group_order"] = hc_group_order
group_lags["mec_group_order"] = mec_group_order
return group_lags
def show_mec_hc_group_average_lag(
ds_dicts,
stats_all,
range_dicts,
condition,
threshold,
random_seed=42,
window_size=15,
colors=["darkblue", "crimson"],
):
np.random.seed(random_seed)
res = get_mec_hc_group_average_lag(
ds_dicts,
stats_all,
range_dicts,
condition=condition,
threshold=threshold,
window_size=window_size,
)
fig, ax = plt.subplots(figsize=(1.5, 0.8))
for group in [0, 1]:
ax.hist(
res[f"group_{group}_peak_frame_lis"],
bins=np.arange(0, window_size * 2),
alpha=0.75,
label=f"group {group+1}",
color=colors[group],
density=True,
)
ax.axvline(window_size, color="k", linestyle="--", linewidth=0.5)
ax.set_xticks([0, window_size, window_size * 2])
ax.tick_params(axis="x", length=1, width=0.5)
ax.set_xticklabels([])
ax.set_yticks([])
ax.spines["bottom"].set_linewidth(0.25) # default is around 1.0
erase_toprightleft_axes(ax)
for group in [0, 1]:
print(f"group {group+1}")
print(f'Mean lag from HPC to MEC is {res[f"group_{group}_mean_lag_ms"]:.1f} ms')
print(f'P-value is {res[f"group_{group}_p_value"]}')
print(f'P-value is {res[f"group_{group}_p_value"]:.10f}')
print(f'T-statistic is {res[f"group_{group}_t_stat"]}')
print(f'T-statistic is {res[f"group_{group}_t_stat"]:.10f}')
print(f'Sample size is {res[f"group_{group}_sample_size"]}')
print()
def compare_strict_and_loose_theta_skipping_cells(
ds_dicts, stats_all, condition, area_name, mask, fig_dir=None
):
strict_theta_skipping_cells = stats_all[condition][area_name]["strict"][
"theta_skipping_cells"
]
only_loose_theta_skipping_cells = np.setdiff1d(
stats_all[condition][area_name]["loose"]["theta_skipping_cells"],
stats_all[condition][area_name]["strict"]["theta_skipping_cells"],
)
fig = show_all_theta_skipping_cells_v2(
ds_dicts[area_name], strict_theta_skipping_cells, mask, n_cols=6
)
fig.suptitle(f"{condition} - {area_name} - Strict theta skipping cells", y=1.05)
if fig_dir is not None:
fig.savefig(
fig_dir / f"{condition}_{area_name}_strict_theta_skipping_cells.png",
dpi=600,
)
plt.close(fig)
fig = show_all_theta_skipping_cells_v2(
ds_dicts[area_name], only_loose_theta_skipping_cells, mask, n_cols=6
)
fig.suptitle(f"{condition} - {area_name} - Only loose theta skipping cells", y=1.05)
if fig_dir is not None:
fig.savefig(
fig_dir / f"{condition}_{area_name}_only_loose_theta_skipping_cells.png",
dpi=600,
)
plt.close(fig)
def show_firing_rate_maps(ds, clusters, n_cols=8, title=False):
n_clusters = len(clusters)
n_rows = (n_clusters + n_cols - 1) // n_cols
fig = plt.figure(figsize=(n_cols * 1, n_rows * 1))
for i, cluster_idx in enumerate(clusters):
ax = fig.add_subplot(n_rows, n_cols, i + 1)
ax.imshow(
ds.firing_rate_map_mask.sel(cluster=cluster_idx).T[::-1],
interpolation="none",
)
ax.axis("off")
if title:
ax.set_title(f"clu {cluster_idx}", fontsize=6)
def get_stats_all(
ds_dicts,
range_dicts,
random_seed=2025,
):
stats_all = {}
for condition in range_dicts.keys():
print(condition)
mask = get_range_mask_from_ds(ds_dicts["hc"], range_dicts[condition])
stats_all[condition] = {
"hc": {"loose": {}, "strict": {}},
"mec": {"loose": {}, "strict": {}},
}
for area_name in ["hc", "mec"]:
print(f" {area_name}")
ds_ = ds_dicts[area_name]
theta_skipping_cell_result = (
get_theta_skipping_cells_with_and_wo_hard_threshold(
ds_, mask, n_samples=1000, hard_threshold_theta_skipping_index=0.5
)
)
for threshold in ["loose", "strict"]:
print(f" {area_name} - {threshold}")
theta_skipping_cells = theta_skipping_cell_result[
f"{threshold}_theta_skipping_cells"
]
affinity_mat = get_affinity_matrix(
ds_, mask, theta_skipping_cells, window_size=40
)
clustering_result = cluster_affinity_matrix(
affinity_mat, cluster_lis=theta_skipping_cells
)
if area_name == "hc":
fmaps = ds_.firing_rate.sel(xy_bin=ds_.track_indices.values).sel(
cluster=theta_skipping_cells
)
elif area_name == "mec":
fmaps = ds_.hd_tuning_curve.sel(cluster=theta_skipping_cells)
corr_map = np.corrcoef(fmaps)
result = (
permutation_test_for_within_across_group_representation_similarity(
corr_map, clustering_result["labels"], random_seed=random_seed
)
)
stats_ = {
"theta_skipping_cells": theta_skipping_cells,
"affinity_mat": affinity_mat,
"fmaps": fmaps,
"clustering_result": clustering_result,
"firing_map_correlation_result": result,
"theta_skipping_cell_result": theta_skipping_cell_result,
"mask": mask,
}
stats_all[condition][area_name][threshold] = stats_
return stats_all
## Process
np.random.seed(2025)
stats_all = get_stats_all(ds_dicts, range_dicts)
stats_supp = get_stats_all(ds_dicts, range_dicts_supp)
## Plotting figures
# for condition in range_dicts.keys():
# mask = get_range_mask_from_ds(ds, range_dicts[condition])
# # for threshold in ['loose', 'strict']:
# for threshold in ['strict']:
# for area_name in ['hc', 'mec']:
# ds_ = ds_dicts[area_name]
# stats_ = stats_all[condition][area_name][threshold]# All theta skipping cells
# fig = show_all_theta_skipping_cells_v2(ds_, stats_['theta_skipping_cells'], mask, n_cols=6)
# fig.suptitle(f'{condition} - {area_name} - {threshold} theta skipping cells', y=1.02)
# labels_dict = stats_['clustering_result']['labels_dict']
# group_order = get_group_direction(ds_, mask, labels_dict, area_name)
# # Theta skipping correlation matrix
# # Fig2g left, Fig3f left, Fig4f left
# fig, ax = plt.subplots(figsize=(3, 3))
# show_correlation_matrix((stats_['affinity_mat'] + stats_['affinity_mat'].T) / 2, stats_['theta_skipping_cells'],
# stats_['clustering_result']['labels'], ax=ax,
# group_order=group_order
# )
# # Summary figure for the clusters
# # Fig2g right, Fig3f right, Fig4f right
# show_group_response(ds_, area_name, mask, stats_['clustering_result']['labels_dict'], group_order=group_order)
# # Time lag
# # Fig2h, Fig3g, Fig4g
# show_mec_hc_group_average_lag(ds_dicts, stats_all, condition, threshold, window_size=20)
### Main figures
condition = "outbound_central_arm"
mask = get_range_mask_from_ds(ds, range_dicts[condition])
threshold = "strict"
for area_name in ["hc", "mec"]:
ds_ = ds_dicts[area_name]
stats_ = stats_all[condition][area_name][threshold] # All theta skipping cells
labels_dict = stats_["clustering_result"]["labels_dict"]
group_order = get_group_direction(ds_, mask, labels_dict, area_name)
# Theta skipping correlation matrix
# Fig2g left, Fig3f left, Fig4f left
show_correlation_matrix(
(stats_["affinity_mat"] + stats_["affinity_mat"].T) / 2,
stats_["theta_skipping_cells"],
stats_["clustering_result"]["labels"],
group_order=group_order,
plot_labels=False,
)
plt.savefig(
fig_dir / f"Fig2-g-left-{area_name}_skipping_cells_affinity_matrix.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
# Summary figure for the clusters
# Fig2g right, Fig3f right, Fig4f right
show_group_response(
ds_,
area_name,
mask,
stats_["clustering_result"]["labels_dict"],
group_order=group_order,
)
plt.savefig(
fig_dir / f"Fig2-g-right-{area_name}_group_response.png",
transparent=True,
bbox_inches="tight",
format="png",
dpi=600,
)
# Time lag
# Fig2h, Fig3g, Fig4g
show_mec_hc_group_average_lag(
ds_dicts, stats_all, range_dicts, condition, threshold, window_size=15
)
plt.savefig(
fig_dir / "Fig2-h-hpc_mec_cycling_cross_correlation_peak_frames.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
# Fig 3
condition = "inbound_from_left_arm"
mask = get_range_mask_from_ds(ds, range_dicts[condition])
threshold = "strict"
for area_name in ["hc", "mec"]:
ds_ = ds_dicts[area_name]
stats_ = stats_all[condition][area_name][threshold] # All theta skipping cells
labels_dict = stats_["clustering_result"]["labels_dict"]
group_order = get_group_direction(ds_, mask, labels_dict, area_name)
# Theta skipping correlation matrix
# Fig2g left, Fig3f left, Fig4f left
show_correlation_matrix(
(stats_["affinity_mat"] + stats_["affinity_mat"].T) / 2,
stats_["theta_skipping_cells"],
stats_["clustering_result"]["labels"],
group_order=group_order,
plot_labels=False,
)
plt.savefig(
fig_dir / f"Fig3-f-left-{area_name}_skipping_cells_affinity_matrix.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
# Summary figure for the clusters
# Fig2g right, Fig3f right, Fig4f right
show_group_response(
ds_,
area_name,
mask,
stats_["clustering_result"]["labels_dict"],
group_order=group_order,
)
plt.savefig(
fig_dir / f"Fig3-f-right-{area_name}_group_response.png",
transparent=True,
bbox_inches="tight",
format="png",
dpi=600,
)
# Time lag
# Fig2h, Fig3g, Fig4g
show_mec_hc_group_average_lag(
ds_dicts, stats_all, range_dicts, condition, threshold, window_size=15
)
plt.savefig(
fig_dir / "Fig3-g-hpc_mec_cycling_cross_correlation_peak_frames.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
# Fig 4
condition = "outbound_topleft_corner"
mask = get_range_mask_from_ds(ds, range_dicts[condition])
threshold = "strict"
for area_name in ["hc", "mec"]:
ds_ = ds_dicts[area_name]
stats_ = stats_all[condition][area_name][threshold] # All theta skipping cells
labels_dict = stats_["clustering_result"]["labels_dict"]
group_order = get_group_direction(ds_, mask, labels_dict, area_name)
# Theta skipping correlation matrix
# Fig2g left, Fig3f left, Fig4f left
show_correlation_matrix(
(stats_["affinity_mat"] + stats_["affinity_mat"].T) / 2,
stats_["theta_skipping_cells"],
stats_["clustering_result"]["labels"],
group_order=group_order,
plot_labels=False,
)
plt.savefig(
fig_dir / f"Fig4-f-left-{area_name}_skipping_cells_affinity_matrix.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
# Summary figure for the clusters
# Fig2g right, Fig3f right, Fig4f right
show_group_response(
ds_,
area_name,
mask,
stats_["clustering_result"]["labels_dict"],
group_order=group_order,
)
plt.savefig(
fig_dir / f"Fig4-f-right-{area_name}_group_response.png",
transparent=True,
bbox_inches="tight",
format="png",
dpi=600,
)
# Time lag
# Fig2h, Fig3g, Fig4g
show_mec_hc_group_average_lag(
ds_dicts, stats_all, range_dicts, condition, threshold, window_size=15
)
plt.savefig(
fig_dir / "Fig4-g-hpc_mec_cycling_cross_correlation_peak_frames.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
### Supp Figures
## Supp Fig2c hpc clustering results, individual rate maps
condition = "outbound_central_arm"
for group in [0, 1]:
cells = stats_all[condition]["hc"]["strict"]["clustering_result"]["labels_dict"][
group
]
show_firing_rate_maps(ds, cells, n_cols=len(cells))
plt.savefig(
fig_dir / f"SuppFig2-c-hpc-group{group+1}-firing-rate-maps.png",
dpi=600,
format="png",
bbox_inches="tight",
transparent=True,
)
## Supp Fig3a hpc clustering results, individual rate maps
condition = "inbound_from_left_arm"
for group in [0, 1]:
cells = stats_all[condition]["hc"]["strict"]["clustering_result"]["labels_dict"][
group
]
show_firing_rate_maps(ds, cells, n_cols=len(cells))
plt.savefig(
fig_dir / f"SuppFig3-a-hpc-group{group+1}-firing-rate-maps.png",
dpi=600,
format="png",
bbox_inches="tight",
transparent=True,
)
# Other corners
supp_figures = {
"SuppFig3-b": "inbound_from_right_arm",
"SuppFig4-b": "inbound_topleft_corner",
"SuppFig4-c": "outbound_topright_corner",
"SuppFig4-d": "inbound_topright_corner",
}
threshold = "strict"
# for condition in range_dicts_supp.keys():
for supp_fig, condition in supp_figures.items():
mask = get_range_mask_from_ds(ds, range_dicts_supp[condition])
for area_name in ["hc", "mec"]:
ds_ = ds_dicts[area_name]
stats_ = stats_supp[condition][area_name][threshold] # All theta skipping cells
labels_dict = stats_["clustering_result"]["labels_dict"]
group_order = get_group_direction(ds_, mask, labels_dict, area_name)
# Theta skipping correlation matrix
# Fig2g left, Fig3f left, Fig4f left
show_correlation_matrix(
(stats_["affinity_mat"] + stats_["affinity_mat"].T) / 2,
stats_["theta_skipping_cells"],
stats_["clustering_result"]["labels"],
group_order=group_order,
plot_labels=False,
)
plt.savefig(
fig_dir
/ f"{supp_fig}-{condition}-{area_name}_skipping_cells_affinity_matrix.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
# Summary figure for the clusters
# Fig2g right, Fig3f right, Fig4f right
show_group_response(
ds_,
area_name,
mask,
stats_["clustering_result"]["labels_dict"],
group_order=group_order,
)
plt.savefig(
fig_dir / f"{supp_fig}-{condition}-{area_name}_group_response.png",
transparent=True,
bbox_inches="tight",
format="png",
dpi=600,
)
# Time lag
# Fig2h, Fig3g, Fig4g
show_mec_hc_group_average_lag(
ds_dicts, stats_supp, range_dicts_supp, condition, threshold, window_size=15
)
plt.savefig(
fig_dir
/ f"{supp_fig}-{condition}-hpc_mec_cycling_cross_correlation_peak_frames.pdf",
dpi=600,
format="pdf",
bbox_inches="tight",
transparent=True,
)
###########
## stats ##
###########
import sys
# open file in write mode (or append 'a')
logfile = open(str(fig_dir / "stats.txt"), "w")
sys.stdout = logfile
#####################################################
### Left-right alternation of population decoding ###
#####################################################
print("#####################################################")
print("### Left-right alternation of population decoding ###")
print("#####################################################")
condition = "outbound_central_arm"
print("")
print("Condition:", condition)
run_shuffling_test_for_latent_position_cycling(
ds,
ML_modes,
bin_labels,
range_dicts[condition],
target_arm_1=1,
target_arm_2=2,
window_size=40,
n_samples=1000,
show_plots=False,
)
# This should output:
# Real skipping index: -0.4034334763948498
# 5th percentile of shuffled skipping index: -0.08027716910775382
# Is cycling between latent positions: True
condition = "inbound_from_left_arm"
print("")
print("Condition:", condition)
run_shuffling_test_for_latent_position_cycling(
ds,
ML_modes,
bin_labels,
range_dicts[condition],
target_arm_1=2,
target_arm_2=0,
window_size=40,
n_samples=1000,
show_plots=False,
)
#####################################################
### Clustering analysis of theta skipping cells ###
#####################################################
window_size = 15
for condition in range_dicts.keys():
for threshold in ["loose", "strict"]:
print("\n-------------------------")
print(f"Condition: {condition}, Threshold: {threshold}\n")
for area_name in ["hc", "mec"]:
stats_ = stats_all[condition][area_name][threshold]
print(f"{area_name} - {threshold}:")
print(
f'Number of theta skipping cells: {len(stats_["theta_skipping_cells"])}'
)
print(
f'Firing map correlation within-across group p-value: {stats_["firing_map_correlation_result"]["p"]:.4f}'
)
res = get_mec_hc_group_average_lag(
ds_dicts,
stats_all,
range_dicts,
condition=condition,
threshold=threshold,
window_size=window_size,
)
for group in [0, 1]:
print(f"group {group+1}")
print(
f'Mean lag from HPC to MEC is {res[f"group_{group}_mean_lag_ms"]:.1f} ms'
)
print(f'P-value is {res[f"group_{group}_p_value"]}')
print(f'P-value is {res[f"group_{group}_p_value"]:.10f}')
print(f'T-statistic is {res[f"group_{group}_t_stat"]}')
print(f'T-statistic is {res[f"group_{group}_t_stat"]:.10f}')
print(f'Sample size is {res[f"group_{group}_sample_size"]}')
print()
for condition in range_dicts_supp.keys():
print("\n-------------------------")
for threshold in ["loose", "strict"]:
print("\n-------------------------")
print(f"Condition: {condition}, Threshold: {threshold}\n")
for area_name in ["hc", "mec"]:
stats_ = stats_supp[condition][area_name][threshold]
print(f"{area_name} - {threshold}:")
print(
f'Number of theta skipping cells: {len(stats_["theta_skipping_cells"])}'
)
print(
f'Firing map correlation within-across group p-value: {stats_["firing_map_correlation_result"]["p"]:.4f}'
)
res = get_mec_hc_group_average_lag(
ds_dicts,
stats_supp,
range_dicts_supp,
condition=condition,
threshold=threshold,
window_size=window_size,
)
for group in [0, 1]:
print(f"group {group+1}")
print(
f'Mean lag from HPC to MEC is {res[f"group_{group}_mean_lag_ms"]:.1f} ms'
)
print(f'P-value is {res[f"group_{group}_p_value"]}')
print(f'P-value is {res[f"group_{group}_p_value"]:.10f}')
print(f'T-statistic is {res[f"group_{group}_t_stat"]}')
print(f'T-statistic is {res[f"group_{group}_t_stat"]:.10f}')
print(f'Sample size is {res[f"group_{group}_sample_size"]}')
print()
logfile.close()

analysis-2.py, under CC-BY-4.0 · at the source

Overview

Authors: Masahiro Nakano¹, Caswell Barry², Claudia Clopath^1,3

ORCID iDs: Masahiro Nakano, Caswell Barry

Institutions: Sainsbury Wellcome Centre (United Kingdom); University College London (United Kingdom); Imperial College London (United Kingdom)

Journal: Hippocampus, volume 36, issue 5, article e70131

Dates: received 13 February 2026; accepted 26 August 2026; published online 18 September 2026; in print September 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1002/hipo.70131 · PMID 42760272 · PMCID PMC13588937 · OpenAlex W7213600936

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: systems (subfield)

Methods: Connectivity, Statistics, Machine learning, Single-unit activity, calcium imaging

Keywords: Bayesian decoding, Hippocampus, medial entorhinal cortex, planning, theta sequences, theta skipping

MeSH: Entorhinal Cortex*, Hippocampus*, Theta Rhythm*, Animals, Bayes Theorem, Decision Making, Male, Maze Learning (* major topic)

Topic: Memory and Neural Mechanisms (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: RCUK | Biotechnology and Biological Sciences Research Council (BBSRC) (BB/Y011775/1); Wellcome Trust (200790/Z/16/Z); EC | EU Framework Programme for Research and Innovation H2020 | H2020 Priority Excellent Science | H2020 European Research Council (H2020 Excellent Science - European Research Council) (MotorAdapt 101169605); Simons Foundation (564408); Engineering and Physical Sciences Research Council (EP/R035806/1)

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

Abstract

The hippocampus is central to memory and spatial planning. A prominent candidate mechanism supporting navigational decision making is hippocampal theta sequences—brief (~120 ms) place‐cell sequences representing future trajectories—which have been interpreted as planning signals based on their alternation between maze arms in T‐maze tasks. However, recent work suggests that intrinsic left–right alternation in medial entorhinal cortex (MEC) theta sequences may underlie this phenomenon. Here, we test this hypothesis using a model of MEC theta dynamics in a T‐maze and show that standard Bayesian decoding yields hippocampal theta sequences that alternate between arms. We then analyzed existing simultaneous MEC‐hippocampal recordings from Vollan et al. (2025). We found tight coupling between MEC and hippocampus, along with coordinated left–right alternation, while the animal was performing an alternation task. Notably, this conjoint alternation occurs not only at choice points but also during inbound and forced‐turn trials. These findings suggest that left–right alternation in hippocampal theta sequences is a ubiquitous phenomenon driven by MEC dynamics, challenging the interpretation that these specifically reflect intentional decision‐making or planning processes only.

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 8 matches between paragraphs and lines of code.

Zenodo 18340856

License: CC-BY-4.0

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Size: 1 file

Software Heritage: not checked

Found in: “Data Availability Statement”

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

Tools: NumPy (17 files), Matplotlib (10 files), JAX (7 files), SciPy (7 files), pandas (5 files), scikit-learn (2 files), xarray (2 files), h5py (1 file), Numba (1 file), OpenCV (1 file)

Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)

26 September 2026: the link answers (HTTP 200)

20 files

zenodo-sciad-share/scripts/analysis-1.py, Python, 1,762 lines, 2 matches
zenodo-sciad-share/scripts/analysis-2.py, Python, 1,664 lines, 2 matches
zenodo-sciad-share/scripts/simulation.py, Python, 635 lines
zenodo-sciad-share/shared_src/__init__.py, Python, 1 line
zenodo-sciad-share/shared_src/decoding_utils.py, Python, 99 lines
zenodo-sciad-share/shared_src/fig_utils.py, Python, 184 lines
zenodo-sciad-share/shared_src/figure/__init__.py, Python, 1 line
zenodo-sciad-share/shared_src/figure/color.py, Python, 86 lines
zenodo-sciad-share/shared_src/figure/ratinabox_utils.py, Python, 113 lines
zenodo-sciad-share/shared_src/figure/trajectory.py, Python, 1,175 lines
zenodo-sciad-share/shared_src/figure/utils.py, Python, 160 lines
zenodo-sciad-share/shared_src/kalmax_utils.py, Python, 714 lines
zenodo-sciad-share/shared_src/kde.py, Python, 570 lines, 1 match
zenodo-sciad-share/shared_src/kde_fast.py, Python, 527 lines, 1 match
zenodo-sciad-share/shared_src/linearlize.py, Python, 287 lines
zenodo-sciad-share/shared_src/ratinabox_utils.py, Python, 314 lines, 2 matches
zenodo-sciad-share/shared_src/spike_autocorrelation.py, Python, 156 lines
zenodo-sciad-share/shared_src/utils.py, Python, 285 lines
zenodo-sciad-share/src/dataset.py, Python, 511 lines
zenodo-sciad-share/src/figure/ds_fig.py, Python, 275 lines

At the source: doi.org/10.5281/zenodo.18340856

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;
20 scripts, each with its path and the digest of its content;
8 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

Datasets cited

doi:10.25493/r5fr-edg, at the source; found in “Data Availability Statement”

Data Availability Statement

Data analyzed in this paper are available at EBRAINS, 10.25493/R5FR‐EDG (https://doi.org/10.25493/R5FR-EDG). Code for reproducing the simulation and analyses in this article will be available at Zenodo 10.5281/zenodo.18340856 (https://doi.org/10.5281/zenodo.18340856).

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 3, 28 September 2026

Publisher: n/a → Wiley

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 6 keywords, 8 MeSH terms, 5 funders, 75 references.

Cite

This paper

Nakano, M., Barry, C., & Clopath, C. (2026). Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences. Hippocampus, 36(5), e70131. https://doi.org/10.1002/hipo.70131

BibTeX

@article{nakano2026decoding, author = {Nakano, Masahiro and Barry, Caswell and Clopath, Claudia}, title = {{Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences}}, journal = {Hippocampus}, year = {2026}, month = sep, volume = {36}, number = {5}, pages = {e70131}, publisher = {Wiley}, issn = {1050-9631}, doi = {10.1002/hipo.70131}, url = {https://doi.org/10.1002/hipo.70131}, pmid = {42760272}, pmcid = {PMC13588937} }

RIS

TY - JOUR AU - Nakano, Masahiro AU - Barry, Caswell AU - Clopath, Claudia TI - Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences T2 - Hippocampus J2 - Hippocampus PY - 2026 DA - 2026/09/01 VL - 36 IS - 5 SP - e70131 SN - 1050-9631 PB - Wiley DO - 10.1002/hipo.70131 UR - https://doi.org/10.1002/hipo.70131 LA - en ER -

CSL-JSON

{ "id": "10.1002/hipo.70131", "type": "article-journal", "title": "Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences", "container-title": "Hippocampus", "author": [ { "family": "Nakano", "given": "Masahiro" }, { "family": "Barry", "given": "Caswell" }, { "family": "Clopath", "given": "Claudia" } ], "container-title-short": "Hippocampus", "volume": "36", "issue": "5", "page": "e70131", "DOI": "10.1002/hipo.70131", "PMID": "42760272", "PMCID": "PMC13588937", "ISSN": "1050-9631", "publisher": "Wiley", "URL": "https://doi.org/10.1002/hipo.70131", "language": "en", "issued": { "date-parts": [ [ 2026, 9, 1 ] ] } }

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.