Simulation of Hodgkin-Hucley model for convolution part

 

%Author: Minster

%Simulation of Hodgkin-Hucley model for convolution part.

 

function rk=rk

%assigned initial values

Vo=-60; %initial volatage

cc=0.005; %Coefficient

t=10;%time-10sec

T=0:cc:t;%X axis

[st en]=size (T);

 

% initial condition for functions

an(1) = An(0);

bn(1) = Bn(0);

%an and bn determine the rate of transfer of K+ throught the n-gate from

%the outside to the inside and from the inside to the outside respectively.

am(1)=Ah(0);

bm(1) = Bm(0);

%am and bm are rate constants that represents the rate of transfer of Na+

%through the m-gate

 

ah(1)=Ah(0);

bh(1)=Bh(0);

%an and bn are rate constants that represent the rate of transfer through

%the h-gate

 

v(1) = Vo;

n(1) = an(1)/(an(1) + bn(1));

m(1)=am(1)/(am(1)+bm(1));

h(1)=ah(1) /(ah(1)+bh(1));

gNa(1) = GNa(m(1),h(1));

%gNa conductance of Na+ through the cell membrane

 

gK(1) = GK(n(1));

%gK conductance of K+ through cell membrane

 

% the 4th order runge kutta method

 

for i=2:en,

p1=(dV(v(i-1),n(i-1),m(i-1),h(i-1)))*cc;

q1=(dN(v(i-1),n(i-1),m(i-1),h(i-1)))*cc;

r1=(dM(v(i-1),n(i-1),m(i-1),h(i-1)))*cc;

s1=(dH(v(i-1),n(i-1),m(i-1),h(i-1)))*cc;

 

p2=(dV(v(i-1)+(p1/2),n(i-1)+(q1/2),m(i-1)+(r1/2),h(i-1)+(s1/2)))*cc;

q2=(dN(v(i-1)+(p1/2),n(i-1)+(q1/2),m(i-1)+(r1/2),h(i-1)+(s1/2)))*cc;

r2=(dM(v(i-1)+(p1/2),n(i-1)+(q1/2),m(i-1)+(r1/2),h(i-1)+(s1/2)))*cc;

s2=(dH(v(i-1)+(p1/2),n(i-1)+(q1/2),m(i-1)+(r1/2),h(i-1)+(s1/2)))*cc;

 

p3=(dV(v(i-1)+(p2/2),n(i-1)+(q2/2),m(i-1)+(r2/2),h(i-1)+(s2/2)))*cc;

q3=(dN(v(i-1)+(p2/2),n(i-1)+(q2/2),m(i-1)+(r2/2),h(i-1)+(s2/2)))*cc;

r3=(dM(v(i-1)+(p2/2),n(i-1)+(q2/2),m(i-1)+(r2/2),h(i-1)+(s2/2)))*cc;

s3=(dH(v(i-1)+(p2/2),n(i-1)+(q2/2),m(i-1)+(r2/2),h(i-1)+(s2/2)))*cc;

 

p4=(dV(v(i-1)+p3,n(i-1)+q3,m(i-1)+r3,h(i-1)+s3))*cc;

q4=(dN(v(i-1)+p3,n(i-1)+q3,m(i-1)+r3,h(i-1)+s3))*cc;

r4=(dM(v(i-1)+p3,n(i-1)+q3,m(i-1)+r3,h(i-1)+s3))*cc;

s4=(dH(v(i-1)+p3,n(i-1)+q3,m(i-1)+r3,h(i-1)+s3))*cc;

 

v(i)=v(i-1) + ((p1+2*p2+2*p3+p4)/6);

n(i) = n(i-1)+((q1+2*q2+2*q3+q4)/6);

m(i) = m(i-1)+((r1+2*r2+2*r3+r4)/6);

h(i) = h(i-1) + ((s1+2*s2+2*s3+s4)/6);

 

gNa(i) = GNa(m(i),h(i));

gK(i) = GK(n(i));

 

end;

 

wo=pi*1;

Dt=0.005;

tt=(0:Dt:10);

L=length(tt);

tp=[2*tt(1):Dt:2*tt(L)];

hh=ustep(t).*(sin(wo*t));%abs((sin(wo*tt)))/abs((wo*tt));

yy=(Dt*conv(t,hh));

%yy(5/Dt)

z=n+ustep(2*tt-6);

%zz=(1.*(ustep(0,tt-1)-2*ustep(0,tt-2.5)+ustep(0,tt-4)));

subplot(4,1,1), plot(tt,gNa),axis( [-2, 12, -2, 1.5]), ylabel(‘x(t)’), grid

subplot(4,1,2), plot(tt,hh),axis( [-2, 12, -4, 4]), ylabel(‘h(t)’), grid

subplot(4,1,3), plot(tt,gNa, ‘:’,tt,hh,’-‘),axis( [-2, 12, -2, 1.5]), ylabel(‘x(t), y(t)’)

subplot(4,1,4), plot(tt,yy), ylabel(‘y=h*x’),axis ([-2,12, -.5, .5]),grid

for i=2:en,

 

x= abs(n(i));

power=(x^2)*0.5;

power=power+power;

end;

end;

%power = abs(sin(wo*t))^2;

fprintf(‘power of periodic signal is %f \n’,power);

%plot(t,x), hold, plot(t,h);

%title(‘potential vs. time’);

%subplot(3,1,1), plot(T, v);

%title(‘n(t)-,h(t):,m(t)– vs. time’)

%subplot(3,1,2), plot(T,n,’-‘, T,h,’:’, T,m,’–‘);

%title(‘gNa, gK vs time’)

%subplot(3,1,3), plot(T, gNa, ‘-‘, T, gK, ‘:’);

 

%function to calculate Conductance of K+ channel

function Kcondition=GK(n)

gK = 36;

Kcondition=gK*(n^4);

 

%function to calculate conductance of Na+ channel

function Nacondition = GNa(m,h)

gNa = 120;

Nacondition = gNa*(m^3)*h;

 

%function to calculate change in voltage

function chdv = dV(V,n,m,h)

 

gL = 0.3;%gL is conductance of K+ through cell membrane

Cm = 1;

VNa=-115;

VK=12;

VL=10.6;

 

chdv=(-GNa(m,h)*(V-VNa)-GK(n)*(V-VK)-gL*(V-VL))/Cm;

 

%function to calculate n gate response

 

function ngate = dN(V,n,m,h)

 

Bn = 0.125*exp(V/80);

An=(0.01 *(V+10))/(exp((V+10)/10)-1);

ngate = (An*(1-n))-Bn*n;

 

%function to calculate m gate response

function mgate = dM(V,n,m,h)

 

Am=(0.1 *(V+25))/(exp((V+25)/10)-1);

Bm=4*exp(V/18);

 

mgate=(Am*(1-m))-Bm*m;

 

%function to calculate h gate response

 

function hgate = dH(V,n,m,h)

 

Bh = 1/(exp((V+30)/10)+1);

Ah = 0.07 * exp(V/20);

 

hgate = (Ah*(1-h))-Bh*h;

 

%function to calculate the constance

function resam = Am(V)

resam = (0.1*(V+25))/(exp((V+25)/10)-1);

function resan = An(V)

resan = (0.01*(V+10))/(exp((V+10)/10)-1);

function resah=Ah(V)

resah=0.07*exp(V/20);

function resbn=Bn(V)

resbn= 0.125*exp(V/80);

function resbm=Bm(V)

resbm=4*exp(V/18);

function resbh=Bh(V)

resbh=1/(exp((V+30)/10)+1);