Transcript Try to save computation time in large

Try to save computation time in largescale neural network
modeling with population density
methods, or just fuhgeddaboudit?
Daniel Tranchina, Felix Apfaltrer & Cheng Ly
New York University
Courant Institute of Mathematical Sciences
Department of Biology
Center for Neural Science
Supported by NSF grant BNS0090159
SUMMARY
•
Can the population density function (PDF) method be made
into a practical time-saving computational tool in large-scale
neural network modeling?
•
Motivation for thinking about PDF methods
•
General theoretical and practical issues
•
The dimension problem in realistic single-neuron models
•
Two dimension reduction methods
1) Moving eigenfunction basis (Knigh, 2000)
 only good news
2) Moment closure method (Cai et al., 2004)
 good news; bad news; worse news; good news
•
Example of a model neuron with a 2-D state space
•
Future directions
Why Consider PDF Methods in Network Modeling?
• Synaptic noise makes neurons behave stochastically: synaptic
failure; random sizes of unitary events; random synaptic delay.
• Important physiological role: mechanism for modulation of
gain/kinetics of population responses to synaptic input;
prevents synchrony in spiking neuron models as in the brain.
• Important to capture somehow the properties of noise in
realistic models.
• Large number of neurons required for modeling physiological
phenomena of interest, e.g. working memory; orientation
tuning in primary visual cortex.
Why Consider PDF Methods (continued)
• Number of neurons is determined by the functional
subunit; e.g. an orientation hypercolumn in V1:
• TYPICAL MODELS: ~ (1000 neurons)/(orientation
hypercolumn) for input layer V1 ( 0.5 X 0.5 mm2 or roughly
0.25 X 0.25 deg2 ).
• REALITY: ~ 34,000 neurons, 75 million synapses
• Many hypercolumns are required to study some problems,
e.g. dependence of spatial integration area on stimulus
contrast.
Why PDF? (continued)
• Tracking by computation the activity of ~103--104 neurons and
~104--106 synapses taxes computational resources: time and
memory.
e.g. 8 X 8 hypercolumn-model for V1 with 64,000 neurons
(Jim Wielaard and collaborators, Columbia):
1 day to simulate 4 seconds real time
• But stunning recent progress by Adi Rangan & David Cai
What to do?
Quest for the Holy Grail: a low-dimensional system of
equations that approximates the behavior of a truly highdimensional system.
Firing rate model (Dyan & Abbott, 2001): system of ODEs or
PDF model (system of PDEs or integro-PDEs)?
The PDF Approach
• Large number of interacting stochastic units suggests a
statistical mechanical approach.
• Lump similar neurons into discrete groups.
• Example: V1 hypercolumn. Course graining over: position,
orientation preference; receptive-field structure (spatialphase); simple--complex; E vs. I may give ~ 50
neurons/population (~ tens OK for PDF methods).
• Each neuron has a set of dynamical variables that
determines its state; e.g.
X  V,Ge ,Gi , for a leaky I&F neuron.
• Track the distribution of neurons over state space and
firing rate for each population.
Rich History of PDF Methods in Computational
Neuroscience
•Wilbur & Rinzel 1983
• Kuramoto 1991
• Abbott & van Vreeswijk 1993
• Gerstner 1995
PDF Methods Recently Espoused and Tested as a Faster
Alternative to Monte Carlo Simulations
• Knight et al., 1996
• Omurtage et al., 2000
• Sirovich et al., 2000
• Casti et al. 2002
• Cai et al., 2004
• Huertas & Smith, 2006
PDF Theory
r
 X k  state variables for population k
r r
r
  x1 , x2 ,..., xn ,t   joint probability density function for the
state variables of populations 1, 2,..., n
r r
r
r r
r
r
r
  x1, x2 ,..., xn ,t dx1    xn  Pr X1, X2 ,..., Xn   at time t





If populations that are not too densely coupled,
 x1 , x2 ,..., xn ,t  factors into 1 x,t 2 x,t    n x,t 
r r
r
r
r
r
 Evolution equations for the individual population density functions,
r r

r
k x,t   gJ k x,t , are coupled via population firing rates.
t
 rk t   total flux of probability across a surface in phase space,
a surface integral of probability flux.
 Can only capture behavior characterized by population firing rates.
Minimal I&F Model: How Many State Variables?
• Most applications of PDF methods as a computational tool
have involved single-neuron models with a 1-D state space:
instantaneous synaptic kinetics; V jumps abruptly up/down
with each unitary excitatory/inhibitory synaptic input event.
• Synaptic kinetics play an enormously important role in
determining neural network dynamics.
• Bite the bullet and include realistic synaptic kinetics.
• Problem with PDF methods: as underlying neuron model is
made more realistic, dimension of the state space increases,
so does the computation time to solve the PDF equations.
• Time saving advantage of PDF over (direct) MC vanishes.
• Minimal I&F model with synaptic kinetics has 3 state variables:
X  V,Ge ,Gi , voltage, excitatory and inhibitory conductances.
Take Baby Steps: Introduce Dimensions One at a
Time and See What We Can Do
 Consider an inetgrate & fire (I&F) neuron receiving excitatory
