Introduce a step-wise increase of intracellular resistivity

Introduce a step-wise increase of intracellular resistivity Ri from 300 Ω*cm to 3000 Ω*cm in the cable center. Construct an appropriate array of intracellular conductance values gix(x). All other parameters are the same as in the uniform case. Calculate Vi(x), Ve(x) and Vm(x). Plot all 3 variables in one figure; plot Vm(x) separately in another figure.

clear all, close all, clc
tic
% Steady-state solution of bidomain equations for Vi, Ve and Vm  
% Intracellular conductivity is function of x: gix=gix(x) 
a = 0.001;          % intracellular domain radius (cm)
b = 0.0015;         % extracellular domain radius (cm)
dx = 0.0010;        % integration step (cm)
L = 1;              % cable length (cm) 
Ri = 0.300;         % specific intracellular resistivity (kohm*cm) 
Re = 0.150;         % specific extracellular resistivity (kohm*cm)
Rm = 5.00;          % specific membrane resistivity (kohm*cm^2)
E = 5;              % shock strength, V/cm
 
ri = Ri/(pi*a^2);            % intracell. resistance per unit length (KOhm/cm)
re = Re/(pi*b^2 - pi*a^2);   % extracell. resistance per unit length KOhm/cm
rm = Rm/(2*pi*a);            % membrane resistance times unit length (KOhm*cm)
lambda = sqrt(rm/(ri + re));
 
gi = (pi*a*a)/(Ri*dx);        % intracell. conductance (mS)
ge = (pi*b^2-pi*a^2)/(Re*dx); % extracell. conductance (mS)
gm = (2*pi*a*dx)/Rm;          % membrane conductance (mS)
 
istim = 1000*E/re;   % stimulation current (uA) at cable edges due to shock field E
 
NX = L/dx+1;         % number of nodes
%  intracellular conductivity is function of x: 
gix = gi*ones(1,NX); % uniform conductivity
 
% -------------------------Bidomain solver---------------------------------
% Solver code is provided by Dr. Pollard
% ----sparse matrix initialization-----
maxnz = 8*NX-4; 
rowvec = zeros(1,maxnz); 
colvec = zeros(1,maxnz); 
coefvec = 0.0001*ones(1,maxnz);
% intracellular assembly
imaxnz = 0; 
rownum = 1;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum;      coefvec(imaxnz)=  gm+gix(rownum);
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum+1;    coefvec(imaxnz)= -gix(rownum);
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum+NX;   coefvec(imaxnz)= -gm;
for rownum = 2:NX-1
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum-1;    coefvec(imaxnz)= -gix(rownum-1);
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum;      coefvec(imaxnz)=  gm+gix(rownum-1)+gix(rownum);
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum+1;    coefvec(imaxnz)= -gix(rownum);
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum+NX;   coefvec(imaxnz)= -gm;
end
rownum = NX;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum-1;    coefvec(imaxnz)= -gix(rownum-1);
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum;      coefvec(imaxnz)= gm+gix(rownum-1);
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum+NX;   coefvec(imaxnz)= -gm;
% extracellular assembly
rownum=NX+1;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum-NX; coefvec(imaxnz)=-gm;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum; coefvec(imaxnz)=gm+ge;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum+1; coefvec(imaxnz)=-ge;
for rownum=NX+2:2*NX-1
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum-NX; coefvec(imaxnz)=-gm;
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum-1; coefvec(imaxnz)=-ge;
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum; coefvec(imaxnz)=gm+2*ge;
    imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum+1; coefvec(imaxnz)=-ge;
end
rownum=2*NX;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum-NX; coefvec(imaxnz)=-gm;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum-1; coefvec(imaxnz)=-ge;
imaxnz=imaxnz+1; rowvec(imaxnz)=rownum; colvec(imaxnz)=rownum; coefvec(imaxnz)=gm+ge;
 
if(imaxnz<maxnz)
    sprintf('lower maxnz from %d to %d',maxnz,imaxnz)
end
S = sparse(rowvec,colvec,coefvec);
% spy(S);
%---- End of sparse matrix initialization ----
 
soln = zeros(2*NX,1);       % vector combining Vi and Ve
rhs = zeros(2*NX,1);        % right-hand side of the equation for stimulation current
rhs(NX+1,1) = +istim;       % current is injected in extracellualr space at one edge
rhs(2*NX,1) = -istim;       % current is withdrawn in extracellualr space at the other edge
 
% solve system of linear equations S*soln = rhs for soln
[soln,flag,relres,iter] = bicg(S,rhs,1.e-6,10000);  
 
if( flag > 0)
    sprintf('failed to converge in %d iterations',iter)
end
 
% extract Vi, Ve; calculate Vm
vi = zeros(1,NX); 
ve = zeros(1,NX); 
vm = zeros(1,NX); 
for i=1:NX
    vi(i) = soln(i);
    ve(i) = soln(NX+i);
end
vm = vi - ve;
%------------------------End of bidomain solver----------------------------
 
% Plot vi, ve, vm. Compare vm with analytical solution
x = -L/2 : dx : L/2;           % coordinates of nodes
vmTheory = 1000*E*lambda*(sinh(x/lambda)/cosh(L/(2*lambda)));
err = vm - vmTheory;                    % numerical error in mV
errrel = abs(100*(err./vmTheory));      % numerical error in %
 
fig1 = figure(1);
set(fig1, 'Position', [50 50 400 600]);
subplot(2,1,1);
plot(x,vi,x,ve,x,vm);
legend('Vi','Ve','Vm');
xlabel('x (cm)');       ylabel('mV');
 
subplot(2,1,2);
plot(x,vmTheory,x,vm,x,err);
text(x(4),vm(4),['Edge: error = ',num2str(err(1),2),' mV or ',num2str(errrel(1),2),'%']);
legend('Analytical','Numerical','Error');
xlabel('x (cm)');       ylabel('Vm (mV)');
 
toc

[darova]

Do you have a question? What is the problem?

[Baylee Sager]

I need to introduce a step-wise increase of intracellular resistivity Ri from 300 Ω*cm to 3000 Ω*cm in the cable center

[darova]

Can you make a sketch?