Bayesian Position Decoding¶

Estimated time: 35-45 minutes

Prerequisites: 11_place_field_analysis.ipynb, 15_simulation_workflows.ipynb

Learning Objectives¶

By the end of this notebook, you will be able to:

• Decode spatial position in one call with decode_session (the golden path)
• Drop down to the manual three-call path when you need custom control
• Build encoding models (place fields) for a population of neurons
• Access and interpret DecodingResult properties (posterior, MAP, mean, posterior entropy)
• Evaluate decoding accuracy with error metrics
• Detect trajectory structure using isotonic/linear regression and Radon transform
• Test significance of decoded sequences using shuffle-based methods

Estimated time: 25-30 minutes

---

Setup¶

First, let's import the necessary libraries:

In [1]:

import matplotlib.pyplot as plt
import numpy as np

from neurospatial import Environment
from neurospatial.decoding import (
    bin_spikes_in_time,
    confusion_matrix,
    decode_position,
    decode_session,
    decoding_correlation,
    decoding_error,
    fit_isotonic_trajectory,
    fit_linear_trajectory,
    median_decoding_error,
)
from neurospatial.encoding import compute_spatial_rate, compute_spatial_rates
from neurospatial.simulation import (
    PlaceCellModel,
    generate_poisson_spikes,
    simulate_trajectory_ou,
)
from neurospatial.stats.shuffle import compute_shuffle_pvalue, shuffle_time_bins

# Set random seed for reproducibility
np.random.seed(42)

# Shared styling (Okabe-Ito palette, consistent figure / font sizes)
import sys
from pathlib import Path

_here = (
    str(Path(__file__).resolve().parent) if "__file__" in globals() else str(Path.cwd())
)
if _here not in sys.path:
    sys.path.insert(0, _here)
from _style import apply_style

apply_style(figsize=(12, 5), font_size=12)

# Colorblind-friendly palette
COLORS = {
    "blue": "#0173B2",
    "orange": "#DE8F05",
    "green": "#029E73",
    "red": "#CC78BC",
    "cyan": "#56B4E9",
}

import matplotlib.pyplot as plt import numpy as np from neurospatial import Environment from neurospatial.decoding import ( bin_spikes_in_time, confusion_matrix, decode_position, decode_session, decoding_correlation, decoding_error, fit_isotonic_trajectory, fit_linear_trajectory, median_decoding_error, ) from neurospatial.encoding import compute_spatial_rate, compute_spatial_rates from neurospatial.simulation import ( PlaceCellModel, generate_poisson_spikes, simulate_trajectory_ou, ) from neurospatial.stats.shuffle import compute_shuffle_pvalue, shuffle_time_bins # Set random seed for reproducibility np.random.seed(42) # Shared styling (Okabe-Ito palette, consistent figure / font sizes) import sys from pathlib import Path _here = ( str(Path(__file__).resolve().parent) if "__file__" in globals() else str(Path.cwd()) ) if _here not in sys.path: sys.path.insert(0, _here) from _style import apply_style apply_style(figsize=(12, 5), font_size=12) # Colorblind-friendly palette COLORS = { "blue": "#0173B2", "orange": "#DE8F05", "green": "#029E73", "red": "#CC78BC", "cyan": "#56B4E9", }

---

Part 1: Generate Synthetic Data¶

We'll create a 1D linear track environment with a population of place cells. This simplified setup makes it easy to visualize and understand decoding.

In [2]:

# Create a 1D linear track (100 cm)
track_length = 100.0  # cm
bin_size = 2.0  # cm per bin

# Create 1D positions along track
positions_1d = np.linspace(0, track_length, 51).reshape(-1, 1)
env = Environment.from_samples(positions_1d, bin_size=bin_size)
env.units = "cm"

print(f"Track length: {track_length} cm")
print(f"Number of spatial bins: {env.n_bins}")
print(f"Bin size: {bin_size} cm")

# Create a 1D linear track (100 cm) track_length = 100.0 # cm bin_size = 2.0 # cm per bin # Create 1D positions along track positions_1d = np.linspace(0, track_length, 51).reshape(-1, 1) env = Environment.from_samples(positions_1d, bin_size=bin_size) env.units = "cm" print(f"Track length: {track_length} cm") print(f"Number of spatial bins: {env.n_bins}") print(f"Bin size: {bin_size} cm")

Track length: 100.0 cm
Number of spatial bins: 51
Bin size: 2.0 cm

Generate Animal Trajectory¶

Simulate a rat running back and forth on the linear track:

In [3]:

# Generate smooth trajectory on track using OU process
duration = 600.0  # 10 minutes
positions, times = simulate_trajectory_ou(
    env,
    duration=duration,
    dt=0.01,  # 10ms timestep
    speed_mean=15.0,  # cm/s
    speed_std=2.0,
    coherence_time=0.5,
    boundary_mode="reflect",
    seed=42,
    speed_units="cm",
)

print(f"Trajectory duration: {times[-1]:.1f} seconds")
print(f"Number of samples: {len(times)}")
print(f"Sampling rate: {1 / (times[1] - times[0]):.0f} Hz")
print(f"Position range: [{positions.min():.1f}, {positions.max():.1f}] cm")

# Generate smooth trajectory on track using OU process duration = 600.0 # 10 minutes positions, times = simulate_trajectory_ou( env, duration=duration, dt=0.01, # 10ms timestep speed_mean=15.0, # cm/s speed_std=2.0, coherence_time=0.5, boundary_mode="reflect", seed=42, speed_units="cm", ) print(f"Trajectory duration: {times[-1]:.1f} seconds") print(f"Number of samples: {len(times)}") print(f"Sampling rate: {1 / (times[1] - times[0]):.0f} Hz") print(f"Position range: [{positions.min():.1f}, {positions.max():.1f}] cm")

Trajectory duration: 600.0 seconds
Number of samples: 60000
Sampling rate: 100 Hz
Position range: [-1.0, 101.0] cm

In [4]:

# Visualize trajectory
fig, ax = plt.subplots(figsize=(14, 4))
ax.plot(times[:5000], positions[:5000, 0], color=COLORS["blue"], linewidth=0.5)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Position (cm)")
ax.set_title("First 50 seconds of trajectory", fontweight="bold")
ax.set_xlim(0, 50)
plt.tight_layout()
plt.show()

