1D cable model

19 views (last 30 days)
Kate Heinzman
Kate Heinzman on 5 Apr 2020
Commented: darova on 6 Apr 2020
When I try to pre-allocate the vectors before the for-loop it keeps giving me an error by saying the size inputs must be integers and I am confused as to why it's giving me this error. Then if I get rid of the preallocation, it throws an error in the line Vm_num(1,j) = Vm_num(2,j) - E * dx; by saying the index in position 1 is invalid. How do I fix this? Also in calculating the steady state equation of the cable what would be the x value used there based on the given model parameters?
% Script to numerically calculate changes of Vm(x,t) produced by
% uniform-field electric shock in a 1D monodomain uniform passive cable
% model
clear all
% Model parameters
a = 0.001; % cell radius in cm
L = 1; % cable length in cm
Rm = 5000; % specific membrane resistance in ohm*cm^2
Cm = 0.001; % specific membrane capacitance in mF/cm^2
Ri = 500; % intracellular resistivity in ohm*cm
E = 5; % shock field strength in V/cm
T = 50; % shock duration in ms
dx = 0.01; % integration step in cm
dt = 20; % integration step in us
Tm = Rm * Cm; % time constant in ms
lamda2 = a*Rm / 2*Ri;
% Number of nodes in the cable is equal to the cable length divided by the
% spatial integration step
N = L / dx;
% Time iterations is the shock duration divided by the time integration
% step
TI = T / dt;
% Preallocate numerical Vm and time vector
Vm_num = zeros(N,TI);
time = zeros(1,TI);
% i is the index corresponding to the nodes of the cable
for i = 1:N
% j is the time index
for j = 1:TI
% Boundary Condition
if i == 1
Vm_num(1,j) = Vm_num(2,j) - E * dx;
% Boundary Condition
elseif i == N
Vm_num(N,j) = Vm_num(N-1,j) + E * dx;
else
Vm_num(i,j+1) = Vm_num(i,j) + (dt/Tm)*((lambda2/dx*dx)*(Vm_num(i+1,j)-2*Vm_num(i,j)+Vm_num(i-1,j))-Vm_num(i,j));
end
end
time(j+1) = time(j) + dt;
end
% Calculate the steady state of the cable equation
lambda = sqrt(lambda2);
Vm_ana = (lambda*E)*(sinh(x???/lambda)/cosh(L/2*lambda));
% Plot the final numerical solution Vm(x) and plot the analytical steady
% state solution of the cable equation for Vm(x) in the same figure
plot(time,Vm_num);
hold on
plot(time,Vm_ana);
title('Comparing Analytical and Numerical Membrane Potentials over Time');
xlabel('Time(ms)');
ylabel('Membrane Potential (mV)');
legend('Numerical Vm','Analytical Vm');
  2 Comments
Tommy
Tommy on 5 Apr 2020
Make sure your units match.
T = 50; % shock duration in ms
...
dt = 20; % integration step in us
so
>> TI = T / dt
TI =
2.5000
But TI should be an integer. With matching units:
T = 50; % shock duration in ms
...
dt = 0.02; % integration step in ms
and
>> TI = T / dt
TI =
2500
darova
darova on 6 Apr 2020
Use button for code inserting

Sign in to comment.

Answers (1)

Star Strider
Star Strider on 5 Apr 2020
You are calculating ‘N’ and ‘TI’. They may appear to be integers, however they are not. The solution to that is to force them to be integers with the fix function:
% Preallocate numerical Vm and time vector
Vm_num = zeros(fix(N),fix(TI));
time = zeros(1,fix(TI));
There are other options such as ceil, floor and round, however fix is likely most appropriate here.
Running your code after that correction reveals this one:
lamda2 = a*Rm / 2*Ri;
and changing that to:
lambda2 = a*Rm / 2*Ri;
solves it.
Then an error occurs with this line:
Vm_ana = (lambda*E)*(sinh(x???/lambda)/cosh(L/2*lambda));
and I have no idea what to do with it.
I will be happy to help with it (as well any other problems that may appear) if your code has problems after you correct the ‘Vm_ana’ calculation.
  5 Comments
