Intro To Computational Neuroscience

Computational Neuroscience H.
Riecke, Northwestern University

Introduction to
Computational Neuroscience
Hermann Riecke
Engineering Sciences and Applied Mathematics
Northwestern University
h-riecke@northwestern.edu
June 8, 2011
c _2011 Hermann Riecke
These notes are based to a large extent on the book Theoretical Neuroscience by P.
Dayan and L.F. Abbott (MIT Press).
To avoid copyright issues this version of the lecture notes does not include any
of the gures.
1
Computational Neuroscience H. Riecke, Northwestern University
Contents
1 Interesting and Fun Links 9
2 Introduction 10
3 Single Neuron 13
3.1 Passive, Single-Compartment Neuron . . . . . . . . . . . . . . . . . . . 13
3.2 Ion Channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.3 The Hodgkin-Huxley Model . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4 Conductance-Based Models: Additional Currents . . . . . . . . . . . . . 32
3.4.1 A-Type Potassium-Current . . . . . . . . . . . . . . . . . . . . . . 33
3.4.2 T-Type Calcium Current . . . . . . . . . . . . . . . . . . . . . . . 34
3.4.3 Sag Current I
h
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.4.4 Calcium-Dependent Potassium Current . . . . . . . . . . . . . . 39
3.5 Two-Dimensional Reduction of Hodgkin-Huxley Model . . . . . . . . . . 41
3.5.1 Phase-Plane Analysis . . . . . . . . . . . . . . . . . . . . . . . . 42
3.6 Integrate-and-Fire Model . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.7 Izhikevichs Simple Model . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4 Cable Equation 53
4.1 Linear Cable Theory: Passive Cable . . . . . . . . . . . . . . . . . . . . 56
4.2 Axons and Active Dendrites . . . . . . . . . . . . . . . . . . . . . . . . . 60
5 Synapses 61
5.1 Gap Junctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2 Chemical Synapses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.2.1 Modeling Chemical Synapses . . . . . . . . . . . . . . . . . . . . 66
5.2.2 Short-Term Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . 73
6 Firing Rates 76
6.1 Poisson Spike Trains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.2 What Makes a Neuron Fire? . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.2.1 Spike-Triggered Average . . . . . . . . . . . . . . . . . . . . . . 84
6.2.2 Linear Response Kernel and Optimal Stimulus . . . . . . . . . . 87
6.2.3 Receptive Fields: Visual System . . . . . . . . . . . . . . . . . . 90
6.3 LNP Cascade Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
2
7 Neural Networks: Rate Models 100
7.1 Firing-Rate Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.2 Feed-Forward Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7.2.1 Hubel-Wiesel Model for Visual Orientation Tuning . . . . . . . . . 104
7.2.2 Neural Coordinate Transformations . . . . . . . . . . . . . . . . . 107
7.3 Recurrent Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
7.3.1 Linear Recurrent Networks . . . . . . . . . . . . . . . . . . . . . 111
7.3.2 Oculomotor Control during Eye Saccades . . . . . . . . . . . . . 113
7.3.3 Limulus Vision. Contrast Enhancement . . . . . . . . . . . . . . 116
7.3.4 Persistent States . . . . . . . . . . . . . . . . . . . . . . . . . . 123
7.3.5 Associative Memory . . . . . . . . . . . . . . . . . . . . . . . . . 126
8 Statistical Approach 131
8.1 Bayes and the Cutaneous Rabbit . . . . . . . . . . . . . . . . . . . . . . 134
8.2 Integration of Multi-Sensory Information . . . . . . . . . . . . . . . . . . 138
A Basic Concepts of Numerical Methods for Differential Equations 151
B Linear Response Kernel using Functional Derivatives 155
3
Recommended Books:
Available online at the library:
P. Dayan and L.F. Abbott, Theoretical Neuroscience, MIT Press, 2001.http://
cognet.mit.edu.turing.library.northwestern.edu/library/books/view?isbn=
0262041995
G.B. Ermentrout and D. Terman, Foundations of Mathematical Neuroscience,
2010. https://wamprod1.itcs.northwestern.edu/amserver/cdcservlet?TARGET=
https%3A%2F%2Fsesame.library.northwestern.edu%3A443%2Fcgi-bin%2FssoEzpCheck.
pl%3Furl%3Dezp.2aHR0cDovL2R4LmRvaS5vcmcvMTAuMTAwNy85NzgtMC0zODctODc3MDgtMg--&Requ
1232338061&MinorVersion=0&ProviderID=http%3A%2F%2Fsesame.library.northwestern.
edu%3A80%2Famagent&MajorVersion=1&IssueInstant=2011-03-28T19%3A47%3A47Z
References
[1] B. Bathellier, D.L. Buhl, R. Accolla, and A. Carleton. Dynamic ensemble odor cod-
ing in the mammalian olfactory bulb: sensory information at different timescales.
Neuron, 57(4):586598, Feb 2008.
[2] R. Ben-Yishai, R. L. Bar-Or, and H. Sompolinsky. Theory of orientation tuning in
visual cortex. Proc. Natl. Acad. Sci., 92:3844, 1995.
[3] P. R. Brotchie, R. A. Andersen, L. H. Snyder, and S. J. Goodman. Head position
signals used by parietal neurons to encode locations of visual-stimuli. Nature,
375(6528):232235, May 1995.
[4] J. J. Bussgang. Crosscorrelation functions of amplitude-distorted gaussian sig-
nals. Technical report, MIT, 1952.
[5] M. V. Chafee and P. S. Goldman-Rakic. Matching patterns of activity in primate
prefrontal area 8a and parietal area 7ip neurons during a spatial working memory
task. J. Neurophysiol., 79(6):29192940, June 1998.
[6] E. J. Chichilnisky. A simple white noise analysis of neuronal light responses.
Network, 12(2):199213, May 2001.
[7] P. S. Churchland and T. J. Sejnowski. The computational brain. MIT, 1999.
[8] J. A. Connor and C. F. Stevens. Voltage clamp studies of a transient outward
membrane current in gastropod neural somata. J. Physiol., 213(1):21, 1971.
[9] J. A. Connor, D. Walter, and R. Mckown. Neural repetitive ring - modications of
hodgkin-huxley axon suggested by experimental results from crustacean axons.
Biophys. J., 18(1):81102, 1977.
4
[10] P. Dayan and L. F. Abbott. Theoretical Neuroscience. MIT Press, 2001.
[11] G. C. DeAngelis, I. Ohzawa, and R. D. Freeman. Receptive-eld dynamics in the
central visual pathways. Trends Neurosci., 18(10):451458, October 1995.
[12] A. Destexhe, D. Contreras, M. Steriade, T. J. Sejnowski, and J. R. Huguenard. In
vivo, in vitro and computational analysis of dendritic calcium currents in thalamic
reticular neurons. J. Neurosci., 16:169, 1996.
[13] A. Destexhe, Z.F. Mainen, and T.J. Sejnowski. Synthesis of models for excitable
membranes, synaptic transmission and neuromodulation using a common kinetic
formalism. J Comput Neurosci, 1(3):195230, Aug 1994.
[14] G. B. Ermentrout and N. Kopell. Parabolic bursting in an excitable system coupled
with a slow oscillation. SIAM J. Appl. Math., 46(2):233253, April 1986.
[15] G. B. Ermentrout and D. H. Terman. Foundations of Mathematical Neuroscience.
2010.
[16] M. O. Ernst and M. S. Banks. Humans integrate visual and haptic information in a
statistically optimal fashion. Nature, 415(6870):429433, Jan 2002.
[17] G. D. Field and E. J. Chichilnisky. Information processing in the primate retina:
Circuitry and coding. Annu. Rev. Neurosci., 30(4):130, 2007.
[18] I. M. Finn, N. J. Priebe, and D. Ferster. The emergence of contrast-invariant
orientation tuning in simple cells of cat visual cortex. Neuron, 54(1):137152,
April 2007.
[19] F. Gabbiani and S. J. Cox. Mathematics for Neuroscientists. 2010.
[20] F. Gabbiani, W. Metzner, R. Wessel, and C. Koch. From stimulus encoding to
feature extraction in weakly electric sh. Nature, 384(6609):564567, Dec 1996.
[21] M. N. Geffen, De Vries S. E. J., and M. Meister. Retinal ganglion cells can rapidly
change polarity from Off to ON. PLoS Biology, 5(3):065, March 2007.
[22] F. A. Geldard and C. E. Sherrick. The cutaneous rabbit: a perceptual illusion.
Science, 178(4057):178, 1972.
[23] D. Goldreich. A bayesian perceptual model replicates the cutaneous rabbit and
other tactile spatiotemporal illusions. PLoS ONE, 3:e333, 2007.
[24] M. S. Graziano, X. T. Hu, and C. G. Gross. Visuospatial properties of ventral
premotor cortex. J. Neurophysiol., 77(5):22682292, May 1997.
[25] H. K. Hartline. Visual receptors and retinal interaction. Science, 164:270, 1969.
[26] H. K. Hartline and F. Ratliff. Inhibitory interaction of receptor units in the eye of
limulus. J. Gen. Physiol., 40(3):357376, 1957.
5
[27] J. Hertz, A. Krogh, and R. G. Palmer. Introduction to the theory of neural compu-
tation. Addison-Wesley, 1991.
[28] A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current
and its applicaiton to conduction and excitation in nerve. J. Physiol., 117:500,
1952.
[29] A. L. Hodgkin, A. F. Huxley, and B. Katz. Measurement of current-voltage relations
in the membrane of the giant axon of loligo. J. Physiol., 116:424, 1952.
[30] A. L. Hodgkin and B. Katz. The effect of sodium ions on the electrical activity of
the giant axon of the squid. J. Physiol., 108:37, 1949.
[31] J. J. Hopeld. Neural networks and physical systems with emergent collective
computational abilities. Proc. N. Acad. U.S.A, 79(8):25542558, 1982.
[32] D. H. Hubel and T. N. Wiesel. Receptive elds, binocular interaction and functional
architecture in cats visual cortex. J. Physiol., 160(1):106, 1962.
[33] J. R. Huguenard and D. A. Mccormick. Simulation of the currents involved in rhyth-
mic oscillations in thalamic relay neurons. J. Neurophysiol., 68(4):13731383,
October 1992.
[34] E. M. Izhikevich. Simple model of spiking neurons. IEEE Trans. Neural Netw.,
14(6):15691572, November 2003.
[35] E. M. Izhikevich. Which model to use for cortical spiking neurons? IEEE Trans.
Neural Netw., 15(5):10631070, September 2004.
[36] E. M. Izhikevich. Dynamical Systems in Neuroscience: the geometry of excitability
and bursting. MIT Press, 2005.
[37] T. Jarsky, A. Roxin, W. L. Kath, and N. Spruston. Conditional dendritic spike prop-
agation following distal synaptic activation of hippocampal ca1 pyramidal neurons.
Nature Neurosc., 8(12):16671676, December 2005.
[38] E. R. Kandel, J. H. Schwartz, and T. M. Jessel. Principles of neural science.
McGraw-Hill, 2000.
[39] J. Keener and J. Sneyd. Mathematical Physiology. Springer, 1998.
[40] B. W. Knight, J. I. Toyoda, and F. A. Dodge. A quantitative description of dynamics
of excitation and inhibition in eye of limulus. J. Gen. Physiol., 56(4):421, 1970.
[41] V. I. Krinskii and IuM Kokoz. Analysis of the equations of excitable membranes.
I. reduction of the Hodgkins-Huxley equations to a 2d order system. Biozika,
18(3):506511, May-Jun 1973.
[42] L. Lapicque. Quantitative investigations of electrical nerve excitation treated as
polarization (translated by N. Brunel and M.C.W. van Rossum in Biol. Cyber. 97
(2007) 341). J. Physiol. Pathol. Generale, 9:620, 1907.
6
[43] H. Markram, Y. Wang, and M. Tsodyks. Differential signaling via the same axon of
neocortical pyramidal neurons. Proc Natl Acad Sci U S A, 95(9):53235328, Apr
1998.
[44] D. A. Mccormick and T. Bal. Sleep and arousal: Thalamocortical mechanisms.
Annu. Rev. Neurosci., 20(13):185215, July 1997.
[45] D. A. McCormick and H. C. Pape. Properties of a hyperpolarization-activated
cation current and its role in rhythmic oscillation in thalamic relay neurones. J
Physiol, 431:291318, Dec 1990.
[46] W. S. Mcculloch and W. Pitts. A logical calculus of the ideas immanent in nervous
activity (reprinted from bulletin of mathematical biophysics, vol 5, pg 115-133,
1943). Bull. Math.Biol., 52(1-2):99115, 1990.
[47] J. L. Mcfarland and A. F. Fuchs. Discharge patterns in nucleus-prepositus-
hypoglossi and adjacent medial vestibular nucleus during horizontal eye-
movement in behaving macaques. J. Neurophysiol., 68(1):319332, July 1992.
[48] A. M. Pastor, R. R. de la Cruz, and R. Baker. Eye position and eye velocity
integrators reside in separate brain-stem nuclei. Proc. N. Acad. U.S.A, 91(2):807
811, January 1994.
[49] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical
recipes, 2007.
[50] R. Rieke, D. Warland, R. de Ruyter van Steveninck, and W. Bialek. Spikes. Ex-
ploring the neural code. MIT, 1999.
[51] E. Salinas and L. F. Abbott. Transfer of coded information from sensory to motor
networks. J. Neurosci., 15(10):64616474, October 1995.
[52] H. S. Seung, D. D. Lee, B. Y. Reis, and D. W. Tank. Stability of the memory of eye
position in a recurrent network of conductance-based model neurons. Neuron,
26(1):259271, April 2000.
[53] M. N. Shadlen and W. T. Newsome. The variable discharge of cortical neurons:
Implications for connectivity, computation, and information coding. J. Neurosci.,
18(10):38703896, May 1998.
[54] G. M. Shepherd. Neurobiology. Oxford University Press, 1988.
[55] S. M. Sherman. Thalamic relays and cortical functioning. In S.M. Sherman
V.A. Casagrande, R.W. Guillery, editor, Cortical function: a view from the thala-
mus, page 107. Elsevier, 2005.
[56] A. P. Sripati, T. Yoshioka, P. Denchev, S. S. Hsiao, and K. O. Johnson. Spatiotem-
poral receptive elds of peripheral afferents and cortical area 3b and 1 neurons
in the primate somatosensory system. J. Neurosci., 26(7):21012114, February
2006.
7
[57] P. Stern, F. A. Edwards, and B. Sakmann. Fast and slow components of unitary
epscs on stellate cells elicited by focal stimulation in slices of rat visual cortex. J
Physiol, 449:247278, Apr 1992.
[58] T. Trappenberg. Fundamentals of Computational Neuroscience. Oxford, 2009.
[59] M. Tsodyks, K. Pawelzik, and H. Markram. Neural networks with dynamic
synapses. Neural Comput., 10(4):821835, May 1998.
[60] A. L. van der Meer, F. R. van der Weel, and D. N. Lee. The functional-signicance
of arm movements in neonates. Science, 267(5198):693695, February 1995.
[61] J. Ma W, J. M Beck, P. E Latham, and A. Pouget. Bayesian inference with proba-
bilistic population codes. Nature Neurosci., 2006.
[62] K. A. Zaghloul, K. Boahen, and J. B. Demb. Contrast adaptation in subthreshold
and spiking responses of mammalian y-type retinal ganglion cells. J. Neurosci.,
25(4):860868, January 2005.
8
1 Interesting and Fun Links
Here is a list of (randomly obtained) links that are of interest for this class or just plain
neural fun.
Flash movie of ions participating in action potential
http://www.brainu.org/action_potential_cartoon.swf
Movies by Rudolph and Destexhe
http://cns.iaf.cnrs-gif.fr/alain_movies.html
Graphical Hodgkin-Huxley Simulator
http://www.cs.cmu.edu/
~
dst/HHsim/
this simulator is also available as Matlab source code.
GUI for Izhikevichs Simple Model [34]
http://vesicle.nsi.edu/users/izhikevich/publications/whichmod.htm#izhikevich
Illusions:
Optical illusions and other visual phenomena
http://www.michaelbach.de/ot/index.html
Animated Necker cube
http://www.dogfeathers.com/java/necker.html
Barber Pole Illusion:
http://www.psychologie.tu-dresden.de/i1/kaw/diversesMaterial/www.
illusionworks.com/html/barber_pole.html
Moving Rhombus:
http://www.cs.huji.ac.il/
~
yweiss/Rhombus/
More optical illusions:
http://www.purveslab.net/seeforyourself/ http://web.mit.edu/persci/
index.html
Best Illusion Contests:
http://illusioncontest.neuralcorrelate.com/
Many interesting aspects of the brain
http://www.youramazingbrain.org.uk/
illusions of various kinds
http://www.youramazingbrain.org.uk/supersenses/default.htm
9
2 Introduction
What is the purpose of the brain?
based on sensory input
mouse sees cat, frog sees a y, ...
and past experience (memory)
learning from experience (narrow escape, tasty), ...
[even bees and ies can be trained: extend proboscis, ying left or right in a
T-channel,...]
the organism must decide what to do and then give generate motor output
jump, eta, run
Figures/ChSe99.f2.1.ps
Figure 1: The multiple scales of the brain [54] (reproduced in [7]) .
Studying brain function requires studies at a wide range of scales:
behavior of the animal
brain regions:
anatomy, connectivity between the regions, activities of the regions related to
behavior (e.g. fMRI measuring blood-oxygen-level-dependent (BOLD) signal)
note extensive feed-back between different regions (not marked in Fig.1)
10
circuits within brain regions:
e.g. mammalian cortex is arranged in columns (half a million, each 0.5mmwide,
2mm thick containing 50, 000 neurons )
columns have 6 layers: input, output to other columns, to other cortical areas, to
spinal cord and other areas outside of cortex
simpler circuits controling breathing, locomotion (worms), stomach grinding (lob-
ster),...
individual neurons:
wide variety of types with different characteristics: excitatory, inhibitory, wide
spread in morphology
neurons consist of
soma with nucleus (DNA and chemical machinery like in other cells that is
needed for metabolism)
axon: output of information to other neurons via an electric pulse (action
potential), currents are carried by ions
dendrite: input of information from other neurons as a result of their action
potential
Figures/webvision.med.utah.edu.visualcortex.ps
Figure 2: Different neuron morphologies in visual cortex. webvision.med.utah.edu/
VisualCortex.html
11
Figures/cajal_cerebellar_cortex.ps
Figure 3: Purkinje cells (A) and basket cells (B) in cerebellar cortex.
dendrites themselves can serve as computational units: inputs can gate other
inputs
electric activity of neuron due to ion channels across the cell membrane
they consist of a combination of a few large molecules, e.g. 4 molecules in
K
+
-channel (tetramer)
ion concentrations are different inside and outside of the cell (e.g. Na
+
much
higher outside)
ion channels typically permeable only to certain types of ions (e.g., Na
+
and
not K
+
, etc.)
ion channels open and close depending on voltage across membrane (voltage-
gated) resulting currents change voltage, which affects in turn ion chan-
nels
neurons are nonlinear devices: excitable by sufciently strong external input
action potential
Figures/hhsim.ps
Figure 4: Action potential simulationhttp://www.cs.cmu.edu/
~
dst/HHsim/.
large number of synapses provide connection between different neurons, mostly
between axons and dendrites ( 2 synapses/mof dendritic length, total O(10, 000),
can be O(100, 000) for Purkinje cell)
neurotransmitters are released frompre-synaptic terminal and diffuse across
the narrow synaptic cleft between the cells
12
neurotransmitters activate numerous types of receptors
similarly to ion channels, receptors are comprised of a combination of a few
large molecules, e.g. AMPA-receptor is a tetramer
binding of the neurotransmitter to the receptor opens an ion channel and
modifes voltage of post-synaptic cell
genetics
different cells express different sub-units of receptors or ion channels giving
them different properties
Computational and theoretical approaches:
detailed modeling at a given level/scale :
individual ion channel, single synapse, single neuron, possibly circuit consisting
of few neurons, ...
extracting sufciently many parameters from experiments is a challenge
effective models to identify essential mechanisms of brain functions:
networks of minimal models of neurons, interaction of neural populations repre-
senting brain areas
e.g. origin of oscillations, principles of memory
statistical approaches and information theory:
encoding of information (e.g. sensory) in spike trains: mean ring rate, timing of
individual spikes, correlation between spikes? Decoding the spike trains.
optimality considerations: natural scenes suggestions for brain function or
architecture (e.g. early visual processing: what type of information should the
retina transmit to the brain?)
3 Single Neuron
3.1 Passive, Single-Compartment Neuron
1
Neurons are bounded by a membrane comprised of a lipid bilayer
hydrophilic heads oriented towards electrolyte (extra- and intra-cellular uid)
hydrophobic tails inside the membrane
electric insulator
1
Dayan & Abbott: DA 5.2
13
Across the membrane a charge can build up: capacitance C
Q = CV C =
Q
V
Units:
charge [Q]=Coulomb, voltage [V ]=Volt, capacitance [C]=Farad
Convention:
The potential outside the cell is dened to be 0
Positive current is dened to be outward:
positive charges ow out of the cell and decrease the potential inside the cell:
V becomes negative
Typical voltages are in the range 90mV to +50mV . These voltages are not far above
the thermal range: at physiological temperatures T the thermal energy corresponds to
a voltage
V
T
=
k
B
T
q
25mV
and thermal uctuations can almost move charges aross these voltages (k
B
is the
Boltzmann constant, q the charge of an electron) expect noisy dynamics (i.e. even
for xed voltage, i.e. in voltage clamp, the current should be expected to uctuate)
A change in the voltage across the layer requires a current I
I
d
dt
Q = C
d
dt
V
The capacitance increases with surface area
C = cA c 1
F
cm
2
= 10
nF
mm
2
Typical cell areas are 0.01 to 0.1 mm
2
:
C 0.1 . . . 1nF
Axial current ow in the intracellular medium enclosed by the double-layer membrane
experiences a resistance R.
For a cylinder of length L and area A one has
R = r
L
A
r 1 . . . 3 k mm = 100 . . . 300 cm
Resistance increases with L (circuit of resistors in series)
14
Resistance decreases with A (circuit of parallel resistors)
unit: [r] = m
Current ow across the membrane possible through ion channels
proteins embedded in membrane
allow ions to diffuse through membrane
proteins can be selective for specic ions
many different types of ion channels
Figures/BeCo96.f2.8a.ps
Figure 5: Ion channels in membrane.
Estimate of channel conductance 1/R assuming that channel pore corresponds to a
cylinder of cellular uid:
L 6nm, A 0.15nm
2
: 1/R = 25pS = 25 10
12
S
unit is Siemens: 1 S = 1
1
the channels cover the membrane: conductance density g
unit: [g] = pS/cm
2
The membrane resistance depends on the density of open channels. It can vary sub-
stantially. For neuron in resting state specic membrane resistance r
m
r
m
1Mmm
2
R
m
=
r
m
A
10 . . . 100M
Unit: [r
m
] = m
2
since channels correspond to resistors in parallel and conductance is
proportional to membrane area.
15
Equivalent Circuit:
Kirchhoffs law:
the ions of the injected current I
e
ow either back out of the cell via the membrane
current I
m
or they increase the charge on the membrane corresponding to a capacitive
current I
C
I
e
= I
m
+ I
C
=
V
R
m
+ C
dV
dt
i.e.
dV
dt
=
1
R
m
C
V +
1
C
I
e
Introducing the membrane time constant
m
= R
m
C one gets
m
dV
dt
= V + R
m
I
e
(1)
Steady state
V
ss
= R
m
I
e
Relaxation towards steady state
2
V (t) = R
m
I
e
_
1 e
t
RmC
_
for V (0) = 0.
Exponential relaxation with time constant
m
= R
m
C.
Typical value:
m
1ms.
3.2 Ion Channels
3
The ion ux through the ion channels depends on the ion concentrations and the volt-
age across the membrane.
The ion concentrations inside and outside the cell are typically different:
[Na
+
] and [Ca
2+
] higher outside the cell
[K
+
] lower outside the cell
These differences are maintained by active ion pumps. Their activity requires metabolic
energy (ATP).
2
The general solution of a linear, inhomogeneous ode like (1) is given by the general solution V
h
of
the homogeneous equation
m
dV/dt = V (here V
h
= Ae
t/m
) and any particular solution V
p
of the
inhomogeneous equation (1) (here V
p
= V
ss
).
3
DA 5.2
16
Figures/sodium_potassium_pump.ps
Figure 6: Na
+
-K
+
-exchange pump.
Movie http://highered.mcgraw-hill.com/olc/dl/120068/bio03.swf
Combined with the semi-permeability of the membrane the concentration differences
lead in the steady state to a difference in the electrical potential across the membrane:
the current across the membrane vanishes for a non-zero voltage called the reversal
potential.
Nernst equation:
Consider only single ion species, e.g. [Na
+
], and a single ion channel for now:
V = 0
number of ions that enter channel from inside is proportional to [Na
+
]
inside
j
out
(V = 0) = k[Na
+
]
inside
with k some constant (related to the mobility of the ions).
number of ions that enter channel from outside is proportional to [Na
+
]
outside
j
in
(V = 0) = k[Na
+
]
outside
for [Na
+
]
outside
> [Na
+
]
inside
ions ux is inward: this ux leads to a net posi-
tive charge inside and a net negative charge outside: V becomes positive
for V > 0
potential is higher inside the cell and only those Na
+
ions can enter the cell
that have an enough energy (above qV ) to overcome the potential energy
barrier
probability distribution for the energy of the ions
p(E)dE = Ne
E
k
B
T
dE
with normalization
_

