SYDE750: Simulating Neurobiological Systems
University of Waterloo
Winter 2010
Chris Eliasmith
MW 11:30 - 12:50
E2-1303E
Course website
Things engineers know that I should look up.
Table of contents (hide)
1. 2010-01-04 (Jan. 4)
RC circuits
- See Wikipedia. Note the curve (similar to leaky integrate and fire neuron).
2. 2010-01-06 (Jan. 6)
Our representation of scalars, vectors, and functions are all basically the same, so if one is not making sense, fall back to a simpler case.
Fourier analysis
Wiener filter (snicker)
In the following equations,
= Input current
= Action potentials
The transformation in a neuron from current to action potentials is as follows:
is the spike generation function
The current,
, is defined by:
is the amount of current based on
, the external stimulus
is the neuron sensitivity (different for each neuron)
is the base-level current (some neurons are always spiking)
Which gives us the following equation for determining the transformation:
There are a number of different spike generation functions observed by neurons. This will be gone into detail later, but for now:
- Rectified linear = just a line
- Leaky integrate and fire = curve similar to RC circuit (equation on page 37)
Note that a response curve is not the same as a tuning curve! The response is a neuron's response to current, and is similar across neurons of the same type. A tuning curve shows a neuron's response to the physical world, and varies across all neurons, due to differences in sensitivity.
For linear decoding, we must use tuning curves.
We'll optimize to find the linear least squares optimal decoder. (Note:
)
Or, in matrix form:
To do this in Matlab, we make a matrix
such that each row is a discretized version of some neuron
's tuning curve:
Then,
(i.e. the dot product between neurons' tuning curves)
(
is a vector, just going from -1 to 1 at granularity
)
Our only unknown is
, which we can obtain by:
But is
invertible? Perhaps not, so we do pseudo-inversion instead. Matlab has a function pinv to do this.
3. 2010-01-11 (Jan. 11)
In neurobiological systems, there is lots of noise, but high efficiency. Sources of noise include:
- Axonal jitter (some randomness (order of nanoseconds) in time down an axon)
- Vesicle release failures
- Thermal noise (minor)
- Ion channel noise (probability of opening or closing)
- Network effects (lots of stuff together -- chaos)
- See "Sources of noise (slide 14)
Have to add
to old equations, so:

This leads to the following characterization of error:
- The first term is the error due to static distortion
- The second term is the error due to noise
A full derivation of this error:
We can group the first three terms as
. Then, since we're doing the average (
is like averaging) we sum over
so we get 0 (since mean of noise is 0).
This is true because:
In matrix form (suitable for Matlab!), we get the following:
Note that adding the diagonal matrix
makes the
matrix invertible, so we don't need pinv!
The question that we're trying to answer is: how do you weight each neuron such that, for all values of
, we can get a good estimate
?
With vector bars, the derivation of our equations is exactly the same as the scalar case.
. This is called the preferred direction vector.- In the scalar case, we just made this on or off (1 or -1)
- In the vector case, it's a vector; e.g. [-0.71, 0.71]

