1D Cable Model

Taking the code below how do I fix the plote of Vm(x) at 11 time points t = 0,2,4,...20? I tried doing it as t = [0:2:20], but it kept giving me an error. Then how do I plot curves of Vm(t) at 11 equally spaced x points along the cable. The code I have for this is also giving me errors. Then once that is resolved how do I make a plot of time constants of Vm(t) as a time to reach (1 - (1/e)) of the final Vm value at 11 points along the cable? I have how to find time constants of Vm(t), but I'm not sure how to plot multiple at 11 points along the cable.

% Model parameters 
a = 0.001;  % cable radius, cm
L = 1;      % cable length, cm
Rm = 5000;  % membrane resistivity, Ohm*cm^2
Cm = 1;     % membrane capacitance, uF/cm^2
Ri = 500;   % intracellular resistivity, Ohm*cm
E = 5;      % shock field, V/cm
T = 50;     % duration of integration, ms
dt = 0.020; % time integration step, ms
dx = 0.0100;% space integration step, cm

tau = Rm*Cm*0.001;              % membrane time constant, ms
lambda = sqrt((a*Rm)/(2*Ri));   % space constant, cm

nt = (T/dt) + 1;            % number of time points
nx = round(L/dx)+1;         % number of nodes
x = -L/2 : dx : L/2;        % coordinates of nodes

mesh = dt*lambda^2/(tau*dx^2);      % stability requirement mesh < 0.5

% Simplifying equation mathematics for loop 
kt = dt/tau;            
kx = lambda^2/dx^2;

% Preallocation 
vm = zeros(nt,nx); 
time = zeros(1,nt);
% Start loop at 2 because of initial condition vm(1,i) = 0
    for j = 2 : nt              % time loop
        for i = 2 : (nx-1)      % space loop for internal nodes; Euler method
            vm(j,i) = vm(j-1,i) + kt*(kx*(vm(j-1,i+1) - 2*vm(j-1,i) + vm(j-1,i-1)) - vm(j-1,i));
        end
        % Updated boundary conditions 
        vm(j,1) = (4*vm(j,2) - vm(j,3) - 2* E*dx)/3;       % vm at left edge from boundary conditions
        vm(j,nx) = (4*vm(j,nx-1) - vm(j,nx-2) + 2*E*dx)/3; % vm at right edge from boundary conditions
        time(j+1) = time(j) +  dt;
    end    
vm = 1000*vm;                   %  convert from V to mV

% Analytical solution
vmAna = 1000*E*lambda*(sinh(x/lambda)/cosh(L/(2*lambda)));


% Plot curves Vm(x) at 11 time points where t = 0,2,4,...20 ms
% Gives an error if try to do t as t = [0:2:20];
t = [2:2:20];
plot(x,vm(t,:));
xlabel('x (cm)');       
ylabel('Vm (mV)');
title('Figure 2: L = 1 cm');
% title('Figure 2: L = 0.1 cm');

% Plot curves Vm(t) at 11 equally spaced x points along the cable 
xx = linspace(-0.5,0.5,11);
plot(time,vm(:,xx));
xlabel('t (ms)');       
ylabel('Vm (mV)');
title('Figure 2: L = 1 cm');
% title('Figure 2: L = 0.1 cm');

% Plot time constants of Vm(t) change as time to reach (1-1/e) of the final
% Vm value at 11 points along the cable


% find time constant of Vm(t) as a time to reach (1-1/e)*Vm_final  
j = 1;
level = 1-1/exp(1);
while (abs(vm(j)) < level*vm(nt))    
    j = j+1;
end;
taunum = j*dt;

[darova]

I dont have any error with this

t = [2:2:20];

Try to clear variables before your script

clear

• % Plot time constants of Vm(t) change as time to reach (1-1/e) of the final. Vm value at 11 points along the cable

I don't understand this. You want to find t=1-1/e? But 1-1/e is bigger than 0.5 (your max time)

>> 1-1/exp(1)
ans =
    0.6321

[Kate Heinzman]

No it gave me an error when I did it as t = [0:2:20];. And to your second comment yes I want to find where t is equal to 0.632 at 11 points along the cable. Did you have any answers about the second question asked about plotting Vm(t) at 11 equally spaced x points along the cable?

[darova]

• No it gave me an error when I did it as t = [0:2:20]

I think it's because of 0 (index starts from 1)

[Kate Heinzman]

Ok thanks I fixed it. Can you please answer my other two questions?

[darova]

If i understood you correctly: you want to find point where t=1-1/e

Solution

imax = find(t>1-1/exp(1),1,'first');    % index

Then to some plots or something

plot(t(imax),vm(imax,:),'or')