Here is Matlab code for 4th order Runge-Kutta solving of the squid axon model. To operate it, you must copy and paste into a Matlab m-file. It is self containing, and does not need any additional functions to work, other than the typical Matlab library of functions. To see how it works, look at my post entitled "Matlab Squid Axon Report".

Code:

clc %clears workspace

clear %clears any variables

%Model Constants:

Cm = 1.0; %(uF/cm^2) Membrane capacitance

ENa = 115; %(mV) Sodium Voltage

EK = -12; %(mV) Potassium Voltage

EL = 10.613; %(mV) Leakage Voltage

g_Na = 120; %(mS/cm^2) Sodium conductance

g_L = 0.3; %(mS/cm^2) Leakage conductance

g_K = 36; %(mS/cm^2) Potassium conductance

I_stim = 200; %(uA/cm^2) Amount of current stimulation for 0.05msec

Simul_time = 20; %(msec) Duration of the simulation

%Initial Conditions:

V(1) = 0; %(mV) Initial membrane voltage condition at time 0

m(1) = ((.1*25)/(exp(25/10)-1))/(((.1*25)/(exp(25/10)-1))+4); %initial value for the gating variable m

h(1) = 0.07/(0.07+1/(exp(30/10)+1)); %initial value for the gating variable h

n(1) = (.1/(exp(1)-1))/(.1/(exp(1)-1)+0.125); %initial value for the gating variable n

%Stepsize

step = 0.01; %(ms) Step size between points where slope is calculated

N = Simul_time/step; %Total number of points of data to be taken

time1(1) = 0; %initialization of time so time at t = 1 is 0 milliseconds

for num = 1:N %for loop initialization that contains calculation of all the data points

for i = 1:4 %embedded for loop to determine the variables necessary to calculate Hodgkin Huxley

if i == 1 %conditional statement for variables at the first slope

variable1 = V(num); %sets variable1 equal to the voltage at point num for the first slope

variable2 = m(num); %sets variable2 equal to the m value at point num for the first slope

variable3 = h(num); %sets variable3 equal to the h value at point num for the first slope

variable4 = n(num); %sets variable4 equal to the n value at point num for the first slope

end %end of the first slope conditional statement

if i<1 || i>4 %conditional statement for variables at the second and third slopes

variable1 = V(num)+(step/2)*K(i-1,1); %determines the voltage for the second and third slopes using the Runge-Kutta 4th order step increments

variable2 = m(num)+(step/2)*K(i-1,2); %determines the m value for the 2nd and 3rd slopes

variable3 = h(num)+(step/2)*K(i-1,3); %determines the h value for the 2nd and 3rd slopes

variable4 = n(num)+(step/2)*K(i-1,4); %determines the n value for the 2nd and 3rd slopes

end %end conditional statement for second and third slope

if i == 4 %conditional statement for variables at the fourth slope

variable1 = V(num)+step*K(i-1,1); %determines voltage for the fourth slope

variable2 = m(num)+step*K(i-1,2); %determines the m value for the fourth slope

variable3 = h(num)+step*K(i-1,3); %determines the h value for the fourth slope

variable4 = n(num)+step*K(i-1,4); %determines the n value for the fourth slope

end %end fourth slope conditional statement

for L = 1:4 %corresponds to slopes of points

if L == 1; %if for the first slope

if time1(num) > 0.05; %conditional statement accounting for the 0.05ms stimulation, so after 0.05ms, the stimulation current is zero

I_stim = 0; %current stimulation is zero after 0.05ms

I_Na(num) = g_Na*variable2^3*variable3*(variable1-ENa); %determines the value of the sodium current at point 1

I_K(num) = g_K*variable4^4*(variable1-EK); %determines the potassium current at point 1

I_L(num) = g_L*(variable1-EL); %determines the leakage current at point 1

K(i,L) = -(I_Na(num)+I_K(num)+I_L(num)-I_stim)/Cm; %determines the value of the dV/dt of the voltage equation

else %condition when time is less than 0.05msec

I_Na(num) = g_Na*variable2^3*variable3*(variable1-ENa); %determines the sodium current at point 1

I_K(num) = g_K*variable4^4*(variable1-EK); %determines the potassium current

I_L(num) = g_L*(variable1-EL); %determines the leakage current

K(i,L) = -(I_Na(num)+I_K(num)+I_L(num)-I_stim)/Cm; %determines the value of dV/dt (slope)

end %end conditional statements

end %end conditional statements

if L == 2; %if for the second slope

Alpha_m = 0.1*(25-variable1)/(exp((25-variable1)/10)-1); %finds the Alpha m H-H gating variable

Beta_m = 4*exp(-variable1/18); %finds the Beta m H-H gating variable

K(i,L) = Alpha_m*(1-variable2)-Beta_m*variable2; %finds the next K value slope

end %end condition

if L == 3; %if for the third slope

Alpha_h = 0.07*exp(-variable1/20); %finds the Alpha h HH gating variable