The derivation is in the appendix of the textbook.
In Matlab, for the vector case...
still ends up as
; you can make the
matrix two dimensional, or you can take row and double the length.
So, in Matlab, you do the same thing, only with more dimensions (d instead of 1).
(note:
is a dot product, which is a cosine.
For Georgopolous decoders, the encoding step is the same, but the decoding is
4. 2010-01-13 (Jan. 13)
Temporal representation in neurons
See slide 4 for LIF neuron graph -- pay close attention to the constants. Note that in real neurons, after the spike, there's a 2-3ms dip below the graph in which no amount of driving current will affect the graph (absolute refractory period).
Slide 3 (should have potassium / sodium ions): sodium is much bigger than potassium.
For temporal representations, we essentially want to model the voltage, as it's the voltage difference across the cell membrane that defines if a neuron spikes or not. So, for the leaky integrate and fire neuron, we have the following derivation (also available in the notes, with extra info):
is the charge
is the capacitance (this seems to be the same across neurons --
(resistance) changes a bit
is current (see diagram)
is the passive leak current
is the current flowing in
is a time constant (see diagram)
So, now we have a first-order differential equation. We want to solve for V. Note that this LIF tuning curve would be from a constant
.
This is a convolution integral. The length you are in the past is exponentially less important (we forget at an exponential rate).
determines the speed of forgetting.
Of course, in Matlab, to do ODEs you can just use
. Deriving it just gives insight (but it is cool insight).
So, to get a tuning curve
isn't a function of time, it's just a constant.
So what is the firing rate?
is the time to threshold
is the firing rate at the time to threshold
is the refractory period
is the threshold current (which we usually set to 1)
Different modelers get information from different sources:
- Rate code: Firing rate over a window (e.g. 70 Hz)
- Timing code: Time between each spike (interspike interval)
Temporal encoder:

We need a temporal decoder,
in
Evaluating gives
.
Fourier analysis defines a domain of functions (see slides 22-24).
Limit the set of dimensions (frequencies) for the set of functions.
5. 2010-01-18 (Jan. 18)
Temporal representation (cont):
- Encoding:
- Decoding:
- Note: * in the decoding refers to a convolution.
Convolution ([1] [2]) and filters ([1] [2]).
How do we generate an x(t) to use in our Matlab simulations? See genSignal.m.
% Inputs: % T = time window (slowest frequency) % dt = timestep % rms = amount of power % bandwidth = which frequencies we care about % rand = random seed for testing and such genSignal(T, dt, rms, bandwidth, rand) % % N = number of samples (in complex space) N = (floor(T ./ 2 * dt)) + 1 % dw (i.e. d-omega) is in radians dw = 2 * pi * (1 ./ T) w = dw * [0 : (N-1)/2] N = length(w) Amps = zeros(1,N) + i * zeros(1,N) index = w <= max(bandwidth) && w > min(bandwidth) % set seed % Next line gives band limited gaussian white noise Amps(index) = randn(1, length(index)) + i * randn(1, length(index)) % conj gives the complex conjugate Amps = [Amps, fliplr(conj(Amps(2:end)))] % iff gives the inverse fourier transform S = real(iff(Amps)) rmsS = sqrt(sum(S.^2) ./ length(S)) normS = rms / rmsS S = S * normS Amps = Amps * normS
Fourier transform of input puts in frequency space (derivation is in text).
The window is there to do sampling.
The optimal decoder compares
(black lines) and
(autocorrelation function).
In the following,
is a filter,
is the spikes.
Note that our optimal filters go into negative time, meaning that they are not causal (as you would have to be able to predict perfectly each spike).
6. 2010-01-20 (Jan. 20)
Complex single cell models
Theta-neuron, saddle-node bifurcation:
Note that
represents where we are in phase space (i.e. the spike).
Also note that there is slightly different notation here. Using our notation:
7. 2010-01-25 (Jan. 25)
Population-temporal representations
Finding the optimal h(t) also gives some weight (in the form of amplitude).
- Implicitly, even the temporal decoders gives a population decoder.
- Same weight for on and off neurons (assignment 2 note: always considering 2 neurons).
N.B. Always normalize h(t) (optimal or not) to area=1 before using decoders.
- For postsynaptic current (PSC) decoding, multiply by the weight given by the optimal signals.
Note that, in this case,
is a weighted function.
As for noise, pre-noising is better for finding optimal filters.
Fluctuations help discern signal from noise over long-distance channels.
is the somatic current? Which is not constant
is the time between spikes (constant for the LIF neuron model)
is the starting time
is constant
- So, for constant input, we get ripples regardless of noise.
8. 2010-01-27 (Jan. 27)
Feedforward transformations
So, to find the transformation for x + y, we find
and
, stick it into
, and then find
.
The recipe on slide 21, using the z = x + y example, boils down to:
- Substitute into
9. 2010-02-01 (Feb. 1)
Nonlinear transformations
We need nonlinearity (though most are small nonlinearities). Still, there are some "big" ones, like in object recognition.
- How common are they? Single cell or network? Generally it is not clear.
Slide on locust visual system: T, not
.
Coincidence detection slide: take the spikes and put little gaussians at each spike.
- See notes for derivation of this formula.
is some constant (for AMPA 10 ms, etc.)
There is not much evidence that coincidence detection happens in mammals. Also, it sucks (since it needs an overconnected structure).
- Of course, we could eventually find some evidence, so let's not knock it too much, kay?
For networks, we do nonlinear computation through layers of linear computations.
As an example, let's try to find an optimal linear decoder to approximate the product. We'll end up with an expression
- where
(instead of
)
See the notes for the full derviation; keep in mind:
is the encoding for the middle layer
means to project
into
space
- The decoders
project back into one-dimensional space
This holds for any non-linearity, not just the product.
Dendritic subunits: embed m on a dendritic spine.
10. 2010-02-03 (Feb. 3)
Function representation
Slide 2: Motor = uniform distribution, ear canal = along one of m dimensions (m indexes dimension)
V1 tuning curve slide: Encoding function is
- where
is the orientation change as a function of orientation
LIP example in textbook.
Think about
. This means that
is parameterized by
.
What is every possible height and position that can be represented by the system? Can do whatever decomposition you want.
If all are Gaussian shapes, (encoding vectors normalized) we can get the function space it represents with SVD (singular value decomposition).
gives the domain for the optimization problem. This is hard, so we'll do random sampling.
contains the coefficients of an orthonormal basis. The only thing that changes is
;
gets integrated out.
Now,
is a function instead of a matrix.
Size of
:
- Scalar:
- Vector:
- Function:
(how you discretize the function space)
However,
is always an
matrix.
If we do SVD on the encoders, we can get the basis functions (
)
- See slide 32 for equations.
Basically, we map into the space of basis function coefficients.
11. 2010-02-08 (Feb. 8)
This be the most math-y of all the lectures!
What are the
s such that
?
is always square and symmetric. So:
For p-inv, if
is 0,
is 0 (otherwise it becomes
because
)
Legendre basis: which functions is it good at computing?
Slide 57: Random and "intelligently designed" neurons do the same!
12. 2010-02-10 (Feb. 10)
Once noise gets close to singular values, we mis-scale our
s. SNR = 1, signal = noise, so it's not helping anymore.
Neural control theory:
Classic control = find a transformation
See notes for dynamics given the post-synaptic currents, and so on.
Note that it's more convenient to work in the Laplace space. This takes the dynamical system to algebra.
Note that convolving with
is the same as mulitplying by
in the Laplace domain.
See the slides (slide 6 on).
We can transfer information about control theory to neural dynamics.
is the thing represented by a population of neurons.
Slide 8: superscript = population, subscript = neuron.
This is principle 3 of the NEF:
- this gives what is on slide 13.
13. 2010-02-22 (Feb. 22)
Communication channel with yourself = neural integrator / memory (because the dynamics attempt to keep the state stable).
with no subscript (from
etc) is
which is defined by the type of neurotransmitter (well, receptor) (e.g. NMDA = 100 ms, AMPA = 5 ms)
R: Drift velocity
Desired (in the Laplace domain):
If we assume that
and
(which are both probably not true) then
;
- which is the desired iff
- note that we're using the same filter for each neuron
Let's take a closer look at the dynamics of time-changing
.
Change to the Laplace domain:
So
tells us how stuff will change.
If we assume
is really big and the others are small (insignificantly small), then:
Less error if we just add a second dimension that holds information on the variable that controls it.