synaptic input only.
 Unitary excitatory-postsynaptic conductance event
 e (t) has a single-exponetial time course:
 t 
 e (t)  exp    , for t  0.
e
 e 
A
dV
1
dGe
1
Ak
  V   r   Ge (t) V   e  ;
  Ge    t  Tk 
dt
m
dt
e
k e
v nullcline
PDF EQUATIONS
v  mebrane voltage
g  excitatory postsynaptic conductance
t  time
 e (t)  total rate ofsynaptic input (from same and/or other populoations)
r

(v, g,t)   gJ (v, g,t)
t
r
J (v, g,t)  JV (v, g,t) , J G (v, g,t)
JV (v, g,t)  
v   r   g v   e  (v, g,t)
m
J G (v, g,t)  
1
1
e
g
%  (g  g ')(v, g ',t) dg '
g  (v, g,t)   e (t)  F
A
e
0

r(t)   JV (vth , g,t) dg. Boundary Condition: JV ( r , g,t)  JV (vth , g,t)
0
r
d r
  L  e (t) , for the discretized problem.
dt
PDF vs. MC and Mean-Field for 2-D Problem.
PDF cpu time is ~ 400 single uncoupled neurons
MC
1000 neurons
PDF
mean-field
100,000 neurons
cpu time: 0.8 s for PDF; 2 s per 1000 neurons for MC.
 e  5ms;  m  20ms;  EPSP  0.3mV;  r  65 mV; vth  55 mV.
Computation Time Comparison: PDF vs. Monte Carlo (MC):
PDF grows linearly; MC grows quadratically
50 neurons per population;1 run; 25% connectivity
• PDF Method is plenty fast for model neurons
with a 2-D state space.
• More realistic models (e.g. with E and I
input) require additional state variables
• Explore dimension reduction methods.
• Use the 2-D problem as a test problem
Dimension Reduction by Moving Eigenvector Basis:
Bowdlerization of Bruce Knight’s (2000) Idea
1) Discretize state-space variables:
r
d r
  L  e (t) .
dt
2) Approximate  e (t) as piecewise constant:
t k  k  t, but t is still continuous;  ek   e (t
 
1
k
2
)
d r
k r
  L  e ,
for t k  t  t k 1 .
dt
r
3) For each fixed  e and t k  t  t k 1 , use eigenvectors  j of L
as a basis set. Eigenvectors of L and LT are bi-orthogonal.
r
4) Equations for the coefficients ak (t) of  j are uncoupled first-order
ODEs, with analytic exponential slutions.
5) Compute only the first couple of coefficients corresponding to
eigenvalues with smallest real parts; throw the rest away.
Dimension Reduction by Moving Eigenvector Basis
Example with 1-D state space, instantaneous synaptic kinetics
• Only 3 eigenvectors
for low, and 7 for
high synaptic input
rates.
• Large time steps
• Eigen-method is 60
times faster than
full 1-D solution
Suggested by Knight, 2000.
Dimension Reduction by Moving Eigenvector Basis
Example with 2-D state space: state variables V & Ge
• Out of 625
eigenvectors: 10
for high, 30 for
medium, and 60
for low synaptic
input rates.
• Large time steps
• Eigen-method is
60 times faster
than full 1-D
solution
Dimension Reduction by Moment Closure

 r(t)   JV (vth , g,t) dg; JV (v, g,t)  
0
 r(t)  
1
m
(v
th
1
m
(v   )  g  (v   )(v, g,t)
r
e
  r )  G|V (vth ,t)  (vth   e ) f0 (vth ,t)
 The firing rate is determined by two functions of t and v only,
evaluated at threshold voltage: fV(0) (v,t) & G|V (v,t).
 No need to know the full PDF, (v,t) .

 

   (v,t) dg  r.h.s involving fV(0) (v,t) & G|V (v,t)
 t

0

 

   g (v,t) dg  r.h.s involving fV(0) (v,t), G|V (v,t) & G 2 |V (v,t)
 t

0
2
 Close the system by ansatz:  G|V
(v,t)   G2 (t)
2
 G 2 |V (v,t)   G2 (t)  G|V
(v,t)
Dimension Reduction by Moment Closure:
2nd Moment
Stimulus
Firing Rate Response
Near perfect agreement between results from dimension
reduction by moment closure, and full 2-D PDF method.
 e  5ms;  m  20ms;  EPSP  0.5mV;  r  65 mV; vth  55 mV.
