Biosequences .. Software .. Molbio soft .. Network News .. FTP

# [Neuroscience] unbounded voltages in active compartmental model

dumb_founded andreajagger_8 at hotmail.com
Fri Dec 9 17:54:30 EST 2005

```I am using matlab to simulate a 10 compartment model with active
Hodgkin Huxley conductances.  The problem I encounter is that the
voltages in all the compartments become infinite.  My alpha and beta
functions are correct as my separate Hodgkin Huxley model for one cell
works correctly.  The code for my compartmental model is below.  I
myself couldn't find anything egregious.

_______________________________________________________________________

clear

g_K=36;
g_Na=120;
g_Cl=.3;

V_Na=115;
V_K=-12;
V_Cl=10.613;

NCMPT = 10;

TMAX = 100;
DT = 0.025;
t = 0:DT:TMAX;

%
% specific resistances and capacitances
%
spRM = 10.0; 	% kOhm-cm^2
spCM = 1.0;	% uF/cm^2
spRA = 0.1;	% kOhm-cm

xc = 0.05 + 0.1*(0:NCMPT-1);

ISTEP = 10 ;

radius = [ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2;
2, 2, 2, 2, 2, 1, 1, 1, 1, 1;
2, 2, 2, 2, 2, 4, 4, 4, 4, 4;
1, 1, 1, 1, 1, 2, 2, 2, 2, 2;
4, 4, 4, 4, 4, 2, 2, 2, 2, 2];

for ca = 1:1

len = .1 * sqrt(spRM/(2*spRA))*sqrt(r);		% make each cmpt 0.1
lambda
sa = 2*pi*r.*len;					% surface area
xa = pi*r.^2;					% cross-sectional area

RM = spRM./sa;					% compute absolute RM, CM, RA
CM = spCM.*sa;
RA = (spRA*len)./xa;

v(1,1:NCMPT) = 0;					% initialization

n(1,1:NCMPT)=alpha_n(v(1,1))/(alpha_n(v(1,1))+beta_n(v(1,1)));
m(1,1:NCMPT)=alpha_m(v(1,1))/(alpha_m(v(1,1))+beta_m(v(1,1)));
h(1,1:NCMPT)=alpha_h(v(1,1))/(alpha_h(v(1,1))+beta_h(v(1,1)));

for i=2:length(t)                               	% MAIN LOOP

% first compartment
vdot(1) = ISTEP - v(i-1,1)/RM(1)...
- (v(i-1,1)-v(i-1,2))/((RA(1)+RA(2))/2)...

-g_K*n(i-1,1)^4*(v(i-1,1)-V_K)-g_Na*m(i-1,1)^3*h(i-1,1)*(v(i-1,1)-V_Na)-g_Cl*(v(i-1,1)-V_Cl);

% middle compartments
for c=2:NCMPT-1;
vdot(c) = -v(i-1,c)/RM(c) ...
-(v(i-1,c)-v(i-1,c-1))./((RA(c)+RA(c-1))/2) ...
-(v(i-1,c)-v(i-1,c+1))./((RA(c)+RA(c+1))/2)...

-g_K*n(i-1,c)^4*(v(i-1,c)-V_K)-g_Na*m(i-1,c)^3*h(i-1,c)*(v(i-1,c)-V_Na)-g_Cl*(v(i-1,c)-V_Cl);

end

% last compartment
vdot(NCMPT) = -v(i-1,NCMPT)/RM(NCMPT) ...

-(v(i-1,NCMPT)-v(i-1,NCMPT-1))/((RA(NCMPT)+RA(NCMPT-1))/2)...

-g_K*n(i-1,NCMPT)^4*(v(i-1,NCMPT)-V_K)-g_Na*m(i-1,NCMPT)^3*h(i-1,NCMPT)*(v(i-1,NCMPT)-V_Na)-g_Cl*(v(i-1,NCMPT)-V_Cl);

v(i,:) = v(i-1,:) + (vdot./CM)*DT;		% Euler integration

n_deriv=alpha_n(v(i,1))*(1-n(i-1,1))-beta_n(v(i,1))*n(i-1,1);
m_deriv=alpha_m(v(i,1))*(1-m(i-1,1))-beta_m(v(i,1))*m(i-1,1);
h_deriv=alpha_h(v(i,1))*(1-h(i-1,1))-beta_n(v(i,1))*h(i-1,1);
n(i,1)=n(i-1,1)+n_deriv*DT;
m(i,1)=m(i-1,1)+m_deriv*DT;
h(i,1)=h(i-1,1)+h_deriv*DT;

for c=2:NCMPT-1
n_deriv=alpha_n(v(i,c))*(1-n(i-1,c))-beta_n(v(i,c))*n(i-1,c);
m_deriv=alpha_m(v(i,c))*(1-m(i-1,c))-beta_m(v(i,c))*m(i-1,c);
h_deriv=alpha_h(v(i,c))*(1-h(i-1,c))-beta_n(v(i,c))*h(i-1,c);
n(i,c)=n(i-1,c)+n_deriv*DT;
m(i,c)=m(i-1,c)+m_deriv*DT;
h(i,c)=h(i-1,c)+h_deriv*DT;
end

n_deriv=alpha_n(v(i,NCMPT))*(1-n(i-1,NCMPT))-beta_n(v(i,NCMPT))*n(i-1,NCMPT);

m_deriv=alpha_m(v(i,NCMPT))*(1-m(i-1,NCMPT))-beta_m(v(i,NCMPT))*m(i-1,NCMPT);

h_deriv=alpha_h(v(i,NCMPT))*(1-h(i-1,NCMPT))-beta_n(v(i,NCMPT))*h(i-1,NCMPT);
n(i,NCMPT)=n(i-1,NCMPT)+n_deriv*DT;
m(i,NCMPT)=m(i-1,NCMPT)+m_deriv*DT;
h(i,NCMPT)=h(i-1,NCMPT)+h_deriv*DT;

end

plot(t, v);

end

```