Skip to content

Latest commit

 

History

1 Commit

Folders and files

Repository files navigation

SynapseDecode: Intracortical BCI Neural Decoding, Point-Process GLMs & ReFIT Kalman Filters

CI License: MIT Python 3.11+ Platform

SynapseDecode is a computational neuroscience and neuroengineering Python library for real-time motor neuroprosthetic decoding from multielectrode intracortical arrays (e.g., 96-channel Utah arrays in primary motor cortex M1).

It implements point-process Poisson Generalized Linear Models (GLMs), closed-form Kinematic Kalman Filtering, continuous-time Laplace point-process state estimation, and the clinical BrainGate2 Recalibrated Feedback Intention-Trained (ReFIT) closed-loop architecture.


Key Features

  • Point-Process Generalized Linear Models (GLM): Penalized Poisson MLE with Georgopoulos directional cosine tuning, modulation depth, and autoregressive post-spike refractory history filters.
  • Kinematic State-Space Kalman Decoder: Closed-form Least-Squares estimation of transition matrix $\mathbf{A}$, process noise $\mathbf{W}$, neural observation matrix $\mathbf{C}$, and observation noise $\mathbf{Q}$.
  • ReFIT Closed-Loop Algorithm: Intention recalibration (velocity rotation towards target, velocity zeroing on target hold, and position-orthogonality enforcement $C_{pos} = 0$) preventing cursor drift and target overshoot.
  • Continuous Point-Process Filter: Eden-Truccolo recursive Bayesian filtering operating on continuous-time discrete spike arrivals via Laplace approximation.
  • Center-Out 8-Target Reaching Paradigm: Flash & Hogan (1985) minimum-jerk kinematics generator with Fitts' law throughput and SPARC movement smoothness metrics.
  • Real-Time ASCII Terminal Visualization: 2D cursor reaching canvas, directional tuning curves, and clinical neuroprosthetic scorecards.

Neuroprosthetic Architecture

flowchart TD
    A["96-Channel Utah Array<br/>Primary Motor Cortex (M1)"] --> B["Spike Detection & Sorting<br/>Multi-Unit Thresholding"]
    B --> C["Time-Binning (Delta t = 20 ms)<br/>Variance-Stabilizing sqrt(N)"]
    C --> D["Point-Process Poisson GLM<br/>log lambda(t) = beta0 + b^T v(t) + gamma^T H(t)"]
    D --> E["Encoding Calibration<br/>Preferred Directions & Modulation"]
    C --> F["Kinematic Kalman Filter<br/>State: x_t = [px, py, vx, vy, 1]^T"]
    F --> G["Prior Prediction<br/>x_{t|t-1} = A x_{t-1}, P_{t|t-1}"]
    F --> H["Measurement Update<br/>K_t = P C^T (C P C^T + Q)^(-1)"]
    H --> I["Decoded Kinematics<br/>Velocity v(t) & Integrated Position p(t)"]
    I --> J["Closed-Loop Visual Feedback<br/>Patient Observes Cursor"]
    J --> K["ReFIT Intention Recalibration<br/>v_intended = ||v|| * (p_target - p) / ||p_target - p||"]
    K --> L["Re-estimate C Matrix with C_pos = 0<br/>Eliminate Drift & Overshoot"]
    L --> M["Clinical Performance Metrics<br/>Fitts' Throughput (bits/s), SPARC, RMSE"]
Loading

Mathematical & Neuroengineering Foundations

1. Point-Process Spiking & GLM Encoding Model

For a recorded motor cortical unit $k \in {1, \dots, K}$, conditional spike intensity $\lambda_k(t \mid H_t, \mathbf{v}(t))$ is modeled as a log-linear Poisson point process:

$$\log \lambda_k(t) = \beta_{0, k} + \mathbf{b}_k^T \mathbf{v}(t) + \sum_{j=1}^J \gamma_{k, j} n_k(t - j \Delta t)$$

where:

  • $\beta_{0, k}$ is the baseline log firing rate.
  • $\mathbf{b}k = [b{x, k}, b_{y, k}]^T$ defines the preferred reaching direction $\theta_k = \text{atan2}(b_{y, k}, b_{x, k})$ and velocity modulation depth $|\mathbf{b}_k|$.
  • $\gamma_{k, j}$ models post-spike refractory suppression and burst dynamics.

Given observed spike counts $y_k(t)$, the penalized log-likelihood is:

