MITOCW9.40-Problem Set 2 Solved

Learning objectives and outcomes: 

 

This problem set investigates two models of how neurons spike in response to current inputs. The leaky integrate and fire (LIF) model is based on the simple RC model of a cell, with an extra rule that the model neuron spikes whenever the membrane potential passes above a threshold. Despite its simplicity, the ideas behind this model constitute the cornerstone of many fundamental computational modelling approaches. In contrast to this simplified model, real neurons produce spikes due to the presence of voltage-dependent ion channels. The Hodgkin-Huxley model is the first model to capture these rich dynamics. It does so by incorporating biophysically accurate descriptions of the sodium and potassium channels. More elaborate biophysically inspired models take the work of  Hodgkin-Huxley as a starting point.

 

These models were covered in lectures 3-5 and recitations 3&4.  

 

The expected learning outcomes for this PSET are:

 

•      Estimate model parameters from different plots of model behavior.

•      Interpret the biophysical role of these parameters.

•      Write code to simulate a model to check your parameter estimations.

•      Produce f-I curves (firing rate versus current injected).   

•      Be able to present your results as a cohesive and well-structured report.

 

MATLAB functions you will need: 

 

There is quite a lot of overlap with PSET 1 in terms of MATLAB functions needed to successfully complete this PSET.

 

For calculating firing rates from spike times, the function  diff will be super helpful.

  

As a reminder the main plotting functions are: figure, plot, subplot, xlim, ylim, title, xlabel, ylabel, legend, hold on, hold off. 

 

For more information on MATLAB functions and commands, check the provided cheat sheet or the more extensive MATLAB documentation.  

 

Problem 1: Estimating the LIF parameters 

 

The provided MATLAB file LIFmodel.p is a scrambled function that simulates a leaky integrate and fire neuron for square pulses of current injection. By scrambled we mean that you can run the function, but you cannot read the code inside the function that generates the outputs.

 

This function is called with the following syntax:

 

>> [V,I,t,spikeTimes] = LIFmodel(amp); 

 

LIFmodel.p takes a single input argument, amp,  which sets the amplitude of the square current pulse in units of nanoAmperes (nA). Upon being called, it returns 4 arguments V, I, t and spikeTimes. The first two outputs are column vectors of the same size containing the following quantities as a function of time:  V, the membrane potential in millivolts (mV); I, the current amplitude in nanoAmperes (nA). The third output is  another vector, t, of the same size containing the time stamps for  V  and I in milliseconds (msec). The last argument, spikeTimes, returns a list of time stamps at  which spikes (action potentials) were generated. Keep in mind that the size of spikeTimes will vary depending on the amplitude of the current injected. In particular, if there are no spikes it will be an empty matrix.  

 

We know that the behavior of an LIF neuron is determined by 5 parameters: the resting membrane potential (EL), the spiking threshold (Vthresh), the reset membrane potential (Vreset), the total membrane resistance (Rm), and the total membrane capacitance (Cm).

 

Your job is to use plots and/or calculations to estimate each of these parameters that are obscured in the provided function. Keep in mind that you will have to run the function with different values of amp, to recover all 5 parameters.

 

Using relevant plots and or calculations complete these tasks. Make sure to justify all your answers properly and provide proper units for your estimates.  

 

1.      Estimate EL  for this neuron.  = -75 mV      set amp = 0

2.      Estimate Vthresh  for this neuron. = -50 mV       set amp =        0.275

3.      Estimate Vreset  for this neuron.   = -70 mv        set amp = 0.275

4.      Estimate R m for this neuron.                   = at Ith = Ie = 0.25 nA, fr = 0                     dV = IR, V= 20mV R=100 M Ohms

5.      Estimate the time constant for this neuron.                                                       = 20.05           ms 

6.      Estimate C m for this neuron.                    = 0. 2005 nF                      tau = RC

7.      Estimate the rheobase     [1] of this neuron.                       = 0.25 nA       

8.      Plot the firing rate   [2] in Hertz (Hz) as a function of the injected current          amplitude. This is called an f-I curve. Describe the different regions of this        curve.      

9.      From the previous plot, estimate           Cm. = 0.2034 nF

                                                              

                 Problem 2: Numerical simulation of an LIF neuron                                                                                                                                                                            

In the previous problem you estimated the 5 parameters that control the behavior of an            LIF neuron. In this problem, your job is to write MATLAB code to implement an LIF   neuron model with those parameters. This will give you an opportunity to verify your  

parameter estimates.

You can use the RCpassive.m from PSET 1 as a starting point. Don’t forget to include a variable to keep track of the spike times. Also, it might be a good idea to wrap this code in a function which similar calling syntax as LIFmodel.p. When you have your code working, complete the following tasks:

 

1.      Plot the f-I curve of this neuron. 

2.      On the previous plot overlay the f-I curve obtained in problem 1 question 8. Comment on any similarities or discrepancies between the two curves. 

 

Problem 3: The role of sodium channel inactivation on the refractory period in the Hodgkin-Huxley model 

 

The gating variables, m,h and n, obey this first order linear differential equation:  

dx

τ x (V )  dt = x ∞ (V )− x,

  

