[Neuroscience] unbounded voltages in active compartmental model

[dumb_founded]

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];
radius = radius*1e-4;	% convert um to cm


for ca = 1:1

    r = radius(ca,:);

    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