AMS/BIO 332: Computational Modeling of Physiological Systems
– Matlab Project 3 –
加拿大Matlab代写 Remember to include: 1. your code (as an appendix) and 2. a statement of collaboration. Simulate tmax = 200 ms
*** Remember to include: 1. your code (as an appendix) and 2. a statement of collaboration. ***
Part I: HodgkinHuxley model. 加拿大Matlab代写
 (Main simulation) Simulate tmax = 200 ms evolution of one HH neuron as defifined in Box 1. Use the forward Euler algorithm with time step dt = 0.01 ms. Stimulate the neuron from time t0 = 40 ms with a constant current of amplitude 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 fifigure but in difffferent subplots:
i) the membrane potential vs. time,
ii) the total membrane current, Im = GL(V − VL) + GNam3h(V − ENa) + GKn4(V − EK), vs. time,
iii) the sodium conductance (=GNam3h) vs. time,
iv) the potassium conductance (=GKn4 ) vs. time, and
v) the current Ie time.
Include the plot in your report. 加拿大Matlab代写
Matlab tip ⇒ Use figure(1) and then start each plot with subplot(5,1,x), where x ranges from 1 to 5.
2.(Spike detection) Endow your code with a mechanism to detect each ‘spike’ (=action potential). To do so, 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. Now, enlarge the previous plot from Exercise 1 to illustrate the dynamics of one single action potential, marking the time at which Vspk was reached. Include the plot in your report.
Matlab tip ⇒ Use xlim([t1 t2]) in each subplot to zoom in all subplots equally; choose as t1 and t2 the time just before and just after the action potential, respectively. ✷

(Action potential threshold) Determine the threshold current to elicit a single action potential. Start from a very small value of the input current (e.g., I0 = 2 pA) and increase it until you fifind the value of the current I1 at which the neuron starts fifiring one action potential only, and then stops fifiring. Report the value I1 in your report.
 (Rheobase current) By increasing the current beyond the value I1, eventually 2, 3 or more action potentials will be generated, but then the neuron stops fifiring again. By increasing the current even more, eventually repetitive fifiring ignites (i.e., the neuron doesn’t stop fifiring). Find Irh, defifined as the smallest value of the current above which repetitive fifiring ignites. Report the value Irh in your report.
 (f − I curve) The previous exercise seems to imply that the HH model will switch from having zero fifiring rate to having a fifinite fifiring rate at a critical value Irh. Run Exercise 1 again, this time for at least 40 values of the current with increasing amplitude, starting from values a bit lower than I = Irh. For each run, estimate the fifiring rate in that run, and then plot the currentfrequency response function once all fifiring rates have been estimated (see Box 2 below for the relevant defifinitions).
Important: Assign zero fifiring rate in the case of no spikes; and assign zero fifiring rate every time the neuron emits a few action potentials and then stops fifiring. In all other cases (repetitive fifiring) the fifiring rate will be a positive number. When plotting the fifiring rate vs. the amplitude of the input current, make sure to include the simulations which resulted in zero fifiring rates. ✷
Include the plot in your report and answer the following questions: 加拿大Matlab代写
i) Describe brieflfly the salient characteristics of this f − I
ii) Via the f − I curve, can you get an independent estimate of Irh? Can you locate Irh on the plot?
iii) Imagine you wanted to use this model neuron to encode a range of input currents into the neuron’s output fifiring rate. Does this f − I curve suggest that this would be an effiffifficient way to do so? Why or why not?
Box 2: Firing rate estimation
The Currentfrequency response function (or ‘f − I curve’ for short) is the curve depicting the fifiring rate as a function of the amplitude of the input current Ie.
You can calculate the fifiring rate in either of two ways:
i) as the number of spikes per second between t0 and tmax (make sure to use a longer simulation time tmax, e.g., 2000 ms, to obtain better estimates of the fifiring rates);
ii) as the inverse of the mean interspike interval (ISI), where the ISI is the interval of time between two successive action potentials.
Use method ii) as it is more accurate (if you try out both methods you will be able to notice the difffference. Can you tell what is the reason behind this difffference?). Note that either method requires the use of the spike detection mechanism developed in Exercise 2 (make sure to count each spike only once!).
Part II: Leaky IntegrateandFire neuron model. 加拿大Matlab代写
Box 3: The Leaky IntergrateandFire model
Consider the leaky integrateandfifire (LIF) neuron model as defifined in Box 3.
 (Main simulation) Repeat Exercise 1 of Part I for the LIF neuron: replace the equations for the HH model with the equations for the LIF neuron (remember to include the boundary conditions for the emission of a spike) and repeat the same exercise. This time, use the forward Euler algorithm with time step dt = 0.1 ms and choose a value for Ie = 1.1 nA. Plot, in the same fifigure but in two difffferent subplots: the membrane potential vs. time and the external current vs. time. Include the plot in your report.
Note The boundary conditions amount to a mechanism of spike initiation followed by clamping the membrane potential to the reset value for some time. In the LIF neuron, the time of spike initiation is also the time of spike detection (compare with Exercise 2 of Part I). If you found a way to detect each spike only once in Exercise 2 of Part I, the same code can be used to implement the refractory period here. But remember that in this case V = Vr during the refractory period.

(f −I curve) Repeat Exercise 5 of Part I for the LIF neuron. Use at least 40 difffferent values for the current. Make sure to use currents which elicit a range of fifiring rates from zero to at least 60 or 70 spikes/s. Plot the f − I curve and include the plot in your report. Answer the following questions in your report: 加拿大Matlab代写
i) What type of f − I curve do you observe in this case? Is this difffferent from the f − I curve of the HH model? In what way?
ii) Is it possible, in the LIF model, to have only a few action potentials and then stop fifiring?
iii) From the f − I plot, locate the rheobase current (analogous to Irh in Part I) and compare it with the theoretical rheobase current given by Irh LIF = GL(Vspk − VL). Do they match?
 (Running time) Answer the following questions:
i) Do your simulations of the HH model work if you use dt = 0.1 ms? And if you use dt = 0.05 ms? Include
one plot in your report.
ii) Do your LIF neuron results change signifificantly when you use dt = 0.2 ms instead of dt = 0.1 ms?
iii) Run both the LIF and the HH model with dt = 0.01 ms and dt = 0.002 ms, and measure how long each simulation takes (see Box 4). Log these values in your report. Is there an appreciable difffference in running times between the HH and the LIF model neurons, at parity of dt? Which one runs faster? Is there a big difffference?
Box 4: Estimate of running time 加拿大Matlab代写
You can estimate how long your simulations take with the tic and toc instructions:
tic (start simulation code) ... (end simulation code) toc
Matlab will output a message such as:
Elapsed time is 2.451381 seconds.
Note: make sure to use the same simulation time tmax with both models.