$$\mathcal{L}(\boldsymbol{\theta}_k) = \sum_{t=1}^T \left[ y_k(t) \log(\lambda_k(t) \Delta t) - \lambda_k(t) \Delta t \right] - \frac{\alpha}{2} |\boldsymbol{\theta}_k|_2^2$$


2. Kinematic State-Space Kalman Decoding

The patient's motor intent is represented by state vector $\mathbf{x}_t = [p_x, p_y, v_x, v_y, 1]^T \in \mathbb{R}^5$:

$$\mathbf{x}_{t+1} = \mathbf{A} \mathbf{x}_t + \mathbf{w}_t, \quad \mathbf{w}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{W})$$

Neural observation vector $\mathbf{y}_t \in \mathbb{R}^K$ (square-root spike counts):

$$\mathbf{y}_t = \mathbf{C} \mathbf{x}_t + \mathbf{q}_t, \quad \mathbf{q}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{Q})$$

The recursive Kalman equations predict and update state belief:

$$\hat{\mathbf{x}}_{t|t-1} = \mathbf{A} \hat{\mathbf{x}}_{t-1|t-1}$$

$$\mathbf{P}_{t|t-1} = \mathbf{A} \mathbf{P}_{t-1|t-1} \mathbf{A}^T + \mathbf{W}$$

$$\mathbf{K}_t = \mathbf{P}_{t|t-1} \mathbf{C}^T \left( \mathbf{C} \mathbf{P}_{t|t-1} \mathbf{C}^T + \mathbf{Q} \right)^{-1}$$

$$\hat{\mathbf{x}}_{t|t} = \hat{\mathbf{x}}_{t|t-1} + \mathbf{K}_t \left( \mathbf{y}_t - \mathbf{C} \hat{\mathbf{x}}_{t|t-1} \right)$$

$$\mathbf{P}_{t|t} = (\mathbf{I} - \mathbf{K}_t \mathbf{C}) \mathbf{P}_{t|t-1}$$

Physical position integration guarantees kinematic consistency:

$$\hat{\mathbf{p}}_t = \hat{\mathbf{p}}_{t-1} + \hat{\mathbf{v}}_t \Delta t$$


3. ReFIT (Recalibrated Feedback Intention-Trained) Architecture

In closed-loop BCI control, open-loop decoders degrade due to visual feedback corrections and target overshoot. ReFIT (Gilja et al. 2012) resolves this through two key clinical principles:

  1. Intention Recalibration: Rotates velocity vectors to point directly toward the target center, and clamps velocity to zero when inside target radius $R_{target}$: $$\mathbf{v}{intended}(t) = \begin{cases} \mathbf{0}, & \text{if } |\mathbf{p}{target} - \mathbf{p}(t)| \le R_{target} \ |\mathbf{v}{decoded}(t)| \frac{\mathbf{p}{target} - \mathbf{p}(t)}{|\mathbf{p}_{target} - \mathbf{p}(t)|}, & \text{otherwise} \end{cases}$$

  2. Position-Velocity Orthogonality: The observation matrix $\mathbf{C} \in \mathbb{R}^{K \times 5}$ is partitioned into $[C_{pos}, C_{vel}, C_{offset}]$. ReFIT explicitly sets: $$\mathbf{C}_{pos} = \mathbf{0}$$ preventing spurious position-dependent neural correlation from causing runaway cursor drift.


4. Fitts' Law & SPARC Smoothness Evaluation

  • Fitts' Law Index of Difficulty (ID): $$\text{ID} = \log_2\left( 1 + \frac{D}{W} \right) \quad [\text{bits}]$$ where $D$ is target distance (e.g. $10\text{ cm}$) and $W$ is target acceptance width ($2\text{ cm}$).
  • Throughput (TP): $$\text{TP} = \frac{\text{ID}}{\text{MT}} \quad [\text{bits/second}]$$
  • Spectral Arc Length (SPARC) Smoothness: $$\eta_{arc} = -\int_0^{\omega_c} \sqrt{\left(\frac{1}{\omega_c}\right)^2 + \left(\frac{d \hat{V}(\omega)}{d\omega}\right)^2} d\omega$$ High jerk and erratic hesitation reduce SPARC, while smooth bell profiles maximize it.

Installation

git clone https://github.com/ZengLLQ/synapse-decode.git
cd synapse-decode
pip install -e ".[dev]"

Quick Start (Python API)