0
p(E)dE = 1 N =
1
k
B
T
17
fraction of ions with energy above E
0
= zqV
0
P(E > E
0
) =
_

E
0
p(E)dE = e
zqV
0
k
B
T
the ux against the positive potential is therefore given by
j
in
(V ) = k [Na
+
]
outside
e
zqV
k
B
T
(z = 1 for Na
+
).
ions leaving the cell need not overcome an energy barrier
j
out
(V ) = j
out
(V = 0) = k[Na
+
]
inside
for a certain voltage, V = E
Na
+, the system is in equilibrium the two uxes bal-
ance
j
in
(V = E
Na
+) = j
out
The reversal potential E
Na
+
for Na
+
is given by the Nernst equation (for Na
+
one has
z = 1)
E
Na
+ =
V
T
z
ln
_
[Na
+
]
outside
[Na
+
]
inside
_
V
T
=
k
B
T
q
.
Notes:
For V < E
Na
+ the electric eld (potential barrier) is not strong enough to keep
the Na
+
from entering the cell more positive charges ow into the cell and
make the potential more positive, i.e. V approaches E
Na
+ and conversely for
V > E
Na
+.
For xed concentrations [Na
+
]
inside,outside
the ion ux pushes the voltage to-
wards the reversal potential E
Na
+.
Squid Na
+
K
+
Cl
Ca
2+
intracellular concentration 50mM 400mM 40mM
extracellular concentration 440mM 20mM 550mM
ratio of concentrations ext/int 8.8 0.05 14
Nernst reversal potential +56mV -77mV -68mV
Birds/Mammals Na
+
K
+
Cl
Ca
2+
intracellular concentration 12mM 155mM 4mM 100nM
extracellular concentration 150mM 4mM 120mM 1.5mM
ratio of concentrations ext/int 12.5 0.03 30 15,000
Nernst reversal potential +67mV -98mV -90mV 130mV
Table 1: Rough values of concentrations. Birds and mammals have much lower con-
centrations than marine invertebrates
Thus, opening Na
+
-channels drive the potential to positive values, while K
+
and
Cl
-channels hyperpolarize cells to negative potentials.

18
In thermodynamic equilibrium it does not matter how ions diffuse through the
channel: result is independent of the mobility k.
If multiple ion species diffuse across the membrane, the steady state does in
general not correspond to thermodynamic equilibrium:
ion uxes can cancel each other:
although charge ux vanishes the individual ion uxes may be non-zero
the ion uxes must be treated explicitly
Goldman-Hodgkin-Katz Equation
4
In general:
multiple ions species diffuse across channel to get the steady state, i.e. the rest
potential, we need the ion ux j
driven by the electric eld E = for each species

Plancks equation
j
= zq
. .
force on ion
D
k
B
T
..
drag coefcient
c = u
z
[z[
c
with
z: valence of ion (zq is total charge of ion)
u: mobility of ion: includes the magnitude of the charge, more strongly charged
ions have larger u
related to Ficks diffusion constant D
u =
[z[q
k
B
T
D
c: concentration of the ion species
Total ion ux
j = D
_
c +
zq
k
B
T
c
_
In general the electric eld depends on the charge distribution (r) through the
Poisson equation
2
=
1
(r) = dielectric constant of the medium

Charge density changes due to the current j coupled equations for and j (Poisson-
Nernst-Planck equations)
4
Keener&Sneyd 2.6.2
19
Figures/channel_sketch.eps
Figure 7: Modeling one-dimensional ion channel.
Simplify for short, narrow channel:
consider the steady state: all quantities are time-independent
assume all quantities depend only on the x-coordinate along the channel:
one-dimensional model:
c
dc
dx
=
d
dx
choose coordinates such that inside of the cell is for x < 0 and outside of the cell
is for x > L:
V = (0) (L)
assume that the electric eld independent of position in the channel:
d
dx
=
(L) (0)
L
=
V
L
This yields the differential equation
dc
dx
c +
j
D
= 0 with =
zq
k
B
T
V
L
with boundary conditions
c(0) = c
inside
c(L) = c
outside
Note:
a rst order differential equation has only a single free parameter: how can we
satisfy both boundary conditions? The boundary conditions will pose an addi-
tional condition on a parameter in the equation: it will determine j.
Using integrating factor e
x
e
x
_
dc
dx
c
_
= e
x
j
D
can be written as
d
dx
_
e
x
c(x)
_
= e
x
j
D
20
integrating from x = 0 to x yields
e
x
c(x) c(0) =
1
j
D
_
e
x
1
_
with c(0) = c
inside
.
To satises the boundary condition at x = L the ux j has to satisfy the condition
j = D
c
inside
e
L
c
outside
1 e
L
To get a current density j
e
we need to multiply j with the charge of the ions. Measuring
c in moles/volume the charge per volume is = zFc with F the Faraday constant
(charge of 1 Mole of single-valence ions).
This yields the Goldman-Hodgkin-Katz current equation
j
e
=
z
2
Fq
k
B
T
P V
c
inside
e
L
c
outside
1 e
L
with L =
zq
k
B
T
V (2)
with membrane permeability
P =
D
L
Notes:
the current-voltage relationship j
e
(V ) is nonlinear ( through )
the membrane permeability in general will depend on the channel and the ion
species
Now we can consider cells with multiple ion species c
i
j
e
=
i
j
i
e
In the steady state the total current density j
e
must vanish.
For simplicity consider single-valence ions, z = 1, with concentrations c
i
0 =
i+
P
i+
c
i+
inside
c
i+
outside
e
q
k
B
T
V
1 e
q
k
B
T
V
+
i
P
i
c
i
inside
c
i
outside
e
q
k
B
T
V
1 e
q
k
B
T
V
rewrite it as
0 =
i+
P
i+
c
i+
inside
c
i+
outside
e
q
k
B
T
V
1 e
q
k
B
T
V
+
i
P
i
c
i
inside
e
q
k
B
T
V
c
i
outside
e
q
k
B
T
V
1
0 =
i+
P
i+
_
c
i+
inside
c
i+
outside
e
q
k
B
T
V
_
i
P
i
_
c
i
inside
e
q
k
B
T
V
c
i
outside
_
21
One can solve this equation for V
V
GHK
=
k
B
T
q
ln
_
i+
P
i+
c
i+
outside
+
i
P
i
c
i
inside
i+
P
i+
c
i+
inside
+
i
P
i
c
i
outside
_
(3)
Notes:
for a single ion species V
GHK
reduces to Nernst reversal potential
5
. For multiple
species the concentrations are weighted with their respective permeability.
at the rest potential V = V
GHK
the electric current across the membrane vanishes
the concentration currents of the various ion species do in general not vanish
the concentration differences across the membrane need to be maintained
by ion pumps to maintain the steady state
the result depends in particular on the assumption of a spatially constant electric
eld along the channel. Depending on the channel type this assumption may or
may not be a good approximation.
unlike the Nernst potential the GHK-result is not universal: P depends on the
mobility of the various ion species in that channel
for the standard three ion species one gets (for z = 1)
V
GHK
=
k
B
T
q
ln
_
P
Na
[Na
+
]
outside
+ P
K
[K
+
]
outside
+ P
Cl
[Cl
]
inside
P
Na
[Na
+
]
inside
+ P
K
[K
+
]
inside
+ P
Cl
[Cl
]
outside
_
For example
Na
+
K
+
Cl
Ca
2+
intracellular concentration 15mM 100mM 13mM
extracellular concentration 150mM 5mM 150mM
ratio of concentrations 10 0.05 11.5
Nernst Potential +62mV -80mV -65mV
P/P
K
0.025 1 0.1
Because of the high membrane permeability for K
+
the resting membrane poten-
tial is dominated by the K
+
ions with a slight depolarization by Na
+
.
3.3 The Hodgkin-Huxley Model
How do neurons function?
5
Note that in the derivation of (3) we have assume z = 1.
22
What makes the action potential?
What is a good mathematical model for it?
Bernstein:
measured action potentials (1968)
based on Nernst equation and knowledge of the concentration difference across
membrane suggested that membrane is semi-permeable to K
+
: V
rest
70mV
(1912)
speculated that during action potential membrane becomes permeable to all
types of ions, which would lead to a break-down of the membrane potential.
Cole and Curtis (1940):
establish that conductance of membrance changes signicantly during action po-
tential
Membrane currents can still not be measured directly until Cole and Marmont develop
space clamp technique (1940s): insert very narrow electrode down the center of the
axon to make voltage spatially uniform and measure associated current.
Figures/KeSn98.f4.1.ps
Figure 8: The giant squid axon is not an axon from the giant squid (from [39])
To measure currents and action potentials in axons they used the giant axon of squid:
100 times thicker than axons in mammals, made experiments feasible.
Hodgkin and Curtis independently show that during action potential the voltage over-
shoots, i.e. crosses V = 0: Bernsteins suggestion of a general permeability of the
membrane cannot be correct.
Hodgkin and Katz (1949):
23
identify Na
+
as essential for action potential generation by reducing NaCl in ex-
tracellular uid by replacing it with some other substance like dextrose or choline
chloride (membrane is impermeable to choline ions)
Figures/HoKa49.f4.ps
Figure 9: Reduction of action potential upon partial replacement of extracellular sea-
water by dextrose, reducing extracellular [Na
+
] . a) 67% dextrose, b) 50%, c) 29%.
Curve 1 and 3 (wash-out) are with seawater [30]
this explains the overshoot, since Na
+
has positive reversal potential
permeability of membrane must change selectively for ions: there must be a
period during which the permeability for Na
+
dominates
Hodgkin and Huxley (1952, 5 papers, one with Katz)
measure currents as a function of membrane voltage:
membrane is rectifying (Fig.10)
hyperpolarization beyond rest potential (V < V
rest
) induces no current
depolarization induces time-dependent current
initially negative (inward) current
the inward current depolarizes the membrane further
later positive (outward) current
24
Figures/HoHu52a.f11.ps
Figure 10: Channels are rectifying: no current for hyperpolarization (to -130mV, top),
transient current for depolarization (to 0mV, bottom gure). Note: sign convention of
current opposite to the one used nowadays. [29]
for channels to open the membrane needs to be sufciently depolarized
shape of transient depends on the voltage step (Fig.11)
different channels open for different voltages and they have different driving
forces (V E)
Figures/HoHu52a.f12.ps
Figure 11: Shape of transient current depends on the membrane voltage.The numbers
marking the curves denote hyper-polarization relative to resting potential: negative
numbers indicate depolarization [29].
extract K
+
and Na
+
-current separately by replacing parts of the extracellular
NaCl with choline chlorid (Fig.12)
25
Na
+
-component is transient despite steady applied voltage
Figures/I_Na_separation.ps
Figure 12: Extraction of Na
+
-current by partial replacement of extracellular Na
+
by
choline. Na
+
-current is transient.
measure conductances: I(V )-curves
procedure:
long depolarizing voltage step to V
1
to open the channel
then second step to V
2
to probe the conductance of the channel at that
voltage V
2
measure current I(V
2
) immediately after second voltage step before the
channel state changes
K
+
(Fig.13)
initial current I(V
2
):
linear I(V
2
)-relationship: channels represent an Ohmic resistance
(Fig.13)
I
K
= g
K
(V E
K
)
current I(V
2
) vanishes at reversal potential E
K
of K
+
slope of I(V
2
)=conductance g
K
(V
2
) depends on V
1
and on duration
of rst stimulus
asymptotic current I(V
1
) (reached after channels adjust to voltage V
1
):
goes to 0 for strongly negative V
1
: channels close
becomes linear for strongly positive V
1
: once channels are maximally
open the channel behaves like an Ohmic resistor
26
Figures/K_I_V_curve_bill.ps
Figure 13: Ohmic resistance of K
+
-channel.
Na
+
:
also Ohmic resistance for given channel state
Evolution of K
+
-conductance (Fig.14)
after step in voltage the conductance g
K
evolves smoothly to asymptotic
value
Figures/HoHu52.f2.ps
Figure 14: Sigmoidal growth and exponential decay of g
k
with voltage [28].
increase from g
K
= 0: sigmoidal (tanh-like) increase
decrease from g
K
> 0: exponential
no oscillations in the evolution
rst-order evolution equation may be sufcient
simplest attempt: linear equation
n
dn
dt
= n
(V ) n (4)
27
n
(V ) asymptotic value reached after a sufciently long time (t

n
)
evolution after step in V :
decay after step from n
(V ) > 0 to n
= 0 is exponential
n = n(0)e
t/n
growth after step from n
= 0 to n
(V ) > 0 is initially linear

n(t) = n
(V )
_
1 e
t/n
_
n
(V )
t
n
not sigmoidal
Hodgkin-Huxley proposed
g
K
= g
K
n
4
(5)
with the exponent 4 resulting from a t to the data (Fig.14)
n
(V ) and
n
(V ) extracted from data by tting
Evolution of Na
+
-conductance
main difference to K
+
- current: the Na
+
-current is transient
Hodgkin-Huxley proposed
g
Na
(V, t) = g
Na
m
3
h
exponent 3 again from tting the data
activation variable m
m
dm
dt
= m
(V ) m
inactivation variable h
h
dh
dt
= h
(V ) h
difference between activation and inactivation variable:
m
(V ) increases with increasing voltage

h
(V ) decreases with increasing voltage

Complete Hodgkin-Huxley model for a single-compartment neuron
C
m
dV
dt
+ g
K
n
4
(V E
K
) + g
Na
m
3
h (V E
Na
) + g
L
(V E
L
) =
I
A
(6)
with
injected current I dened positive if inward (opposite to the denition of the cur-
rents through the ion channels. A is area of cell membrane
28
leak current g
L
(V E
L
) contains ungated channels (e.g. for Cl
and other ions)

Hodgkin and Huxley solved the coupled differential equations numerically with a me-
chanical calculator.
Using ts to the voltage dependencies of m
(V ),
m
(V ) etc. that they measured ex-
perimentally for the squid giant axon they obtained very good agreement with the ex-
perimentally observed temporal evolution of the currents after applied voltage steps
(Fig.16) (Nobel Prize in 1963, jointly with Eccles).
The ts obtained by Hodgkin and Huxley for the squid axon [28] are
n
=
0.01 (V + 55)
1 e
0.1(V +55)

n
= 0.125 e
0.0125(V +65)
m
=
0.1 (V + 40)
1 e
0.1(V +40)

m
= 4e
0.0556(V +65)
h
= 0.07e
0.05(V +65)
h
=
1
1 + e
0.1(V +35)
where the evolution equation (4) for n is written as
dn
dt
=
n
(V ) (1 n)
n
(V )n (7)
with
n
(V ) =
1
n
(V ) +
n
(V )
n
(V ) =

n
(V )
n
(V ) +
n
(V )
and analogously for m and h.
Here voltages are in mV and times are in ms.
Figures/DaAb01.f5.10.ps
Figure 15: Steady-state levels m
, h
, n
and the time scales

m,h,n
for the activation
and inactivation variables of the Na
+
- and K
+
-channels of the Hodgkin-Huxley model.
29
Figures/HoHu52.f11.ps
Figure 16: Comparison of current trace from numerical simulation of HH-model (left)
with the experimental curves (right). [28].
Simulations (use, e.g., hhsim, cf. Fig.4)
mechanism underlying the generation of the action potential:
depolarizing stimulation
activation variable m of Na
+
increases rapidly (cf. Fig.15)
Na
+
-channels open Na
+
rushes into the cell leading to further rapid
depolarization (positive feed-back)
inactivation variable h of Na
+
and activation variable of K
+
evolve more
slowly but speed up as cell becomes more depolarized: Na
+
channels close
and K
+
channels open Na
+
-ux is stopped and K
+
ows out of the cell
voltage decreases, cell repolarizes again
upon reaching V
r
some K
+
-channels are still open
overshoot of V beyond V
r
inhibition can also trigger action potential (rebound)
hyperpolarization beyond V
r
by injecting negative current increases inactiva-
tion variable h of Na
+
(Na
+
-channels become deinactivated) and decreases
activation n variable of K
+
(K
+
-channels become more inactivated)
at off-set of negative current mincreases rapidly to its equilibrium value while
h decreases only slowly and n also increases only slowly increased Na
+
-
current that is not balanced by increased K
+
-current depolarization of the
cell that can lead to an action potential
steady current injection periodic spike train (homework: does frequency change
signiciantly when current is reduced below the value needed for stimulating ac-
tion potential? )
stimulus during action potential does not trigger second action potential: refrac-
tory period (see homework)
30
Notes:
Hodgkin-Huxley modeled squid axon:
the same Na
+
- and K
+
- currents arise across all types of species (evolutionarily
conserved), with somewhat modied parameters HH-model forms basis for
almost all other neuron models
modeling of the currents:
in HH-model I(V )-dependence is assumed linear
GHK-current (2) depends nonlinearly on V
so far theory not good enough to be predictive from rst principles: theory
provides functional form of dependence. Parameters need to be determined
by tting to experiments.
Usually no drastic difference, given the variability between cells. (Fig.17)
Figures/KeSn.f3.1.ps
Figure 17: a) Comparison of the Ohmic resistance in the Hodgkin-Huxley model with
the current given by the Goldman-Huxley-Katz model (2) [39].
b) Voltage at which no current ows when channels for Na
+
and K
+
are both open
as a function of the K
+
concentration for GHK-model and for a combination of Ohmic
resistances for the Na
+
- and the K
+
-conductance, with E
Na
+ and E
K
+ determined
separately.
other neurons can exhibit more complex behavior due to additional channel types
Interpretation of the activation/inactivation equations:
Write (4) again as
dn
dt
=
n
(V ) (1 n)
n
(V )n (8)
with
n
(V ) =
1
n
(V ) +
n
(V )
n
(V ) =

n
(V )
n
(V ) +
n
(V )
31
or equivalently
n
=
n
n
=
1 n
n
Eq.(7) can be read as the evolution equation for the probability n of a gate to be open
n(t + t) = n(t)
..
probability for gate to be open
+
+ t
n
(V )
. .
probability for a closed gate to open
(1 n(t))
. .
probability of gate to be closed
t
n
(V )
. .
probability for an open gate to close
n(t)
..
probability for a gate to be open
alternatively n can also be read as the mean number of open gates.
If
l gates are required to be open for the channel to be open and
the gates open and close independently of each other
then the probability for the channel to be open (= fraction of open channels) is given by
g
K
g
K
= n
l
Thus, (5) suggests that the K
+
-channel consists of 4 independently operating gates
that all have to be open for the K
+
-channel to allow ions to pass.
Microsopically, 4 equivalent gates have been identied only very recently by MacK-
innon in the form of 4 proteins, which can switch between 2 different conformations
(geometric arrangements) one of which blocks the channel (Nobel Prize in 2003).
Analogously, for the Na-channel one can think of m as the probability of one of three
gates at the channel pore to be open and of h as the probability that an additional
blocking particle is not sitting in the pore of the channel. The probability for the pore to
let charge pass is then m
3
h. [Note: this picture is not to be taken too literally, although
the blocking particle can be thought of as a tethered plug formed by parts of the -unit
protein that forms the core of the channel.]
3.4 Conductance-Based Models: Additional Currents
6
In Hodgkin-Huxley model only two active currents are included
6
DA Chap 6.2
32
transient Na
+
-current
delayed rectifying K
+
-current (kicks in more slowly than Na
+
)
In the brain there is a zoo of neurons, which exhibit a variety of additional currents
giving the neurons specic properties.
Very commonly they are modeled similar to the currents in the HH-model in terms of a
reversal potential and a
Ohmic resistance with a conductance that may be controled by
an activation variable
an inactivation variable
activation and inactivation variables are modeled by rst-order differential equa-
tions
3.4.1 A-Type Potassium-Current
The A-type K
+
- current was introduced by Connor & Stevens [8, 9] to model repetitive
ring in walking-leg neurons of crustaceans. They can re repetitive spikes at a very
low frequency. Not possible in HH-model.
In contrast to the delayed rectifying K
+
-current, the A-type K
+
-current I
A
inactivates.
The Connor-Stevens model includes both K
+
-currents.
I
A
= g
A
a
3
b (V E
A
)
with
activation variable a
a
da
dt
= a
(V ) a
inactivation variable b
b
db
dt
= b
(V ) b
33
Figures/connor_stevens_activation.eps
Figure 18: Voltage dependence of the parameters in the Connor-Stevens model. Note
that the parameters for I
Na
and I
K
are also different than in the HH-model [9].
One nds (see homework)
I
A
delays spiking after depolarizing voltage step
ring rate grows from 0 when spiking threshold is crossed, allows arbitrarily slow
repetitive ring: Type-I neuron (HH-neuron is Type-II neuron)
3.4.2 T-Type Calcium Current
Thalamus functions as a relay fromsensory organs to cortex and also between different
cortical areas.
Figures/brain_facts_p5b.ps
Figure 19: Brain structures: thalamus (above spinal cord and hindbrain) functions as
relay for sensory information.
34
in awake states and during REM sleep the relay neurons re single spikes: signal
can be transmitted faithfully
during slow-wave sleep (EEG shows slow waves) the thalamo-cortical relay neu-
rons re repetitive bursts
T-type Ca
2+
-current is involved in the repetitive bursting [44, 33].
I
CaT
= g
CaT
M
2
H (V E
Ca
)
with M and H satisfying
M
dM
dt
= M
(V ) M
H
dH
dt
= H
(V ) H
Figures/HuMc92.f1.AB.ps
Figures/HuMc92.f1.D.ps
Figure 20: a) Steady-state activation and inactivation M
and H
for I
T
, c) time con-
stant
M
of activation, d) time constant
H
of inactivation. Note the very slow recovery
from inactivation for V < 70mV [33].
Notes:
35
extracellular Ca
2+
-concentration is much higher than the intracellular one
E
Ca
= +120mV , I
CaT
is depolarizing similar to the Na
+
-current
I
CaT
is transient (like Na
+
-current);
almost no overlap of activation and inactivation curves for no voltage signi-
cant steady current
I
CaT
is much slower than Na
+
-current; deinactivation of I
CaT
is very slow(Fig.20d)
I
CaT
can induce spiking even without I
Na
: Ca-spikes (broader than Na
+
-spikes)
Figure 21: Burst due to I
CaT
. Delay of the action potential due to I
A
. Note decrease in
frequency across burst [10].
Burst after hyperpolarization:
injection of hyperpolarizing current I
CaT
deinactivates, i.e. inactivation is re-
moved: H 1
stop current injection increase in voltage activates I
CaT
further depolariza-
tion
strong depolarization triggers multiple Na-spikes, which occur on a much faster
time scale than the variation in I
CaT
activation
eventually I
CaT
becomes inactivated (H 0) and bursting stops
Tonic depolarizing current injection does not deinactivate I
CaT
: single Na-spikes, no
bursting
What could be the function of these two modes?
tonic spiking regime: single spikes, relay cells can transmit sensory input faithfully
in a graded manner
bursting regime: immediately above the threshold for spike output the spike fre-
quency is high during the burst (could be like a wake-up call).
36
Figures/Sh05a.f2.ps
Figure 22: Response of a thalamic neuron (in LGN) in the tonic and the bursting regime.
Same current injection induces tonic ring from a holding potential of V = 59mV (a)
but to a transient burst for V = 70mV (b). Graded vs. discontinuous response of the
neuron in the two regimes (c). In the tonic regime the ring rate encodes the sinusoidal
input very well (d), in the bursting regime only the onset of each wave is encoded (e)
[55]
Both modes arise in a given neuron in the awake and the sleeping animal. But the tonic
mode increases in prevalence with increased alertness/wakefulness of the animal.
3.4.3 Sag Current I
h
In thalamic neurons I
CaT
can exhibit not only single bursts but also rhythmic bursting
(1-2Hz).
Figures/McPa90.f14.ps
Figure 23: Rhythmic bursting of thalamic relay neurons. [45]
37
This is mediated by a hyperpolarization-activated cation current, I
h
, [45]
permeable to Na
+
and K
+
reversible potential E
h
40mV
I
h
is activated below V 70mV
activation and deactivation is slow O(100ms 1, 000ms)
Figures/McPa90.f1.eps
Figure 24: Sag current I
h
. Current and voltage clamp results. The sag refers to
the slow decrease in hyperpolarization after strong hyperpolarizing current injection
(arrows). [45]
Figures/McPa90.f2a.eps Figures/McPa90.f2b.eps
Figure 25: a) Activation curve of I
h
(obtained from 7 different neurons, solid symbols
1 neuron). b) Time scales for activation and deactivation, c) enlargement of the tail
current after repolarization to V = 49mV (between arrows in b). [45]
38
Figures/McPa90.f3.eps
Figure 26: Voltage dependence of activation and deactivation times. Fits to transients
using single exponentials are quite good. [45]
Rhythmic bursting when the cell receives input (or leak or other current) that tends to
hyperpolarize it:
the I
h
-current can depolarize the cell from hyperpolarized states
activates I
CaT
inducing a (rebound) burst
during the depolarized state I
h
deactivates
when I
CaT
inactivates the cell becomes hyperpolarized again due to the inputs
(or other conductances)
I
h
becomes activated again ....
3.4.4 Calcium-Dependent Potassium Current
Central pattern generators control locomotion (swimming of lamprey, walking, running)
or chewing and digestive rhythms. Often the control is in the form of periodic bursting.
For instance, California spiny lobster (and other crustaceans) grinds its food in the
stomach. Controlled by stomatogastric ganglion cells that re repetitive bursts with
low frequency in the absence of any input. The frequency depends mainly on a Ca
2+
-
dependent K
+
-current I
KCa
.
I
KCa
is mainly activated by Ca
2+
, although activation may also depend somewhat on
voltage.
Model it as
I
KCa
= g
KCa
c
4
(V E
K
)
39
with c satisfying
c
dc
dt
= c
(V, [Ca
2+
]) c
Activation variable c depends on voltage V and the Ca
2+
-concentration [Ca
2+
].
Evolution of [Ca
2+
]:
Ca
d[Ca
2+
]
dt
=
Ca
I
Ca
[Ca
2+
]
Ca
2+
inux I
Ca
through Ca
2+
-channels (I
CaT
), note I
Ca
< 0
slow buffering by Ca
2+
-buffers (remove Ca
2+
from the solution) with time constant
Ca
Figure 27: Periodic bursting in a model for a stomatogastric ganglion cell.
Burst Cycle:
I
CaT
induces Ca-spike with Na-spikes riding on top of it (burst)
rise in [Ca
2+
] activates I
KCa
burst turned off by inactivation of I
CaT
and by activation of I
KCa
after burst [Ca
2+
] decreases (slow buffering) and deactivates I
KCa
when I
CaT
is deinactivated (very slow process; Fig.20) and I
KCa
is deactivated
next burst can arise.
Note:
over multiple spikes I
KCa
builds up an after-hyperpolarization conductance (AHP)
40
it increases time until next spike
frequency of spikes decreases despite constant stimulation: spike-frequency adap-
tation
Thus:
even individual, single-compartment (point) neurons can exhibit a wide range of
behavior (see also Sec.3.7)
3.5 Two-Dimensional Reduction of Hodgkin-Huxley Model
Reconsider HH-neuron. Even standard HH-model with only fast Na-conductance I
Na
and delayed rectier I
K
is a four-dimensional dynamical system: hard to visualize.
Can the equations be simplied to get at least a semi-quantitative intuitive understand-
ing?
Observations:
activation variable m of Na
+
-current is very fast
inactivation variable h of Na
+
and activation variable n of K
+
are strongly corre-
lated
h 0.89 1.1n
replace
m(t) m
(V )
h(t) 0.89 1.1n(t)
Figures/hh_h-n_V-n.eps
Figure 28: Strong correlation between n and h in standard Hodgkin-Huxley model.
41
Figures/hh_adiabatic_approx.eps
Figure 29: Deviations of the activation variables m, n, h from their steady-state values
m
, n
, h
during regular ring. The large values of m/m
and h/h
are due to m
and h
becoming very small during deactivation and inactivation, respectively.