# Visualize trajectory fig, ax = plt.subplots(figsize=(14, 4)) ax.plot(times[:5000], positions[:5000, 0], color=COLORS["blue"], linewidth=0.5) ax.set_xlabel("Time (s)") ax.set_ylabel("Position (cm)") ax.set_title("First 50 seconds of trajectory", fontweight="bold") ax.set_xlim(0, 50) plt.tight_layout() plt.show()

/tmp/claude-501/ipykernel_99754/3936418870.py:9: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

Create Place Cell Population¶

Generate a population of place cells with fields distributed along the track:

In [5]:

# Create place cells with fields distributed along the track
n_neurons = 30
field_centers = np.linspace(5, track_length - 5, n_neurons)  # Spread across track
field_widths = np.random.uniform(8, 15, n_neurons)  # Variable field widths
peak_rates = np.random.uniform(10, 30, n_neurons)  # Variable peak firing rates

# Create PlaceCellModel for each neuron and generate spikes
spike_times_list = []
place_cells = []

for i in range(n_neurons):
    cell = PlaceCellModel(
        env,
        center=np.array([field_centers[i]]),
        width=field_widths[i],
        max_rate=peak_rates[i],
        baseline_rate=0.1,
    )
    place_cells.append(cell)

    # Generate spike train
    rates = cell.firing_rate(positions, times)
    spikes = generate_poisson_spikes(rates, times, refractory_period=0.002, seed=42 + i)
    spike_times_list.append(spikes)

print(f"Created {n_neurons} place cells")
print(f"Total spikes: {sum(len(s) for s in spike_times_list)}")
print(f"Mean spikes per neuron: {np.mean([len(s) for s in spike_times_list]):.0f}")

# Create place cells with fields distributed along the track n_neurons = 30 field_centers = np.linspace(5, track_length - 5, n_neurons) # Spread across track field_widths = np.random.uniform(8, 15, n_neurons) # Variable field widths peak_rates = np.random.uniform(10, 30, n_neurons) # Variable peak firing rates # Create PlaceCellModel for each neuron and generate spikes spike_times_list = [] place_cells = [] for i in range(n_neurons): cell = PlaceCellModel( env, center=np.array([field_centers[i]]), width=field_widths[i], max_rate=peak_rates[i], baseline_rate=0.1, ) place_cells.append(cell) # Generate spike train rates = cell.firing_rate(positions, times) spikes = generate_poisson_spikes(rates, times, refractory_period=0.002, seed=42 + i) spike_times_list.append(spikes) print(f"Created {n_neurons} place cells") print(f"Total spikes: {sum(len(s) for s in spike_times_list)}") print(f"Mean spikes per neuron: {np.mean([len(s) for s in spike_times_list]):.0f}")

Created 30 place cells
Total spikes: 78049
Mean spikes per neuron: 2602

In [6]:

# Visualize a few place fields
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