Dimension Reduction by Moment Closure:
3rd Moment
Response to Square-Wave Modulation of Synaptic Input Rate
ZOOM
3rd-moment closure performs better than 2nd at high input rates.
 e  5ms;  m  20ms;  EPSP  0.5mV;  r  65 mV; vth  55 mV.
Trouble with Moment Closure and Troubleshooting
• Dynamical solutions “breakdown” when synaptic input
rates drop below ~ 1240 Hz, where actual firing rate
(determined by MC and full 2-D solution) ~ 60 spikes/s.
• Numerical problem or theoretical problem?
• Is moment closure problem ill-posed for some
physiological parameters?
• Examine the more tractable steady-state problem
Steady-State Moment Closure Problem: Existence Study
 A coupled pair of ODEs for fV(0) (v) and G|V (v) can be reduced to
a single ODE with a boundary condition (B.C.).
 Firing rate, r  0  
1
m
(v   )  
r
(0)
(v)

(


v)

f

G|V
e
V (v)  0
v  r
 G|V (v) 
e  v
v  r
q(v)  G|V (v) 
 0;
e  v
 mr
f (v) 
q(v)  ( e  v)
(0)
V
   e e  r  
v  r  
  G 
 q 1 





v


v
 
 

m
e
e
q m  

e 
 e  v



dq

dv
q 2   G2


+ B.C.
Phase Plane Analysis of Steady-State Moment Closure
Problem to Study Existence/Nonexistence of Solutions
 Introduce an auxiliary variable, s, and study the dynamical system:
dq
 m 1    e  e   r  
v   r  
 q
q 1 

   G 

ds
 e ( e  v)    m  e  v  
 e  v  
dv
 q 2   G2
ds
 Requirements for a v(s), q(s) trajectory to be a solution:
1) q   G at any point along trajectory.
v  r
2) q  G 
at least once.
e  v
3) q(v) satisfies its boundary condition.
Phase Plane and Solution at High Synaptic Input Rate
 e  1341Hz;  e  5ms;  m  20ms;  EPSP  0.5mV
solution trajectory
must intersect
must not intersect
Steady-State Solution Doesn’t Exist for Low Synaptic Input Rate
trajectory 1
must intersect
must not intersect
trajectory 2
Promise of a New Reduced Kinetic Theory with
Wider Applicability, Using Moment Closure
(Cai, Tao, Shelley, McLaughlin, 2004)
A numerical
method on a
fixed voltage
grid that
introduces a
boundary layer
with numerical
diffusion finds
solutions in
good
agreement with
direct
simulations.
SUMMARY
• PDF methods show promise
• Small population size OK, but connectivity cannot be dense
• Realistic synaptic kinetics introduce state-space variables
• Time saving benefit lost when state space dimension is high
• Dimension reduction methods could maintain efficiency:
• Moving eigenvector basis speeds up 2-D PDF method 60 X
• Moment closure method (unmodified) has existence problems
• Numerical implementations suggest moment closure can
work well
• Challenge is to find methods that work for >= 3 dimensions
THANKS
• Bruce Knight
• Charles Peskin
• David McLaughlin
• David Cai
• Adi Rangan
• Louis Tao
• E. Shea-Brown
• B. Doiron
• Larry Sirovich
Edges of parameter space:
e
Minimal
input rate:
Minimal e
r from 2D PDM
5 ms
1804.11 Hz
51.35 Hz
2 ms
2377.84 Hz
74.22 Hz
1 ms
2584.88 Hz
85.82 Hz
0.5 ms
2760.47 Hz
109.16 Hz
0.2 ms
539.05 Hz
<0.01 Hz
0.1 ms
535.19 Hz
<0.01 Hz
 r  70mV ,  EPSP  0.5mV
Minimal EPSP
with fixed mean G:
e
Min EPSP
r from 2D PDM
5 ms
7.47 mV
61.l5 Hz
2 ms
4.89 mV
65.70 Hz
1 ms
3.40 mV
69.24 Hz
 r  70mV , fix G  e  A at mean-field threshold,
increase EPSP (  A) until solution exists
Parameter Values
 Voltage values are: v [ r , vth ].
 With enough excitatory input, a neuron’s voltage will cross
threshold voltage, vth .
 Upon crossing threshold, a spike is said to occur,
and the membrane voltage is reset (instantaneously for sake of
simplicity) to  r .
 Typical parameter values (our values):
 r  65mV
vth  55mV
 m  20ms
 e  5ms
 A  2.44  10 4 s, giving  EPSP  0.5mV


0
0
fV (k ) (v,t)   g k (v, g,t) dg   g k fG|V (v, g,t) fV (0) (v,t) dg
 G k |V (v,t) fV (0) (v,t)