These assumptions yield reduction of HH-model to a two-dimensional system
C
dV
dt
= g
K
n
4
(v E
K
) g
Na
m
3
(V ) (0.89 1.1n) (V E
Na
) (9)
g
L
(V E
L
) + I
e
F
V
(V, n)
n
(V )
dn
dt
= n
(V ) n F
n
(V, n) (10)
Note:
this type of reduction underlies also the FitzHugh-Nagumo model (cf. [39] [41])
3.5.1 Phase-Plane Analysis
7
Many aspects of two-dimensional dynamical systems can be understood using a phase-
plane analysis
7
For introductions to nonlinear dynamics see, e.g., S.H. Strogatz Nonlinear Dynamics and Chaos or
the lecture notes for Interdisciplinary Nonlinear Dynamics 438 http:// people.esam.northwestern.
edu/
~
riecke/new_ www/ lectures.htm on the web
42
Figures/Iz05.f5.20.ps
Figure 30: a) stable xed point, large perturbation can lead to excitation. b) unstable
xed point, periodic spiking.
To get V
rest
need xed points of (9,10). They are given by the intersections of the two
nullclines
F
V
(V, n) = 0 (11)
F
n
(V, n) = 0 (12)
The nullcline F
V
separates regions with
dV
dt
> 0 from those with
dV
dt
< 0, and analo-
gously for F
n
.
Any xed point (V
0
, n
0
) satises (11,12) simultaneously.
What happens in the vicinity of the xed point?
Consider small deviations (V, n)
V = V
0
+ V n = n
0
+ n
and linearize (9,10) around the xed point (V
0
, n
0
)
d
dt
(V
0
+ V ) = F
V
(V
0
+ V, n
0
+ n
0
) F
V
(V
0
, n
0
)
. .
=0
+
F
V
V
(V
0
n
0
)
V +
F
V
n
(V
0
n
0
)
n
d
dt
(n
0
+ n) = F
n
(V
0
+ V, n
0
+ n
0
) F
n
(V
0
, n
0
)
. .
=0
+
F
n
V
(V
0
n
0
)
V +
F
n
n
(V
0
n
0
)
n
d
dt
V =
F
V
V
(V
0
n
0
)
V +
F
V
n
(V
0
n
0
)
n (13)
d
dt
n =
F
n
V
(V
0
n
0
)
V +
F
n
n
(V
0
n
0
)
n (14)
43
We have a linear system of homogeneous differential equations with constant coef-
cients exponential ansatz
V = V
1
e
t
n = n
1
e
t
_
v
1
n
1
_
e
t
=
_
F
V
V
(V
0
n
0
)
F
V
n
(V
0
n
0
)
Fn
V
(V
0
n
0
)
Fn
n
(V
0
n
0
)
_
_
v
1
n
1
_
e
t
L
_
v
1
n
1
_
e
t
with L the Jacobian of (9,10) at the xed point (V
0
, n
0
)
has to be an eigenvalue of L and
_
v
1
n
1
_
is the associated eigenvector
8
The number of eigenvalues is equal to the number N of time-derivatives of the system.
The general solution of (13,14) is a linear superposition of the eigenmodes with N
unknown constants to satisfy arbitrary initial conditions
_
V (t)
n(t)
_
= A
1
_
v
(1)
1
n
(1)
1
_
e
1
t
+ A
2
_
v
(2)
1
n
(2)
1
_
e
2
t
The eigenvalues determine the qualitative character of the ow in the phase space
near the xed point:
real :
saddle (
1
2
< 0), stable (
1,2
< 0) or unstable (
1,2
> 0) node
complex z
1,2
i:
L is real
2
=
2
and v
(2)
= v
(1)
V (t) = e
t
_
A
1
v
1
e
it
+ A
1
v
1
e
it
_
stable ( < 0) and unstable ( > 0) spiral
xed point is linearly stable all eigenvalues have negative real parts
8
To calculate the eigenvalues you determine in general for which values of the determinant of the
matrix L I vanishes. Here I is the identity matrix. In the 2x2 case at hand here one can also solve
the rst row of the equation for n
1
,
n
1
=
_

F
V
V
(V0n0)
_
v
1
1
FV
n
(V0n0)
. (15)
Inserting this expression in the second row yields
_

F
n
n
(V0n0)
__

F
V
V
(V0n0)
_
v
1
1
FV
n
(V0n0)
=
F
n
V
(V0n0)
v
1
.
For v
1
,= 0 this results in a quadratic equation for the eigenvalues
(1,2)
and (15) yields an equation for
the eigenvectors
_
v
(1,2)
1
n
(1,2)
1
_
associated with these eigenvalues.
44
Figures/page215a.eps
Figure 31: Generic (typical) trajectories in the phase plane close to a xed point. Eigen-
vectors need not be orthogonal.
Note:
linear stability: stable with respect to innitesimal perturbations
even a linearly stable xed point can be unstable to nite-size perturbations
Need to get insight into global properties of the phase plane
Often the two variables evolve on different time scales (time-scale separation):
fast-slow analysis allows approximation to action potential, ring frequency, ...
If V is much faster than n: fast-slow analysis
V reaches its nullcline F
V
(V, n) = 0 quickly without n changing much
most of the time the system follows the V -nullcline F
V
(V, n) = 0 V = V
nc
(n)
slow one-dimensional dynamics driven by n(t)
dn
dt
= F
n
(V
nc
(n), n)
at turning points (V
TP
, n
TP
) the system cannot follow the V -nullcline any more:
fast one-dimensional evolution ( jump) to another branch of the V -nullcline
without much change in n
dV
dt
= F
V
(V, n
TP
)
Can construct an approximation to the complete action potential from simpler pieces.
45
3.6 Integrate-and-Fire Model
9
The Hodgkin-Huxley model is computationally slow
4 differential equations: V , m, n, h
even in the two-dimensional reduction the equations are stiff:
resolving the voltage spike requires very small time steps
time between spikes can be long compared to the spike width
single neuron: adaptive time step to step quickly between spikes (implicit
scheme for stability)
many coupled neurons: likely that some neuron spikes for any given time
small time steps all the time or complex code that updates only neurons
undergoing spike
Consider simplied models that capture essential qualitative features of the neuron
The shape of the voltage spike is often not important:
do not model the dynamics of the spiking dynamics explicitly:
do not resolve the dynamics of activation and inactivation/deactivation of the fast
Na
+
-current and the delayed K
+
-rectier.
Assume the conductances of the Na
+
- and K
+
currents do not depend on volt-
age: lump all reversal potentials and conductances together into
m
dV
dt
= E
L
V + R
m
I
e
(16)
the input current I
e
(injection or from other neurons) can drive the voltage to the
spiking threshold V
th

the voltage is reset explicitly to a reset voltage V
reset
V (t
spike
) = V
th
V (t
+
spike
) = V
reset
(17)
spike-like output, which serves as input to other neurons, is generated by
hand :
V
out
(t) = V
sp
(t t
spike
)
or
V
out
(t) = V
sp
_
e
1
e
2
_
with
2
<
1
This leaky Integrate-and-Fire model was proposed by Lapicque already in 1907 by
Lapicque [42] http://www.springerlink.com/content/x03533370x281257/.
Note:
9
DA 5.4
46
for small uctuations of V around resting potential the assumption of constant
Na
+
and K
+
- conductance is reasonable; further away from V
rest
it is not a
good assumption
Firing Rate:
IF-model allows simple analytical solution for I
e
= const.
m
d
dt
_
e
t/m
V
_
= e
t/m
(E
L
+ R
m
I
e
)
V (t) = (E
L
+ R
m
I
e
)
_
1 e
t/m
_
+ V
0
e
t/m
with V
0
V (t = 0). Rewrite using
V
E
L
+ R
m
I
e
as
V (t) = V
(V
V
0
) e
t/m
Without the spike triggering mechanism one would have for t
V (t) V
For small I
e
, more precisely for V
< V
th
, no spike is triggered.
For
V
> V
th
i.e. I
e
>
1
R
m
(V
th
E
L
)
a spike is triggered at time t
spike
, which is obtained by solving V (t
spike
) = V
th
for t
spike
t
spike
(V
0
) =
m
ln
_
V
V
th
V
V
0
_
The spike time t
spike
(V
0
) depends on the initial condition V
0
.
For xed I
e
the neuron res periodically:
the initial condition for the voltage after the n
th
-spike is given by V (t
(n)+
spike
) = V
reset
the inter-spike interval (ISI) T
ISI
t
(n+1)
spike
t
(n)+
spike
is given by
T
ISI
=
m
ln
_
V
V
reset
V
V
th
_
=
m
ln
_
1
Vreset
V
1
V
th
V
_
for V
V
th
> V
reset
and the ring rate by
r
ISI
=
_
m
ln
_
V
V
reset
V
V
th
__
1
47
Notes:
as I
e
is increased the spiking sets in rst with vanishing ring rate:
in that respect the Type-I neuron is like the neuron of the Connor-Stevens model
(cf. the standard HH-neuron sets in with a xed nite frequency: Type-II neuron)
for large I
e
and large V
one has
T
ISI

m
ln
_
(1
V
reset
V
)(1 +
V
th
V
)
_

m
V
th
V
reset
V
Thus, for large I

e
the ring rate grows linearly with I
e
,
r
ISI

1
m
R
m
V
th
V
reset
I
e
in the IF-model (16,17) sufciently large inputs I
e
trigger spikes with arbitrarily
small ISI (no absolute refractory period, see below)
time-dependence of the voltage in the IF-model does not show the up-swing that
is characteristic of the regenerative action potential; instead the voltage has a
negative second derivative (concave downward) before the spike
Refractory Period
In the Hodgkin-Huxley model the Na
+
-current is still inactivated shortly after an action
potential: generating another action potential is more difcult
To get an action potential the depolarizing Na
+
-current must overcome the hyperpolar-
izing K
+
-current. This requires
g
K
n
4
(V E
K
) < g
Na
m
3
h (V E
Na
)
V < V
bal

g
K
n
4
E
K
+ g
Na
m
3
hE
Na
g
K
n
4
+ g
Na
m
3
h
V
bal
determined by reversal potentials that are weighted by the degree that the respec-
tive channels are open.
For action potential we need roughly V
bal
V
peak
When Na
+
-current is inactivated, h 1, V
bal
E
K
m reaches quickly m
(V ) 1, while n and h vary only slowly.

Shortly after action potential:
h very small
48
V
bal
E
k
< V
peak
even for m = 1
even arbitrarily strong inputs cannot trigger action potential, i.e. no regener-
ative event that is driven by the activation of Na is possible: absolute refractory
period
Later after an action potential
larger h
if V
bal
> V
peak
for m > m
(V
min
)
for inputs that drive the cell to V = V
min
the activation variable can quickly
grow to m
(V
min
) and push V
bal
above V
peak
action potential triggered by super-threshold inputs
as long as h has not recovered, triggering an action potential requires larger
inputs than the inputs to trigger an action potential from the resting state: relative
refractory period
ring threshold decreases with time
Modeling in IF-model:
absolute refractory period: explicitly disallow ring for a duration after any spike
relative refractory period:
V
reset
< V
rest
: effective ring threshold decreases as V relaxes from V
reset
to
V
rest
introduce a decaying hyper-polarizing conductance (cf. I
KCa
) or time-dependent
(relaxing) threshold V
th
(cf. [10]) or some generalized recovery variable (cf.
Sec.3.7)
3.7 Izhikevichs Simple Model
10
Aim: model a wide spectrum of types of neurons phenomenologically with a single,
computationally efcient model
Extend the IF-model to include additional variable that may capture effects of I
A
, I
KCa
,
I
CaT
, ...
If the spike-shape is not essential: focus on subthreshold behavior
in reduced HH-model we the dynamics of V and n
10
[39, 35, 36]
49
extensions of HH (I
A
, I
CaT
, I
KCa
, ...) would introduce additional conductances
with dynamics that are slow compared to the spike generating Na
replace n by a generalized recovery variable u
Approximate the nullclines of V and u near the minimum of the V -nullcline
V -nullcline: u = u
min
+ P (V V
min
)
2
u-nullcline: u = s (V V
0
)
thus
V
dV
dt
= P (V V
min
)
2
(u u
min
)
u
du
dt
= s (V V
0
) u
Adding the injected I current, Izhikevich writes the model in the form [35]
dv
dt
= 0.04v
2
+ 5v + 140 u + I (18)
du
dt
= a (bv u) (19)
Note:
numbers chosen such that time is measured in ms, v in mV
v can diverge in nite time
dv
dt
= v
2
_
v
v
0
dv
v
2
=
_
t
0
dt
_
1
v

1
v
0
_
= t
v =
v
0
1 v
0
t
for t
1
v
0
> 0
escape from the region near the left branch of the V -nullcline: action potential
the nonlinearity gives an up-swing before the action potential
Supplement (18,19) with a reset condition like in the integrate-and-re model
v(t
sp
) = 30mV
_
v(t
+
sp
) = c
u(t
+
sp
) = u + d
(20)
Notes:
50
a is the relaxation time of u
b characterizes the dependence of the steady-state value of the recovery variable
u on v, e.g. deinactivation of I
CaT
by hyperpolarization
c is a reset voltage, in general not the resting potential; it could, for instance,
include an after-hyperpolarization (for V
reset
< V
rest
)
d characterizes the change in the recovery variable during the spike excursion
to the other branch of the nullcline, e.g. Ca
2+
-inux controling I
KCa
or the inacti-
vation of I
CaT
due to the high voltage during the spike.
size of d provides a second time scale to the slow evolution of u
Figures/Iz03.f2a.ps
Figure 32: Parameters of Izhikevichs simple model.
Note:
the simple model is computationally efcient and captures qualitative aspects of
many neuronal types (Fig.34)
the model can be thought of as an extension of the quadratic integrate-and-re
model (QIF)
dv
dt
= b + v
2
with (21)
v(t
sp
) = v
peak
v(t
+
sp
) = v
reset
(22)
(21) is the normal form for a saddle-node bifurcation
the saddle-node bifurcation leads to a periodic orbit due to the reset, which im-
plements - by hand - the existence of an attracting invariant circle in phase space.
51
Figures/page155.eps
Figure 33: Sketch of a pair of xed points on an invariant circle (parametrized by
the phase of the oscillation during regular ring). The mutual annihilation of the
xed points in a saddle-node bifurcation generates an innite-period limit cycle (SNIC-
bifurcation leading to Type-I oscillation).
the normal form for the transition to spiking in a saddle-node bifurcation on an
invariant circle (SNIC) is the -neuron by Ermentrout-Kopell [14]
d
dt
= (1 cos ) + (1 + cos ) g (23)
where describes the phase of the neuron (oscillator) and = corresponds
to a spike, and g represents input to the neuron. The rest state is near = 0. For
small and weak forcing it reduces to (21),
d
dt
=
2
+ g +O(
4
, g
2
)
For parameters close to the saddle-node bifurcation at g = 0 the neuron spends
almost all of the time near the ghost of the xed point the period T can be
approximated by the time it takes the neuron to go from = to = +
_
+
2
+ g
=
_ T
2
T
2
dt
i.e.
T
1
g
g arctan

=
1
thus the ring rate goes to zero like r

g, with g measuring the distance from
the saddle-node bifurcation.
52
Figures/Iz03.f2b.ps
Figure 34: Some of the neuronal types captured by Izhikevichs simple model [35].
Compare behavior labeled thalamo-cortical with that discussed for I
CaT
-current in
Sec.3.4.2.
4 Cable Equation
11
So far we assumed that the neuron is characterized by a single value of the voltage:
point neuron.
Many neurons are, however, very long and branched.
what determines how good the single-compartment model is?
how to describe spatial structure of neurons
11
DA 6.3
53
Figures/neuron_morph_kath.ps
Figure 35: Neurons are not pointlike. Their morphology varies widely.
Figure 36: Voltage attenuation along a dendrite. A) current injection in soma triggers
large action potential there but delayed and reduced action potential in dendrite. B)
converse arrangement with stimulation (via synapses) of dendrite [10].
Voltage can be strongly attenuated further away from the stimulation site
Most parts of a neuron are long and thin: cable of radius a
assume voltage depends only on the axial position, but not on the transverse
coordinate: one-dimensional theory, radius may depend on position as well: a =
a(x)
two types of currents:
54
longitudinal current I
L
transverse current across membrane I
m
Kirchhoffs law for charge conservation
Figure 37: Sketch of currents in an axon [10]
Longitudinal current I
L
:
Consider Ohms law for a section of length x of the cable
I
L
=
V (x + x) V (x)
R
L
=
V
R
L
current driven by the voltage drop V across the section
consider I
L
positive if it ows in the positive x-direction
the current (positive charges) ows from the more positive potential (voltage) to
the more negative potential
in terms of the intra-cellular resistivity r
L
R
L
=
x
a
2
r
L
For small x
I
L
=
a
2
r
L
V
x
becomes I
L
=
a
2
r
L
V
x
for x 0
similar to Ficks law of diffusion: current driven by gradient in voltage
55
Transverse current I
m
:
In terms of the membrane current density i
m
I
m
= 2axi
m
Similarly, write the current injected by an electrode in terms of the density i
e
I
e
= 2axi
e
Kirchhoffs law for charge conservation for a cable segment of length x
2a(x)xc
m
V
t
. .
capacitive current
=
_
a(x)
2
r
L
V
x
_
left
. .
I
L
into the segment
+
_
a(x)
2
r
L
V
x
_
right
. .
I
L
out of the segment
2a(x)xi
m
. .
im out of the segment
+2a(x)xi
e
. .
ie into the segment
Dividing by x and using
1
x
_
_
V
x
_
right
_
V
x
_
left
_

x
_
V
x
_
for x 0
one obtains
c
m
V
t
=
1
2ar
L
x
_
a(x)
2
V
x
_
i
m
+ i
e
(24)
diffusive partial differential equation for the voltage
in general i
m
contains all the ion conductances discussed for the point neuron
(single compartment neuron): in a Hodgkin-Huxley framework the PDE is coupled
to ODEs for the activation and inactivation variables (they do not contain any
spatial derivatives)
4.1 Linear Cable Theory: Passive Cable
To make analytical progress consider passive cable with only leak current
i
m
= g
L
(V E
L
)
For passive neuron
g
L
=
1
r
m
E
L
= V
rest
56
Rewrite in terms of deviation v from V
rest
: v = V V
rest
For simplicity assume constant radius a = const.
c
m
v
t
=
a
2r
L
2
v
x
2