import numpy as np
from synapse_decode import (
    CenterOutReachingTask,
    NeuralPopulation,
    InhomogeneousPoissonSpikeGenerator,
    KinematicKalmanDecoder,
    BCIMetricsEvaluator,
    AsciiBCIVisualizer,
)

# 1. Setup 8-target center-out task with 64 cortical units
dt = 0.020  # 50 Hz control rate
task = CenterOutReachingTask(radius_m=0.10, dt=dt)
population = NeuralPopulation(num_neurons=64, seed=42)
spike_gen = InhomogeneousPoissonSpikeGenerator(dt=dt)

# 2. Synthesize reaching training data
trials = task.generate_all_trials(hold_time_s=0.2)
kin_train = np.vstack([np.column_stack([tr.positions, tr.velocities]) for tr in trials])
spikes_train = np.vstack([
    np.sqrt(spike_gen.generate_spike_counts(population.encode_velocities(tr.velocities), seed=i))
    for i, tr in enumerate(trials)
])

# 3. Fit Kinematic Kalman Filter
decoder = KinematicKalmanDecoder(dt=dt, integrate_position=True)
decoder.fit(kin_train, spikes_train)

# 4. Decode test reaching movement
test_trial = trials[1]  # Target 45 deg
test_spikes = np.sqrt(spike_gen.generate_spike_counts(population.encode_velocities(test_trial.velocities), seed=999))
decoded_trajectory = decoder.decode_sequence(test_spikes)

# 5. Render ASCII reaching canvas
print(AsciiBCIVisualizer.render_center_out_trajectory(
    task.target_coords, decoded_trajectory[:, 0:2], radius_m=task.radius_m
))

CLI Guide

The package provides a synapse-decode CLI tool:

1. Run Complete Center-Out BCI Decoding Demo

synapse-decode demo --neurons 64 --trials 3

Sample output:

Initializing Intracortical BCI Simulation with 64 M1 Units...
Kalman filter trained in 1.25 ms.
Decoded 48 steps in 3.80 ms.

+---------------------------------------------------+
|  CENTER-OUT BCI CURSOR TRAJECTORY CANVAS          |
+---------------------------------------------------+
|                           **                      |
|                        *3*                        |
|                        ** ***                     |
|           4            *    **        2           |
|       *               *     *                     |
|                 *    *     **                     |
|                   **      @ *                     |
|      5       ** * ***   O                  1      |
|           6                           8           |
|                         7                         |
+---------------------------------------------------+
| Legend: [O] Origin  [1-8] Targets  [*] Path  [@] End Cursor |

+---------------------------------------------------------+
|      NEUROPROSTHETIC BCI PERFORMANCE SCORECARD          |
+---------------------------------------------------------+
| Decoder Architecture         : Kinematic Kalman Filter  |
| Velocity Correlation Y (ry)  :    0.789                 |
| Velocity RMSE                :    0.058 m/s             |
| Position RMSE                :     4.18 cm              |
| Fitts' Throughput            :     2.69 bits/s           |
| Smoothness (SPARC)           :    -2.01                  |
+---------------------------------------------------------+

2. Benchmark ReFIT Closed-Loop Algorithm

synapse-decode refit-benchmark
+------------------------------+-------------------+-------------------+
| Performance Metric           | Standard Kalman   | ReFIT Kalman      |
+------------------------------+-------------------+-------------------+
| Velocity Correlation (r_mean)|             0.879 |             0.920 |
| Velocity RMSE (m/s)          |             0.046 |             0.038 |
| Position RMSE (cm)           |              2.08 |              1.82 |
| Throughput (bits/s)          |              2.58 |              2.85 |
| Smoothness (SPARC)           |             -4.89 |             -3.42 |
+------------------------------+-------------------+-------------------+

3. Fit Poisson Point-Process GLM to Spike Data

synapse-decode glm-fit

Testing & Continuous Integration

SynapseDecode includes 15 automated unit tests covering Poisson GLM parameter recovery, Kalman matrices, point-process Laplace filtering, ReFIT intention rotation, and Fitts' law metrics:

pytest -v

CI runs automatically on GitHub Actions across:

  • Operating Systems: Ubuntu (ubuntu-latest), macOS (macos-latest), Windows (windows-latest)
  • Python Versions: 3.11, 3.12

License

MIT License. Copyright (c) 2026 ZengLLQ.

Releases

Packages

Contributors

Languages