

%
%  cab3.m
%  Steve Cox, CAAM 335
%
%  Compute and plot the potentials in an N compartment
%  resistive fibre with resting potential resulting 
%  from a DC current stimulus at node number inj
%
%  usage:    cab3(N,inj)
%  example:  cab3(16,8)
%

function cab3(N,inj)

ell = 1;	% fiber length (cm)
a = 0.025;	% fiber radius (cm)
rho_i = 30;	% cytoplasmic resistivity (Ohm cm)
rho_m = 1000;	% membrane resistivity (Ohm cm^2)
E_m = -70;	% resting potential mV

R_i = rho_i*(ell/N)/(pi*a*a);  	% axial resistance (Ohm)
R_m = rho_m/(2*pi*a*ell/N);	% lateral resistance (Ohm)

i_0 = .001;	% current stimulus (mA)

m = 2*N;	% number of edges
n = N+1;	% number of nodes

A = zeros(m,n);	% initialize incidence matrix

for j=1:N,	% build A one compartment at a time

    A(2*j-1:2*j,j:j+1) = [-1 1;0 -1];

end

b = zeros(m,1);
b(2:2:m) = -E_m;

G = zeros(m,1);	% initialize and assemble G
G(1:2:m) = 1/R_i;
G(2:2:m) = 1/R_m;
G = diag(G);

S = A'*G*A;	% assemble main matrix

f = zeros(n,1); % assemble the stimulus vector
f(inj) = i_0;
f = A'*G*b + f;	

x = S \ f;	% solve for the potentials

z = 0:ell/N:ell;	% nodal positions for plotting

plot(z,x,'x-')	% plot and label potential vs. distance
box off
xlabel('z  (cm)','fontsize',16)
ylabel('x  (mV)','fontsize',16)