for idx, ax in enumerate(axes):
    cell_idx = idx * (n_neurons // 3)
    result_cell = compute_spatial_rate(
        env,
        spike_times_list[cell_idx],
        times,
        positions,
        smoothing_method="diffusion_kde",
        bandwidth=5.0,
    )
    field = result_cell.firing_rate

    # Plot field
    x_positions = env.bin_centers[:, 0]
    ax.bar(x_positions, field, width=bin_size * 0.8, color=COLORS["blue"], alpha=0.7)
    ax.axvline(
        field_centers[cell_idx],
        color=COLORS["red"],
        linestyle="--",
        linewidth=2,
        label="True center",
    )
    ax.set_xlabel("Position (cm)")
    ax.set_ylabel("Firing rate (Hz)")
    ax.set_title(f"Neuron {cell_idx + 1}", fontweight="bold")
    ax.legend()

plt.suptitle("Example Place Fields", fontsize=14, fontweight="bold")
plt.tight_layout()
plt.show()

# Visualize a few place fields fig, axes = plt.subplots(1, 3, figsize=(15, 4)) for idx, ax in enumerate(axes): cell_idx = idx * (n_neurons // 3) result_cell = compute_spatial_rate( env, spike_times_list[cell_idx], times, positions, smoothing_method="diffusion_kde", bandwidth=5.0, ) field = result_cell.firing_rate # Plot field x_positions = env.bin_centers[:, 0] ax.bar(x_positions, field, width=bin_size * 0.8, color=COLORS["blue"], alpha=0.7) ax.axvline( field_centers[cell_idx], color=COLORS["red"], linestyle="--", linewidth=2, label="True center", ) ax.set_xlabel("Position (cm)") ax.set_ylabel("Firing rate (Hz)") ax.set_title(f"Neuron {cell_idx + 1}", fontweight="bold") ax.legend() plt.suptitle("Example Place Fields", fontsize=14, fontweight="bold") plt.tight_layout() plt.show()

/Users/edeno/Documents/GitHub/neurospatial/.venv/lib/python3.13/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py:606: SparseEfficiencyWarning: splu converted its input to CSC format
  return splu(A).solve
/Users/edeno/Documents/GitHub/neurospatial/.venv/lib/python3.13/site-packages/scipy/sparse/linalg/_matfuncs.py:707: SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
  return spsolve(Q, P)
/tmp/claude-501/ipykernel_99754/2985946172.py:33: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

---

Part 2: Decode in One Call (the golden path)¶

For most analyses you do not need to wire up the pieces by hand. decode_session is the one-call golden path: give it the environment, the population spike trains, the timestamps, and the positions, and it builds the encoding models (place fields), bins the spikes onto a regular time grid, and runs the Bayesian decoder for you — returning a DecodingResult.

In [7]:

# Decode position from spikes + behavior in a single call.
result = decode_session(
    env,
    spike_times_list,
    times,
    positions,
    dt=0.1,  # 100 ms time bins
    bandwidth=5.0,  # place-field smoothing bandwidth (cm)
    min_occupancy=0.5,  # match the manual path below so the two agree exactly
)

print("Decoding complete (one call)!")
print(f"Posterior shape: {result.posterior.shape}  (n_time_bins, n_bins)")
print(f"MAP position shape: {result.map_position.shape}")

# Evaluate: the decoded MAP position should track the true trajectory.
actual_track = np.interp(result.times, times, positions[:, 0]).reshape(-1, 1)
median_err = median_decoding_error(result.map_position, actual_track)
print(
    f"Median decoding error: {median_err:.1f} cm (track length {track_length:.0f} cm)"
)

# Decode position from spikes + behavior in a single call. result = decode_session( env, spike_times_list, times, positions, dt=0.1, # 100 ms time bins bandwidth=5.0, # place-field smoothing bandwidth (cm) min_occupancy=0.5, # match the manual path below so the two agree exactly ) print("Decoding complete (one call)!") print(f"Posterior shape: {result.posterior.shape} (n_time_bins, n_bins)") print(f"MAP position shape: {result.map_position.shape}") # Evaluate: the decoded MAP position should track the true trajectory. actual_track = np.interp(result.times, times, positions[:, 0]).reshape(-1, 1) median_err = median_decoding_error(result.map_position, actual_track) print( f"Median decoding error: {median_err:.1f} cm (track length {track_length:.0f} cm)" )

Decoding complete (one call)!
Posterior shape: (5999, 51)  (n_time_bins, n_bins)
MAP position shape: (5999, 1)
Median decoding error: 2.5 cm (track length 100 cm)

In [8]:

# Plot the decoded posterior with the true trajectory overlaid.
fig, ax = plt.subplots(figsize=(14, 5))
n_show = 500
result.plot(ax=ax, show_map=True, colorbar=True)
ax.plot(
    result.times[:n_show],
    actual_track[:n_show, 0],
    color=COLORS["cyan"],
    linewidth=2,
    linestyle="--",
    label="Actual position",
)
ax.set_xlim(0, n_show * 0.1)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Position (cm)")
ax.set_title("decode_session: posterior, MAP (white), actual (cyan)", fontweight="bold")
ax.legend(loc="upper right")
plt.tight_layout()
plt.show()

# Plot the decoded posterior with the true trajectory overlaid. fig, ax = plt.subplots(figsize=(14, 5)) n_show = 500 result.plot(ax=ax, show_map=True, colorbar=True) ax.plot( result.times[:n_show], actual_track[:n_show, 0], color=COLORS["cyan"], linewidth=2, linestyle="--", label="Actual position", ) ax.set_xlim(0, n_show * 0.1) ax.set_xlabel("Time (s)") ax.set_ylabel("Position (cm)") ax.set_title("decode_session: posterior, MAP (white), actual (cyan)", fontweight="bold") ax.legend(loc="upper right") plt.tight_layout() plt.show()

/tmp/claude-501/ipykernel_99754/3616812306.py:19: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

That is the entire encode → bin → decode pipeline in one call. The rest of this notebook opens up the same pipeline so you can see (and customize) each stage — and then builds on it for trajectory analysis and replay significance testing.

---

Advanced: manual three-call path (custom control)¶

decode_session covers the common case. Reach for the manual three-call path (compute_spatial_rates → bin_spikes_in_time → decode_position) when you need to:

• pass custom encoding models (e.g. fields fit on a different session, or from a non-default encoder),
• reuse fitted place fields across many decodes without recomputing them,
• inspect intermediates — the occupancy map, the smoothed place fields, or the binned spike counts — for debugging or figures.

The next four parts walk through those stages explicitly. They reproduce what decode_session did above, so the decoded result matches.

Build Encoding Models¶

Encoding models are place fields that describe how each neuron's firing rate varies with spatial position. These are the "tuning curves" we use for decoding.

In [9]:

# Compute place fields for all neurons (encoding models) with the batch API.
#
# ``compute_spatial_rates`` precomputes the shared occupancy map and
# diffusion kernel once, then processes all neurons in a single pass.
# Passing ``fill_value=0.0`` replaces low-occupancy bins (NaN under the
# default ``fill_value=None``) with an explicit zero firing rate, so the
# encoding models compose directly with ``decode_position`` — no manual
# ``np.nan_to_num`` scrubbing required.
batch_result = compute_spatial_rates(
    env,
    spike_times_list,
    times,
    positions,
    smoothing_method="diffusion_kde",
    bandwidth=5.0,
    min_occupancy=0.5,
    fill_value=0.0,
)
encoding_models = batch_result.firing_rates  # shape: (n_neurons, n_bins)

print(f"Encoding models shape: {encoding_models.shape}")
print(f"  (n_neurons, n_bins) = ({n_neurons}, {env.n_bins})")
print(f"Max firing rate: {np.nanmax(encoding_models):.2f} Hz")
print(f"NaN values: {np.isnan(encoding_models).sum()}")

# Compute place fields for all neurons (encoding models) with the batch API. # # ``compute_spatial_rates`` precomputes the shared occupancy map and # diffusion kernel once, then processes all neurons in a single pass. # Passing ``fill_value=0.0`` replaces low-occupancy bins (NaN under the # default ``fill_value=None``) with an explicit zero firing rate, so the # encoding models compose directly with ``decode_position`` — no manual # ``np.nan_to_num`` scrubbing required. batch_result = compute_spatial_rates( env, spike_times_list, times, positions, smoothing_method="diffusion_kde", bandwidth=5.0, min_occupancy=0.5, fill_value=0.0, ) encoding_models = batch_result.firing_rates # shape: (n_neurons, n_bins) print(f"Encoding models shape: {encoding_models.shape}") print(f" (n_neurons, n_bins) = ({n_neurons}, {env.n_bins})") print(f"Max firing rate: {np.nanmax(encoding_models):.2f} Hz") print(f"NaN values: {np.isnan(encoding_models).sum()}")

Encoding models shape: (30, 51)
  (n_neurons, n_bins) = (30, 51)
Max firing rate: 21.62 Hz
NaN values: 0

In [10]:

# Visualize all encoding models as a heatmap
fig, ax = plt.subplots(figsize=(12, 6))

# Sort neurons by peak position for better visualization
peak_positions = np.argmax(encoding_models, axis=1)
sorted_idx = np.argsort(peak_positions)
sorted_models = encoding_models[sorted_idx]

im = ax.imshow(
    sorted_models,
    aspect="auto",
    cmap="hot",
    extent=[0, track_length, n_neurons, 0],
    vmin=0,
    vmax=np.percentile(sorted_models, 95),
)
ax.set_xlabel("Position (cm)")
ax.set_ylabel("Neuron (sorted by peak position)")
ax.set_title("Population Encoding Models (Place Fields)", fontweight="bold")
plt.colorbar(im, label="Firing rate (Hz)")
plt.tight_layout()
plt.show()

# Visualize all encoding models as a heatmap fig, ax = plt.subplots(figsize=(12, 6)) # Sort neurons by peak position for better visualization peak_positions = np.argmax(encoding_models, axis=1) sorted_idx = np.argsort(peak_positions) sorted_models = encoding_models[sorted_idx] im = ax.imshow( sorted_models, aspect="auto", cmap="hot", extent=[0, track_length, n_neurons, 0], vmin=0, vmax=np.percentile(sorted_models, 95), ) ax.set_xlabel("Position (cm)") ax.set_ylabel("Neuron (sorted by peak position)") ax.set_title("Population Encoding Models (Place Fields)", fontweight="bold") plt.colorbar(im, label="Firing rate (Hz)") plt.tight_layout() plt.show()

/tmp/claude-501/ipykernel_99754/1584133039.py:22: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

Bin Spikes for Decoding¶

Bayesian decoding works on spike counts in discrete time bins. We need to convert spike times to spike counts:

In [11]:

# Time bin parameters
dt = 0.1  # 100 ms time bins (typical for spatial decoding)

# Bin spikes for all neurons at once using the canonical helper.
# ``bin_spikes_in_time`` owns the time-grid definition so binned counts and
# bin centers are always consistent with the same grid used by decode_position.
spike_counts, time_bin_centers = bin_spikes_in_time(
    spike_times_list,
    dt,
    t_start=float(times[0]),
    t_stop=float(times[-1]),
    orient="time_x_neuron",  # shape: (n_time_bins, n_neurons) — what decode_position expects
)
n_time_bins = len(time_bin_centers)

# Get actual positions at each time bin center
actual_bin_indices = np.searchsorted(times, time_bin_centers) - 1
actual_bin_indices = np.clip(actual_bin_indices, 0, len(positions) - 1)
actual_positions = positions[actual_bin_indices]

print(f"Time bin width: {dt * 1000:.0f} ms")
print(f"Number of time bins: {n_time_bins}")
print(f"Spike counts shape: {spike_counts.shape}")
print(f"Total spikes in binned data: {spike_counts.sum()}")
print(f"Mean spikes per time bin: {spike_counts.sum(axis=1).mean():.2f}")

# Time bin parameters dt = 0.1 # 100 ms time bins (typical for spatial decoding) # Bin spikes for all neurons at once using the canonical helper. # ``bin_spikes_in_time`` owns the time-grid definition so binned counts and # bin centers are always consistent with the same grid used by decode_position. spike_counts, time_bin_centers = bin_spikes_in_time( spike_times_list, dt, t_start=float(times[0]), t_stop=float(times[-1]), orient="time_x_neuron", # shape: (n_time_bins, n_neurons) — what decode_position expects ) n_time_bins = len(time_bin_centers) # Get actual positions at each time bin center actual_bin_indices = np.searchsorted(times, time_bin_centers) - 1 actual_bin_indices = np.clip(actual_bin_indices, 0, len(positions) - 1) actual_positions = positions[actual_bin_indices] print(f"Time bin width: {dt * 1000:.0f} ms") print(f"Number of time bins: {n_time_bins}") print(f"Spike counts shape: {spike_counts.shape}") print(f"Total spikes in binned data: {spike_counts.sum()}") print(f"Mean spikes per time bin: {spike_counts.sum(axis=1).mean():.2f}")

Time bin width: 100 ms
Number of time bins: 5999
Spike counts shape: (5999, 30)
Total spikes in binned data: 78036
Mean spikes per time bin: 13.01

Decode Position (manual)¶

Now we can decode position using Bayesian inference. The decode_position() function computes the posterior probability distribution over positions for each time bin. This reproduces the decode_session result from Part 2 — now with the encoding models and binned spike counts exposed as explicit intermediates you can swap or inspect.

In [12]:

# Decode position from spike counts
result = decode_position(
    env,
    spike_counts,
    encoding_models,
    dt,
    prior=None,  # Uniform prior
    times=time_bin_centers,
)

print("Decoding complete!")
print(f"Result type: {type(result).__name__}")
print(f"Posterior shape: {result.posterior.shape}")
print(f"  (n_time_bins, n_bins) = ({result.n_time_bins}, {env.n_bins})")

# Decode position from spike counts result = decode_position( env, spike_counts, encoding_models, dt, prior=None, # Uniform prior times=time_bin_centers, ) print("Decoding complete!") print(f"Result type: {type(result).__name__}") print(f"Posterior shape: {result.posterior.shape}") print(f" (n_time_bins, n_bins) = ({result.n_time_bins}, {env.n_bins})")

Decoding complete!
Result type: DecodingResult
Posterior shape: (5999, 51)
  (n_time_bins, n_bins) = (5999, 51)

DecodingResult Properties¶

The DecodingResult container provides several useful properties (computed lazily):

In [13]:

# Access result properties
print("DecodingResult Properties:")
print(f"  posterior shape: {result.posterior.shape}")
print(f"  map_estimate shape: {result.map_estimate.shape} (bin indices)")
print(f"  map_position shape: {result.map_position.shape} (coordinates)")
print(f"  mean_position shape: {result.mean_position.shape} (coordinates)")
print(f"  uncertainty shape: {result.posterior_entropy.shape} (entropy in bits)")
print(f"  times shape: {result.times.shape if result.times is not None else 'None'}")

print(
    f"\nPosterior sum check (should be ~1.0): {result.posterior.sum(axis=1).mean():.6f}"
)
print(f"Mean uncertainty: {result.posterior_entropy.mean():.2f} bits")
print(f"Max uncertainty (uniform): {np.log2(env.n_bins):.2f} bits")

# Access result properties print("DecodingResult Properties:") print(f" posterior shape: {result.posterior.shape}") print(f" map_estimate shape: {result.map_estimate.shape} (bin indices)") print(f" map_position shape: {result.map_position.shape} (coordinates)") print(f" mean_position shape: {result.mean_position.shape} (coordinates)") print(f" uncertainty shape: {result.posterior_entropy.shape} (entropy in bits)") print(f" times shape: {result.times.shape if result.times is not None else 'None'}") print( f"\nPosterior sum check (should be ~1.0): {result.posterior.sum(axis=1).mean():.6f}" ) print(f"Mean uncertainty: {result.posterior_entropy.mean():.2f} bits") print(f"Max uncertainty (uniform): {np.log2(env.n_bins):.2f} bits")

DecodingResult Properties:
  posterior shape: (5999, 51)
  map_estimate shape: (5999,) (bin indices)
  map_position shape: (5999, 1) (coordinates)
  mean_position shape: (5999, 1) (coordinates)
  uncertainty shape: (5999,) (entropy in bits)
  times shape: (5999,)

Posterior sum check (should be ~1.0): 1.000000
Mean uncertainty: 2.94 bits
Max uncertainty (uniform): 5.67 bits

Visualize Decoding Results¶

In [14]:

# Plot posterior probability as heatmap (first 100 time bins)
fig, ax = plt.subplots(figsize=(14, 5))

n_show = 500  # Number of time bins to show
result.plot(ax=ax, show_map=True, colorbar=True)
ax.set_xlim(0, n_show * dt)

# Overlay actual position
ax.plot(
    time_bin_centers[:n_show],
    actual_positions[:n_show, 0],
    color=COLORS["cyan"],
    linewidth=2,
    linestyle="--",
    label="Actual position",
)

ax.set_xlabel("Time (s)")
ax.set_ylabel("Position (cm)")
ax.set_title(
    "Decoded Posterior with MAP Estimate (white) and Actual Position (cyan)",
    fontweight="bold",
)
ax.legend(loc="upper right")
plt.tight_layout()
plt.show()

# Plot posterior probability as heatmap (first 100 time bins) fig, ax = plt.subplots(figsize=(14, 5)) n_show = 500 # Number of time bins to show result.plot(ax=ax, show_map=True, colorbar=True) ax.set_xlim(0, n_show * dt) # Overlay actual position ax.plot( time_bin_centers[:n_show], actual_positions[:n_show, 0], color=COLORS["cyan"], linewidth=2, linestyle="--", label="Actual position", ) ax.set_xlabel("Time (s)") ax.set_ylabel("Position (cm)") ax.set_title( "Decoded Posterior with MAP Estimate (white) and Actual Position (cyan)", fontweight="bold", ) ax.legend(loc="upper right") plt.tight_layout() plt.show()

/tmp/claude-501/ipykernel_99754/1944021853.py:26: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

In [15]:

# Compare decoded vs actual position
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Time series comparison
ax = axes[0]
ax.plot(
    time_bin_centers[:n_show],
    actual_positions[:n_show, 0],
    color=COLORS["blue"],
    linewidth=1,
    alpha=0.7,
    label="Actual",
)
ax.plot(
    time_bin_centers[:n_show],
    result.map_position[:n_show, 0],
    color=COLORS["orange"],
    linewidth=1,
    alpha=0.7,
    label="Decoded (MAP)",
)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Position (cm)")
ax.set_title("Decoded vs Actual Position Over Time", fontweight="bold")
ax.legend()

# Scatter plot
ax = axes[1]
ax.scatter(
    actual_positions[:, 0], result.map_position[:, 0], alpha=0.3, s=3, c=COLORS["blue"]
)
ax.plot(
    [0, track_length], [0, track_length], "k--", linewidth=2, label="Perfect decoding"
)
ax.set_xlabel("Actual position (cm)")
ax.set_ylabel("Decoded position (cm)")
ax.set_title("Decoded vs Actual Position", fontweight="bold")
ax.legend()
ax.set_aspect("equal")

plt.tight_layout()
plt.show()

# Compare decoded vs actual position fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # Time series comparison ax = axes[0] ax.plot( time_bin_centers[:n_show], actual_positions[:n_show, 0], color=COLORS["blue"], linewidth=1, alpha=0.7, label="Actual", ) ax.plot( time_bin_centers[:n_show], result.map_position[:n_show, 0], color=COLORS["orange"], linewidth=1, alpha=0.7, label="Decoded (MAP)", ) ax.set_xlabel("Time (s)") ax.set_ylabel("Position (cm)") ax.set_title("Decoded vs Actual Position Over Time", fontweight="bold") ax.legend() # Scatter plot ax = axes[1] ax.scatter( actual_positions[:, 0], result.map_position[:, 0], alpha=0.3, s=3, c=COLORS["blue"] ) ax.plot( [0, track_length], [0, track_length], "k--", linewidth=2, label="Perfect decoding" ) ax.set_xlabel("Actual position (cm)") ax.set_ylabel("Decoded position (cm)") ax.set_title("Decoded vs Actual Position", fontweight="bold") ax.legend() ax.set_aspect("equal") plt.tight_layout() plt.show()

/tmp/claude-501/ipykernel_99754/2816290781.py:42: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

---

Part 3: Evaluate Decoding Accuracy¶

Let's quantify how well the decoder performs using error metrics:

In [16]:

# Compute decoding error
errors = decoding_error(result.map_position, actual_positions)
median_err = median_decoding_error(result.map_position, actual_positions)

# Also compute error for mean position estimate
mean_errors = decoding_error(result.mean_position, actual_positions)
median_mean_err = median_decoding_error(result.mean_position, actual_positions)

print("Decoding Error Summary:")
print("\nMAP estimate:")
print(f"  Median error: {median_err:.2f} cm")
print(f"  Mean error: {np.nanmean(errors):.2f} cm")
print(f"  Std error: {np.nanstd(errors):.2f} cm")

print("\nMean position estimate:")
print(f"  Median error: {median_mean_err:.2f} cm")
print(f"  Mean error: {np.nanmean(mean_errors):.2f} cm")
print(f"  Std error: {np.nanstd(mean_errors):.2f} cm")

# Compute decoding error errors = decoding_error(result.map_position, actual_positions) median_err = median_decoding_error(result.map_position, actual_positions) # Also compute error for mean position estimate mean_errors = decoding_error(result.mean_position, actual_positions) median_mean_err = median_decoding_error(result.mean_position, actual_positions) print("Decoding Error Summary:") print("\nMAP estimate:") print(f" Median error: {median_err:.2f} cm") print(f" Mean error: {np.nanmean(errors):.2f} cm") print(f" Std error: {np.nanstd(errors):.2f} cm") print("\nMean position estimate:") print(f" Median error: {median_mean_err:.2f} cm") print(f" Mean error: {np.nanmean(mean_errors):.2f} cm") print(f" Std error: {np.nanstd(mean_errors):.2f} cm")

Decoding Error Summary:

MAP estimate:
  Median error: 2.47 cm
  Mean error: 3.06 cm
  Std error: 2.77 cm

Mean position estimate:
  Median error: 2.24 cm
  Mean error: 2.64 cm
  Std error: 2.17 cm

In [17]:

# Plot error distribution
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Histogram of errors
ax = axes[0]
ax.hist(errors, bins=50, color=COLORS["blue"], alpha=0.7, edgecolor="white")
ax.axvline(
    median_err,
    color=COLORS["red"],
    linewidth=2,
    linestyle="--",
    label=f"Median = {median_err:.2f} cm",
)
ax.set_xlabel("Decoding error (cm)")
ax.set_ylabel("Count")
ax.set_title("Error Distribution (MAP estimate)", fontweight="bold")
ax.legend()

# Error vs uncertainty
ax = axes[1]
ax.scatter(result.posterior_entropy, errors, alpha=0.3, s=5, c=COLORS["blue"])
ax.set_xlabel("Uncertainty (bits)")
ax.set_ylabel("Decoding error (cm)")
ax.set_title("Error vs Posterior Uncertainty", fontweight="bold")

plt.tight_layout()
plt.show()

# Compute correlation
corr = decoding_correlation(result.map_position, actual_positions)
print(f"\nDecoding correlation: {corr:.3f}")

# Plot error distribution fig, axes = plt.subplots(1, 2, figsize=(12, 4)) # Histogram of errors ax = axes[0] ax.hist(errors, bins=50, color=COLORS["blue"], alpha=0.7, edgecolor="white") ax.axvline( median_err, color=COLORS["red"], linewidth=2, linestyle="--", label=f"Median = {median_err:.2f} cm", ) ax.set_xlabel("Decoding error (cm)") ax.set_ylabel("Count") ax.set_title("Error Distribution (MAP estimate)", fontweight="bold") ax.legend() # Error vs uncertainty ax = axes[1] ax.scatter(result.posterior_entropy, errors, alpha=0.3, s=5, c=COLORS["blue"]) ax.set_xlabel("Uncertainty (bits)") ax.set_ylabel("Decoding error (cm)") ax.set_title("Error vs Posterior Uncertainty", fontweight="bold") plt.tight_layout() plt.show() # Compute correlation corr = decoding_correlation(result.map_position, actual_positions) print(f"\nDecoding correlation: {corr:.3f}")

Decoding correlation: 0.993

/tmp/claude-501/ipykernel_99754/2870634866.py:27: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

Confusion Matrix¶

A confusion matrix shows which positions are confused with each other:

In [18]:

# Compute confusion matrix
actual_bins = env.bin_at(actual_positions)
cm = confusion_matrix(env, result.posterior, actual_bins, method="map")

# Plot
fig, ax = plt.subplots(figsize=(8, 7))
im = ax.imshow(cm, cmap="Blues", aspect="auto")
ax.set_xlabel("Decoded bin")
ax.set_ylabel("Actual bin")
ax.set_title("Confusion Matrix", fontweight="bold")
plt.colorbar(im, label="Count")
plt.tight_layout()
plt.show()

print(f"Diagonal (correct) proportion: {np.diag(cm).sum() / cm.sum():.1%}")

# Compute confusion matrix actual_bins = env.bin_at(actual_positions) cm = confusion_matrix(env, result.posterior, actual_bins, method="map") # Plot fig, ax = plt.subplots(figsize=(8, 7)) im = ax.imshow(cm, cmap="Blues", aspect="auto") ax.set_xlabel("Decoded bin") ax.set_ylabel("Actual bin") ax.set_title("Confusion Matrix", fontweight="bold") plt.colorbar(im, label="Count") plt.tight_layout() plt.show() print(f"Diagonal (correct) proportion: {np.diag(cm).sum() / cm.sum():.1%}")

Diagonal (correct) proportion: 24.2%

/tmp/claude-501/ipykernel_99754/1378071997.py:13: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

---

Part 4: Trajectory Analysis (for Replay Detection)¶

For replay detection, we often want to detect whether decoded positions follow a coherent trajectory. Let's analyze a short segment that might resemble replay.

In [19]:

# Select a short segment (simulating a replay event)
start_idx = 100
end_idx = 150  # 50 time bins = 5 seconds

segment_posterior = result.posterior[start_idx:end_idx]
segment_times = time_bin_centers[start_idx:end_idx]
segment_actual = actual_positions[start_idx:end_idx, 0]

print(f"Segment: {segment_times[0]:.1f}s to {segment_times[-1]:.1f}s")
print(f"Duration: {segment_times[-1] - segment_times[0]:.1f}s")
print(f"Number of time bins: {len(segment_times)}")

# Select a short segment (simulating a replay event) start_idx = 100 end_idx = 150 # 50 time bins = 5 seconds segment_posterior = result.posterior[start_idx:end_idx] segment_times = time_bin_centers[start_idx:end_idx] segment_actual = actual_positions[start_idx:end_idx, 0] print(f"Segment: {segment_times[0]:.1f}s to {segment_times[-1]:.1f}s") print(f"Duration: {segment_times[-1] - segment_times[0]:.1f}s") print(f"Number of time bins: {len(segment_times)}")

Segment: 10.1s to 15.0s
Duration: 4.9s
Number of time bins: 50

Isotonic Regression¶

Fit a monotonic (increasing or decreasing) trajectory to the posterior:

In [20]:

# Fit isotonic trajectory
iso_result = fit_isotonic_trajectory(
    segment_posterior,
    segment_times,
    method="expected",  # Use posterior mean
    increasing=None,  # Auto-detect direction
)

print("Isotonic Regression Results:")
print(f"  Direction: {iso_result.direction}")
print(f"  R-squared: {iso_result.r_squared:.4f}")
print(
    f"  Fitted positions range: [{iso_result.fitted_positions.min():.1f}, {iso_result.fitted_positions.max():.1f}] bins"
)

# Fit isotonic trajectory iso_result = fit_isotonic_trajectory( segment_posterior, segment_times, method="expected", # Use posterior mean increasing=None, # Auto-detect direction ) print("Isotonic Regression Results:") print(f" Direction: {iso_result.direction}") print(f" R-squared: {iso_result.r_squared:.4f}") print( f" Fitted positions range: [{iso_result.fitted_positions.min():.1f}, {iso_result.fitted_positions.max():.1f}] bins" )

Isotonic Regression Results:
  Direction: increasing
  R-squared: 0.2023
  Fitted positions range: [1.6, 2.9] bins

In [21]:

# Visualize isotonic fit
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Posterior heatmap with isotonic fit
ax = axes[0]
extent = [segment_times[0], segment_times[-1], 0, env.n_bins]
ax.imshow(
    segment_posterior.T, aspect="auto", origin="lower", cmap="viridis", extent=extent
)
ax.plot(
    segment_times,
    iso_result.fitted_positions,
    color="white",
    linewidth=2,
    label="Isotonic fit",
)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Position (bin index)")
ax.set_title(f"Isotonic Fit (R² = {iso_result.r_squared:.4f})", fontweight="bold")
ax.legend()

# Residuals
ax = axes[1]
ax.bar(
    range(len(iso_result.residuals)),
    iso_result.residuals,
    color=COLORS["blue"],
    alpha=0.7,
)
ax.axhline(0, color="black", linestyle="-", linewidth=1)
ax.set_xlabel("Time bin")
ax.set_ylabel("Residual (bins)")
ax.set_title("Isotonic Regression Residuals", fontweight="bold")

plt.tight_layout()
plt.show()

# Visualize isotonic fit fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # Posterior heatmap with isotonic fit ax = axes[0] extent = [segment_times[0], segment_times[-1], 0, env.n_bins] ax.imshow( segment_posterior.T, aspect="auto", origin="lower", cmap="viridis", extent=extent ) ax.plot( segment_times, iso_result.fitted_positions, color="white", linewidth=2, label="Isotonic fit", ) ax.set_xlabel("Time (s)") ax.set_ylabel("Position (bin index)") ax.set_title(f"Isotonic Fit (R² = {iso_result.r_squared:.4f})", fontweight="bold") ax.legend() # Residuals ax = axes[1] ax.bar( range(len(iso_result.residuals)), iso_result.residuals, color=COLORS["blue"], alpha=0.7, ) ax.axhline(0, color="black", linestyle="-", linewidth=1) ax.set_xlabel("Time bin") ax.set_ylabel("Residual (bins)") ax.set_title("Isotonic Regression Residuals", fontweight="bold") plt.tight_layout() plt.show()

/tmp/claude-501/ipykernel_99754/3992581432.py:36: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

Linear Regression with Uncertainty¶

Fit a linear trajectory using Monte Carlo sampling to account for posterior uncertainty:

In [22]:

# Fit linear trajectory with sampling
lin_result = fit_linear_trajectory(
    env,
    segment_posterior,
    segment_times,
    method="sample",  # Monte Carlo sampling
    n_samples=1000,
    rng=42,
)

print("Linear Regression Results (sampling method):")
print(f"  Slope: {lin_result.slope:.2f} bins/s")
print(f"  Slope std: {lin_result.slope_std:.2f} bins/s")
print(f"  Intercept: {lin_result.intercept:.2f} bins")
print(f"  R-squared: {lin_result.r_squared:.4f}")

# Convert slope to cm/s
slope_cm_per_s = lin_result.slope * bin_size
print(f"\n  Speed: {slope_cm_per_s:.1f} cm/s")

# Fit linear trajectory with sampling lin_result = fit_linear_trajectory( env, segment_posterior, segment_times, method="sample", # Monte Carlo sampling n_samples=1000, rng=42, ) print("Linear Regression Results (sampling method):") print(f" Slope: {lin_result.slope:.2f} bins/s") print(f" Slope std: {lin_result.slope_std:.2f} bins/s") print(f" Intercept: {lin_result.intercept:.2f} bins") print(f" R-squared: {lin_result.r_squared:.4f}") # Convert slope to cm/s slope_cm_per_s = lin_result.slope * bin_size print(f"\n Speed: {slope_cm_per_s:.1f} cm/s")

Linear Regression Results (sampling method):
  Slope: 0.17 bins/s
  Slope std: 0.19 bins/s
  Intercept: 0.08 bins
  R-squared: 0.0774

  Speed: 0.3 cm/s

---

Part 5: Shuffle-Based Significance Testing¶

To determine if a decoded sequence is significant, we compare it to a null distribution generated by shuffling. This is essential for replay detection.

In [23]:

# Extract the segment of spike counts for shuffling
segment_spikes = spike_counts[start_idx:end_idx]

# Use R² from isotonic fit as our sequence score
observed_score = iso_result.r_squared

# Generate null distribution by shuffling time bins
n_shuffles = 500  # Use fewer for demo (typically 1000+)
null_scores = []

print(f"Running {n_shuffles} shuffles...")
for shuffled_spikes in shuffle_time_bins(segment_spikes, n_shuffles=n_shuffles, rng=42):
    # Decode shuffled spikes
    shuffled_result = decode_position(env, shuffled_spikes, encoding_models, dt)
    # Fit isotonic trajectory
    shuffled_fit = fit_isotonic_trajectory(
        shuffled_result.posterior, segment_times, method="expected"
    )
    null_scores.append(shuffled_fit.r_squared)

null_scores = np.array(null_scores)
print("Done!")

# Extract the segment of spike counts for shuffling segment_spikes = spike_counts[start_idx:end_idx] # Use R² from isotonic fit as our sequence score observed_score = iso_result.r_squared # Generate null distribution by shuffling time bins n_shuffles = 500 # Use fewer for demo (typically 1000+) null_scores = [] print(f"Running {n_shuffles} shuffles...") for shuffled_spikes in shuffle_time_bins(segment_spikes, n_shuffles=n_shuffles, rng=42): # Decode shuffled spikes shuffled_result = decode_position(env, shuffled_spikes, encoding_models, dt) # Fit isotonic trajectory shuffled_fit = fit_isotonic_trajectory( shuffled_result.posterior, segment_times, method="expected" ) null_scores.append(shuffled_fit.r_squared) null_scores = np.array(null_scores) print("Done!")

Running 500 shuffles...

Done!

In [24]:

# Compute p-value
p_value = compute_shuffle_pvalue(observed_score, null_scores, tail="greater")

print("Significance Test Results:")
print(f"  Observed R²: {observed_score:.4f}")
print(f"  Null mean: {null_scores.mean():.4f}")
print(f"  Null std: {null_scores.std():.4f}")
print(f"  p-value: {p_value:.4f}")
print(f"  Significant (p < 0.05): {'Yes' if p_value < 0.05 else 'No'}")

# Compute p-value p_value = compute_shuffle_pvalue(observed_score, null_scores, tail="greater") print("Significance Test Results:") print(f" Observed R²: {observed_score:.4f}") print(f" Null mean: {null_scores.mean():.4f}") print(f" Null std: {null_scores.std():.4f}") print(f" p-value: {p_value:.4f}") print(f" Significant (p < 0.05): {'Yes' if p_value < 0.05 else 'No'}")

Significance Test Results:
  Observed R²: 0.2023
  Null mean: 0.1239
  Null std: 0.1156
  p-value: 0.1457
  Significant (p < 0.05): No

In [25]:

# Visualize null distribution
fig, ax = plt.subplots(figsize=(10, 5))

ax.hist(
    null_scores,
    bins=50,
    color=COLORS["blue"],
    alpha=0.7,
    edgecolor="white",
    label="Null distribution",
)
ax.axvline(
    observed_score,
    color=COLORS["red"],
    linewidth=3,
    linestyle="--",
    label=f"Observed (R² = {observed_score:.4f})",
)
ax.axvline(
    np.percentile(null_scores, 95),
    color=COLORS["orange"],
    linewidth=2,
    linestyle=":",
    label="95th percentile",
)

ax.set_xlabel("Isotonic R²")
ax.set_ylabel("Count")
ax.set_title(f"Shuffle Test: p = {p_value:.4f}", fontweight="bold")
ax.legend()

plt.tight_layout()
plt.show()

# Visualize null distribution fig, ax = plt.subplots(figsize=(10, 5)) ax.hist( null_scores, bins=50, color=COLORS["blue"], alpha=0.7, edgecolor="white", label="Null distribution", ) ax.axvline( observed_score, color=COLORS["red"], linewidth=3, linestyle="--", label=f"Observed (R² = {observed_score:.4f})", ) ax.axvline( np.percentile(null_scores, 95), color=COLORS["orange"], linewidth=2, linestyle=":", label="95th percentile", ) ax.set_xlabel("Isotonic R²") ax.set_ylabel("Count") ax.set_title(f"Shuffle Test: p = {p_value:.4f}", fontweight="bold") ax.legend() plt.tight_layout() plt.show()

/tmp/claude-501/ipykernel_99754/3187274689.py:33: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

---

Part 6: Export Results¶

The DecodingResult can be exported to a pandas DataFrame for further analysis:

In [26]:

# Export to DataFrame
df = result.to_dataframe()

print(f"DataFrame shape: {df.shape}")
print(f"\nColumns: {list(df.columns)}")
print("\nFirst 5 rows:")
df.head()

# Export to DataFrame df = result.to_dataframe() print(f"DataFrame shape: {df.shape}") print(f"\nColumns: {list(df.columns)}") print("\nFirst 5 rows:") df.head()

DataFrame shape: (5999, 5)

Columns: ['time', 'map_bin', 'map_x', 'mean_x', 'posterior_entropy']

First 5 rows:

Out[26]:

	time	map_bin	map_x	mean_x	posterior_entropy
0	0.05	1	2.0	5.805973	3.048502
1	0.15	4	8.0	7.927778	3.280566
2	0.25	4	8.0	7.136182	3.068832
3	0.35	0	0.0	3.627466	2.560201
4	0.45	3	6.0	6.388784	3.102327

---

Summary¶

In this notebook, you learned:

The Golden Path (one call)¶

• decode_session(env, spike_times, times, positions, dt=...) runs the whole encode → bin → decode pipeline and returns a DecodingResult
• Reach for it first; it covers the common case in a single line

The Manual Three-Call Path (custom control)¶

• Compute place fields for all neurons with compute_spatial_rates() (access .firing_rates) — result shape (n_neurons, n_bins), ready for decode_position()
• Bin spikes with bin_spikes_in_time(spike_trains, dt, t_start, t_stop) → shape (n_time_bins, n_neurons)
• Decode with decode_position() to get the posterior distribution
• Use this when you need custom encoding models, to reuse fitted fields, or to inspect intermediates
• Access DecodingResult properties: posterior, map_position, mean_position, posterior_entropy

Error Metrics¶

• decoding_error() - Per-time-bin position error
• median_decoding_error() - Summary statistic
• decoding_correlation() - Weighted Pearson correlation
• confusion_matrix() - Spatial confusion analysis

Trajectory Analysis¶

• fit_isotonic_trajectory() - Monotonic trajectory fitting
• fit_linear_trajectory() - Linear fit with uncertainty estimation
• R² measures trajectory coherence

Significance Testing¶

• shuffle_time_bins() - Generate null distribution
• compute_shuffle_pvalue() - Monte Carlo p-value
• Essential for determining if decoded sequences are significant

Next Steps¶

• Try different shuffle methods (shuffle_cell_identity, shuffle_place_fields_circular)
• Explore the Radon transform for trajectory detection (detect_trajectory_radon)
• Apply to real neural data
• Experiment with custom priors

---

References¶

1. Zhang, K., Ginzburg, I., McNaughton, B. L., & Sejnowski, T. J. (1998). "Interpreting neuronal population activity by reconstruction: Unified framework with application to hippocampal place cells." Journal of Neurophysiology.

2. Davidson, T. J., Kloosterman, F., & Wilson, M. A. (2009). "Hippocampal replay of extended experience." Neuron.

3. Karlsson, M. P., & Frank, L. M. (2009). "Awake replay of remote experiences in the hippocampus." Nature Neuroscience.