Finding membrane time constant from a graph

I wrote a script to numerically calculate membrane potential produced by intracellular current injection in a cell model using euler integration method to compare it to the analytical membrane potential. I'm not sure how to calculate the membrane time constant from the numerical Vm(t)as a time to reach ~ 0.632 of the maximal value.

% Model parameters

a = 0.001; % cell radius in cm

Rm = 5000; % specific membrane resistance in ohm*cm^2

Cm = 0.001; % specific membrane capacitance in mF/cm^2

I0 = 0.1; % current strength in mA

T = 50; % integration duration in ms

t = 0.02; % integration step in ms

Tm = Rm * Cm; % time constant in ms

num_iter = T / t; % number of interations of time loop

% Preallocate Vm_num (numerical Vm), Vm_ana (analytical Vm), and time

% vectors for the for loop

Vm_num = zeros(1,num_iter);

Vm_ana = zeros(1,num_iter);

time = zeros(1,num_iter);

% When calculating Rm * I0 in the Euler method, Rm is in terms of the area

% of the cell so need to multiply it by the surface area to get rid of the

% cm^2 in the units. So Rm * I0 is in terms of mV instead of mV * cm^2

% Use surface area bc that is where the shock is being applied

% Surface area of a circle is 4*pi*r^2

RM = Rm / (4*pi*a^2);

for j = 1:num_iter

% initial condition: Vm(1) = 0 so start at 2nd index

Vm_num(j+1) = Vm_num(j) + (t/Tm)*((RM*I0)-Vm_num(j));

Vm_ana(j+1) = (RM*I0)*(1 - exp(-time(j)/Tm));

time(j+1) = time(j) + t;

end

% Plot numerical and analytical Vm on the same plot

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');

Answer 1

[Star Strider]

Try adding this line to your code (after the loop):

tau_Vm_num = interp1(Vm_num, time, 0.632*max(Vm_num));

so the entire code is now:

% Model parameters
a = 0.001; % cell radius in cm
Rm = 5000; % specific membrane resistance in ohm*cm^2
Cm = 0.001; % specific membrane capacitance in mF/cm^2
I0 = 0.1; % current strength in mA
T = 50; % integration duration in ms
t = 0.02; % integration step in ms 
Tm = Rm * Cm; % time constant in ms 
num_iter = T / t; % number of interations of time loop
% Preallocate Vm_num (numerical Vm), Vm_ana (analytical Vm), and time 
% vectors for the for loop 
Vm_num = zeros(1,num_iter);
Vm_ana = zeros(1,num_iter);
time = zeros(1,num_iter);
% When calculating Rm * I0 in the Euler method, Rm is in terms of the area
% of the cell so need to multiply it by the surface area to get rid of the 
% cm^2 in the units. So Rm * I0 is in terms of mV instead of mV * cm^2
% Use surface area bc that is where the shock is being applied 
% Surface area of a circle is 4*pi*r^2
RM = Rm / (4*pi*a^2);
 
for j = 1:num_iter
    % initial condition: Vm(1) = 0 so start at 2nd index 
    Vm_num(j+1) = Vm_num(j) + (t/Tm)*((RM*I0)-Vm_num(j));
    Vm_ana(j+1) = (RM*I0)*(1 - exp(-time(j)/Tm));
    time(j+1) = time(j) + t;
end
tau_Vm_num = interp1(Vm_num, time, 0.632*max(Vm_num));
% Plot numerical and analytical Vm on the same plot 
plot(time,Vm_num);
hold on
plot(time,Vm_ana);
plot(tau_Vm_num, 0.632*max(Vm_num), 'sr')
title('Comparing Analytical and Numerical Membrane Potentials over Time');
xlabel('Time(ms)');
ylabel('Membrane Potential (mV)'); 
legend('Numerical Vm','Analytical Vm','\tau');

with ‘tau_Vm_num’ being the time constant, and the square marker in the plot denoting the time constant as you have defined it.

[Kate Heinzman]

Thank you so much!!!

[Star Strider]

As always, my pleasure!