Brain machine interfaces modeling strategies for neural signal processing
Download
1 / 68

Brain Machine Interfaces: Modeling Strategies for Neural Signal Processing - PowerPoint PPT Presentation


  • 86 Views
  • Uploaded on

Brain Machine Interfaces: Modeling Strategies for Neural Signal Processing. Jose C. Principe, Ph.D. Distinguished Professor ECE, BME Computational NeuroEngineering Laboratory Electrical and Computer Engineering Department University of Florida www.cnel.ufl.edu [email protected]

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about ' Brain Machine Interfaces: Modeling Strategies for Neural Signal Processing' - arden-carney


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Brain machine interfaces modeling strategies for neural signal processing

Brain Machine Interfaces: Modeling Strategies for Neural Signal Processing

Jose C. Principe, Ph.D.

Distinguished Professor ECE, BME

Computational NeuroEngineering Laboratory

Electrical and Computer Engineering Department

University of Florida

www.cnel.ufl.edu

[email protected]


Brain machine interfaces bmi
Brain Machine Interfaces (BMI) Signal Processing

  • A man made device that either substitutes a sensory input to the brain, repairs functional communication between brain regions or translates intention of movement.


Types of bmis
Types of BMIs Signal Processing

  • Sensory (Input BMI): Providing sensory input to form percepts when natural systems are damaged.

    • Ex: Visual, Auditory Prosthesis

  • Motor (Output BMI): Converting motor intent to a command output (physical device, damaged limbs)

    • Ex: Prosthetic Arm Control

  • Cognitive BMI: Interpret internal neuronal state to deliever feedback to the neural population.

    • Ex: Epilepsy, DBS Prosthesis

      Computational Neuroscience and Technology developments are playing a larger role in the development of each of these areas.


General Architecture Signal Processing

BCI (BMI)bypasses the brain’s normal pathways of peripheral nerves (and muscles)

J.R. Wolpaw et al. 2002


The fundamental concept
The Fundamental Concept Signal Processing

BRAIN MACHINE

INTENT

ACTION

Decoding

PERCEPT

STIMULUS

Coding

Neural Interface Physical Interface

Need to understand how brain processes information.


Levels of abstraction for neurotechnology
Levels of Abstraction for Neurotechnology Signal Processing

  • Brain is an extremely complex system

  • 1012 neurons

  • 1015 synapses

  • Specific interconnectivity


Tapping into the nervous system
Tapping into the Nervous System Signal Processing

  • The choice and availability of brain signals and recording methods can greatly influence the ultimate performance of the BMI.

    The level of BMI performance may be attributed to selection of electrode technology, choice of model, and methods for extracting rate, frequency, or timing codes.


Coarse(mm) Signal Processing

http://ida.first.fhg.de/projects/bci/bbci_official/



Spatial resolution of recordings
Spatial Resolution of Recordings Signal Processing

Moran


Florida multiscale signal acquisition
Florida Multiscale Signal Acquisition Signal Processing

Develop a experimental paradigm with a nested hierarchy for studying neural population dynamics.

Least Invasive

EEG

NRG IRB Approval for Human Studies

ECoG

NRG IACUC Approval for Animal Studies

Microelectrodes

Highest Resolution