where x is an index to reflect any of those gating variables. It is important to notice that the time constant and steady state value depend on the membrane potential:

                                                                                    1                                                               α x (V )       

                                   τ x (V ) =  α x (V )+ β x (V )                   x ∞ =  α x (V ) + β x (V )

 The voltage dependency is through the rate constants   α x  and   β x which have the following mathematical expressions:

           

αh = 0.07exp(−0.05(V + 70)))        βh = 1 + exp(−0.1(1 V + 40))

0.1(V + 45)

                                    αm  =               β m = 4exp(−0.0556(V + 70))           

1−exp(−0.1(V + 45))

0.01(V + 60)

αn  =                β n = 0.125exp(−0.0125(V + 70))

                                                  1−exp(−0.1(V + 60))

 

In these formulae the membrane potential is in mV and the rate constants in ms-1. In the folder HHmodel you will find a set of six functions that implements each of these 6 rate constants.  

           
The file HH.m in the same folder implements a simulation of a Hodgkin-Huxley neuron under a current clamp experiment. Current clamp means that we control the external current applied and measure the membrane potential in response to our manipulation. Note that this is different from voltage clamp in which we control the membrane potential and measure the ionic current.  

In particular, this script is set to apply a depolarizing external current of 10 µA/cm2  for 10 milliseconds. If you run this file you will see three plots: the membrane potential as a function of time, the external current applied as a function of time, and the gating variables as a function of time.  

Your job is to use those files and make modifications as instructed to study the role of sodium channel inactivation and its effect on the refractory period.

           

1.      Plot the steady state of the inactivation gating variables (h      ∞) as a function q1. h vs V.jpg file has the solution    of membrane potential. Restrict your plot to the range -100 to 100 mV.             

2.      Make a copy HH.m and call it               HHpaired.m . Modify this file so that it                                                   applies two current pulses, each of amplitude 3                          µA/cm        2 and a duration of                                                    10 milliseconds. The pulses should have an inter-pulse interval of 10                                                                                                      milliseconds. This is they are separated by 10 milliseconds measured                                                           from end of first pulse to the start of second pulse. Look at the plots and                                                                                                                                                                                comment on what you see.                                                   

3.      Modify the code so that the inter-pulse interval is now 2                                                  milliseconds but                                             keeping the amplitude at 3 µA/cm2 . Look at the output plots, what has                                                                                                    changed from the previous case? What is the role of the inactivation                                                                                                           gating variable in this behavior?                                      

4.      Modify the file alphah.m            so that the rate is double the value of the                                               original                                                                       function for all values of membrane potential. In a single plot                                                                 overlay the                                    h∞ versus membrane potential curves for the original and                              modified                                                        

        version. What is the effect of this change in    ah   on h∞ ? q3hv.fig file has the graphs from this question. The 0.6                                                                                                                                                                       

      peaking graph is original and 0.8 peaking is modified     5. Redo question 3 but using the modified    alphah.m   function. Comment on      the changes you see. Relate this to the sodium channel disorders   discussed

in lecture. The neuron is not able to spike for the 2nd pulse when rate is original. If the rate is  

doubled then the neuron can spike but there is a delay wrt the time of pulse. Hence                                                               the body won't be able to react and the temporary paralysis discussion from lecture                                                                                                                      

          Problem 4: Hyperpolarization can excite neurons                                                                                                                                
                                                                                                                               
In this problem we want to study the effects of prolonged hyperpolarization on the          behavior of the Hodgkin-Huxley model. In order to do so, make a copy of     HH.m and name it HHhyper.m. Modify this code such as to deliver a square pulse of current that starts a 10 ms and finishes at a 100 ms. Set the amplitude to -10 microamperes/cm     2. Simulate for 150 ms total.

                                                                                                                                                                      Inactivation is almost                                   

          After running your simulation answer the following questions:                                                                                                                                                                                                                  none during -10 nA. 

                                                                                                                                                                               and Na/K  gates are

1.      Plot V, m, h and n as functions of time. You will see the neuron spiking at the

                    end of the simulation                                                                                                         almost            closed. The  

                                                                                                                                                                         time it takes for                                                     

2.      In a few sentences provide a rationale for why this is happening. Construct

                               your argument in terms of the gating variables m, h and n. Hint: It will be    inactivation to occur

                                                                                                                           

                      helpful to consider the time scales involved in the activation/inactivation of         allows sodium/

                                                                 

                    the gating variables (problem 3).                                                                             potassium gates to

3.      In a few sentences discuss whether (and why or why not) a leaky integrate open which leads to a and fire neuron can exhibit this behavior. spike towards the end

LIF cannot exhibit this behavior as this is the property of the gating variables. Negative     current cannot be injected in a LIF as it doesn't involve this opening/closing of  gates. There is

nothing to account for more hyperpolarization caused by negative current                                                          4

MIT OpenCourseWare 

https://ocw.mit.edu/ 

 

9.40 Introduction to Neural Computation 

 

For information about citing these materials or our Terms of Use, visit: https://ocw.mit.edu/terms.  

 

[1] The current amplitude such that the steady-state is equal to V                  thresh    

 
[2] We are defining the rate the same way as we did in class: as the reciprocal of the time elapsed between two consecutive spikes.