Star Strider
Star Strider on 6 Apr 2020
The loop is not necessary. Use the vectorisation abilities MATLAB provides, as I did to calculate ‘Vm_ana’.
This revised code makes all the vector and matrix sizes work correctly, however there are obviously still problems that you need to sort:
% Model parameters
a = 0.001; % cell radius in cm
L = 1; % cable length in cm
Rm = 5000; % specific membrane resistance in ohm*cm^2
Cm = 0.001; % specific membrane capacitance in mF/cm^2
Ri = 500; % intracellular resistivity in ohm*cm
E = 5; % shock field strength in V/cm
T = 50; % shock duration in ms
dx = 0.01; % integration step in cm
dt = 20; % integration step in us
Tm = Rm * Cm; % time constant in ms
lambda2 = a*Rm / 2*Ri;
% Number of nodes in the cable is equal to the cable length divided by the
% spatial integration step
N = L / dx;
% Time iterations is the shock duration divided by the time integration
% step
TI = T / dt;
% Preallocate numerical Vm and time vector
Vm_num = zeros(fix(N),fix(TI));
time = zeros(1,fix(N));
% i is the index corresponding to the nodes of the cable
for i = 1:N
% j is the time index
for j = 1:TI
% Boundary Condition
if i == 1
Vm_num(1,j) = Vm_num(2,j) - E * dx;
% Boundary Condition
elseif i == N
Vm_num(N,j) = Vm_num(N-1,j) + E * dx;
else
Vm_num(i,j+1) = Vm_num(i,j) + (dt/Tm)*((lambda2/dx*dx)*(Vm_num(i+1,j)-2*Vm_num(i,j)+Vm_num(i-1,j))-Vm_num(i,j));
end
end
time(i) = (i-1) + dt;
end
% Calculate the steady state of the cable equation
lambda = sqrt(lambda2);
x = linspace(-L/2, L/2, numel(time));
Vm_ana = (lambda*E)*(sinh(x/lambda)/cosh(L/2*lambda));
% Plot the final numerical solution Vm(x) and plot the analytical steady
% state solution of the cable equation for Vm(x) in the same figure
plot(time,Vm_num);
hold on
plot(time,Vm_ana);
title('Comparing Analytical and Numerical Membrane Potentials over Time');
xlabel('Time(ms)');
ylabel('Membrane Potential (mV)');
legend('Numerical Vm','Analytical Vm');
We’ll work on ‘lambda’ when I understand how you want to define it. Just now, that’s not obvious.
Tommy
Tommy on 6 Apr 2020
Edited: Tommy on 6 Apr 2020
A few notes which I hope are helpful:
1. Don't forget about the units. Keeping the fix() is a good idea, but with a total duration of 50 milliseconds and a step duration of 20 microseconds, you should have 2500 time steps, not 2. See my comment above.
2. "Plot the final numerical solution Vm(x) and plot the analytical steady state solution of the cable equation for Vm(x) in the same figure" makes me think you are supposed to plot Vm_ana over x, as well as the final Vm_num over x. Meaning
plot(x,Vm_num(:,end));
hold on
plot(x,Vm_ana);
instead of
plot(time,Vm_num);
hold on
plot(time,Vm_ana);
3. Currently to fill Vm_num, you pick a given node (i) and loop through all times (j) for that node, then move to the next node. But the value of Vm_num(i,j+1) depends on Vm_num(i+1,j), meaning you need to know the next node's Vm_num (from the previous time). I think it makes more sense to loop through all the nodes at a given time, then move to the next time - in other words, switch the order of the loops.

Sign in to comment.

Categories

Find more on Contour Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!