Quantitative Methods in Neuroscience (neu 366M) Homework 10

Quantitative Methods in Neuroscience
(neu 366M)
Homework 10
Due: Thursday Nov 20, by the beginning of class
In this assignment, we will explore first-order linear differential equations and neural
models whose subthreshold response is characterized by this dynamics.
1) The leaky integrate-and-fire (LIF) model neuron Consider a simple spiking neuron
model, the leaky integrate-and-fire neuron (LIF). The LIF model neuron is identical
to the RC circuit model for subthreshold voltage that we studied in class:
Cm
dV
= grest (Vrest − V ) + Iapp (t),
dt
with the addition of a threshold condition for spiking: once the voltage reaches a
threshold Vth , it is assumed to have fired a spike, and the voltage is reset to some
value Vreset .
The LIF model differs from the Hodgkin-Huxley neuron model in two key ways. In
both models, the voltage dynamics well below threshold approximate an RC circuit.
However, in Hodgkin-Huxley, the different ionic conductances that produce the resting
potential are explicitly modeled while in LIF neurons, all the different conductances
are replaced by a single “rest” conductance grest and its associated battery Vrest .
Next, in both models, neurons fire an action potential once their voltage is sufficiently
depolarized. In the Hodgkin-Huxley model, the ionic conductances that produce the
stereotyped spike shape are modeled in detail, but in the LIF model, the spike shape
is unimportant and spikes are put in by hand. We saw the LIF model in action in
class. Now it’s your turn to play with it.
a. Simple subthreshold RC dynamics: Download IFsim.m and run the RC circuit numerical simulation for subthreshold voltage dynamics (no spike-generation
mechanism), trying a few different constant levels of input current Iapp . Plot the
resulting voltage as a function of time. Add to this plot the analytical curve
for voltage we derived in class, using the values of Cm , grest , Vrest etc. that are
used in the simulation. Remember that to get a τ on the left side, you need to
divide the entire equation by grest , which means that Iapp (t) in your analytical
expression should be scaled appropriately.
b. Let Iapp be a time-varying random Gaussian noise (uncomment the line for random Iapp ). Plot the resulting numerically obtained voltage response V (t). Also
plot on top of this curve the analytical result, that V (t) is the convolution of the
exponential filter τ1m e−t/τm with Iapp , summed with Vrest . Evaluate the convolution itself in matlab, using conv. Again, agreement should be excellent!
c. Vary the time-step ∆t over a large range and explore when the numerical simulation result starts to degrade (by comparing with the analytical result from
[b.] above). Plot the squared error (sum across time-steps of squared difference
between the analytical and numerical result, divided by total number of timesteps) as a function of ∆t. As you should see, ∆t should be much smaller than τ ,
the RC time-constant of the membrane voltage, for good numerical performance
(optional: verify that if τ is larger, then ∆t can be made larger without hurting
accuracy).
The numerical integration method we are using, known as “forward Euler”, is
one of the less efficient/robust numerical integration techniques, but it is conceptually simple and easy to implement. In applications in the real world, you
would consider using Runge-Kutta and other “higher-order” numerical integration methods, which can produce more accurate results for a given ∆t.
d. The f-I curve. Now run the spiking model. Set the reset voltage Vreset = Vrest .
Systematically vary Iapp to be at different constant levels, and determine the
corresponding firing rate of the LIF neuron (use the spks variable in IFsim.m to
compute firing rate, by dividing number of spikes by elapsed time). Plot firing
rate versus Iapp . This is the f-I curve of the LIF neuron. Change Vreset so that it
is at a more hyperpolarized value than Vrest . This resetting of voltage post-spike
to a hyper- polarized value is effectively a “soft” refractory period. Replot the
f-I curve. How has it changed?
e. In the spiking model, allow Iapp to be random. Run the simulation for a longer
duration, for instance 5 seconds. How does the autocorrelation of the spike train
(spks) look? What happens to this autocorrelation if Vreset is more hyperpolarized than Vrest ? How about when it is more depolarized?
2) The first-order linear differential equation as a low-pass filter. In class, we
showed that first-order linear differential equations of the form
τ
dx
= −x + b(t)
dt
(1)
have solutions of the form x(t) = conv(b(t), h(t)), where h(t) is an exponentially
decaying filter with time-constant τ : h(t) = τ1 e−t/τ . This filter is called a linear lowpass filter, and Equation (1) is called a linear low-pass filtering equation for the input
b(t). Here, we’ll see why it’s called this.
a. Derive the result of the convolution conv(b(t), h(t), where h(t) is the exponential filter given above and b(t) is the cosine function cos(ωt), where ω is the
frequency. [Hint: It is much easier to evaluate the integral if you write write
ix
−ix
cos(x) = e +e
]. Also, in the end, express your result after removing any
2
complex numbers in the denominator by multiplying the top and the bottom
by the complex-conjugate of the denominator. What does convolution with the
exponential filter do to cosine waves of different frequencies?
b. Write Matlab code to numerically integrate Equation (1), with b(t) set to a cosine
wave of unit amplitude and some frequency ω. Let x(0) = 0 and make a choice
for τ . Plot the input b(t) and the output x(t) on the same plot. What has
changed and what has remained the same? Use max to get the amplitude of x(t).
Obtain the amplitude of x for different input frequencies ω, and plot amplitude
as a function of frequency ω (make sure you use frequencies much larger than and
much smaller than 1/τ ). Add to this plot the amplitude part of the expression
you derived analytically in [a.]. Are they a good match (they should be!). Do
you see why the filter is called a low-pass filter?
c. In your numerical integration code, let b(t) step from zero to some non-zero value
starting at t = 200 ms. Plot b(t) and x(t). What has the low-pass filter done to
the sharp step? Any function can be written as a sum of sines and cosines, and
fast changes like the sharp edge of the step onset, is composed of high-frequency
elements. Interpret the response x(t) to the step in light of this fact and your
observations in [a.-b.] above.

Download Report