Beta_h = 1/(exp((30-variable1)/10)+1); %finds the Beta h HH gating variable

K(i,L) = Alpha_h*(1-variable3)-Beta_h*variable3; %finds the next slope

end %end condition

if L == 4; %4th slope

Alpha_n = 0.01*(10-variable1)/(exp((10-variable1)/10)-1); %HH gating variable

Beta_n = 0.125*exp(-variable1/80); %HH gating variable

K(i,L) = Alpha_n*(1-variable4)-Beta_n*variable4; %4th slope

end %end condition

end %end for loop

if i == 4; %conditional statement

for P = 1:4; %for loop to determine the values of k needed for finding the next points

kbar(P) = (1/6)*(K(1,P)+2*K(2,P)+2*K(3,P)+K(4,P)); %finds the weighted average of the slope voltage values

end %end for loop

V(num+1) = V(num)+step*kbar(1); %determines the value of the next voltage point

m(num+1) = m(num)+step*kbar(2); %determines the value of the next m

h(num+1) = h(num)+step*kbar(3); %determines the value of the next h

n(num+1) = n(num)+step*kbar(4); %determines the value of the next n

time1(num+1) = time1(num)+step; %increments the time so that it continues on

end %end conditional statement

end %end for loop

end %end initial for loop

I_Na(num+1) = g_Na*m(num+1)^3*h(num)*(V(num+1)-ENa); %sodium ionic current for the next point

I_K(num+1) = g_K*n(num+1)^4*(V(num+1)-EK); %potassium ionic current for the next point

I_L(num+1) = g_L*(V(num+1)-EL); %leakage current for the next point

Voltage = -100:.1:150; %voltage span for plotting

Alpha_m = 0.1.*(25-Voltage)./(exp((25-Voltage)./10)-1); %alpha m gating variable

Beta_m = 4.*exp(-Voltage./18); %beta m gating variable

Tao_m = 1./(Alpha_m+Beta_m); %Tm time constant

Alpha_h = 0.07.*exp(-Voltage./20); %alpha h gating variable

Beta_h = 1./(exp((30-Voltage)./10)+1); %beta h gating variable

Tao_h = 1./(Alpha_h+Beta_h); %Th time constant

Alpha_n = 0.01.*(10-Voltage)./(exp((10-Voltage)./10)-1); %alpha n gating variable

Beta_n = 0.125.*exp(-Voltage./80); %beta n gating variable

Tao_n = 1./(Alpha_n+Beta_n); %Tn time constant

m_bar = Alpha_m./(Alpha_m+Beta_m); %sodium current m value

h_bar = Alpha_h./(Alpha_h+Beta_h); %sodium current h value

n_bar = Alpha_n./(Alpha_n+Beta_n); %potassium current n value

figure(1) %produces the first figure

plot(time1,V) %plots voltage versus time, an action potential

xlabel('Time (ms)') %labels the x axis with time

ylabel('Membrane Voltage (mV)')%labels the y axis with voltage

title('Action Potential') %provides a title for the plot

figure(2) %produces the second figure

plot(Voltage,m_bar,Voltage,h_bar,Voltage,n_bar) %plots the voltage versus gating variables, all three on one plot

xlabel('Voltage (mV)') %labels the x axis with voltage

ylabel('Gating variables') %labels the y axis with gating variable values

title('Gating Variables vs Voltage') %provides a title for the plot

legend('m value','h value','n value') %produces a legend to label the plot curves

figure(3) %produces the third figure

plot(Voltage,Tao_m,Voltage,Tao_h,Voltage,Tao_n) %plots voltage versus time constants for all three time constants

xlabel('Membrane Voltage (mV)') %labels the x axis with voltage

ylabel('Time Constant Value') %labels the y axis with time constant values

title('Time Constants as a function of Membrane Voltage') %provides a title to the plot

legend('Tau m','Tau h','Tau n') %provides a legend to label the curves

figure(4) %produces the fourth figure

plot(time1,I_Na,time1,I_K,time1,I_L) %plots the ionic currents versus time

title('Ionic Currents vs Time') %adds a title to the plot

xlabel('Time (ms)') %labels the x axis with time

ylabel('Ionic Current (uA)') %labels the y axis with current

legend('Sodium Current','Potassium Current','Leakage Current') %provides a legend to label the plots

figure(5) %produces the final figure

plot(time1,m,time1,h,time1,n) %plots the gating variables versus time

xlabel('Time (ms)') %labels the x axis

ylabel('Gating variables') %labels the y axis

title('Gating Variables vs Time') %adds a title to the plot

legend('m value','h value','n value') %provides a legend to label the plots

Can u post a code for Myelinated Nerve as the code posted on your blog is seems to be of unmyelinated one.

ReplyDeleteso i have a question... i'm trying to understand your code, as we have the exact same project we're working on now, and i seem to understand everything except for why you said i<1 or i<4 intuitively you would think it should be the opposite. i just can't wrap my brain around how i is being incremented. if you could explain this that would be great! :D

ReplyDelete