1
r
m
v + i
e
Using
m
= r
m
c
m
the linear cable equation can be written as
m
v
t
=
2
2
v
x
2
v + r
m
i
e
(25)
with
=
_
ar
m
2r
L
has the dimension of a length: this electrotonic length characterizes the length
scale on which the voltage varies in the cable
with increasing length x of a segment
the membrane resistance goes down
R
m
(x) =
1
2ax
r
m
the longitudinal resistance of the intracellular medium goes up
R
L
(x) =
x
a
2
r
L
at what length is the two resistances equal?
R
m
() =
r
m
2a
= R
L
() =
r
L
a
2
2
=
ar
m
2r
L
=
2
Thus: the electrotonic length is that length for which the membrane resistance
and the longitudinal resistance are equal
R
=
r
m
2a
=
r
L
a
2
Consider injecting current at a point. Flowing away from the injection site the current
ows
along the cable
across the membrane
57
The situation is somewhat similar to the ow of water in a soaker house (i.e. a leaky
pipe).
How the current is distributed between these two paths depends on the relative resis-
tances:
if R
L
< R
m
the current predominantly ows into the cable and spreads away from
the injection site
if R
L
> R
m
the current predominantly ows across the membrane and does not
spread away from the injection site
the current spreads to a length at which R
L
() = R
m
(), i.e. = .
Explicit solution for the steady state v = v(x) for a point injection current
i
e
(x) =
I
e
2a
(x)
with (x) being the Dirac -function
12
2
d
2
v
dx
2
= v r
m
I
e
2a
(x) (26)
Consider the two domains separately
x > 0 : the general solution is given by
v
+
(x) = A
+
e
+x/
+ B
+
e
x/
x < 0 : the general solution is given by
v
(x) = A
e
+x/
+ B
e
x/
boundary conditions:
for x the voltage must remain nite:
v(x) =
_
B
+
e
x/
for x > 0
A
e
x/
for x < 0
At x = 0 we need to match the two solutions:
d
2
v
dx
2
must be innite to balance (x), i.e. the rst derivative
dv
dx
makes a jump
d
2
v
dx
2
= lim
x0
dv
dx
x+
1
2
x
dv
dx
x
1
2
x
x
12
Denition of the Dirac--function: (x) = 0 for x ,= 0 and
_
+
(x) dx = 1. The -function has

dimension 1/length.
58
dv
dx
is well-dened on both sides of x = 0
v(x) is continuous
13
A
= B
+
v
0
Integrate (26) across x = 0
2
_
+
d
2
v
dx
2
dx =
_
+
v(x) r
m
I
e
2a
(x) dx
In the limit 0 the integral
_
+
v(x) dx goes to 0 since the integrand is nite
2
_
dv
dx
dv
dx
_
= r
m
I
e
2a
2
2
v
o
= r
m
I
e
2a
v
0
=
r
m
I
e
4a
=
1
2
R
I
e
Thus
v(x) =
1
2
R
I
e
e
|x|/
As expected, the injected current spreads into the cable to a distance .
Notes:
the point-neuron model is a good approximation as long as the spatial extend of
the neuron is smaller than

a for thin axons or thin dendrites, which have relatively more surface
area than cross section, is small and the space-dependence of the voltage has
to be taken into account to capture, e.g., the propagation of the action potential
along the axon to the next neuron.
The cable equation (25) is a diffusive equation: a pulse-like current injection diffuses
away from the source like heat (cf. Fig.38).
13
can be derived by taking antiderivative twice of (26). Alternatively, if v(x) was not continuous, its rst
derivative would be a -function and its second derivative would be the derivative of the -function; but
the equation contains only a -function to balance.
59
Figure 38: Voltage v(x) for steady current injection (A) and for short current pulse (B)
[10].
4.2 Axons and Active Dendrites
For almost all neurons the linear cable theory is insufcient
all axons support action potentials
most dendrites have voltage-gated channels or Ca
2+
-dependent channels, some
even support action potentials
i
m
in (24) includes currents like I
Na
, I
K
, I
CaT
, I
KCa
, ...
coupled, nonlinear PDE on a complex geometry: cannot expect to solve it analyti-
cally
Numerical solution:
discretize dendrite and/or axon in space compartments (small cylinders)
ODEs in each compartment for V and activation/inactivation variables
coupling between the compartments via the current owing into and out of the
compartments
Numerical simulation package/language:
NEURON developed by M. Hines, J.W. Moore, and T. Carnevale: NEURON web site
http://www.neuron.yale.edu has free downloads for any relevant operating system,
documentation, tutorials. In addition, it has a large database called ModelDB, which
contains NEURON codes (models) that have been developed and used for neuro-
science publications. (For example see Fig. 39)
60
Figures/DeCo96.ps
Figure 39: Burst in thalamic reticular neuron driven by low-threshold Ttype Ca
2+
-
current which is slower than the T-current in thalamocortical relay cells (cf. Sec.3.4.2)
[12]. For movies see
Thalamic reticular neurons http://cns.iaf.cnrs-gif.fr/alain_movies.html
Figures/gating_roxin.ps
Figure 40: Gating of dendritic inputs by other dendritic inputs in pyramidal cells in
CA1 in hippocampus [37]. Movie at http://people.esam.northwestern.edu/
~
kath/
gating.html.
5 Synapses
Connections between neurons are provided by synapses
electrical synapses: gap junctions
chemical synapses
5.1 Gap Junctions
Gap junctions provide a direct electrical connection between neurons
61
Figures/gap_junction_bill.ps
Figure 41: Sketch of gap junctions.
channel (hole) formed by the protein connexin, of which there exist different
types
hemi-channels (connexons) in the two cell membranes connect to each other
and form the gap junction
gap junction channels are permeable for any ion (unspecic)
permeability can be modulated by neuromodulators
Ohmic resistance
I
gap
= g
gap
(V
1
V
2
) with V
i
voltage of cell i
Notes:
direct electrical coupling tends to make coupled neurons to behave similarly, e.g.
synchronizes them
gap junctions are bi-directional, coupled cells affect each other reciprocally
5.2 Chemical Synapses
Chemical synapses do not provide a direct electrical coupling between cells.
The information is transmitted by a neurotransmitter
62
Figures/synapse_bill.ps
Figure 42: Electromicrograph of synapse displaying many vesicles and the active zone
(darker area at synaptic cleft).
Mechanism of synaptic transmission
neurotransmitter is stored in vesicles (small bags made of a lipid bilayer) in the
pre-synaptic terminal
depolarization of the pre-synaptic neuron by an incoming action potential induces
Ca
2+
-inux into the pre-synaptic cell
increased [Ca
2+
] cause vesicles to dock at the cell membrane and merge with it.
This releases the neurotransmitter into the synaptic cleft (narrow space) between
the pre-synaptic terminal and the post-synaptic cell
the neurotransmitter diffuses across the cleft ( 20 40nm)
binding of neurotransmitter to the receptors activates them
ionotropic receptor
binding opens ion channel directly
ion ux through channel creates
excitatory postsynaptic current (EPSC) (Na
+
and K
+
, also Ca
2+
for
NMDA) induces excitatory postsynaptic potential (EPSP)
inhbitory postsynaptic current (IPSC) (Cl
) induces a IPSP
response is fast and decays quickly (in millisecond range)
metabotropic receptor
binding activates a G-protein in the membrane, which triggers a signal-
ing pathway
eventually opens channels
can induce other changes (e.g. protein synthesis)
63
response is slow and can last very long (up to hours)
neurotransmitter is removed from the cleft through re-uptake by the pre-synaptic
cell
new vesicles are formed in pre-synaptic terminal by budding
Figures/synapse_mousedb_erlangen.ps
Figure 43: Sketch of chemical synapse(www.biochem.uni-erlangen.de/MouseDB)
Notes:
transmission is uni-directional from the pre-synaptic to the post-synaptic cell
pre-synaptic cell is not affected at all by post-synaptic cell
different neurotransmitters bind to different receptors, which open different chan-
nels
the post-synaptic cell can become excited or inhibited by the pre-synaptic
cell
synapses can amplify signals: single vesicle release can open thousands of ion
channels
Dales principle:
essentially all neurons express only a single neurotransmitter type a neuron
can generate either excitatory or inhibitory output:
it is either an excitatory neuron or an inhibitory neuron
neurons can express excitatory and inhibitory receptors:
they can receive excitatory and inhibitory inputs simultaneously
64
Excitatory synapses
neurotransmitters: mostly glutamate and also kainate
two major types of glutamate receptors:
AMPA (responds also to AMPA
14
):
activates very fast and deactivates quickly ( 5ms)
permeable to Na
+
and K
+
(mixed kation channel)
reversal potential E
AMPA
0mV
Ohmic I(V )-dependence
NMDA (responds also to NMDA
15
):
activates more slowly ( 2ms) and deactivates much more slowly (
150ms) than AMPA
permeable to Na
+
, K
+
, and to Ca
2+
reversal potential E
NMDA
0mV
nonlinear I(V )-dependence:
extra-cellular Mg
2+
binds to the pore in the open channel and blocks
it, depolarization of the postsynaptic cell expells the Mg
2+
Inhibitory synapses
neurotransmitters: mostly GABA and glycine
two types of GABA receptors:
GABA
A
ionotropic receptor
permeable to Cl
with reversal potential E

Cl
80mV
GABA
B
metabotropic receptor
slower response
can induce long-lasting change inside the post-synaptic neuron
permeable to K
+
Synapses respond probabilistically
14
AMPA=alpha-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid
15
NMDA=N-methyl-D-aspartic acid
65
pre-synaptic neuron: docking of vesicles at the membrane is a random process
docking energy is comparable to thermal energy
process characterized by probability of release
vesicles have a typical size: quantal release consisting of 1, 2, ... vesicles
EPSCs and IPSCs are roughly quantized (multiples of a minimal unit)
post-synaptic neuron: binding with post-synaptic receptor is also a random pro-
cess
Thus:
not each action potential needs to lead to synaptic transmission
even without an action potential one can have spontaneous release of neuro-
transmitter
different types of synapses have different levels of reliability
if many synapses connect the same neurons the random uctuations average out
Note:
by contrast, gap junctions are not probabilistic.
5.2.1 Modeling Chemical Synapses
16
Ionotropic receptors are like other ion channels
I
s
= g
s
(V E
syn
)
with E
syn
the reversal potential of the synapse, which determines whether it is excitatory
or inhibitory.
The conductance is determined by the number of open synaptic channels:
g
s
= g
s
P(t)
where
P = Prob(postsynaptic channel open)
= Prob(release occurs)
. .
pre-synaptic
Prob(channel open [ transmitter released)
. .
post-synaptic
= P
rel
P
s
(27)
16
DA 5.8
66
Evolution of P
s
:
Model the probability P
s
(t) of a channel to be open
dP
s
dt
=
s
(1 P
s
)
s
P
s
(28)
opening rate
s
depends on the concentration of the neurotransmitter:
for larger concentration the probability that a neurotransmitter molecule is close
to the receptor and therefore can bind to the receptor is larger
closing rate
s
is typically assumed to be independent of the concentration
channel closes when the neurotransmitter molecule unbinds from the receptor
Need to model the transmitter concentration in the cleft.
Very simple square-wave model
triggered by the pre-synaptic action potential at t = 0 the concentration jumps to
a positive value and remains constant for a duration T
s

s
0 < t T
after the duration T the concentration jumps to 0 again
s
= 0 T < t
Evolution of P
s
P
s
(t) =
_
1 + (P
s
(0) 1) e
st
0 < t T
P
s
(T)e
s(tT)
T < t
(29)
Figure 44: Fit of (29) to experimentally obtained EPSC (averaged over many stimuli)
[13].
Note:
67
release of multiple vesicles will not change the concentration signicantly com-
pared to teh release of single vesicles
if the synaptic cleft is much narrower than the distance between release
sites, since then concentration at receptor right across from release site will
be much larger than at neighboring receptors
if re-uptake is very quick, since then the neurotransmitter will be removed
before the next vesicle is released.
then one has essentially independent release sites and times, each of which
contributes P
rel
P
s
to the overall current.
Fast Synapses
For fast synapses with large
s
one need not resolve P
s
(t) during the short rise time:
replace the evolution by a jump.
For P
s
(0) = 0 the jump is given by P
s
= 1 e
sT
.
For general P
s
(0) one can rewrite P
s
(T) as given by (29) as
P
s
(T) = P
s
(0) + P
s
(1 P
s
(0))
with
P
s
= 1 e
sT
.
Model the synapse then by
_
P
s
P
s
+ P
s
(1 P
s
) when synaptic event occurs
dPs
dt
=
s
P
s
between synaptic events
Note:
Since P
s
1 one has P
s
1 after the spike
due to the slow decay of P
s
subsequent synaptic currents can overlap
Slow Synapses
Need to resolve also the rising phase of the synaptic current. The rise time introduces a
delay in the information transmission. This can have signicant effects on the dynamics
of the neuronal network.
Could use (28) with solution (29). The switching from one solution to the other at t = T
is awkward: one would need to keep track of the spike and the time switching time T
later.
For simplicity, often two other approaches are taken (keeping in mind that the square-
wave model is not really realistic either)
68
1. Difference of two exponentials
For P
s
(0) = 0
P
s
(t) = P
s
N
_
e
t/
1
e
t/
2
_
(30)
with
1
>
2
Decay time:
for large t
P
s
(t) = P
s
N e
t/
1
decay
=
1
Maximal conductance is reached at
t
max
= ln
_
2
_

1
2
ln
_
2
_
rise
N normalizes P
s
(t) such that P
s
(t
max
) = P
s
,
N =
_
_
1
_

2
1
_

1
2
_
1
For general P
s
(0) implement the evolution of P
s
(t) as
P
s
(t) = N (A(t) B(t)) (31)
where between spikes
1
dA
dt
= A
2
dB
dt
= B
and at spike the amplitudes are reinitiated (somewhat analogous to treat-
ment of fast synapse)
A(t
+
) A(t
) + P
s
_
1 P
s
(t
)
_
(32)
B(t
+
) B(t
) + P
s
_
1 P
s
(t
)
_
(33)
such that P
s
(t) is continuous at spike:
A(t
+
) B(t
+
) = A(t
) B(t
)
2. -function
for P
s
(0) = 0
P
s
(t) = P
s
t
s
e
t/s
has decay time
s
and reaches maximal conductance at
t
max
=
s
69
Example:
Consider 2 IF-neurons with neuron 1 giving excitatory synaptic input to neuron 2 and
vice versa
dV
1
dt
= V
1
+ I
e
(t) g
1
P
2
(t) (V
1
E
syn
) (34)
dV
2
dt
= V
2
+ I
e
(t) g
2
P
1
(t) (V
2
E
syn
) (35)
with V
threshold
= 30 and V
reset
= 55. P
i
(t) is given by P
s
(t) triggered by spike in neuron
i (assuming P
rel
= 1).
a)
Figures/2cells_IF_membrane_nosummation.eps
b)
Figures/2cells_IF_membrane_summation.eps
Figure 45: Integration of multiple inputs by slow membrane. IF-modell (34,35) with
= 1,
1
= 0.05, g
1
= 0,
2
= 0.025, g
2
= 50, E
syn
= 50 a) I
e
= 50 b) I
e
= 250.
Fig.45a:
weak synaptic coupling: single spike of neuron 1 does not trigger a spike in
neuron 2
inter-spike interval (ISI) large enough for most channels to close between synap-
tic events
slow membrane dynamics (large
m
):
some temporal summation of post-synaptic currents by the membrane
Fig.45b:
higher ring frequency (smaller ISI): more spikes during membrane relaxation
stronger integration in time second neuron spikes.
multiple synaptic inputs needed to trigger spike in neuron 2
70
a)
Figures/2cells_IF_synapse_nosummation.eps
b)
Figures/2cells_IF_synapse_summation.eps
Figure 46: Integration of multiple inputs by slow synapse. IF-modell (34,35) with =
0.1,
1
= 0.5, g
1
= 0,
2
= 0.25, g
2
= 10, E
syn
= 50, a) I
e
= 302, b) I
e
= 400.
Fig.46:
synapses slower than membrane:
channels do not close during ISI temporal summation of open synaptic chan-
nels
large fraction of channels are open persistently
Synaptic Conductance g
s
:
For AMPA and GABA receptors g
s
is independent of voltage.
For NMDA receptors the conductance depends on the voltage of the post-synaptic
neuron
g
s
= g
s
(V
post
)
near the resting potential the NMDA receptors are blocked by Mg
2+
ions
the Mg
2+
-block is reduced/removed as voltage is increased
g
NMDA
=

g
NMDA
1
1 + a [Mg
2+
] e
Vpost/V
block
71
Figures/StEd92.f10.eps
Figure 47: Different time scales and different voltage dependence of synaptic currents
evoked by AMPA and by NMDA receptors. CNQX blocks AMPA receptors, APV blocks
NMDA receptors. [57]
Figure 48: Dependence of NMDA-conductance on [Mg
2+
] and voltage [10].
Notes:
NMDA-mediated excitatory currents require
pre-synaptic action potential
post-synaptic depolarization
NMDA-receptors can function as coincidence detectors
NMDA-receptors are permeable to Ca
2+
72
Ca
2+
plays important role in modifying the strength of synapses (synaptic plastic-
ity) and learning
NMDA-receptors are likely a neural substrate for Hebbian learning : neurons
that re together are wired together (see Sec.7.3.5)
5.2.2 Short-Term Plasticity
The magnitude and probability of synaptically evoked currents can depend on the his-
tory of the activity at the synapes
short-term plasticity: modications that last from milliseconds to seconds
long-term plasticity: modications that last much longer, possibly as long as the
animal lives
Many mechanisms can contribute to plasticity. They can be pre-synaptic and post-
synaptic.
We focus here on the short-term plasticity through changes in the probability of release
P
rel
(cf. (27)).
depression: P
rel
is reduced by action potential
e.g. the pool of releasable vesicles can become depleted if it does not become
replenished sufciently fast after the release
facilitation: P
rel
is increased by action potential
e.g. the presynaptic [Ca
2+
] and with it the release probability can rise through
multiple action potentials
Simple phenomenological model:
Between action potentials P
rel
relaxes to its usual value P
rel
dP
rel
dt
= P
P
rel
(36)
relaxation to the steady-state value could be that of the number of releasable
vesicles (vesicle fusion, transmitter re-uptake) or of [Ca
2+
] (action of Ca
2+
pumps).
Depression:
P
rel
f
D
P
rel
, f
D
< 1 at t
+
spike
immediately after spike
Facilitation:
P
rel
P
rel
+ f
F
(1 P
rel
) at t
+
spike
immediately after spike
73
Consider the steady state that is reached when transmitting a regular train of action
potentials with period T
relaxation between spikes must be balanced by increments/decrements at spikes
Facilitation:
After a spike at t = 0
+
P
rel
gets incremented
P
rel
(0
+
) = P
rel
(0
) + f
F
_
1 P
rel
(0
)
_
(37)
Then P
rel
relaxes according to (36)
P
rel
(T
) = P
+
_
P
rel
(0
+
) P
_
e
T/
rel
For periodic solution we need
P
rel
(T
) = P
rel
(0
)
Inserting P
rel
(0
+
) into P
rel
(T
) yields
P
+
_
P
rel
(0
) + f
F
_
1 P
rel
(0
)
_
P
_
e
T/
rel
= P
rel
(0
)
Thus
P
rel
(0
) = P
1
_
1
f
F
P
_
e
T/
rel
1 (1 f
F
) e
T/
rel
P
For low ring rate (T

rel
) the exponential is small and the denominator can be
expanded
P
rel
(0
) = P
+ f
F
(1 P
) e
T/
rel
+O(e
2T/
rel
)
Facilitation decays between spikes, P
rel
only slightly increased
For high ring rate (T
rel
)
P
rel
(0
)
f
F
+ (P
f
F
)
T
rel
f
F
+ (1 f
F
)
T
rel
= 1
1 P
f
F
T
rel
+O
_
_
T
rel
_
2
_
Facilitation raises P
rel
to values close 1
Depression:
(37) is replaced by
P
rel
(0
+
) = f
D
P
rel
(0
)
yielding for the periodicity condition
P
+
_
f
D
P
rel
(0
) P
_
e
T/
rel
= P
rel
(0
)
74
P
rel
(0
) = P
1 e
T/
rel
1 f
D
e
T/
rel
P
For low ring rate (T

rel
)
P
rel
(0
) = P
_
1 (1 + f
D
)e
T/
rel
_
+O
_
e
2T/
rel
_
For high ring rate (T
rel
)
P
rel
(0
) = P
1
1 f
D
T
rel
+O
_
_
T
rel
_
2
_
depression severely suppresses release.
How many spikes are actually transmitted by the synapse?
Per period of the presynaptic train of action potentials 1 P
rel
spikes are transmitted
the transmission rate of spikes is given by
R = P
rel
1
T
For high ring rates the transmission rate of the depressing synapse saturates
R =
P
1 f
D
1
rel
+ . . .
Note:
For high ring rates this depressing synapse does not transmit any information
about the value of the input ring rate 1/T.
R is the output ring rate of an array of parallel synapses
For high ring rates consider step change in input ring rate r
1
T
r r + r
Before the step
P
rel
=
P
(1 f
D
)
rel
1
r
R
=
P
(1 f
D
)
rel
Immediately after the step P
rel
is unchanged
R
+
= P
rel
(r + r) = R
r + r
r
and
R
+
R
=
r
r
where for long plateaus R
is close to the asymptotic value R =

P
1f
D
1
rel
.
Thus
75
For high ring rates this depressing synapse responds only to changes in input
ring rate
The change in output ring rate is given by the relative change of the input ring
rate, independent of its absolute value (scale-invariant)
synapse performs a computation on the input
Figure 49: Response of depressing synapse to varying input ring rates [10]
Note:
for natural conditions irregular input spiking is a better approximation
for randomly distributed spike times (Poisson train) the analogous computation
yields an equivalent result (cf. Dayan & Abbott Ch.5.8)
depression and facilitation can occur at the same synapse on different time scales.
For a more detailed model of short-term plasticity see [59, 43].
6 Firing Rates
17
Purpose of neurons is to transmit information
from sensory organs
17
DA 1.2,1.4
76
to cortex (merging of sensory information, storage in memory, comparison with
memory, decisions)
to motor control
to muscles
In most parts of the brain the information is transmitted via action potentials. How do
they transmit the information?
the shape of the action potential is stereo-typical
presumably it contains little information
can replace it by an elementary event
timing of the spike
rate at which spikes are red
correlations between spikes of pairs of neurons
synchrony
given time difference between spikes
...
correlations between spikes of neuron triplets
ring sequences
...
Figures/visual_receptive_fields_HW.ps
Figure 50: Receptive elds of neurons in primary visual cortex. Hubel & Wiesel video.
Long movie http://people.esam.northwestern.edu/riecke/Vorlesungen/Comp_
Neuro/Notes/Movies/VisualCortex.mov
There are a number of short movies on www.youtube.com.
77
Figures/BaBu08.f2A.ps
Figure 51: Spike responses to odor inhalation of a mouse. Bars in horizontal lines
denote spikes in different trials under the same conditions [1].
Figure 52: Neuron in primary motor cortex. a) raster plots during reaching movements
in different directions. b) average ring rate as a function of the reaching direction [10].
Notes:
The rate at which spikes are red seems to contain signicant information.
There is a signicant variation in the individual spike times across repeated trials.
To assess the signicance of correlations between spikes of different neurons
one needs to measure multiple/many neurons simultaneously
78
barn owl: in echo-location the coincidence of spikes from neurons from the
left and the right ear provides the information about the location of the object
in general, role of spike synchrony and correlations is still being debated
very actively.
oscillations and rhythms are observed in coarse measures (e.g. local eld
potential, EEG) of many large ensembles of neurons: function?
difculties:
data acquisition
voltage-sensitive dies and Ca
2+
-imaging are slow
multi-electrode arrays require spike sorting (neuron identication)
...
data analysis: binomial explosion of correlations (which neurons are
read?)
Firing Rate
There are different ways to dene a ring rate of a neuron.
Consider the neural reponse function
(t) =
n
i=1
(t t
i
)
where t
i
, i = 1 . . . n, are the spike times and (t) is the Dirac -function.
Spike-count rate
r =
n
T
=
1
T
_
T
0
(t)dt
counts the total number of spikes over the whole duration T of the trial
does not give any temporal resolution
for uncorrelated spikes the variability in r depends on the total number of
spikes in the interval:
variability becomes small for rT 1
Time-dependent ring rate
to get temporal resolution need an integration window t that is short com-
pared to the time scale over which the ring rate changes
uctuations in ring rate are small if
if the ring rate is high and
if it varies only on a slow time scale allowing t = O():
the number of spikes during must be large: r 1.
79
if t is very short it will contain very few spikes
to get a reliable measure average over multiple trials of the same experiment
r(t) =
_
1
t
_
t+t
t
(t
)dt
_
trials
=
1
t
_
t+t
t
(t
)) dt
where . . .)
trials
denotes the average over trials
the size of t limits the temporal resolution
t can be chosen small if many trials are available
in the limit t 0 at most 1 spike occurs during t
r(t)t is the fraction of trials during which a spike occurred in the interval
[t, t + t]
thus, r(t)t is the probability of ring a spike during that interval
Figure 53: Spike train and ring rates obtained by binning in xed windows (b), in sliding
rectangular window (c), in sliding Gaussian window (d), causal window (e). [10].
Population ring rate
in many situations an animal does not have the privilege to average over
many trials
if multiple, statistically independent neurons convey the same information
the animal can use a population average . . .)
pop
summing up the spikes of
those neurons during an interval t
r
pop
(t) =
1
t
_
t+t
t
(t
))
pop
dt
6.1 Poisson Spike Trains