Common bmi bci methods
Common BMI-BCI Methods Signal Processing

  • BMIs --- Invasive, work with intention of movement

    • Spike trains, field potentials, ECoG

    • Very specific, potentially better performance

  • BCIs --- Noninvasive, subjects must learn how to control their brain activity

    • EEG

    • Very small bandwidth


  • Computational neuroscience
    Computational NeuroScience Signal Processing

    • Integration of probabilistic models of information processing with the neurophysiological reality of brain anatomy, physiology and purpose.

    • Need to abstract the details of the “wetware” and ask what is the purpose of the function. Then quantify it in mathematical terms.

    • Difficult but very promising. One issue is that biological evolution is a legacy system!

    • BMI research is an example of a computational neuroscience approach.


    How to put it together
    How to put it together? Signal Processing

    • NeoCortical Brain Areas Related to Movement

    Posterior Parietal (PP) – Visual to motor transformation

    Premotor (PM) and Dorsal Premotor (PMD) -

    Planning and guidance (visual inputs)

    Primary Motor (M1) – Initiates muscle contraction


    Ensemble correlations local in time are averaged with global models
    Ensemble Correlations – Signal ProcessingLocal in Time – are Averaged with Global Models


    Computational models of neural intent
    Computational Models of Neural Intent Signal Processing

    • Two different levels of neurophysiology realism

      • Black Box models – no realism, function relation between input desired response

      • Generative Models – minimal realism, state space models using neuroscience elements


    Signal processing approaches with black box modeling
    Signal Processing Approaches with Black Box Modeling Signal Processing

    Accessing 2 types of signals (cortical activity and behavior) leads us to a general class of I/O models.

    Data for these models are rate codes obtained by binning spikes on 100 msec windows.

    • Optimal FIR Filter – linear, feedforward

    • TDNN – nonlinear, feedforward

    • Multiple FIR filters – mixture of experts

    • RMLP – nonlinear, dynamic


    Linear model wiener hopf solution
    Linear Model (Wiener-Hopf solution) Signal Processing

    Consider a set of spike counts from M neurons, and a hand position vector dC (C is the output dimension, C = 2 or 3). The spike count of each neuron is embedded by an L-tap discrete time-delay line. Then, the input vector for a linear model at a given time instance n is composed as x(n) = [x1(n), x1(n-1) … x1(n-L+1), x2(n) … xM(n-L+1)]T, xLM, where xi(n-j) denotes the spike count of neuron i at a time instance n-j.

    A linear model estimating hand position at time instance n from the embedded spike counts can be described as

    where yc is the c-coordinate of the estimated hand position by the model, wji is a weight on the connection from xi(n-j) to yc, and bc is a bias for the c-coordinate.


    Linear model wiener hopf solution1

    x Signal Processing1(n)

    yx(n)

    z-1

    z-1

    z-1

    z-1

    yy(n)

    xM(n)

    yz(n)

    Linear Model (Wiener-Hopf solution)

    In a matrix form, we can rewrite the previous equation as

    where y is a C-dimensional output vector, and W is a weight matrix of dimension (LM+1)C. Each column of W consists of [w10c, w11c, w12c…, w1L-1c, w20c, w21c…, wM0c, …, wML-1c]T.


    Linear model wiener hopf solution2
    Linear Model (Wiener-Hopf solution) Signal Processing

    For the MIMO case, the weight matrix in the Wiener filter system is estimated by

    R is the correlation matrix of neural spike inputs with the dimension of (LM)(LM),

    where rij is the LL cross-correlation matrix between neurons i and j (i ≠ j), and rii is the LL autocorrelation matrix of neuron i.

    P is the (LM)C cross-correlation matrix between the neuronal bin count and hand position, where pic is the cross-correlation vector between neuron i and the c-coordinate of hand position. The estimated weights WWiener are optimal based on the assumption that the error is drawn from white Gaussian distribution and the data are stationary.


    Linear model wiener hopf solution3
    Linear Model (Wiener-Hopf solution) Signal Processing

    The predictor WWiener minimizes the mean square error (MSE) cost function,

    Each sub-block matrix rij can be further decomposed as

    where rij() represents the correlation between neurons i and j with time lag . Assuming that the random process xi(k) is ergodic for all i, we can utilize the time average operator to estimate the correlation function. In this case, the estimate of correlation between two neurons, rij(m-k), can be obtained by


    Linear model wiener hopf solution4
    Linear Model (Wiener-Hopf solution) Signal Processing

    The cross-correlation vector pic can be decomposed and estimated in the same way, substituting xj by the desired signal cj.

    From the equations, it can be seen that rij(m-k) is equal to rji(k-m). Since these two correlation estimates are positioned at the opposite side of the diagonal entries of R, the equality leads to a symmetric R.

    The symmetric matrix R, then, can be inverted effectively by using the Cholesky factorization. This factorization reduces the computational complexity for the inverse of R from O(N3) using Gaussian elimination to O(N2) where N is the number of parameters.


    Optimal linear model
    Optimal Linear Model Signal Processing

    • Normalized LMS with weight decay is a simple starting point.

    • Four multiplies, one divide and two adds per weight update

    • Ten tap embedding with 105 neurons

    • For 1-D topology contains 1,050 parameters (3,150)

    • Alternatively, the Wiener solution


    Time delay neural network tdnn
    Time-Delay Neural Network (TDNN) Signal Processing

    • The first layer is a bank of linear filters followed by a nonlinearity.

    • The number of delays to span I second

    • y(n)= Σ wf(Σwx(n))

    • Trained with backpropagation

    • Topology contains a ten tap embedding and five hidden PEs– 5,255 weights (1-D)

    Principe, UF


    Multiple switching local models
    Multiple Switching Local Models Signal Processing

    • Multiple adaptive filters that compete to win the modeling of a signal segment.

    • Structure is trained all together with normalized LMS/weight decay

    • Needs to be adapted for input-output modeling.

    • We selected 10 FIR experts of order 10 (105 input channels)

    d(n)


    Recurrent multilayer perceptron rmlp nonlinear black box
    Recurrent Multilayer Perceptron (RMLP) – Nonlinear Signal Processing“Black Box”

    • Spatially recurrent dynamical systems

    • Memory is created by feeding back the states of the hidden PEs.

    • Feedback allows for continuous representations on multiple timescales.

    • If unfolded into a TDNN it can be shown to be a universal mapper in Rn

    • Trained with backpropagation through time


    Motor tasks performed
    Motor Tasks Performed Signal Processing

    Data

    Task 1

    • 2 Owl monkeys – Belle, Carmen

    • 2 Rhesus monkeys – Aurora, Ivy

    • 54-192 sorted cells

    • Cortices sampled: PP, M1, PMd, S1, SMA

    • Neuronal activity rate and behavior is time synchronized and downsampled to 10Hz

    Task 2


    Model building techniques
    Model Building Techniques Signal Processing

    • Train the adaptive system with neuronal firing rates (100 msec) as the input and hand position as the desired signal.

    • Training - 20,000 samples (~33 minutes of neuronal firing)

    • Freeze weights and present novel neuronal data.

    • Testing - 3,000 samples – (5 minutes of neuronal firing)


    Results belle

    Signal to error ratio (dB) Signal Processing

    Correlation Coefficient

    (average)

    (max)

    (average)

    (max)

    LMS

    0.8706

    7.5097

    0.6373

    0.9528

    Kalman

    0.8987

    8.8942

    0.6137

    0.9442

    TDNN

    1.1270

    3.6090

    0.4723

    0.8525

    Local Linear

    1.4489

    23.0830

    0.7443

    0.9748

    RNN

    1.6101

    32.3934

    0.6483

    0.9852

    Results (Belle)

    • Based on 5 minutes of test data, computed over 4 sec windows (training on 30 minutes)


    Physiologic interpretation
    Physiologic Interpretation Signal Processing

    • When the fitting error is above chance, a sensitivity analysis can be performed by computing the Jacobian of the output vector with respect to each neuronal input i

    • This calculation indicates which inputs (neurons) are most important for modulating the output/trajectory of the model.


    Computing sensitivities through the models
    Computing Sensitivities Through the Models Signal Processing

    Identify the neurons that affect the output the most.

    Feedforward Linear Eq.

    Feedforward RMLP Eqs.

    General form of Linear Sensitivity

    General form of RMLP Sensitivity


    Data analysis the effect of sensitive neurons on performance
    Data Analysis : The Effect of Sensitive Neurons on Performance

    Decay trend appears in all animals and behavioral paradigms


    Directional tuning vs sensitivity of ranked cells
    Directional Tuning vs. Sensitivity of ranked cells Performance

    Tuning

    Sensitivity

    Significance: Sensitivity analysis through trained models automatically delivers deeply tuned cells that span the space.


    Reaching movement segmentation
    Reaching Movement Segmentation Performance

    How does each cortical area contribute to the reconstruction of this movement?


    Cortical contributions belle day 2
    Cortical Contributions Belle Day 2 Performance

    Train 15 separate RMLPs with every combination of cortical input.


    Is there enough information in spike trains for modeling movement
    Is there enough information in spike trains for modeling movement?

    • Analysis is based on the time embedded model

      • Correlation with desired is based on a linear filter output for each neuron

    • Utilize a non-stationary tracking algorithm

      • Parameters are updated by LMS

    • Build a spatial filter

      • Adaptive in real time

      • Sparse structure based on regularization for enables selection


    Architecture

    z movement?-1

    z-1

    z-1

    z-1

    Architecture

    x1(n)

    w11

    y1(n)

    //

    c1

    w1L

    y2(n)

    c2

    xM(n)

    cM

    wM1

    yM(n)

    //

    wML

    Adapted by on-line LAR

    (Kim et. al., MLSP, 2004)

    Adapted by LMS


    Training algorithms

    x movement?1

    k

    . . .

    y

    r= y-X = y

    Find argmaxi |xiTr|

    i=0

    xj

    xj

    xk

    r= y-X = y-xjj

    Adjustj s.t.

    k, |xkTr|=|xiTr|

    r= y-(xjj+ xkk)

    Adjustj & k s.t.

    q, |xqTr|=|xkTr|=|xiTr|

    j

    j

    Training Algorithms

    • Tap weights for every time lag is updated by LMS

    • Then, the spatial filter coefficients are obtained by on-line version of least angle regression (LAR) (Efron et. al. 2004)



    Application to bmi data neuronal subset selection
    Application to BMI Data – Neuronal Subset Selection movement?

    Hand Trajectory

    (z)

    Early

    Part

    Neuronal Channel Index

    Late

    Part


    Generative models for bmis
    Generative Models for BMIs movement?

    • Use partial information about the physiological system, normally in the form of states.

    • They can be either applied to binned data or to spike trains directly.

    • Here we will only cover the spike train implementations.

      Difficulty of spike train Analysis:

      Spike trains are point processes, i.e. all the information is contained in the timing of events, not in the amplitude fo the signals!


    Goal movement?

    Build an adaptive signal processing framework for BMI decoding in the spike domain.

    • Features of Spike domain analysis

      • Binning window size is not a concern

      • Preserve the randomness of the neuron behavior.

      • Provide more understanding of neuron physiology (tuning) and interactions at the cell assembly level

      • Infer kinematics online

      • Deal with nonstationary

      • More computation with millisecond time resolution


    Recursive bayesian approach

    Time-series model movement?

    State

    cont. observ.

    Prediction

    Updating

    P(state|observation)

    Recursive Bayesian Approach


    Recursive bayesian approach1
    Recursive Bayesian approach movement?

    • State space representation

      First equation (system model) defines a first order Markov process.

      Second equation (observation model) defines the likelihood of the observations p(zt|xt) . The problem is completely defined by the prior distribution p(x0).

      Although the posterior distribution p(x0:t|u1:t,z1:t) constitutes the complete solution, the filtering density p(xt|u1:t, z1:t) is normally used for on-line problems.

      The general solution methodology is to integrate over the unknown variables (marginalization).


    Recursive bayesian approach2
    Recursive Bayesian approach movement?

    • There are two stages to update the filtering density:

      • Prediction (Chapman Kolmogorov)

        System model p(xt|xt-1) propagates into the future the posterior density

      • Update

        Uses Bayes rule to update the filtering density. The following equations are needed in the solution.


    Kalman filter for bmi decoding
    Kalman filter for BMI decoding movement?

    For Gaussian noises and linear prediction and observation models, there

    is an analytic solution called the Kalman Filter.

    Continuous Observation

    Linear

    Neuron tuning function

    Kinematic State

    Firing rate

    Linear

    Gaussian

    Prediction

    P(state|observation)

    Updating

    [Wu et al. 2006]


    Particle filter for bmi decoding
    Particle Filter for BMI decoding movement?

    In general the integrals need to be approximated by sums using Monte Carlo integration with a set of samples drawn from the posterior distribution of the model parameters.

    Continuous Observation

    Exponential

    Neuron tuning function

    Kinematic State

    Firing rate

    Linear

    nonGaussian

    Prediction

    P(state|observation)

    Updating

    [Brockwell et al. 2004]


    State estimation framework for bmi decoding in spike domain

    v movement?

    x

    k-1

    k-1

    x

    k

    F

    k

    (

    )

    =

    ,

    =

    (

    )

    ,

    n

    z

    H

    x

    k

    k

    k

    k

    State estimation framework for BMI decoding in spike domain

    Tuning function

    Neural Tuning function

    Kinematics

    state

    Multi-spike trains observation

    Decoding

    Kinematic dynamic model

    Key Idea: work with the probability of spike firing which is a

    continuous random variable


    Adaptive algorithm for point processes
    Adaptive algorithm for point processes movement?

    Poisson Model

    nonlinear

    Point process

    Neuron tuning function

    Kinematic State

    spike train

    Linear

    Gaussian

    Prediction

    P(state|observation)

    Updating

    [Brown et al. 2001]


    Monte carlo sequential estimation for point process
    Monte Carlo Sequential estimation for point process movement?

    nonlinear

    Point process

    Neuron tuning function

    Kinematic State

    spike train

    nonLinear

    nonGaussian

    Prediction

    P(state|observation)

    Updating

    Sequential Estimate PDF

    [Wang et al. 2006]


    Monte carlo sequential estimation framework for bmi decoding in spike domain
    Monte Carlo sequential estimation framework for BMI decoding in spike domain

    • STEP 1. Preprocessing

      1. Generate spike trains from stored spike times 10ms interval, (99.62% binary train)

      2. Synchronize all the kinetics with the spike trains.

      3. Assign the kinematic vector to reconstruct.

      X=[position velocity acceleration]’

      (more information, instantaneous state avoid error accumulation,

      less computation)


    Step 2 neural tuning analysis

    kinematics in spike domain

    Neural spike trains

    STEP 2- Neural tuning analysis

    Encoding

    (Tuning)

    • A example of a tuned neuron

    • Metric: Tuning depth:

    • how differently does a neuron fire across directions?

    • D=(max-min)/std (firing rate)

    Neuron 72: Tuning Depth 1


    Step 2 information theoretic metric of tuning
    Step 2- Information Theoretic Metric of Tuning in spike domain

    kinematics direction angle

    Information

    neural spikes


    Step 2 i nformation theoretic tuning depths for 3 kinds of kinematics log axis
    Step 2- I in spike domainnformation theoretic Tuning depths for 3 kinds of kinematics (log axis)


    Step 2 tuning function estimation

    velocity in spike domain

    spikes

    Linear filter

    nonlinear f

    Poisson model

    Step 2- Tuning Function Estimation

    • Neural firing Model

    • Assumption :

      generation of the spikes depends only on the kinematic vector we choose.


    Step 2 linear filter estimation
    Step 2- Linear Filter Estimation in spike domain

    • Spike Triggered Average (STA)

    • Geometry interpretation

    2nd Principal component

    1st Principal component


    Step 2 nonlinear f estimation
    Step 2- Nonlinear in spike domainf estimation


    Step 2 diversity of neural nonlinear properties
    Step 2- Diversity of neural nonlinear properties in spike domain

    Ref: Paradoxical cold

    [Hensel et al. 1959]



    Step 3 sequential estimation algorithm for point process filtering
    Step 3: Sequential Estimation Algorithm for Point Process Filtering

    • Consider the neuron as an inhomogenous Poisson point process

    • Observing N(t) spikes in an interval DT, the posterior of the spike model is

    • The probability of observing an event in Dt is

    • And the one step prediction density (Chapman-Kolmogorov)

    • The posterior of the state vector, given an observation DN


    Step 3 sequential estimation algorithm for point process filtering1
    Step 3: Sequential Estimation Algorithm for Point Process Filtering

    • Monte Carlo Methods are used to estimate the integral. Let represent a random measure on the posterior density, and represent the proposed density by

    • The posterior density can then be approximated by

    • Generating samples from using the principle of Importance sampling

    • By MLE we can find the maximum or use direct estimation with kernels of mean and variance




    Step 3: Information Estimated Delays Filtering

    Figure 3-14 Mutual information as function of time delay for 5 neurons.

    For 185 neurons, average delay is 220.108 ms


    Step 4 monte carlo sequential kinematics estimation

    Neural Tuning function Filtering

    Kinematic State

    spike trains

    Prediction

    Updating

    P(state|observation)

    Step 4: Monte Carlo sequential kinematics estimation

    NonGaussian



    Results comparison Filtering

    Table 3-2 Correlation Coefficients between the Desired Kinematics and the Reconstructions

    Table 3-3 Correlation Coefficient Evaluated by the Sliding Window

    [Sanchez, 2004]


    Conclusion
    Conclusion Filtering

    • Our results and those from other laboratories show it is possible to extract intent of movement for trajectories from multielectrode array data.

    • The current results are very promising, but the setups have limited difficulty, and the performance seems to have reached a ceiling at an uncomfortable CC < 0.9

    • Recently, spike based methods are being developed in the hope of improving performance. But difficulties in these models are many.

    • Experimental paradigms to move the field from the present level need to address issues of:

      • Training (no desired response in paraplegic)

      • How to cope with coarse sampling of the neural population

      • How to include more neurophysiology knowledge in the design


    ad