Finding membrane time constant from a graph

8 views (last 30 days)
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');

Accepted Answer

Star Strider
Star Strider on 5 Apr 2020
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.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!