Averaging over trials relies on the variability of the spiking response across trials (cf.Fig.51).
If the spikes always occurred at precisely the same times averaging would not make
any smoother.
80
Consider random spike trains
The trains are thought of as being generated by a stochastic process. We are inter-
ested in
P[t
1
, t
2
, . . . , t
n
] = Probability of spikes occuring at the specied times
Spikes are events: a stochastic process that generates a sequence of discrete events
is called a point process.
Events could depend on all previous events. Simpler processes are
renewal process: spiking probability depends only on the previous spike time
Poisson process: spiking probability is completely independent of all previous
spike times
Note:
if only the ring rate is relevant a suitable model for the spike train is a Poisson
process
constant ring rate r: homogeneous Poisson process
time-varying ring rate r(t): inhomogeneous Poisson process
Homogeneous Poisson process:
Determine the probability that during interval T exactly n spikes are red (at arbitrary
times): P
T
[n]
Divide T into M = T/t bins with M eventually
P
T
[n] = lim
t0
(rt)
n
. .
spike in n bins
(1 rt)
Mn
. .
no spikes in Mn bins
M!
(M n!) n!
. .
choose n from M
Rewrite second factor using = rt
lim
0
(1 + )
T
t
n
= lim
0
(1 + )
Tr
n
= lim
0
_
_
(1 + )
1
_
rT
(1 + )
n
_
= e
rT
Using Stirling formula for ln N for large N
ln N! = N ln N N +O(ln N)
for large M at xed n we get then
ln
M!
(M n)!
M ln M M
_
(M n) ln
_
M (1
n
M
)
__
+ M n
= M ln M (M n) ln M +O(n)
nln M
81
yields the limit
M!
(M n)!
M
n
for M
Thus
P
T
[n] = lim
t0
(rt)
n
e
rT
1
n!
_
T
t
_
n
Thus, the Poisson process is described by the Poisson distribution
P
T
[n] =
1
n!
(rT)
n
e
rT
Figure 54: a) Probability of having exactly n spikes as a function of rescaled ring rate
rT. b) For not too small rT the Poisson distribution is well approximated by a Gaussian
distribution with a variance that is equal to its mean [10].
Compute mean and variance
n) =
n=0
nP
T
[n]
2
n
n
2
) n)
2
=
n=0
(n n))
2
P
T
[n]
To compute these it is useful to introduce what is called the moment-generating function
g() =
n=0
e
n
P
T
[n]
It allows to calculate the moments
n
k
_
very easily,
n
k
) =
d
k
d
k
g()
=0
Now
g() =
n=0
e
n
1
n!
(rT)
n
e
rT
= e
rT
n=0
1
n!
(e
rT)
n
= e
rT
e
e
rT
82
using the Taylor series of the exponential
e
x
=
n=0
1
n!
x
n
dg
d
=0
= e
rT
e
rTe
rTe
=0
= rT
d
2
g
d
2
=0
= e
rT
e
rTe
(rTe
)
2
+ e
rT
e
rTe
rTe
=0
= (rT)
2
+ rT
Thus:
n) = rT
2
n
= rT
Note:
the Fano factor is dened as
F =

2
n
n)
for the Poisson distribution the mean and the variance are equal, independent of
T
F = 1
Figures/ShNe98.f1.ps
Figure 55: High ring variability of cortical neurons (area MT of an alert monkey) reect-
ing balanced excitatory and inhibitory input from a large number of neurons. Variance
grows (super-linearly) with mean ring rate [53].
83
Figures/FiCh07.f4_top.ps
Figure 56: Low ring variability of retinal ganglion cells (sub-Poisson). [17]
Figures/SrYo06.f11.ps
Figure 57: Mean ring rate and variability of ring rate of peripheral (SA1 and RA) and
cortical neurons (3b and 1) in response to tactile stimulation of macaque monkey [56].
6.2 What Makes a Neuron Fire?
18
Tuning curves like Fig.52 show how a neuron responds for a range of selected stimuli,
e.g. different orientations of a given grating pattern.
Can one assess what kind of stimulus would elicit maximal response from the neuron
without restricting the stimulus to a specic type (e.g. grating pattern with a given
wavelength)?
6.2.1 Spike-Triggered Average
Consider a neuron responding to a stimulus s(t) with a spike train t
i
, i = 1 . . . n
18
DA 1.3 +2.2
84
Figures/DaAb01.f1.9_left.ps
Figure 58: Time-dependent stimulus and resulting spiking activity of electrosensory
neuron in weakly electric sh [10].
The stimulus provides a time-varying input current to the neuron.
The neuron accumulates the input and spikes when the threshold voltage is sur-
passed (cf. Fig.45).
what is the average stimulus history that leads to a spike?
For each spike at time t
i
in the experimental interval [0, T] measure the stimulus at a
time t
i
contributing to the spike
To get spike-triggered average stimulus average the spike trains for many presentations
(trials) of the same stimulus s(t)
C() =
_
1
n
n
i=1
s(t
i
)
_
Figure 59: Procedure to compute the spike-triggered average stimulus [10].
Assuming small uctuations in the number of spikes across the trials ( 1/n)
1/2
) to
pull the factor
1
n
out of the average over trials and using the neural response function
=
n
i=1
(t t
i
)
85
C()
1
n)
_
n
i=1
s(t
i
)
. .
R
T
0
(t
t
i
) s(t
)dt
_
=
=
1
n)
__
T
0
(t
) s(t
) dt
_
=
1
n)
_
T
0
(t
)) s(t
) dt
=
1
n)
_
T
0
r(t
) s(t
) dt
Consider the correlation function between stimulus and ring rate

Q
rs
() =
1
T
_
T
0
r(t
) s(t
+ ) dt
Using n) = rT we get
C() =
Q
rs
()
r
Note:
the spike-triggered average stimulus is proportional to the reverse correlation (be-
cause of the minus-sign)
C() depends on the stimulus set used:
to cover all possible inputs: use random stimulus ensemble (e.g. white
noise)
statistically stationary stimulus
the response properties of a neuron could depend on its mean ring
rate (e.g. by facilitation or depression); statistically nonstationary stimuli
would average over different states of the neuron.
. . .) is the average over many realizations from that ensemble
To characterize stimulus ensemble: auto-correlation function of s(t)
Q
ss
() =
1
T
_
T
0
s(t
)s(t
+ )dt
Use white noise without any correlations between different times

Q
ss
() =
2
s
()
86
Figure 60: Spike-triggered average stimulus for a pyramidal neuron in the electrosen-
sory lateral line-lobe, which is the rst central nucleus of the electrosensory pathway in
weakly electric sh [20] (in [10]).
6.2.2 Linear Response Kernel and Optimal Stimulus
The spike-triggered average gives information about the stimulus that contributed to
triggering a spike, i.e. to information about the stimulus given a spike.
Can one conversely estimate/predict how a neuron will respond to a given stimulus?
Assume a linear response of the neuron
r
est
(t) = r
0
+
_

0
D()s(t )d (38)
r
0
is the spontaneous ring rate for s = 0
effect of mean stimulus s
0

_
T
0
s(t
)dt
= 0 can be absorbed in r
0
:
assume s
0
= 0
the kernel D() measures the degree to which the stimulus at the earlier time
t contributes to the ring rate at time t
Notes:
neurons are nonlinear
assumption of linear response is non-trivial
approximation expected to be good if the spontaneous ring rate r
0
is only
slightly modulated by the stimulus
87
one can go beyond the linear response using a Volterra expansion (or a
Wiener expansion), which is an analog of a Taylor expansion)
ring-rate framework requires that ring rate and stimulus vary sufciently slowly
Determine D() by minimizing the error in estimating the ring rate r
est
compared to
the measured ring rate r
E =
1
T
_
T
0
(r
est
(t
) r(t
))
2
dt
Using the linear response assumption

ED() =
1
T
_
T
0
_
r
0
+
_

0
D()s(t
)d r(t
)
_
2
dt
Goal:
nd D() such that E is minimal, i.e. arbitrary small changes () of D() do not
change E (cf.
df
dx
= 0 at the minimum of the function f(x)), i.e. require
E ED() + () ED() = 0 +O(
2
)
E =
1
T
_
T
0
_
r
0
r(t
) +
_

0
D()s(t
)d +
_

0
(
)s(t
)d
_
2
dt
1
T
_
T
0
_
r
0
r(t
) +
_

0
D()s(t
)d
_
2
dt
= 2
1
T
_
T
0
__
r
0
r(t
) +
_

0
D()s(t
)d
__

0
(
)s(t
)d
_
dt
+O(
2
)
=
2
T
_

0
(
)
__
T
0
s(t
)
_
r
0
r(t
) +
_

0
D()
_

0
s(t
)d
_
dt
_
d
+O(
2
)
For D() to be optimal the term proportional to must vanish for arbitrary (
). This
requires the term inside to vanish
19
_
T
0
(r
0
r(t
)) s(t
) dt
+
_

0
D()
__
T
0
s(t
)s(t
)dt
_
d = 0
Thus
_

0
D()
_
1
T
_
T
0
s(t
)s(t
) dt
_
d =
1
T
_
T
0
(r(t
) r
0
) s(t
) dt
19
In Appendix B this is derived using functional derivatives.
88
Shifting integration variable t
= t
1
T
_
T
0
s(t
)s(t
) dt
=
1
T
_
T
s(t
)s(t
) dt
= Q
ss
(
)
Using
_
T
0
s(t
)dt
= 0 and interchanging and
we get
_

0
D(
)Q
ss
(
)d
= Q
rs
()
For a white-noise stimulus Q
ss
(
) =
2
s
(
)
D() =
Q
rs
()
2
s
=
r)
2
s
C() (39)
Notes:
for white-noise stimuli the spike-triggered average gives the optimal kernel to
predict the ring rate
r
est
is obtained by projecting the stimulus onto D (see (38))
via (39) the maximal response
20
is obtained by a stimulus with the shape of
the spike-triggered average stimulus
Figure 61: Response prediction using the optimal kernel (39) for the H1 neuron in the
visual system of the y. The H1 is sensitive to horizontal motion of large portions of the
visual eld.
a) image velocity used as input. b) two spike sequences obtained in response to that
input. c) measured (dashed) and predicted (solid) ring rate [10].
20
at xed energy. See DA Appendix 2.9B
89
6.2.3 Receptive Fields: Visual System
Neurons may respond not only to a single stimulus, but to many different stimuli. The
response depends then not only on the temporal characteristics of the stimulus but on
the combination of the various stimuli.
Example: spatial information in the visual system
Sketch of initial stages:
retina:
rods and cones, temporal processing by bipolar cells, lateral connections via hor-
izontal cells
optical chiasm: left and right elds of vision combined for both eyes
lateral geniculate nucleus in the thalamus: some integration with other informa-
tion
primary visual cortex V1
...
a)
Figures/DaAb01.f2.4a.ps
b)
Figures/KaSc00.f26-6.ps
Figure 62: Structure of the retina. a) Drawing by Cajal (1911). b) Sketch of cone and
rod pathway in the retina [38].
Figure 63: Visual pathway from the retina to primary visual cortex [10].
90
Each stage - even the retina - involves extensive neural computation.
Notes:
The mapping between the stages is retinotopic: neighboring neurons respond to
neighboring pixels in the retina.
In each stage there are different types of neurons with different response proper-
ties
The neural response is characterized by the receptive eld D(x, y, t), i.e. the spatio-
temporal stimulus that elicits optimal response.
Focus here on neuron types for which the linear response theory appears to be ade-
quate.
In analogy to (39) one then has
D(x, y, t) =
r)
2
s
C(x, y, t)
with the spatio-temporal spike-triggered averaged stimulus C(x, y, t).
Retinal ganglion cells/LGN cells:
The receptive eld of many cells has a center-surround structure in space
D(x, y) =
A
center
2
2
center
e
x
2
+y
2
2
2
center

A
surround
2
2
surround
e
x
2
+y
2
2
2
surround
on-center cell: A
center
> 0, A
surround
> 0
off-center cell: A
center
< 0, A
surround
< 0
a)
Figures/DeOh95.f1.ps
b)
Figures/DeOh95.f3.ps
Figure 64: a) Center-surround receptive eld of an on-center LGN-neuron in cat. b)
Temporal structure of receptive eld of LGN-neurons in cat. Horizontal axis: x, vertical
axis: t, green: excitatory, red: inhibitory. A) Excitatory response at light onset. B)
Excitatory response delayed, but stronger than immediate inhibitory response. [11]
91
Typically the receptive elds have also a temporal structure,
A
center,surround
= A
center,surround
(t),
as shown in Fig.64b.
6.3 LNP Cascade Models
21
So far:
To characterize response properties of neuronswe used spike-triggered average
and the linear response.
Linear response only expected to be successful if the ring rate varies by small
amounts about a non-zero mean rate.
Linear response does not capture
spiking threshold
saturation of ring rate
Rate model does not include individual spikes, which can occur probabilistically
(Poisson spiking)
Beyond linear response:
Volterra series: t the input-ouput relation with higher-order polynomial in the stimulus
s
r
est
(t) = r
0
+
_

0
D
1
()s(t )d +
_

0
_

0
D
2
(,
)s(t )s(t
)dd
+ . . .
The kernels D
i
(,
, . . .) are to be determined to give a good approximation for a wide

range of stimuli: use white-noise stimulus with
Q
ss
(,
) =
2
s
(
)
Then we had
D
1
() =
r)
2
s
C()
with C() the spike-triggered average
C() =
_
n
i=1
s(t
i
)
_
trials
The higher-order kernels can be determined via higher powers of the stimulus with the
output (spiking).
Notes:
21
cf. [6]
92
as in tting data points with polynomials the coefcients of the polynomials change
when the order of the polynomial changes (the t is not a Taylor expansion)
The interdependence between the various kernels can be removed by switching to a
Wiener series[50]
r
est
(t) = g
0
+
_

0
g
1
()s(t)d+
__

0
_

0
g
2
(,
)s(t )s(t
)dd
2
s
_

0
g
2
(, )d
_
+. . .
with the lters given by
g
0
= r(t))
g
1
() =
1
2
s
r(t)s(t )) = D
1
()
g
2
(,
) =
1
4
s
r(t)s(t )s(t
))
. . .
But:
Sigmoidal dependence of the ring rate on the input is poorly approximated by
polynomials
Volterra and Wiener series require high orders to give good results im-
practicable, would require huge amounts of data
Can one make use of the preferred stimulus (spike-triggered average) to construct a
simpler nonlinear model?
Cascade Model:
spike-triggered average generates from the input history a scalar quantity r
lin
(t)
r
lin
(t) is the input to a static nonlinearity F which gives a spiking probability
spikes are generated from a Poisson distribution based on the time-dependent
spiking probability F(r
lin
(t))
r
est
(t) = F
__

0
g()s(t )d
_
How to determine F and g?
This model has only a single lter, g()
Require:
93
the rst Wiener lter of the true data should agree with the rst Wiener lter of the
estimated data
r(t) s(t )) = r(t)
est
s(t ))
Note:
since the Wiener lters are independent of each other this condition does not
require any knowledge of higher Wiener or Volterra lters
Slightly easier to consider discrete time series (binned data) s
i
, i = 1 . . . N, for the N
past values of the stimulus, which result in the current value of the output r and its
estimate r
est
. The model then reads
r
est
= F (g s)
with the discrete lter g.
Consider rst spherically symmetric probability distributions for the stimuli
P(s) = P([s[)
r
est
s) =
_
r
est
sP(s)ds
1
. . . ds
N
=
_

P([s[)
. .
spherical
_
F (g s) s + F
_
g s
+
_
s
+
ds
1
. . . ds
N
where s
+
is chosen such that g s = g s
+
, i.e.
s
+
= 2
s g
[g[
2
g s
and the sum
_
is adjusted to avoid double counting.

Then
r
est
s) =
_
_
_

P([s[)
. .
spherical
F (g s)
_
_
_
s +s
+
_
. .
g
ds
1
. . . ds
N
g
Since we require r
est
s) = r s) the lter g satises
g = r s)
The nonlinearity is then obtained by plotting r vs the linear lter output g s.
94
For Gaussian, not necessarily spherical, distributions one can show the equivalent
result (Bussgang theorem[4, 10])
r
est
s) =
_
r
est
(s) se
1
2
s
t
C
1
s
ds
1
. . . ds
N
=
_
r
est
(s) C
s
_
e
1
2
s
t
C
1
s
_
ds
1
. . . ds
N
=
..
integration by parts
C
_
e
1
2
s
t
C
1
s
s
r
est
(s)ds
1
. . . ds
N
= C
s
r
est
(s))
= C
s
F (g s))
= Cg F
(g s))
Thus, we have
g C
1
r
est
s) = C
1
r s)
Notes:
for Gaussian white noise: the lter g is proportional to the spike-triggered average
for correlated Gaussian noise: the spike-triggered average needs in principle to
be decorrelated.
however: when the data sets are not large enough the average r s) often con-
tains still noise. The decorrelation also amplies this noise and may not improve
the lter.
the LN-model approach can also be used to characterize the membrane voltage
response of neuron in response to a stimulus
95
a)
Figures/ZaBo05.f3a1.ps
b)
Figures/ZaBo05.f3a2.ps
Figure 65: LN-model for OFF Y-type retinal ganglion cell for low and high contrast of
the stimulus. a) Membrane voltage response. b) Spike rate response. [62]
Figures/ZaBo05.f2b.ps
Figure 66: Comparison of the voltage prediction of LN-model for the OFF Y-type retinal
ganglion cell shown in Fig.65 for low and high contrast of the stimulus (stimulus in top
row) [62]..
LN-models can also be used to characterize quantitatively certain aspects of the neu-
ronal response, like the gain.
Typically the gain depends on the contrast of the stimulus.
For the cell shown in Fig.65 the amplitude of the linear lter can be rescaled such that
the nonlinearity is the same for both contrast conditions
r
est
= F (
low,high
g s) with
high
<
low
Note:
96
The gain is reduced for high contrast. This helps avoiding saturation of the re-
sponse.
a)
Figures/ZaBo05.f3b1.ps
b)
Figures/ZaBo05.f3b2.ps
Figure 67: LN-model for the same OFF Y-type retinal ganglion cell as in Fig.65. After
rescaling of the lter the nonlinearities become the same for low and high contrast of
the stimulus (grey=low contrast, black=high contrast). a) Membrane voltage response.
For high contrast gain is reduced to 0.83 of the gain for low contrast. b) Spike rate
response. Gain reduced to 0.64 [62]
97
a)
Figures/ZaBo05.f3c1.ps
b)
Figures/ZaBo05.f3c2.ps
Figure 68: LN-model for ON Y-type retinal ganglion cell. After rescaling of the lter
the nonlinearities become the same for low and high contrast of the stimulus (grey=low
contrast, black=high contrast). a) Membrane voltage response. Gain redcued to 0.74.
b) Spike rate response. Gain reduced to 0.54. [62]
The description of some cells require multiple lters and nonlinearities.
Shifting backgrounds in visual scenes can modify the response of ganglion cells to
inputs in their receptive eld. This could be relevant during saccades.
98
a)
Figures/GeDe07.f1a.eps
b)
Figures/GeDe07.f2.eps
Figure 69: Input outside the classical receptive eld can modify a ganglion cells re-
sponse. a) grating exciting the periphery is shifted every 9s. b) The response of this
ganglion cell is transiently changed from OFF-center to ON-center briey after a pe-
ripheral shift. [21]
Figures/GeDe07.f4.eps
Figure 70: Extraction of two temporal lters using a principal component analysis (PCA)
of the spike-triggered covariance matrix [21]
99
7 Neural Networks: Rate Models
The power of the brain stems to a large extent from the interaction between different
neurons that communicate with each other.
Characteristics of most functional networks in the brain
comprised of many neurons
each neuron receives many inputs
neurons are connected only to a fraction of the other neurons:
sparse connectivity
Modeling with detailed spiking neurons very challenging
time scales:
channel opening/closing: < 1ms
action potential: 2ms
...
period of -rhythm 100ms
persistent activity in working memory > 1s
...
multiple variables V , m, n, h,...
many compartments for branching dendrites, axons
complex nonlinearities
many parameters, of which many are poorly known
For large networks ring-rate models are useful
do not resolve short time scales of action potential
fewer variables
simpler nonlinearities
few parameters
cannot capture aspects of spike synchrony
Note:
In the Blue Brain Project a single neocortical column consisting of 10,000 neurons
is being simulated with morphological details of neurons retained. It uses a 8192-
processor Blue Gene computer.
http://bluebrain.epfl.ch/page18699.html
100
7.1 Firing-Rate Models
22
In a spiking model the activity of the network is characterized by the neural response
function
j
(t) =
i
(t t
i
) of each neuron j.
In a ring-rate model
j
(t) is replaced by its trial-average, the ring rate r(t) = (t)).
Necessary for ring-rate approximation:
network behavior does not depend on specic timing of the spikes
The dynamics of a neuron depend on its total input.
The response of a neuron is insensitive to individual spike times of its inputs if
large number N of uncorrelated spikes contribute to the response
many synapses
slow synapses
slow membrane properties (
m
= RC large)
(cf. Fig.45b)
input uctuations due to individual spikes
N, mean input N.
We need
synaptic input current I
s
to neuron (measured at soma) due to the changes in all
the synaptic conductances as a function of pre-synaptic ring rate u
output ring rate v as a function of I
s
Synaptic Current I
s
Consider a neuron with N
u
synaptic inputs
I
s
(t) =
Nu
i=1
w
i
t
j
<t
K
s
(t t
(i)
j
) =
Nu
i=1
w
i
_
t
K
s
(t )
i
()d
with
synaptic kernel K
s
effect at the soma of the synaptic input (spike times given by
i
(t))
includes synaptic dynamics and propagation along dendritic cable
22
DA 7.1
101
assume K
s
equal for all synapse on the neuron
need not be the case: the impact of a synapse may well depend on its
distance from the soma
normalized:
_
t
K
s
(t )d =
_
0
K
s
(t
)dt
= 1
w
i
strength of synapse i
assume no short-term plasticity, w
i
does not change from spike to spike,
does not depend on inter-spike interval
w
i
> 0: excitatory synapse
w
i
< 0: inhibitory synapse
assume synaptic inputs sum up linearly
however, active dendrites can sum up their input sub-linearly or super-linearly
Firing-rate assumption
the synaptic current I
s
can be approximated by its trial average
I
s
I
s
(t)) =
Nu
i=1
w
i
_
t
K
s
(t )
i
())
. .
u
i
()
d (40)
Assume a single-exponential kernel (fast rise time)
K
s
(t) =
1
s
e
t/s
H(t)
with the Heaviside function
H(t) =
_
1 for t > 0
0 for t < 0
H(t) = 1 for t > 0 and H(t) = 0 for t < 0.
Then one can express (40) in terms of a differential equation
Using the derivative of (40)
dI
s
dt
=
Nu
i=1
w
i
_
1
s
u
i
(t)
_
t
e
(t)/s
u
i
()d
_
yields
s
dI
s
dt
= I
s
+
Nu
i=1
w
i
u
i
I
s
+w u (41)
Firing Rate
For steady input current I
s
the output ring rate of the neuron is given by
v
= F(I
s
)
Typical choices
102
threshold linear function
F(I
s
) = [I
s
I
threshold
]
+

_
I
s
I
threshold
for I
s
I
threshold
> 0
0 for I
s
I
threshold
0
sigmoid function for smoothness and for saturation
For time-dependent I
s
(t) a simple assumption is that v relaxes exponentially to v
on
a time scale
r
.
This yields the coupled system
r
dv
dt
= v + F (I
s
(t)) (42)
s
dI
s
dt
= I
s
+w u (43)
i)
r

s
: v relaxes very quickly to v
v = F(I
s
)
s
dI
s
dt
= I
s
+w u (44)
ii)
r

