Spike Train to Spatial Field Conversion

This guide covers the primitives for converting spike trains into occupancy-normalized spatial fields (firing rate maps), a foundational operation in place field analysis and spatial neuroscience.

Overview

The neurospatial encoding API provides one canonical entry point for spike-to-field conversion:

• compute_spatial_rate() - Converts spike trains to occupancy-normalized spatial firing-rate maps and returns a SpatialRateResult

The function follows the neurospatial API convention: environment comes first in the parameter order.

Why Occupancy Normalization Matters

Firing rate fields must be normalized by the time spent in each spatial bin (occupancy) to produce meaningful estimates of spatial tuning:

\[ \text{firing rate}_i = \frac{\text{spike count}_i}{\text{occupancy}_i} \]

Without occupancy normalization, bins visited more frequently would appear to have higher firing rates simply due to increased sampling, not actual spatial preference. This normalization is standard practice in place field analysis and is used universally in neuroscience literature (O'Keefe & Dostrovsky, 1971; Muller et al., 1987).

What is Occupancy?

Occupancy is the time-weighted measure of how long an animal spent in each spatial bin. Neurospatial computes occupancy by:

1. Binning the position trajectory into spatial bins
2. Computing the duration of each time interval (difference between consecutive timestamps)
3. Allocating these durations to the bins visited

By default, env.occupancy() returns occupancy in seconds (time-weighted). You can optionally get interval counts (unweighted) using return_seconds=False, but for firing rate computation, always use the default return_seconds=True.

Basic Usage

Converting Spike Trains to Firing Rate Fields

The compute_spatial_rate() function converts spike train data into spatial firing rate maps:

import numpy as np
from neurospatial import Environment
from neurospatial.encoding import compute_spatial_rate

# Create environment from trajectory data
env = Environment.from_samples(positions, bin_size=2.5)

# Spike times (seconds)
spike_times = np.array([1.2, 3.5, 5.1, 7.8, ...])

# Trajectory timestamps and positions
times = np.array([0.0, 0.033, 0.067, ...])  # 30 Hz sampling
positions = np.array([[x1, y1], [x2, y2], ...])  # Shape (n_timepoints, 2)

# Compute firing rate field
result = compute_spatial_rate(
    env,
    spike_times,
    times,
    positions,
    smoothing_method="binned",
    min_occupancy=0.5,  # Exclude bins with < 0.5 seconds occupancy
)
firing_rate = result.firing_rate

# Result: firing_rate.shape = (env.n_bins,)
# Units: spikes/second
# Bins with insufficient occupancy are set to NaN

Key points:

• Parameter order: environment comes first, followed by spike data and trajectory data
• Returns firing rate in spikes/second
• Bins with occupancy less than min_occupancy are set to NaN
• Default min_occupancy=0.0 includes all bins

One-Liner with Smoothing: compute_spatial_rate()

For typical place field analysis workflows, use compute_spatial_rate() with a smoothing method:

from neurospatial.encoding import compute_spatial_rate

# Compute place field with Gaussian smoothing
result = compute_spatial_rate(
    env,
    spike_times,
    times,
    positions,
    min_occupancy=0.5,
    smoothing_method="gaussian_kde",
    bandwidth=5.0,  # Gaussian kernel bandwidth (cm)
)
place_field = result.firing_rate

Use smoothing_method="binned" when you want the unsmoothed occupancy-normalized rate map.

Parameter Order: Environment First

Important: All neurospatial functions follow a consistent API pattern where the environment comes first:

# Correct parameter order
compute_spatial_rate(env, spike_times, times, positions, ...)

# This matches the existing neurospatial API:
env.occupancy(times, positions)
env.smooth(field, bandwidth)
map_points_to_bins(points, env)
distance_field(env.connectivity, sources)

This design decision ensures:

• Consistency across the entire neurospatial API
• Readability - the environment context is always explicit and comes first
• Composability - functions can be easily chained and combined

Min Occupancy Threshold: Best Practices

The min_occupancy parameter controls reliability filtering:

# Include all bins (default)
firing_rate = compute_spatial_rate(
    env, spike_times, times, positions, min_occupancy=0.0
).firing_rate

# Exclude bins with < 0.5 seconds occupancy (typical for place fields)
firing_rate = compute_spatial_rate(
    env, spike_times, times, positions, min_occupancy=0.5
).firing_rate

# More conservative (1 second threshold)
firing_rate = compute_spatial_rate(
    env, spike_times, times, positions, min_occupancy=1.0
).firing_rate

Recommendations:

• Default (0.0): Use when you need complete spatial coverage or will apply additional filtering
• 0.5 seconds: Standard for place field analysis - excludes bins with too few samples for reliable rate estimates
• 1.0+ seconds: Conservative threshold for high-confidence place field detection

Bins with occupancy below the threshold are set to NaN (not zero), which:

• Allows you to distinguish between "silent" bins (visited but no spikes) and "unsampled" bins
• Prevents these bins from contaminating spatial statistics
• Works correctly with env.smooth() which handles NaN values appropriately

Why Not Use Default Filtering?

The default min_occupancy=0.0 (no filtering) is intentional:

1. Flexibility: Users can inspect all bins and make informed filtering decisions
2. Transparency: All data is preserved, no hidden thresholding
3. Composability: You can apply custom filters after computing firing rates

For typical place field analysis, explicitly set min_occupancy=0.5 to match standard neuroscience practices.

Edge Cases and Validation

Empty Spike Trains

Empty spike arrays (no spikes) produce fields of zeros:

spike_times = np.array([])  # No spikes

firing_rate = compute_spatial_rate(env, spike_times, times, positions).firing_rate
# Result: all zeros (or NaN where occupancy < min_occupancy)

Out-of-Bounds Spikes

Spikes occurring outside the trajectory time range are filtered with a warning:

spike_times = np.array([1.0, 100.0, 200.0])  # Some spikes beyond times.max()
times = np.array([0.0, 0.033, ..., 50.0])  # Max time = 50.0

firing_rate = compute_spatial_rate(env, spike_times, times, positions).firing_rate
# UserWarning: "X spikes fall outside trajectory time range [0.0, 50.0] and were excluded"

Similarly, spikes whose interpolated positions fall outside the environment bounds are filtered with a warning.

1D Trajectories

The function handles 1D trajectories (linear tracks, circular tracks) correctly:

# 1D positions - either shape (n_timepoints, 1) or (n_timepoints,)
positions_1d = np.array([0.0, 2.5, 5.0, 7.5, ...])  # Shape (n_timepoints,)

firing_rate = compute_spatial_rate(env, spike_times, times, positions_1d).firing_rate
# Automatically handles 1D interpolation

Input Validation

The function validates inputs and raises informative errors:

# Mismatched lengths
times = np.array([0.0, 0.033, 0.067])
positions = np.array([[0, 0], [1, 1]])  # Different length!
# Raises ValueError: "times and positions must have the same length"

# Negative min_occupancy
compute_spatial_rate(env, spike_times, times, positions, min_occupancy=-1.0)
# Raises ValueError: "min_occupancy must be non-negative (got -1.0)"

NaN Handling in Smoothing

When using compute_spatial_rate() with smoothing, NaN values from low-occupancy bins are handled by the smoothing implementation:

1. NaN bins are temporarily filled with 0
2. Gaussian smoothing is applied
3. NaN values are restored to their original locations

Important caveat: This approach can artificially reduce firing rates near unvisited regions due to the zero-filling. For scientific applications requiring high precision near boundaries, use the two-step workflow:

# Manual workflow with custom NaN handling
firing_rate = compute_spatial_rate(
    env,
    spike_times,
    times,
    positions,
    smoothing_method="binned",
    min_occupancy=0.5,
).firing_rate

# Option 1: Smooth only valid bins (NaN stays NaN)
smoothed = env.smooth(firing_rate, bandwidth=5.0)

# Option 2: Inpaint NaN before smoothing (advanced)
from scipy.ndimage import distance_transform_edt
# ... custom NaN handling ...

For most use cases, the default NaN handling in compute_spatial_rate() is appropriate and matches common neuroscience practice.

Visualizing Results

Typical visualization workflow:

import matplotlib.pyplot as plt

# Compute occupancy and firing rate
occupancy = env.occupancy(times, positions, return_seconds=True)
firing_rate = compute_spatial_rate(
    env,
    spike_times,
    times,
    positions,
    smoothing_method="binned",
    min_occupancy=0.5,
).firing_rate
place_field = compute_spatial_rate(
    env,
    spike_times,
    times,
    positions,
    min_occupancy=0.5,
    smoothing_method="gaussian_kde",
    bandwidth=5.0,
).firing_rate

# Visualize using environment plotting
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

env.plot(occupancy, ax=axes[0], cmap='viridis')
axes[0].set_title('Occupancy (seconds)')

env.plot(firing_rate, ax=axes[1], cmap='hot')
axes[1].set_title('Firing Rate (raw)')

env.plot(place_field, ax=axes[2], cmap='hot')
axes[2].set_title('Place Field (smoothed)')

plt.tight_layout()
plt.show()

Tips:

• Use cmap='viridis' for occupancy (sequential)
• Use cmap='hot' or cmap='inferno' for firing rates (perceptually uniform, peaks visible)
• The env.plot() method automatically handles NaN values (shown as white/transparent)

Complete Example

Here's a complete workflow from trajectory data to smoothed place field:

import numpy as np
from neurospatial import Environment
from neurospatial.encoding import compute_spatial_rate

# 1. Load or generate trajectory data
times = np.linspace(0, 60, 1800)  # 1 minute at 30 Hz
positions = np.column_stack([
    10 * np.sin(0.5 * times),  # X coordinate
    10 * np.cos(0.5 * times)   # Y coordinate
])  # Circular trajectory

# 2. Generate spike train (simulated place cell)
preferred_location = np.array([5.0, 5.0])
distances = np.linalg.norm(positions - preferred_location, axis=1)
firing_prob = 0.5 * np.exp(-distances**2 / (2 * 3**2))  # Gaussian tuning
spike_times = times[np.random.rand(len(times)) < firing_prob * 0.033]

# 3. Create environment
env = Environment.from_samples(positions, bin_size=2.5)
env.units = "cm"

# 4. Compute place field
result = compute_spatial_rate(
    env,
    spike_times,
    times,
    positions,
    min_occupancy=0.5,
    smoothing_method="gaussian_kde",
    bandwidth=5.0,
)
place_field = result.firing_rate

# 5. Visualize
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 6))
env.plot(place_field, ax=ax, cmap='hot')
ax.set_title('Place Field')
plt.show()

print(f"Environment: {env.n_bins} bins")
print(f"Peak firing rate: {np.nanmax(place_field):.2f} Hz")
print(f"Spatial information: {np.sum(place_field > 0.5 * np.nanmax(place_field))} bins > 50% peak")

API Reference

For complete parameter descriptions and examples, see:

• compute_spatial_rate() - Place field computation
• bin_spike_train() - Spike binning helper
• Environment.occupancy() - Occupancy computation
• Environment.smooth() - Spatial field smoothing

Related Topics

• Spatial Analysis - Broader spatial analysis workflows
• RL Primitives - Reward field generation for reinforcement learning
• Workflows - Complete analysis pipelines