AMS/BIO 332: Computational Modeling of Physiological Systems - Matlab Project 3 - Part I: Hodgkin-Huxley model. Consider the following Hodgkin-Huxley (HH) neuron model, dV = GL (V VL ) GN a m3 h(V EN a ) GK n4 (V EK ) + Ie , (1) dt where V is the membrane potential, Ie is an external current and the gating variables m, n and h obey the rst-order ODEs dm = m (V )(1 m) m (V )m, (2) dt dh = h (V )(1 h) h (V )h, (3) dt dn = n (V )(1 n) n (V )n, (4) dt with transition rates (V ) and (V ) given by C m (V ) = 0.1(V + 40) ; 1 e0.1(V +40) m (V ) = 4e0.0556(V +65) , (5) h (V ) = 0.07e0.05(V +65) ; h (V ) = 1 , 1 + e0.1(V +35) (6) n (V ) = 0.01(V + 55) ; 1 e0.1(V +55) n (V ) = 0.125e0.0125(V +65) . (7) The values of the remaining parameters are: GN a = 400 nS, GK = 200 nS, GL = 2 nS, EN a = 99 mV, EK = 85 mV, VL = 65 mV, C = 2 pF and Ie as specied below. Set the initial conditions as V (0) = VL , m(0) = m (V (0)), h(0) = h (V (0)), and n(0) = n (V (0)), where m (V ) = m (V ) , m (V ) + m (V ) (8) and analogous for h (V ) and n (V ). 1. Simulate tmax = 200 ms evolution of one HH neuron dened as above. Use the forward Euler algorithm with time step dt = 0.01ms. Stimulate the neuron from time t0 = 40 ms with a constant current of magnitude I0 = 200 pA and lasting until tmax (i.e., Ie = 0 for t < 40 ms and Ie = I0 for t > 40 ms). Plot in the same gure but in dierent subplots: i) the membrane potential vs. time, ii )the total membrane current Im = GL (V VL ) + GN a m3 h(V EN a ) + GK n4 (V EK ), iii) the sodium conductance (=GN a m3 h) vs. time, iv) the potassium conductance (=GK n4 ) vs. time, and v) the current Ie vs. time (use figure(1) and then start each plot with subplot(5,1,x), where x ranges from 1 to 5). Include the plot in your report. 2. i) Decrease the current from I0 = 200 pA until you nd the value of the current I1 at which the neuron starts ring one action potential only, and then stops ring. Enlarge the plot to illustrate the dynamics of one single action potential, marking the temporal point at which Vspk was reached. (Use xlim([t1 t2 ]) in each subplot to zoom out all subplots equally; choose as t1 and t2 the time just before and just after the action potential, respectively). Include the plot in your report. ii) By increasing this current slightly, 2, 3 or more action potentials will be generated, before the neuron stops ring again. By increasing the current even more, eventually repetitive ring ignites (i.e., the neuron doesn't stop ring). Find the smallest value of the current I2 , above which repetitive ring ignites. Note that to carry out this exercise, you must decide on a critical value Vspk at which you might say that a spike has been emitted (values between 30 and 0 mV should be good choices. Make sure to count each spike only once!). 1 3. The previous exercise seems to imply that the HH model will switch from having zero ring rate to having a nite ring rate at a critical value I2 (found in the previous exercise). Run Exercise 1 again, this time for 30 to 40 values of the current with increasing magnitude, starting from values a bit lower than I = I2 . Include in your code a calculation of the ring rate (=the number of spikes per second between t0 and tmax ) for each value of the current. i) Plot the current-frequency response function (the 'f I curve' for short), dened as the ring rate as a function of the magnitude of the input current Ie . Include values of the current at which the neuron doesn't re; assign zero ring rate to those cases, as well as to those cases where no repetitive ring is observed (make sure to count each spike only once!). Include the plot in your report and describe briey the salient characteristics of this f I curve. ii) With this method, can you get a better estimate of I2 compared to the value established in Exercise 2 by sampling many points around I2 ? (You may want to increase the simulation time tmax to obtain better estimates of the ring rates). iii) Imagine you wanted to use this model neuron to encode a range of input currents into the neuron's output ring rate. Does this f I curve suggest that this would be an ecient way to do so? Why or why not? Part II: Leaky Integrate-and-Fire neuron model. Consider the leaky integrate-and-re (LIF) neuron model, whose subthreshold membrane potential V obeys dV = GL (V VL ) + Ie , (9) dt where Ie is an external current (values specied below), C = 1 nF, GL = 50 nS, and VL = 65 mV. This model is completed by the following boundary conditions: a spike is said to occur whenever V hits a threshold = 45 mV, after which it is reset to Vr = 65 mV for an absolute refractory period arp = 2 ms. Afterwards, the dynamics of Eq. 9 resumes. Set the initial condition as V (0) = VL . C 1. In your code for Part I, Exercise 3, replace the equations for the HH model with the equations for the LIF neuron as dened above (remember to include the boundary conditions for the emission of a spike; this is the trickiest part) and repeat the same exercise. This time, use the forward Euler algorithm with time step dt = 0.1 ms and use at least 40 dierent values for the current. Make sure to use currents which elicit a range of ring rates from zero to at least 60 or 70 spikes/s. i) Plot the f I curve and include the plot in your report. From the f I plot, try to locate the rheobase current (=the current at which the rst spikes appear, analogous to I2 in Part I) and compare it with the theoretical rheobase current given by ILIF = GL ( VL ). Do they match? 2 ii) What type of f I curve do you observe in this case? Is this dierent from the f I curve of the HH model? In what way? Is it possible, now, to have only a few action potentials and then stop ring? LIF iii) Next, choose a value for Ie larger than I2 and plot, in the same gure but in dierent subplots: the current vs. time and the membrane potential vs. time. Include the plot in your report. 2. Play around with dt: i) what happens if you use dt = 0.05 ms in the HH model? Do you have the same problem with the LIF neuron? ii) Try also dt = 0.1 ms and dt = 0.002 ms with both models. Explain what you observe and why. iii) Do your LIF neuron results change signicantly when you use dt = 0.2 ms instead of dt = 0.1 ms? iv) Use 'tic' and 'toc' to measure how long each simulation takes (in your code, invoke tic at the beginning of your simulation and toc at the end). Log these values in your report. Is there an appreciable dierence in simulation times between the HH and the LIF model neurons at parity of dt? Note The exercises of Part II are almost identical to those in Part I, apart from the equations to integrate and the appropriate boundary con ditions. Once the code for Part I has been written, minimally modify it to replace the model equations with the new ones, and the new results should follow by using almost exactly the same code (don't forget to reset the membrane to Vr after a spike in the LIF neuron). Note also that nding the time at which a spike is emitted is trivial in the LIF neuron: this is the time at which the threshold is reached. As a nal word of caution, remember that now you have to reset V after a spike and wait until the refractory period is over: implementation of this requires a bit of thinking and a few lines of new code compared to the simulations in Part I. 2 2