s
: I
s
(t) relaxes quickly
r
dv
dt
= v + F(w u) (45)
Note:
if the neurons providing the input are in turn described by a ring-rate model one
gets for (44)
s
dI
s
dt
= I
s
+
Nu
i=1
w
i
F (I
si
(t)) (46)
with I
si
the input current of the input neurons.
comparing (46) with (45) shows:
the two models differ then in the position of the nonlinearity F.
the position of the nonlinearity, which is non-negative, is also related to the fact
that I
s
can be negative but v cannot.
7.2 Feed-Forward Networks
The brain can roughly be thought of as a modularly structured network
103
feed-forward connections from one module to the next (e.g. sensory thalamus
initial processing in cortex ... higher level processing .... motor control
motor output)
recurrent connections within modules
top-down (feed-back) input: most regions that receive feed-forward input from
lower levels also project back to those levels
General form within the frame-work of ring-rate models:
Feed-forward
dv
dt
= v +F(Wu) (47)
with W the matrix of synaptic weights of the inputs.
In components
dv
i
dt
= v
i
+ F
_
j
W
ij
u
j
_
Recurrent
dv
dt
= v +F(Wu +Mv) (48)
with M the matrix of synaptic weights of the recurrent connections.
Consider rst feed-forward networks. Although simpler to treat than recurrent ones,
they can still perform very powerful computations.
Note:
without input u feedforward networks are silent, v = 0
(unless the individual, uncoupled neuron re spontaneously, F(0) ,= 0.)
7.2.1 Hubel-Wiesel Model for Visual Orientation Tuning
Observations of visual receptive elds:
retinal ganglion cells and cells in lateral geniculate nucleus (LGN) in thalamus
have center-surround receptive elds (cf.Sec.6.2.3)
cells in the input layer (layer IV) of visual cortex V1 have more complex receptive
elds:
for many cells the receptive eld is based on an oriented stripe pattern
different orientations
different wavelengths
104
stationary or moving at a certain speed
Figures/HuWi62.f2.ps
Figure 71: Receptive elds: A) ON-center in LGN, B) OFF-center in LGN, C)-G) simple-
cell receptive elds in V1 [32]
How can orientation tuning in V1 arise from center-surround tuning of the ganglion cells
in the retina?
Figure 72: Feed-forward network for simple cell [32].
Simple cell responds most strongly to
bright (dark) stripe in the center
dark (bright) stripes next to the bright stripe
selective for orientation, wavelength and position.
Complex cells are not selective for position (phase) of the stripes, only for their orien-
tation
105
Figure 73: Receptive eld of a complex cell. A)-D) same orientation but different lo-
cation of the bar within the receptive eld. E)-G) different orientations. H) bar with
orientation as in A)-D) moved rapidly across the receptive eld. The horizontal bars
above the record indicate when the stimulus is on [32].
Figure 74: Possible connectivity of a complex cell: input from multiple simple cells with
different preferred locations of the bar [32].
For any phase (position) of the grating some simple cell in the receptive eld of the
complex cell responds, while the others are quiet: no cancellation of the input from the
active cells
Thus:
the complex cell is not sensitive to the position of the stripe:
it encodes the generalized information some stripe with this orientation
Note:
for linearly responding cells with high spontaneous ring rate the reduced input
from the simple cells that are responding only weakly would compensate the
enhanced input from the strongly responding simple cells.
106
Notes:
it has been argued that the feed-forward architecture is insufcient to explain
sharpness of selectivity
contrast invariance of the selectivity
selectivity is as sharp for high-contrast pattern as for low-contrast pattern (cf.
ice-berg effect)
recurrent network architectures have been proposed [2]
it seems, however, that if data and models are analyzed in more detail the feed-
forward network may well be sufcient [18]
7.2.2 Neural Coordinate Transformations
23
Receptive elds of LGN cells and of cells in V1 are in terms of positions on the retina
(retinotopic map), which depends on eye and head direction (gaze).
To manipulate an object its position needs to be known in a body-centered coordinate
system.
brain needs to transform from the retinal to the body coordinate system
One-dimensional example:
s = angle in retina coordinates g = angle of the gaze s + g = body-centered angle
Figure 75: Coordinate transformation between body-centered coordinates and retina-
centered coordinates.
23
DA 7.3
107
Figures/GrHu97.f13.ps
Figure 76: Response of a bimodal (visual-tactile) neuron in the ventral pre-motor cortex
of a macaque monkey. Its tactile receptive eld is near the left cheek. The visual
response to an object approaching along paths marked I-V does not depend on the
eye position (A1, B1, C1) (or the arm position (B2)), but it does depend on the head
position (B3), i.e. it depends on the position of the object relative to its tactile receptive
eld (similar data for tactile eld on the arm) [24].
a)
b)
Figure 77: Tuning curves of face+visual bimodal neurons in ventral pre-motor cortex
are in head-centered coordinates. a) The receptive eld is shifted when the head is
turned. b) The receptive eld is not shifted if the gaze, but not the head is rotated [24].
Figures/DaAb01.f7.6a.ps
Figure 78: Tuning curve of neurons in the posterior parietal cortex. The preferred
location is given in retinal coordinates, but the magnitude of the response, i.e. the
height of the tuning curve, is modulated by the gaze direction [3].
Consider two model networks at the interface between sensory function and motor
control [51]
108
i) Sensory Network (posterior parietal cortex PP):
The response of bimodal neurons in PP depends on the retinal coordinate s and the
gaze coordinate g. Different neurons have different preferred retinal angles and dif-
ferent characeristic gaze angle .
The neurons have tuning curves in retina coordinates, but the gaze enters via the gain,
u
(s, g) = f
u
(s , g ) = Ae
(s)
2
2
2
s
[M (g )]
M
+
Experimentally it is found that the gain varies rougly linearly with the gaze, here mod-
eled with a threshold and saturation.
Depending on the neuron the gain can increase or decrease with the gaze angle.
Assuming that the second network reads both types of neurons one can lump them
together into an effective population described by
f
u
(s , g ) = Ae
(s)
2
2
2
s
[M [g []
M
+
which effectively has a preferred gaze angle.
For simplicity assume a continuum of neurons characterized by and .
ii) Motor Control Network (ventral pre-motor cortex PMv):
Neurons in PMv receive input from the neurons of the sensory network
Can the input weights be chosen such that the neurons in the motor control network
have tuning curves in body-coordinates?
Steady-state ring rate for one such neuron
v
(s, g) = F (w u)
Assuming a continuum for the preferred angles and the sum can be replaced by
integrals and one obtains
v
(s, g) = F
_
_
_
_
+
_
+
w(, )
. .
connection from neuron labeled (,)
f
u
(s , g )
dd
_
_
_
where
d and
d are the numbers of neurons with preferred location and gaze in

the ranges [, +d] and [, +d], respectively. They are assumed to be independent
of and and can then be pulled out of the integral.
The tuning curve for this neuron should depend only on the angle in body coordinates
v
= v
(s + g)
To bring s and g together shift the variables
s = s + g
g =
109
i.e.
=
g =
+ g
then
v
(s, g) = F
_
_
+
_
+
w(
g,
+ g) f
u
(s + g
) d
_
To make the rst integrand independent of g choose
w(, ) = w( + )
then
v
(s, g) = F
_
_
+
_
+
w(
) f
u
(s + g
) d
_
Figures/coordinate_transform_1.eps Figures/coordinate_transform_2.eps
Figure 79: a) The neuron in PMv integrates the contributions fromall neurons in PP with
a weight w(, ) depending on their preferred location (, ). Their receptive elds are
represented by circles. b) For a different stimulus (s
, g
) the each of the PP-neurons

responds differently. But for each neuron at (, ) there is a neuron at (
) with
= + that responds like the neuron (, ) did for stimulus (s, g). The output of
the PMv-neuron depends only on s + g if the weight for the neuron (
) is equal to
the weight of neuron (, ), i.e., if the weight w(, ) depends only on + .
Note:
In [51] it is discussed that the required connectivity w( +) can be obtained in a
natural learning procedure: baby follows its own random hand motion [60].
7.3 Recurrent Networks
24
Recurrent networks are harder to analyze than feed-forward networks. Consider rst
linear neural dynamics, F(u) = u.
24
DA 7.4
110
7.3.1 Linear Recurrent Networks
Need to consider N coupled odes for the components of v
r
dv
dt
= v + Wu
..
input h
+Mv (49)
Solve in terms of the eigenvectors e
of M
Me
with eigenvalues
.
For general matrices M the eigenvalues and the eigenvectors can be complex, C.
Assume: M is symmetric, M
ij
= M
ji
eigenvalues
R
eigenvectors e
mutually orthogonal
e
(50)
all vectors v can be written as a linear combination of the eigenvectors
v(t) =
Nv
=1
v
(t) e
(51)
for any xed t with N
v
the number of components of v. v
(t) is the time-dependent

amplitude of the eigenmode e
.
Insert expansion (51) into (49)
dv
dt
e
(1
) v
+h
Use orthogonality (50) to project the equation onto e
r
dv
dt
= (1
) v
+e
h (52)
Now we have only a single ode. Its solution is given by
v
(t) =
e
h
1
+ A
1
r
t
Using the initial condition
v
(0) = e
v(0)
111
we get
A
= e
v(0)
e
h
1
The full solution is then given by the sum over all the eigenmodes with their respective
amplitudes
v(t) =
N
=1
_
e
h
1
+ A
1
r
t
_
e
(53)
The evolution of the components v
of v depends on the associated eigenvalue
For
< 1:
a steady state is reached at large times
v
(t) v
h
1
for t
for
< 0 the output component v
is reduced compared to the input com-

ponent
for 0 <
< 1 the output component v
is larger than the input component

e
h: selective amplication of components of h

for
1 the approach to the steady state becomes innitely slow
=

1
(critical slowing down)

For
> 1:
no steady state is reached
v
diverges exponentially within this linear approximation

nonlinearity F may lead to saturation
For
= 1:
it is easier to determine solution directly from (52)
r
dv
dt
= e
h
i.e.
v
(t) = e
h
t
no saturation, linear growth

112
more generally, for h = h(t)
v
(t) =
1
r
_
t
t
0
e
h(t
)dt
network functions as an temporal integrator of its input

Notes:
temporal integration is needed, e.g., to obtain positions from a velocities
to compensate head orientation on visual input need to obtain head orienta-
tion from rotation measured in vestibular system.
> 0 implies excitatory recurrent inputs: self-excitation
7.3.2 Oculomotor Control during Eye Saccades
In many animals eye position is xed for a period of time and then a sudden jump
(saccade) to a new position is made.
Even with closed eyes the eye position relative to the object remains xed during head
rotations.
How is the eye position controled?
Figure 80: Sketch indicating location of medulla and other parts of the brain stem [38].
Different types of neurons are involved (located in parts of the medulla at the end of
the spinal cord)
neurons ring tonically: ring rate proportional to eye position
neurons ring bursts during saccades
113
both neurons feed into the motor neurons:
motor neurons have pulse-step behavior
Note:
If position neurons fail (by lesion, e.g.), animal can still make a saccade to a new
eye position but the eye relaxes back to its rest position afterwards
a)
Figures/McFu92.f4.ps
b)
Figures/McFu92.f8.ps
Figure 81: a) Burst-position neuron (A,C) and position neuron (B,D) in the monkey
medulla (HE=horizontal eye position, HT=horizontal target position, FR=ring rate). b)
The duration of the bursts is strongly correlated with the duration of the saccades and
usually precedes them by 10ms (solid bars in inset for burst-position cells, open bars
for position cells) [47].
114
Figures/PaCr94.f2.ps
Figure 82: Neurons in goldsh medulla. coding for eye position (A,B) or eye velocity
(D). A) spontaneous eye movement. Firing rate (FR) correlated with left eye position
(LE). B) Sinusoidal head rotation (in the dark) with compensating eye movement. Firing
rate correlated with eye position, not eye velocity. D) Sinusoidal head rotation. Firing
rate correlated with eye velocity. (LE=left eye position,

LE=left eye velocity,

H=head
velocity) [48].
Position neurons need to integrate velocity information, e.g. from vestibular system
about head rotation
Model for an integrator:
single neuron with excitatory autapse (synapse of neuron onto itself)
network with excitatory recurrent connections
Integrator needs to receive input corresponding to the motor commands (efference
copy
25
) to update the position.
25
Efference copies inform the brain also about expected sensory input (nger touching hand) to dis-
tinguish sensory signals created by the animal itself rather than by external inuence. Efference copies
presumably underly the observation that one cannot tickle oneself.
115
Figures/SeLe00a.f1.ps
Figures/SeLe00a.f3.ps
Figure 83: a) Network model for integrator in okulomotor control via recurrent exci-
tation. The integrator network receives efference copies of the motor control signals
(excitatory and inhibitory bursts) to update the stored eye position. b) Results from the
integrator model of [52].
For these integrators exact integration requires
= 1
For
< 1 the eye position would slowly relax back to center position
For
> 1 the eye position would diverge to the sides

The relaxation/divergence occurs with the time-constant
=

r
1
Experimentally one nds that the eye position relaxes with a time constant
relax
= O(10s)
(in goldsh [52]) and
r
= 5ms . . . 100ms
For this type of integrator to be sufcient one would need the relevant eigenvalue to be
precisely tuned,
[
1[ =

r
relax
=
1
2000
. . .
1
100
.
So far, it is not clear whether there are mechanisms that would allow such a precise
tuning of the synaptic strengths.
7.3.3 Limulus Vision. Contrast Enhancement
Output of the retina is not simply a pixelated picture of the visual scene:
The receptive elds of retinal ganglion cells have a spatial center-surround structure
116
and in addition temporal structure (Fig.64)
non-trivial computations are performed already in the retina
Simple example: contrast enhancement by lateral inhibition
Figures/mach_band_david_anson_bw_no_separator.ps
Figures/mach_band_david_anson_bw_with_separator.ps
Figure 84: Mach bands. Contrast enhancement by lateral inhibition.
The vertebrate (and mammalian) retina is very complex (Fig.85).
Figure 85: Many different types of cells contribute to the function of the retina, among
them horizontal cells that provide lateral inhibition [38].
Consider simpler visual system: Limulus (horse-shoe crab)
117
Figures/horseshoe-crab.jpg Figures/Horseshoe-crab-eye-detail.jpg
Figure 86: Horseshoe crab and its compound eye
The compound eye consists of O(1000) ommatidia. Each ommatidium has an individual
lens and is functioning as 1 pixel.
Vision in Limulus was in particular studied by H.K. Hartline, starting in the 1930s (Nobel
Prize 1967 [25]).
118
Figures/Ha69.f1.ps
Figure 87: Electrical activity in single optic nerve. Stimulation with 3 different intensities
[25].
Figures/HaRa57.f1.ps
Figures/HaRa57.f2.ps
Figure 88: Stimulated ommatidia inhibit each other reciprocally with threshold [26].
The output from a given ommatidium depends on the stimulation of the ommatidia in
its surround:
lateral inhibition
linear in the ring rate of the inhibiting ommatidium
threshold: minimal ring rate needed for inhibition
Note:
The inhibition is not driven directly by the stimulus to the other ommatidium but
by the output of the other ommatidium: recurrent network
119
Simple Rate Model:
r
dv
i
dt
= v
i
+ F(I
s,i
(t))
s
dI
s,i
dt
= I
s,i
+ h
i
..
light stimulus
j
w
ij
[v
j
v
thresh
]
+
Note:
It is found that the synaptic strength w
ij
depends on the distance between the
ommatidia.
Focus on stationary state:
v
i
= F(I
s
) I
s
= h
i
j
w
ij
_
v
j
v
thresh
+
thus
v
i
= F
_
h
i
j
w
ij
_
v
j
v
thresh
+
_
Simplications
one-dimensional continuum of neurons: v
i
v(x)
to avoid boundaries assume periodic domain < x +
observation: suppression of ring is linear in the ring rate of the suppressing
neuron
equations become linear
distance dependence of synaptic strength:
w
ij
g M(x x
) =
_
g x x
x +
0 otherwise
Thus
v
(x) = h(x) g
_
+
M(x x
)v
(x
)dx
(54)
To solve this integral equation note
equation is linear
120
corresponds to matrix equation (49)
solve using eigenvectors, i.e. eigenfunctions, satisfying
26
_
x+
x
M(x x
)e
n
(x
)dx
=
n
e
n
(x)
equation is translation invariant: M(x, x
) = M(x x
)
eigenfunctions are Fourier modes
e
(c)
n
(x) = cos nx e
(s)
n
= sin nx
Determine eigenvalues:
_
x+
x
M(x x
) cos nx
dx
=
_
x+
x
cos nx
dx
=
=
1
n
(sin n(x + ) sin n(x ))
2
n
cos nx sin n
_
x+
x
M(x x
) sin nx
dx
=
_
x+
x
sin nx
dx
=
2
n
sin nx sin n
Thus
0
= 2 for eigenvectore
(c)
0
= 1
n
=
2
n
sin n for e
(c)
n
= cos nx and fore
(s)
n
= sin nx n > 0
Expand v
(x) and h(x) into Fourier series

v
(x) =
n=0
v
(c)
n
cos nx +
n=1
v
(s)
n
sin nx h(x) =
n
h
(c)
n
cos nx + h
(s)
n
sin nx
and insert into (54)
_
x+
x
M(x x
)v
(x
) dx
n=1
_
x+
x
M(x x
)
_
v
(c)
n
cos nx
+ v
(s)
n
sin nx
_
dx
=
=
n=1
n
v
(c)
n
cos nx +
n
v
(s)
n
sin nx
Using the orthogonality of different trigonometric functions (e.g.
_
+
cos nx cos mxdx =
nm
)
one can collect all cos nx and sin nx separately to get
v
(c,s)
n
=
1
1 + g
n
h
(c,s)
n
26
cf. DA Chap. 7.4
121
Consider small step in luminosity
h(x) =
_
1 x < 0
1 + x > 0
with expansion
h
(s)
2n+1
=
4
(2n + 1)
n = 0, 1, 2 . . .
h
(s)
2n
= 0 h
(c)
n
= 0 n = 0, 1, 2, . . .
h
(c)
0
= 1
Thus
v
(x) =
1
1 + 2g
+
n=0
1
1 + g
2
2n+1
sin(2n + 1)
4
(2n + 1)
sin (2n + 1) x (55)
a)
Figures/mach_plot.eps
b)
Figures/visual_lateral.1.eps
Figure 89: Mach band: lateral inhibition enhances small differences in luminosity (g =
1). a) solution (55) b) qualitative sketch.
Note:
differences between pixels are enhanced by lateral inhibition
Mach bands
contrast enhancement
a)
Figures/Ha69.f13a.ps
b)
Figures/KnTo70.f6.ps
Figure 90: a) Sharpening of on-transient by delayed lateral inhibition. Upper curve:
response if only a single ommatidium is illuminated. Lower non-monotonic curve: il-
lumination covers also neighbors, which provide delayed inhibition [25]. b) Temporal
shape of lateral interaction: short excitation followed by longer inhibition [40].
122
Dynamical behavior:
Lateral inhibition arises with a delay
initially strong response of all cells, then inhibition kicks in ring of cells de-
crease inhibition decreases ring rate increases again after a delay
non-monotonic approach to steady state
In addition self-inhibition:
rises very quickly, decays then more slowly than lateral inhibition
Could model lateral interaction with a double-exponential as for the slow synapse (cf.
(30))
rapidly decaying excitation
slowly decaying inhibition
I
s,i
= h
i
+ I
(e)
s,u
+ I
(i)
s,i
e
dI
(e)
s,i
dt
= I
(e)
s,i
+
j
w
(e)
ij
_
v
j
(t) v
(e)
thresh
_
+
i
dI
(i)
s,i
dt
= I
(i)
s,i

j
w
(i)
ij
_
v
j
(t) v
(i)
thresh
_
+
7.3.4 Persistent States
Consider again nonlinear recurrent network (17)
dv
dt
= v +F(h +Mv)
So far we did not investigate any possible role for the nonlinearity F
Consider for simplicity minimal model:
single neuron with sigmoid ring-rate function with threshold v
dv
dt
= v + F(h + mv) F(h + mv) =
1
1 + e
(h+m(vv
))
Can combine threshold and feed-forward input: v
th
= (v
h) /m
dv
dt
= v +
1
1 + e
m(vv
th
)
(56)
123
Steady states v = v
satisfy the transcendental equation

v
=
1
1 + e
m(vv
th
)
Cannot be solved in terms of simple functions: use graphical approach
small m:
straight line intersections sigmoid in exactly 1 point: single steady state (xed
point)
large m:
large positive v
th
: single xed point with v
0
large negative v
th
: single xed point with v
1
intermediate values of v
th
: three xed points simultaneously
can show
outer xed points v
(1)
0 and v
(3)
1 are linearly stable

intermediate xed point v
(2)
is linearly unstable: solution exists, but any

small perturbation will push v away from v
(2)
to either v
(1)
or v
(3)
Easier visualization of dynamics using a potential
dv
dt
=
dU
dv
with U(v) =
1
2
v
2
_
v
1
1 + e
m(v
v
th
)
dv
Consider large m and approximate F by a step function

F(v) =
_
1 for v > v
th
0 for v < v
th
then
U(v) =
_
1
2
v
2
(v v
th
) for v > v
th
1
2
v
2
for v < v
th
Minima at
v = 1 for v
th
< 1
v = 0 for v
th
> 0
Thus:
v
th
< 0 single minimum at v
(1)
= 1
0 < v
th
< 1 two minima at v
(1)
= 0 v
(3)
= 1
1 < v
th
single minimum at v
(3)
= 0
124
Figures/bistable_potential.eps
Figure 91: Potential U(v) for three values of v
th
.
During the evolution the value of the potential U(v) cannot increase
dU(v(t))
dt
=
dU(v)
dv
dv
dt
=
_
dU
dv
_
2
0
Particle with position v is sliding down the potential landscape U
Thus:
eventually the neuron will always go into a steady state
neuron can be active v
1 even without any input h: persistent activity

neuron can be bistable:
depending on initial conditions the neuron will go into the steady state v
(1)
or v
(3)
neuron can function as a memory

Working memory used to remember items for a brief time (e.g. phone numbers,...)
Figures/delay1.eps Figures/delay2.eps Figures/delay3.eps
Figure 92: Delayed memory task. The cue (red circle) tells the monkey which of the
two blue buttons it is to press after the delay, i.e. when the red go-signal comes on.
125
Figures/ChGo98.f8b.ps
Figure 93: Firing rate of a population of neurons in prefrontal cortex in rhesus macaque
monkeys. Neurons re persistently during the delay period while the stimulating cue is
off. The neurons stop ring when the monkey makes a saccade to the goal indicated
by the cue [5].
7.3.5 Associative Memory
27
We would like to store multiple patterns in a network.
Simplify the systemand its dynamics yet further by considering a network of N McCulloch-
Pitts neurons [46]
neurons are ON (v = +1) or OFF (v = 1)
input: sum of weighted inputs
step function F(x) = sgn(x)
discrete time steps: iteration v
i
v
i
v
i
= sgn
_
j
M
ij
v
j
_
(57)
Figures/hyper_cube.eps
Figure 94: State space of a network consisting of 3 McCulloch-Pitts neurons.
Note:
27
[27] Ch.2
126
the space of possible states are the corners of an N-dimensional hypercube
the dynamics consists in hopping from one corner to another
Goal:
nd M such that given k patterns v
(s)
, s = 1 . . . k, are stable xed points of the
iteration (57)
Consider symmetric M, M
ij
= M
ji
, with vanishing diagonal
28
, M
ii
= 0,
Introduce the energy function
E =
1
2
j
M
ij
v
i
v
j
The energy is non-increasing under the iteration (57):
when some component v
l
changes to
v
l
= sgn
_
j
M
ij
v
j
_
v
l
+ v
l
the energy change is (using M
ii
= 0)
E = v
l
j
M
lj
v
j
Then
E = sgn
_
r
M
lr
v
r
_

r
M
lr
v
r
+ v
l
_
r
M
lr
v
r
_
=
r
M
lr
v
r
+ v
l
..
=1
r
M
lr
v
r
Thus
E =
_
2 [
r
M
lr
v
r
[ forsgn(
r
M
lr
v
r
) = sgn(v
l
) v
l
did change
0 for sgn(
r
M
lr
v
r
) = sgn(v
l
) no change in v
l
Moreover, E is bounded from below since v
l
= 1
Note:
for any such M and any initial conditions the dynamics always leads to a xed
point
28
Diagonal elements need not be chosen to vanish, but the performance is somewhat better that way
[27].
127
discrete motion in a discrete energy landscape
no energy exists in general if M is not symmetric
How to choose the coupling matrix M to store k specic patterns, i.e. to have specic
patterns as stable xed points?
Single pattern: Consider v = ( v
1
, . . . , v
N
),
sgn
_
j
M
lj
v
j
_
= v
i
Try
M
ij
=
_
v
i
v
j
for i ,= j
0 for i = j
Then each term in the energy sum has the same sign for that pattern v making the
energy extremal for that pattern
E( v) =
1
2
j
M
ij
v
i
v
j
=
1
2
ij
v
2
i
v
2
j
=
..
M
ii
=0
1
2
_
N
2
N
_
Consider the iteration (57)
j
M
ij
v
j
= v
i
N
i=j=1
v
2
j
= (N 1) v
i
Choose = 1/ (N 1). Then v is a xed point of the iteration.
What happens if some of the components of the initial conditions v do not agree with
v? E.g. v agrees with v only in the rst m components,
v = ( v
1
, v
2
, . . . v
m
, v
m+1
, . . . , v
N
)
Then
29
j
M
lj
v
j
=
m
j=1
M
lj
v
j
j=m+1
M
lj
v
j
=
1
N
v
l
(m(N m)) =
1
N
(2mN) v
l
Thus, for m > N/2 the sign of this sum is the same as for v
v is also mapped to v.
Notes:
v is an attractor
29
Strictly speaking, since M
ll
= 0 one needs to distinguish between l m and l > m. But it makes
only a small difference. For l m one gets (N 1)
1
(2m N 1) whereas for l > m one gets
(N 1)
1
(2mN + 1).
128
the mapping M rst projects
30
the vector v onto v and then scales up the vector
again to full size.
for initial conditions within the basin of attraction the iteration corrects incorrect
components of the initial condition
initial conditions that contain too many errors do not converge to the attractor,
they are outside the basin of attraction of this attractor and in the basin of attrac-
tion of another attractor
following the work of John Hopeld [31] this type of associative memory has been
studied extensively from a computational perspective.
Thus:
The network serves as associative memory (content-addressable memory):
the memorized state v is retrieved (reached) even if it is only partially remem-
bered in the initial condition
Multiple patterns. Consider v
()
, 1 k, and try
M
ij
=
1
N
v
()
i
v
()
j
Is each pattern indeed a xed point of the iteration?
Start the iteration at one of the patterns v
()
j
M
ij
v
()
j
=
1
N
v
()
i
j
v
()
j
v
()
j
= v
()
i
+
1
N
j
v
()
i
v
()
j
v
()
j
. .
cross-talk
Thus,
analogous to the case of a single stored memory M rst projects v onto the
subspace spanned by the patterns v
()
and then scales the vector up again via
the nonlinearity
if
1
N
v
()
i
v
()
j
v
()
j
< 1
the sign of
j
M
ij
v
()
j
is the same as that of v
()
j
the pattern v
()
is correctly retrieved.
30
M can be written as M = v v
t
, with v
t
being a row vector and v as column vector.
129
again the iteration can complete the pattern if the initial condition v
(init)
differs not
too much from v
()
because of the cross-talk the error tolerance is reduced
expect that the cross-talk becomes more signicant for larger numbers of stored
patterns
the network has only a nite capacity k
max
for storing patterns, which depends
on the error rate accepted. For random patterns one can show
k
max
0.14 N
the network has often also spurious memories, i.e. attractors that do not corre-
spond to the k stored patterns
because of symmetry v
i
v
i
inverted pattern is always also an attractor
some mixtures of patterns like v
i
= v
(
1
)
i
v
(
2
)
v
(
3
)
i
are also attrac-
tors (mixtures composed of even number of patterns could add up to 0 and
cannot be in the state space.
if the patterns to be stored are orthogonal to each other,

j
v
()
j
v
()
j
=
, the
cross-talk vanishes: no error in the iteration
but: for k = N one has, using matrix notation (V)
i
v
i
,
cross-talk: V
t
V = I V
t
= V
1
weight matrix: VV
t
= VV
1
= I
thus the weight matrix is diagonal and the network preserves any input: no
memory at all.
i.e. since retrieval corresponds to projection onto the subspace spanned
by the stored patterns: N orthogonal vectors span the whole space the
projection becomes the identity.
for k < N the performance does improve if the patterns are orthogonal
often preprocessing is employed to make the patterns of interest more or-
thogonal
Note:
If the network learns one pattern at a time
M
(+1)
ij
= M
()
ij
+
1
N
v
(+1)
i
v
(+1)
j
the weight update resembles Hebbs rule (neurons that re together wire to-
gether):
130
weight is strengthend if pre- and post-synaptic neuron are active
weight is weakened if only one of them is active
but:
weight is also strengthened when neither neuron is active
weight can be positive or negative: same synapse can switch from exci-
tatory to inhibitory or vice versa
model can be modied to make it physiologically more reasonable
8 Statistical Approach
In particular in cortex
variability of spiking (cf. Fig.55)
deterministic variability or noise?
relevant for function or nuisance?
during short intervals only small number of spikes, which uctuates strongly be-
tween trials
imprecise and unreliable transmission of information
need to merge multiple sources/types of sensory inputs:
multi-sensory integration: how to weigh different sources of information
Gain insight into working of brain by considering illusions, i.e. situations in which the
brain seemingly draws incorrect conclusions
Examples:
unreliable information: moving rhombi at low contrast
http://www.cs.huji.ac.il/
~
yweiss/Rhombus/
merging information for disambiguation: Necker cube with additional information
http://www.cut-the-knot.org/Curriculum/Geometry/Necker.shtml
http://www.dogfeathers.com/java/necker.html
barber pole illusion
http://www.psychologie.tu-dresden.de/i1/kaw/diverses%20Material/www.illusionworks
com/html/barber_pole.html
In particular at low contrast: imprecise information due to low spike rate:
interpret the spiking of the neuron as information about the probability of a certain
sensory input, say.
131
Consider simple example: neuron sensing the position s of an object (prey, say)
Goal for the brain area that is reading that unreliable information:
What is the probability for the prey to be at a certain position s under the condition that
a certain neuron spikes during a brief time t (r = 1/t),
P(s[r)
What information does the neuron provide directly?
the neuron has a tuning curve, i.e. r = f(s)
Assuming a Poisson distribution for the spiking of the neuron as a function of its mean
ring rate r one obtains for the probability of the neuron to spike during a brief intervat
t
P(r[s) =
1
1!
( rT)
1
e
rt
rt = f(s) t
How do we convert between these two probabilities?
Consider joint probability
P(r, s) = P(r[s)P(s) = P(s[r)P(r)
Which implies Bayes Theorem
31
P(s[r)
. .
Posterior
=
1
P(r)
P(r[s)
. .
Likelihood
P(s)
..
Prior
Notes:
Prior: probabiliy distribution of stimuli independent of whether the neuron spikes
or not
The prior reects the knowledge (expectation) that the system (animal) has a
priori about its environment.
It may have learned the prior on different time scales
under evolutionary control (long time scales)
through previous input, e.g.
memories stored from previous exposures
previous spikes of the same neuron (affecting neural spiking via facilia-
tion and/or depression) e.g. ambient light
Posterior: probability distribution of the stimulus incorporating the fact that the
neuron spiked
Best guess for the position of the prey presumably at the maximum of the pos-
terior (MAP=maximum a posteriori estimate)
dP(s[r)
ds
smax
= 0
subsequent brain areas may use this best guess to base decisions on
31
Thomas Bayes, mathematician and Presbyterian minister (1702-1761).
132
P(r) does not affect the dependence on the stimulus: we can usually ignore it
Examples:
i) Gaussian tuning curve f(s) e
(ss
0
)
2
/2
2
ii) Gaussian tuning curve with linear prior: e.g., prey is approaching from the right
shift in s
max
iii) successive measurement by the same neuron: for interpretation of second spike
use the posteriori obtained from the rst spike as prior sharpening of distribution.
a)
Figures/bayes_tuning.eps
b)
Figures/bayes_two_spikes.eps
Figure 95: a) Posteriori for Gaussian tuning curve and linear prior. b) Flat initial prior,
posterior of rst and second spike used as prior for second and third spike, respectively.
iv) two neurons:
if the spiking of the two neurons is given by two independent Poisson processes
P(s[r
1
, r
2
) =
1
P(r
1
, r
2
)
P(r
1
, r
2
[s)P(s)
1
P(r
1
, r
2
)
P(r
1
[s)P(r
2
[s)P(s)
a)
Figures/bayes_two_neurons_close.eps
b)
Figures/bayes_two_neurons_further.eps
Figure 96: Two neurons with Gaussian tuning curve and different preferred location. a)
equal widths. b) unequal widths.
When combining the information from two neurons, the peak is shifted towards the pre-
ferred location of neuron with narrower likelihood (which has more precise information).
Note:
In the case of the Necker cube the likelihood for both congurations are exactly
the same as is the prior
both congurations have also equal posteriors and the percept switches.
integration of additional information necessary to disambiguate the percept.
133
8.1 Bayes and the Cutaneous Rabbit
A tactile illusion:
stimulate mechanically the skin with brief pulses (2 ms) sequentially at two locations
(cf. Fig.97a)
Perception:
...There is a smooth progression of jumps up the arm, as if a tiny rabbit were hopping
from wrist to elbow...[22]
...In one experiment it was possible to induce a vivid hopping that traveled up one
arm, across the neck, and down the other arm by the use of ve contacts, two widely
separated on each arm and one on the nape of the neck...[22]
The touch can even be perceived at locations that are anesthetized.
a)
Figures/Go07a.f1A.ps
b)
Figures/Go07a.f1B.ps
c)
Figures/Go07a.f1C.ps
Figure 97: Sketch of stimulation protocol (positions in cm, time in seconds). a) Rabbit
can be interpreted as length contraction. b) -effect: the more rapidly traversed of
two equal distances is perceived as shorter. c) Two-arm comparison -effect: different
lengths are perceived as the same if theu are traversed in different times [23]
Can these illusions be understood as arising from a rational model for the perception
of the brain?
Need a high-level model:
consider a Bayesian observer, i.e. model probabilities for response and percept:
describe perceived trajectory as
x = b + mt
with perceived endpoints
x(0) = b x(t) = b + mt
for simplicity consider two neurons with Gaussian receptive elds centered at the
true endpoints of the trajectory, x
1
and x
2
, respectively,
f
i
(x)
1
s
e
(xx
i
)
2
/2
2
s
, i = 1, 2
134
Likelihood of the activation pattern D of these neurons if the trajectory was given
by x = b + mt and the time interval was t,
32
P(D[m, b)
1
s
e
(bx
1
)
2
/2
2
s
1
s
e
(b+mtx
2
)
2
/2
2
s
the values m and b can be considered as hypotheses about the trajectory
assume: eventual percept of the trajectory is given by the maximum (mode)
of the probability for the trajectory given the stimulus-evoked neural activation
pattern D
P(m, b[D).
m
max
and b
max
would constitute the most likely hypotheses about the trajectory.
i.e. use maximum a posteriori (MAP) estimate.
to obtain P(m, b[D) in terms of P(D[m, b) use Bayes theorem
P(m, b[D) P(D[m, b) P(m, b)
we need a prior P(m, b):
assume:
the initial point could be anywhere (at prior for b)
small velocities are more likely than large velocities
P(m, b) 1
1
v
e
m
2
/2
2
v
combined
P(m, b[D)
1
s
e
(bx
1
)
2
/2
2
s
1
s
e
(b+mtx
2
)
2
/2
2
s
1
1
v
e
m
2
/2
2
v
Figure 98: Probability distributions of the model in terms of initial position and velocity
(a) or initial and end position (b). Factual values inidicated by cross [23].
32
naive Bayes: assume that the conditional probability of the response of the two sensory modalities
are independent of each other, i.e. can factorize the joint conditional probability. noise of the two
modality is uncorrelated
135
If the brain operates using the MAP estimate then the perceived trajectory should cor-
respond to m
max
and b
max
that maximize P(m, b[D)
Determine maximum by minimizing exponent
b
ln (P(m, b[D)) =
1
2
s
((b
max
x
1
) + (b
max
+ m
max
t x
2
)) = 0
b
max
=
1
2
(x
1
+ x
2
m
max
t)
m
ln (P(m, b[D)) =
1
2
s
(b
max
+ m
max
t x
2
) t
1
2
v
m
max
= 0
b
max
=

2
s
t
2
v
m
max
m
max
t + x
2
yielding
m
max
=
x
t
1
1 +
2
(t)
2
and
b
max
= x
1
+ x
1
2 + (t)
2
with
=

v
s
characterizing the reliability of the position information and the degree of bias towards
low velocities.
Maximum of posterior is at a smaller velocity
Consequently the perceived lengths is contracted (since times are kept xed)
x
m
max
t = x
1
1 +
2
(t)
2
i.e.
x
x
=
1
1 +
2
(t)
2
The initial point is perceived as shifted towards the nal point by
b
max
x
1
= x
1
2 + (t)
2
This model has a single t parameter :
large : high spatial acuity and/or low expectation for slow motion
136
small : low spatial acuity and/or strong expectation for slow motion
a)
Figures/Go07a.f3A.ps
b)
Figure 99: a) Perceived length contraction as a function of duration of trajectory. True
distance between taps is 10cm, perceived distance is l
. b) Two-arm -effect: ratio of

lengths that yield the same perceived length as a function of the ratio of the durations
of the trajectories (t
2
= 1s t
1
). Value for was t separately for a) and b). [23]
Two-arm -effect: different lengths are perceived as equal if the durations are different
x
1
x
2
= 1 =
x
1
1
1+
2
(t
1
)
2
x
2
1
1+
2
(t
2
)
2
In the experiment a specic protocol was used (t
max
= 1s)
t
2
= t
max
t
1
Rewrite length ratio in terms of
=
t
2
t
1
to yield
x
1
x
2
=
1 +
2(1+)
2
2
t
2
max
1 +
2(1+)
2
2
t
2
max
2
Notes:
model provides coherent rationalization of
137
perception of length contraction
-effect (dependence of perceived distance on time between stimuli)
-effect (dependence of perceived time between stimuli on distance between
them)
merges the expectation of slow movements with the observed imprecisely deter-
mined positions
reasonable success suggests that
the brain may be using such probabilities in the information processing
the brain may be employing optimal merging of unreliable information in
Bayes fashion
how could probabilities be encoded?
one can show that for Poisson-like spike statistics [61]
neural populations can encode probabilities
implementing Bayes theorem, which requires multiplying two probability dis-
tributions, can be achieved by adding the activities of the neurons encoding
the two distributions
this gives a function to the high variability of cortical spiking activity
8.2 Integration of Multi-Sensory Information
33
How should possibly conicting information from different sensory inputs be com-
bined?
How does the brain do it?
Experimental study of the integration of haptic and visual information
visual information about the height of a bar via three-dimensional image:
to modify the reliability of that information the images are corrupted by vari-
ous degrees of noise
haptic (touch) information provided by simulated force feed-back
the haptic and visual information disagree with each other by O(10%)
33
following [16]
138
a)
Figures/ErBa02.f2a.ps
b)
Figures/ErBa02.f2b.ps
Figure 100: a) Set-up: three-dimensional image of a bar and simulated force feed-back
from a bar of possibly different height. b) Stereoscopic image of bar without and with
noise [16].
Theoretical treatment:
Assume independent Gaussian likelihoods for the visual and the haptic information
P
v
(r
v
[S)
1
v
e
(SSv)
2
2
2
v
P
h
(r
h
[S)
1
h
e
(SS
h
)
2
2
2
h
thus
P
vh
(r
v
, r
h
[S)
1
h
e
(SSv)
2
2
2
v
e
(SS
h
)
2
2
2
h
Assume a constant prior
P(S) = P
0
Posterior
P
vh
(S[r
v
, r
h
) P
vh
(r
v
, r
h
[S) P
0
Note:
Since the prior is constant:
maximum a posteriori estimate (MAP) is equivalent to a maximum likelihood es-
timate (MLE)
Rewrite in terms of a single Gaussian
(S S
v
)
2
2
v
+
(S S
h
)
2
2
h
=
1
2
vh
(S S
vh
)
2
+ c
0
with
1
2
vh
=
1
2
v
+
1
2
h

2
vh
=

2
v
2
h
2
v
+
2
h
(58)
S
vh
= w
v
S
v
+ w
h
S
h
with w
v
=

2
vh
2
v
w
h
=

2
vh
2
h
(59)
Note:
139
if the uncertainty in the two responses is not equal, the response of the input that
is less uncertain, i.e. with
2
smaller, is to be weighted more.
To test whether humans use MLE in this context one needs to know S
vh
as well as S
v
and S
h
and the associated uncertainties,
2
v
,
2
h
, and
2
vh
Measure uncertainty via psychometric function:
have subjects repeatedly compare the height of the test stimulus with that of
comparison stimuli, which have various heights
psychometric function (S): relative frequency of the subjects ranking the com-
parison stimulus (variable height) to be higher than that of a test stimulus (height
55mm)
point of subjective equality (PSE) = 50% level:
(S
PSE
) = 0.5
discrimination threshold T
v,h
: (S
PSE
+ T
v,h
) = 0.84
Figures/ErBa02.f1.ps
Figure 101: Sketch of the likelihoods (top) and of the psychometric function (bottom)
for two values of the ratio
2
h
/
2
v
with xed discrepancy between the visual and haptic
bar height. Top: solid line gives estimated height based on combined information.
Bottom: The psychometric function gives the fraction of trials in which the comparison
stimulus (variable height) is perceived as higher than the test stimulus (xed average
height
1
2
(S
v
+ S
h
) = 5.5cm) [16].
the cumulative Gaussian (error function) gives a good t of (Fig.102). For that
function one has
T
v,h
=
v,h
this determines the individual
2
h
and
2
v
for the different noise levels in the visual
input.
140
Figures/ErBa02.f3a.ps
Figure 102: Measured psychometric function for estimates based on haptic or visual
information alone (visual information corrupted by noise) [16]
Figures/ErBa02.f3b.ps
Figure 103: Noise dependence of the combined visual-haptic height estimate [16]
the measured
2
v,h
give the theoretically optimal weights w
v,h
for the integration
w
v
=

2
vh
2
v
w
h
=

2
vh
2
h
experimentally determined weights from (59) using w
h
= 1 w
v
w
v
=
S
vh
S
h
S
v
S
h
quite good agreement with the optimal weights (Fig.104)
Note:
the weights must adjust within the 1 second of the presentation, since noise levels
vary random from one presentation to the next.
141
Figures/ErBa02.f3c.ps
Figure 104: Measured and MLE-optimal weights entering the combined estimate as a
function of noise level [16]
integration of visual and haptic information reduces the variance (58), which low-
ers the threshold for discrimination (T
vh
=
vh
)
Figures/ErBa02.f3d.ps
Figure 105: The threshold for combined visual-haptic height discrimination is lower
than each individual threshold [16].
Not quite class-appropriate example of an illusion arising from the merging of unreliable
auditory with precise, but incorrect visual information
http://www.youtube.com/watch?v=7-ZnPE3G_YY
142
List of Figures
1 The multiple scales of the brain [54] (reproduced in [7]) . . . . . . . . . . 10
2 Different neuron morphologies in visual cortex. webvision.med.utah.
edu/VisualCortex.html . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3 Purkinje cells (A) and basket cells (B) in cerebellar cortex. . . . . . . . . 12
4 Action potential simulation
http://www.cs.cmu.edu/
~
dst/HHsim/. . . . . . . . . . . . . . . . . . . 12
5 Ion channels in membrane. . . . . . . . . . . . . . . . . . . . . . . . . . 15
6 Na
+
-K
+
-exchange pump.
Movie http://highered.mcgraw-hill.com/olc/dl/120068/bio03.swf 17
7 Modeling one-dimensional ion channel. . . . . . . . . . . . . . . . . . . 20
8 The giant squid axon is not an axon from the giant squid (from [39]) . . 23
9 Reduction of action potential upon partial replacement of extracellular
seawater by dextrose, reducing extracellular [Na
+
] . a) 67% dextrose, b)
50%, c) 29%. Curve 1 and 3 (wash-out) are with seawater [30] . . . . . 24
10 Channels are rectifying: no current for hyperpolarization (to -130mV,
top), transient current for depolarization (to 0mV, bottom gure). Note:
sign convention of current opposite to the one used nowadays. [29] . . . 25
11 Shape of transient current depends on the membrane voltage.The num-
bers marking the curves denote hyper-polarization relative to resting po-
tential: negative numbers indicate depolarization [29]. . . . . . . . . . . 25
12 Extraction of Na
+
-current by partial replacement of extracellular Na
+
by
choline. Na
+
-current is transient. . . . . . . . . . . . . . . . . . . . . . 26
13 Ohmic resistance of K
+
-channel. . . . . . . . . . . . . . . . . . . . . . . 27
14 Sigmoidal growth and exponential decay of g
k
with voltage [28]. . . . . 27
15 Steady-state levels m
, h
, n
and the time scales

m,h,n
for the ac-
tivation and inactivation variables of the Na
+
- and K
+
-channels of the
Hodgkin-Huxley model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
16 Comparison of current trace from numerical simulation of HH-model
(left) with the experimental curves (right). [28]. . . . . . . . . . . . . . . 30
17 a) Comparison of the Ohmic resistance in the Hodgkin-Huxley model
with the current given by the Goldman-Huxley-Katz model (2) [39].
b) Voltage at which no current ows when channels for Na
+
and K
+
are
both open as a function of the K
+
concentration for GHK-model and for a
combination of Ohmic resistances for the Na
+
- and the K
+
-conductance,
with E
Na
+ and E
K
+ determined separately. . . . . . . . . . . . . . . . . 31
143
18 Voltage dependence of the parameters in the Connor-Stevens model.
Note that the parameters for I
Na
and I
K
are also different than in the
HH-model [9]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
19 Brain structures: thalamus (above spinal cord and hindbrain) functions
as relay for sensory information. . . . . . . . . . . . . . . . . . . . . . . 34
20 a) Steady-state activation and inactivation M
and H
for I
T
, c) time
constant
M
of activation, d) time constant
H
of inactivation. Note the
very slow recovery from inactivation for V < 70mV [33]. . . . . . . . . 35
21 Burst due to I
CaT
. Delay of the action potential due to I
A
. Note decrease
in frequency across burst [10]. . . . . . . . . . . . . . . . . . . . . . . . 36
22 Response of a thalamic neuron (in LGN) in the tonic and the bursting
regime. Same current injection induces tonic ring from a holding po-
tential of V = 59mV (a) but to a transient burst for V = 70mV (b).
Graded vs. discontinuous response of the neuron in the two regimes (c).
In the tonic regime the ring rate encodes the sinusoidal input very well
(d), in the bursting regime only the onset of each wave is encoded (e) [55] 37
23 Rhythmic bursting of thalamic relay neurons. [45] . . . . . . . . . . . . . 37
24 Sag current I
h
. Current and voltage clamp results. The sag refers to the
slow decrease in hyperpolarization after strong hyperpolarizing current
injection (arrows). [45] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
25 a) Activation curve of I
h
(obtained from 7 different neurons, solid sym-
bols 1 neuron). b) Time scales for activation and deactivation, c) en-
largement of the tail current after repolarization to V = 49mV (between
arrows in b). [45] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
26 Voltage dependence of activation and deactivation times. Fits to tran-
sients using single exponentials are quite good. [45] . . . . . . . . . . . 39
27 Periodic bursting in a model for a stomatogastric ganglion cell. . . . . . 40
28 Strong correlation between n and h in standard Hodgkin-Huxley model. 41
29 Deviations of the activation variables m, n, h from their steady-state val-
ues m
, n
, h
during regular ring. The large values of m/m
and
h/h
are due to m
and h
becoming very small during deactivation

and inactivation, respectively. . . . . . . . . . . . . . . . . . . . . . . . . 42
30 a) stable xed point, large perturbation can lead to excitation. b) unstable
xed point, periodic spiking. . . . . . . . . . . . . . . . . . . . . . . . . 43
31 Generic (typical) trajectories in the phase plane close to a xed point.
Eigenvectors need not be orthogonal. . . . . . . . . . . . . . . . . . . . 45
32 Parameters of Izhikevichs simple model. . . . . . . . . . . . . . . . . . . 51
33 Sketch of a pair of xed points on an invariant circle (parametrized by the
phase of the oscillation during regular ring). The mutual annihilation of
the xed points in a saddle-node bifurcation generates an innite-period
limit cycle (SNIC-bifurcation leading to Type-I oscillation). . . . . . . . . 52
144
34 Some of the neuronal types captured by Izhikevichs simple model [35].
Compare behavior labeled thalamo-cortical with that discussed for I
CaT
-
current in Sec.3.4.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
35 Neurons are not pointlike. Their morphology varies widely. . . . . . . . . 54
36 Voltage attenuation along a dendrite. A) current injection in soma trig-
gers large action potential there but delayed and reduced action potential
in dendrite. B) converse arrangement with stimulation (via synapses) of
dendrite [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
37 Sketch of currents in an axon [10] . . . . . . . . . . . . . . . . . . . . . 55
38 Voltage v(x) for steady current injection (A) and for short current pulse
(B) [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
39 Burst in thalamic reticular neuron driven by low-threshold Ttype Ca
2+
-
current which is slower than the T-current in thalamocortical relay cells
(cf. Sec.3.4.2) [12]. For movies see
Thalamic reticular neurons http://cns.iaf.cnrs-gif.fr/alain_movies.
html . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
40 Gating of dendritic inputs by other dendritic inputs in pyramidal cells in
CA1 in hippocampus [37]. Movie at http://people.esam.northwestern.
edu/
~
kath/gating.html. . . . . . . . . . . . . . . . . . . . . . . . . . . 61
41 Sketch of gap junctions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
42 Electromicrograph of synapse displaying many vesicles and the active
zone (darker area at synaptic cleft). . . . . . . . . . . . . . . . . . . . . 63
43 Sketch of chemical synapse
(www.biochem.uni-erlangen.de/MouseDB) . . . . . . . . . . . . . . . . . 64
44 Fit of (29) to experimentally obtained EPSC (averaged over many stim-
uli) [13]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
45 Integration of multiple inputs by slow membrane. IF-modell (34,35) with
= 1,
1
= 0.05, g
1
= 0,
2
= 0.025, g
2
= 50, E
syn
= 50 a) I
e
= 50 b)
I
e
= 250. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
46 Integration of multiple inputs by slow synapse. IF-modell (34,35) with
= 0.1,
1
= 0.5, g
1
= 0,
2
= 0.25, g
2
= 10, E
syn
= 50, a) I
e
= 302, b)
I
e
= 400. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
47 Different time scales and different voltage dependence of synaptic cur-
rents evoked by AMPA and by NMDA receptors. CNQX blocks AMPA
receptors, APV blocks NMDA receptors. [57] . . . . . . . . . . . . . . . 72
48 Dependence of NMDA-conductance on [Mg
2+
] and voltage [10]. . . . . 72
49 Response of depressing synapse to varying input ring rates [10] . . . . 76
145
50 Receptive elds of neurons in primary visual cortex. Hubel & Wiesel
video.
Long movie http://people.esam.northwestern.edu/riecke/Vorlesungen/
Comp_Neuro/Notes/Movies/VisualCortex.mov
There are a number of short movies on www.youtube.com. . . . . . . . 77
51 Spike responses to odor inhalation of a mouse. Bars in horizontal lines
denote spikes in different trials under the same conditions [1]. . . . . . 78
52 Neuron in primary motor cortex. a) raster plots during reaching move-
ments in different directions. b) average ring rate as a function of the
reaching direction [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
53 Spike train and ring rates obtained by binning in xed windows (b), in
sliding rectangular window (c), in sliding Gaussian window (d), causal
window (e). [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
54 a) Probability of having exactly n spikes as a function of rescaled ring
rate rT. b) For not too small rT the Poisson distribution is well approxi-
mated by a Gaussian distribution with a variance that is equal to its mean
[10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
55 High ring variability of cortical neurons (area MT of an alert monkey)
reecting balanced excitatory and inhibitory input from a large number
of neurons. Variance grows (super-linearly) with mean ring rate [53]. . 83
56 Low ring variability of retinal ganglion cells (sub-Poisson). [17] . . . . . 84
57 Mean ring rate and variability of ring rate of peripheral (SA1 and RA)
and cortical neurons (3b and 1) in response to tactile stimulation of
macaque monkey [56]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
58 Time-dependent stimulus and resulting spiking activity of electrosensory
neuron in weakly electric sh [10]. . . . . . . . . . . . . . . . . . . . . . 85
59 Procedure to compute the spike-triggered average stimulus [10]. . . . . 85
60 Spike-triggered average stimulus for a pyramidal neuron in the elec-
trosensory lateral line-lobe, which is the rst central nucleus of the elec-
trosensory pathway in weakly electric sh [20] (in [10]). . . . . . . . . . 87
61 Response prediction using the optimal kernel (39) for the H1 neuron in
the visual system of the y. The H1 is sensitive to horizontal motion of
large portions of the visual eld.
a) image velocity used as input. b) two spike sequences obtained in
response to that input. c) measured (dashed) and predicted (solid) ring
rate [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
62 Structure of the retina. a) Drawing by Cajal (1911). b) Sketch of cone
and rod pathway in the retina [38]. . . . . . . . . . . . . . . . . . . . . . 90
63 Visual pathway from the retina to primary visual cortex [10]. . . . . . . . 90
146
64 a) Center-surround receptive eld of an on-center LGN-neuron in cat. b)
Temporal structure of receptive eld of LGN-neurons in cat. Horizontal
axis: x, vertical axis: t, green: excitatory, red: inhibitory. A) Excitatory
response at light onset. B) Excitatory response delayed, but stronger
than immediate inhibitory response. [11] . . . . . . . . . . . . . . . . . 91
65 LN-model for OFF Y-type retinal ganglion cell for low and high contrast
of the stimulus. a) Membrane voltage response. b) Spike rate response.
[62] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
66 Comparison of the voltage prediction of LN-model for the OFF Y-type
retinal ganglion cell shown in Fig.65 for low and high contrast of the
stimulus (stimulus in top row) [62].. . . . . . . . . . . . . . . . . . . . . . 96
67 LN-model for the same OFF Y-type retinal ganglion cell as in Fig.65. Af-
ter rescaling of the lter the nonlinearities become the same for low and
high contrast of the stimulus (grey=low contrast, black=high contrast). a)
Membrane voltage response. For high contrast gain is reduced to 0.83
of the gain for low contrast. b) Spike rate response. Gain reduced to
0.64 [62] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
68 LN-model for ON Y-type retinal ganglion cell. After rescaling of the l-
ter the nonlinearities become the same for low and high contrast of the
stimulus (grey=low contrast, black=high contrast). a) Membrane voltage
response. Gain redcued to 0.74. b) Spike rate response. Gain reduced
to 0.54. [62] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
69 Input outside the classical receptive eld can modify a ganglion cells
response. a) grating exciting the periphery is shifted every 9s. b) The
response of this ganglion cell is transiently changed from OFF-center to
ON-center briey after a peripheral shift. [21] . . . . . . . . . . . . . . . 99
70 Extraction of two temporal lters using a principal component analysis
(PCA) of the spike-triggered covariance matrix [21] . . . . . . . . . . . . 99
71 Receptive elds: A) ON-center in LGN, B) OFF-center in LGN, C)-G)
simple-cell receptive elds in V1 [32] . . . . . . . . . . . . . . . . . . . . 105
72 Feed-forward network for simple cell [32]. . . . . . . . . . . . . . . . . . 105
73 Receptive eld of a complex cell. A)-D) same orientation but different
location of the bar within the receptive eld. E)-G) different orientations.
H) bar with orientation as in A)-D) moved rapidly across the receptive
eld. The horizontal bars above the record indicate when the stimulus is
on [32]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
74 Possible connectivity of a complex cell: input from multiple simple cells
with different preferred locations of the bar [32]. . . . . . . . . . . . . . . 106
75 Coordinate transformation between body-centered coordinates and retina-
centered coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
147
76 Response of a bimodal (visual-tactile) neuron in the ventral pre-motor
cortex of a macaque monkey. Its tactile receptive eld is near the left
cheek. The visual response to an object approaching along paths marked
I-V does not depend on the eye position (A1, B1, C1) (or the arm posi-
tion (B2)), but it does depend on the head position (B3), i.e. it depends
on the position of the object relative to its tactile receptive eld (similar
data for tactile eld on the arm) [24]. . . . . . . . . . . . . . . . . . . . . 108
77 Tuning curves of face+visual bimodal neurons in ventral pre-motor cortex
are in head-centered coordinates. a) The receptive eld is shifted when
the head is turned. b) The receptive eld is not shifted if the gaze, but
not the head is rotated [24]. . . . . . . . . . . . . . . . . . . . . . . . . . 108
78 Tuning curve of neurons in the posterior parietal cortex. The preferred lo-
cation is given in retinal coordinates, but the magnitude of the response,
i.e. the height of the tuning curve, is modulated by the gaze direction [3]. 108
79 a) The neuron in PMv integrates the contributions from all neurons in PP
with a weight w(, ) depending on their preferred location (, ). Their
receptive elds are represented by circles. b) For a different stimulus
(s
, g
) the each of the PP-neurons responds differently. But for each

neuron at (, ) there is a neuron at (
) with
= + that
responds like the neuron (, ) did for stimulus (s, g). The output of the
PMv-neuron depends only on s + g if the weight for the neuron (
)
is equal to the weight of neuron (, ), i.e., if the weight w(, ) depends
only on + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
80 Sketch indicating location of medulla and other parts of the brain stem
[38]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
81 a) Burst-position neuron (A,C) and position neuron (B,D) in the mon-
key medulla (HE=horizontal eye position, HT=horizontal target position,
FR=ring rate). b) The duration of the bursts is strongly correlated with
the duration of the saccades and usually precedes them by 10ms
(solid bars in inset for burst-position cells, open bars for position cells)
[47]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
82 Neurons in goldsh medulla. coding for eye position (A,B) or eye veloc-
ity (D). A) spontaneous eye movement. Firing rate (FR) correlated with
left eye position (LE). B) Sinusoidal head rotation (in the dark) with com-
pensating eye movement. Firing rate correlated with eye position, not
eye velocity. D) Sinusoidal head rotation. Firing rate correlated with eye
velocity. (LE=left eye position,

LE=left eye velocity,

H=head velocity) [48].115
83 a) Network model for integrator in okulomotor control via recurrent ex-
citation. The integrator network receives efference copies of the motor
control signals (excitatory and inhibitory bursts) to update the stored eye
position. b) Results from the integrator model of [52]. . . . . . . . . . . . 116
84 Mach bands. Contrast enhancement by lateral inhibition. . . . . . . . . . 117
148
85 Many different types of cells contribute to the function of the retina,
among them horizontal cells that provide lateral inhibition [38]. . . . . . 117
86 Horseshoe crab and its compound eye . . . . . . . . . . . . . . . . . . . 118
87 Electrical activity in single optic nerve. Stimulation with 3 different inten-
sities [25]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
88 Stimulated ommatidia inhibit each other reciprocally with threshold [26]. 119
89 Mach band: lateral inhibition enhances small differences in luminosity
(g = 1). a) solution (55) b) qualitative sketch. . . . . . . . . . . . . . . 122
90 a) Sharpening of on-transient by delayed lateral inhibition. Upper curve:
response if only a single ommatidium is illuminated. Lower non-monotonic
curve: illumination covers also neighbors, which provide delayed inhibi-
tion [25]. b) Temporal shape of lateral interaction: short excitation fol-
lowed by longer inhibition [40]. . . . . . . . . . . . . . . . . . . . . . . . 122
91 Potential U(v) for three values of v
th
. . . . . . . . . . . . . . . . . . . . 125
92 Delayed memory task. The cue (red circle) tells the monkey which of the
two blue buttons it is to press after the delay, i.e. when the red go-signal
comes on. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
93 Firing rate of a population of neurons in prefrontal cortex in rhesus macaque
monkeys. Neurons re persistently during the delay period while the
stimulating cue is off. The neurons stop ring when the monkey makes
a saccade to the goal indicated by the cue [5]. . . . . . . . . . . . . . . 126
94 State space of a network consisting of 3 McCulloch-Pitts neurons. . . . 126
95 a) Posteriori for Gaussian tuning curve and linear prior. b) Flat initial
prior, posterior of rst and second spike used as prior for second and
third spike, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
96 Two neurons with Gaussian tuning curve and different preferred location.
a) equal widths. b) unequal widths. . . . . . . . . . . . . . . . . . . . . 133
97 Sketch of stimulation protocol (positions in cm, time in seconds). a)
Rabbit can be interpreted as length contraction. b) -effect: the more
rapidly traversed of two equal distances is perceived as shorter. c) Two-
arm comparison -effect: different lengths are perceived as the same if
theu are traversed in different times [23] . . . . . . . . . . . . . . . . . . 134
98 Probability distributions of the model in terms of initial position and ve-
locity (a) or initial and end position (b). Factual values inidicated by cross
[23]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
99 a) Perceived length contraction as a function of duration of trajectory.
True distance between taps is 10cm, perceived distance is l
. b) Two-
arm -effect: ratio of lengths that yield the same perceived length as a
function of the ratio of the durations of the trajectories (t
2
= 1s t
1
).
Value for was t separately for a) and b). [23] . . . . . . . . . . . . . . 137
149
100 a) Set-up: three-dimensional image of a bar and simulated force feed-
back from a bar of possibly different height. b) Stereoscopic image of
bar without and with noise [16]. . . . . . . . . . . . . . . . . . . . . . . . 139
101 Sketch of the likelihoods (top) and of the psychometric function (bot-
tom) for two values of the ratio
2
h
/
2
v
with xed discrepancy between
the visual and haptic bar height. Top: solid line gives estimated height
based on combined information. Bottom: The psychometric function
gives the fraction of trials in which the comparison stimulus (variable
height) is perceived as higher than the test stimulus (xed average height
1
2
(S
v
+ S
h
) = 5.5cm) [16]. . . . . . . . . . . . . . . . . . . . . . . . . . . 140
102 Measured psychometric function for estimates based on haptic or visual
information alone (visual information corrupted by noise) [16] . . . . . . 141
103 Noise dependence of the combined visual-haptic height estimate [16] . 141
104 Measured and MLE-optimal weights entering the combined estimate as
a function of noise level [16] . . . . . . . . . . . . . . . . . . . . . . . . 142
105 The threshold for combined visual-haptic height discrimination is lower
than each individual threshold [16]. . . . . . . . . . . . . . . . . . . . . 142
150
A Basic Concepts of Numerical Methods for Differen-
tial Equations
Consider differential equations of the type
dy
dt
= F(t, y) y(t = 0) = y
0
In general y is a vector: coupled equations.
Note:
all higher-order equations can be written as systems.
E.g.
d
2
y
dt
2
= F(t, y,
dy
dt
)
can be written as
dy
dt
= v
dv
dt
= F(t, y, v)
For simplicity we consider scalar equations
dy
dt
= F(t, y)
For numerical approximation discretize time
y
n
y(t
n
)
We need then
For small time steps t we can perform a Taylor expansion
y
n+1
= y(t
n
+ t) =
= y(t
n
) + t
dy
dt
tn
+
1
2
t
2
d
2
y
dt
2
tn
+ ...
We cannot keep all terms in the expansion: truncate
y
n+1
= y
n
+ t F(t
n
, y
n
) +c
local
(t, t
n
) (60)
where the error c
local
(t, t
n
) contains all the omitted terms.
This is the forward Euler
34
scheme.
Expect
c
local
(t, t
n
) = O(t
2
) = a(t
n
)t
2
+ b(t
n
)t
3
+ . . .
Note:
34
Leonhard Euler (1707-1783)
151
in each time step the scheme makes an error O(t
2
): this is the local error
to obtain y(t
max
) = y(Nt) we need to make N steps
in these N steps the error could accumulate to a global error
c
global
=
N
n=1
c
local
(t, t
n
)
N
n=1
a(t
n
)t
2
N at
2
= at
max
t
expect that the error decreases linearly when t is increased.
if c
global
0 for t 0 a scheme is called convergent
The Euler scheme is called a rst-order accurate scheme.
To obtain accurate solutions the Euler scheme requires very small t: slow.
Higher-order schemes with c
global
= O(t
n
) with n = 2, 3, 4, .. can be much more
efcient (faster) to obtain accurate solutions. Examples are Runge-Kutta, Crank-
Nicholson, predictor-corrector,... [49]
Example:
dy
dt
= y
with exact solution
y(t) = y(0)e
t
Note:
for < 0 the solution decays to 0 for any initial condition
Euler scheme
y
n+1
= y
n
+ t y
n
= (1 + t) y
n
= (1 + t)
2
y
n1
= . . .
i.e.
y
n
= (1 + t)
n
y(0)
Does this solution decay for any t if < 0?
[y
n
[ = [ (1 + t)
n
[ [y(0)[ = [ (1 + t) [
n
[y(0)[
Thus: the solution decays only if [1 + t[ < 1,
1 < 1 + t < +1
For < 0 this is only the case if
t < t
max

1
[[
The forward Euler scheme has a stability limit.
For t > t
max
:
152
solution grows exponentially in magnitude with the number of time steps (al-
though exact solution decays to 0)
in each time step the solution switches sign
the scheme is unstable for these large time steps: the numerical solution has
nothing to do with the true solution of the differential equation.
Thus:
the numerical solution of an equation can only be trusted if it converges as t 0
i.e. one has to solve the differential equation for sequentially smaller values of t
and the solution is to change less and less as t becomes smaller and smaller
c
global
t
p
p = 1 for Euler
ln c
global
p ln t
p is given by slope on log-log-plot. Always check whether the scheme obtains the
correct slope over at least a decade in t
unstable schemes typically generate solutions that oscillate with each time step
independent of the size of the step:
halng the time step doubles the frequency of the oscillations
any oscillation that corresponds to a true solution of the differential equation will
have a period that is
independent of t for small t
signicantly larger than t (to resolve an oscillation one needs on the order
of 10 time steps/period).
Choice of time step:
accuracy: time step needs to be small enough to resolve fastest change in the
true solution
stability: scheme needs to remain stable
Stiff problems:
true solution evolves on a slow time scale: would allow large t
system has very rapid modes that decay and are not visible in the nal solution:
requires small t
153
Example (see homework):
dy
dt
= y + Asin t < 0 [[
after a rapid decay during a time 1/[[ the solution evolves slowly with frequency .
In the Taylor formula expand around t
n+1
y
n
= y
n+1
t
dy
dt
t
n+1
+O(t
2
)
yielding the scheme
y
n+1
= y
n
+ t F(y
n+1
, t
n+1
) (61)
This is the backward Euler scheme, which is also rst-order accurate.
In our simple linear example we get:
y
n+1
= y
n
+ t y
n+1
which is easily solved for y
n+1
y
n+1
=
1
1 t
y
n
Consider magnitude of the solution
[y
n
[ =
1
[1 t[
n
[y(0)[
thus, [y
n
[ 0 with n for arbitrarily large values t if < 0.
The backward Euler scheme is unconditionally stable.
Note:
of course, for very large t the solution may not be very accurate, but no insta-
bility will arise.
choice of time step purely based on the desired accuracy of the solution.
But:
the scheme (61) is implicit :
the unknown y
n+1
appears also on the right-hand side inside the function F(y, t)
the equation is an implicit equation for y
n+1
, which may be difcult to solve if
F(y, t) is nonlinear in y.
implicit schemes typically require additional effort (e.g. Newton iteration method
for solving the nonlinear equation)
For more details about numerical schemes see, e.g. Numerical Recipes
35
[49] or the
lecture notes for ESAM 346 Modeling and Computation in Science and Engineering
36
.
35
Numerical Recipes http: // www. nr. com
36
ESAM 346 Modeling and Computation in Science and Engineering http://people.esam.
northwestern.edu/
~
riecke/Vorlesungen/346/Notes/notes.346.pdf
154
B Linear Response Kernel using Functional Derivatives
In a more fancy fashion the optimal kernel describing the linear linear response of a
neuron can be derived using functional derivatives (cf. Sec.6.2.2).
E is a functional of D()
ED() =
1
T
_
T
0
_
r
0
+
_

0
D()s(t
)d r(t
)
_
2
dt
For functions minima are characterized by a vanishing derivative. Analogously, minima

of functionals can be found elegantly by using the functional derivative. One can avoid
the explicit use of functional derivatives by converting the integral into a sum (see
Appendix 2.9A in [10])
At a minimum the functional derivative vanishes
ED()
D()
= 0
Basic property of functional derivatives
D()
D(
)
= (
) Dirac function
In short, this is analogous to
x
x
= 1
y
x
= 0
or for a function that depends on the components v
i
of a vector one has
v
i
v
j
=
ij
Using this basic property
ED()
D(
)
=
1
T
_
T
0
D(
)
_
r
0
+
_

0
D()s(t
)d r(t
)
_
2
dt
=
1
T
_
T
0
_
2
_
r
0
+
_

0
D()s(t
)d r(t
)
_
_

0
_
_
_
_
_
_

D(
)
D()
_
. .
(
)
s(t
)
_
_
_
_
_
d
_
_
dt
=
2
T
_
T
0
_
r
0
+
_

0
D()s(t
)d r(t
)
_
s(t
)dt
!
..
= 0 at a minimum
155
NOTES:
discuss also some of the networks stuff in Rolls&Treves (good for homework sim-
ulations and to get a feeling for these kind of networks)
check out book by Trappenberg (Tr09)[58]
look at book by Cox (GaCo10)[19]
also Ermentrout and Terman [15]
include also population coding: Fisher information etc.
other simple spiking models in addition to Izhikevich (e.g. Gerstner: Badel 2008).
Put something like that in the homework? Fit things.
LNP and GLM (i.e. LNP with feedback of spike history for refractory effects). cf.
Alonso paper DeJi10
156

Intro To Computational Neuroscience

Uploaded by

Copyright:

Available Formats

Intro To Computational Neuroscience

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Intro To Computational Neuroscience

Uploaded by

Copyright:

Available Formats

Computational Neuroscience H.

Riecke, Northwestern University

-channels hyperpolarize cells to negative potentials.

driven by the electric eld E = for each species

(r) = dielectric constant of the medium

(V ) asymptotic value reached after a sufciently long time (t

(V ) > 0 is initially linear

(V ) increases with increasing voltage

(V ) decreases with increasing voltage

and other ions)

and the time scales

during regular ring. The large values of m/m

becoming very small during deactivation and inactivation, respectively.

Thus, for large I

(V ) 1, while n and h vary only slowly.

thus the ring rate goes to zero like r

(x) dx = 1. The -function has

v(x) dx goes to 0 since the integrand is nite

with reversal potential E

For low ring rate (T

For low ring rate (T

is close to the asymptotic value R =

6.1 Poisson Spike Trains

Consider the correlation function between stimulus and ring rate

Use white noise without any correlations between different times

Using the linear response assumption

= 0 and interchanging and

, . . .) are to be determined to give a good approximation for a wide

is adjusted to avoid double counting.

d are the numbers of neurons with preferred location and gaze in

) the each of the PP-neurons

(t) is the time-dependent

of v depends on the associated eigenvalue

< 0 the output component v

is reduced compared to the input com-

< 1 the output component v

is larger than the input component

h: selective amplication of components of h

1 the approach to the steady state becomes innitely slow

(critical slowing down)

diverges exponentially within this linear approximation

no saturation, linear growth

network functions as an temporal integrator of its input

> 1 the eye position would diverge to the sides

(x) and h(x) into Fourier series

cos nx cos mxdx =

satisfy the transcendental equation

1 are linearly stable

is linearly unstable: solution exists, but any

Easier visualization of dynamics using a potential

Consider large m and approximate F by a step function

1 even without any input h: persistent activity

neuron can function as a memory

. b) Two-arm -effect: ratio of

and the time scales

during regular ring. The large values of m/m

becoming very small during deactivation

) the each of the PP-neurons responds differently. But for each

For functions minima are characterized by a vanishing derivative. Analogously, minima