Effect of Conductance Percentage On Voltage Over Time

Effect of Voltage Shift On Voltage Over Time

MATLAB Code

function [t, vv] = LR3(q,shift,final_time)
%% Emilia Entcheva
%% Matlab template for class project in BME301
%% This code describes a model of an excitable cell according to Luo-Rudy
%% 1991
%%
%% Choose your own input parameters, currently only one is specified.
%% ‘final_time’ is the total simulation time in [ms], for
%% example 1500ms (1sec)

%%**************** model constants ********************/
R = 8.314;
F = 96490;
T = (273+37);
nernst = 26.71;
%% ***********Step One: Initialize!!! ***********************

Vm_rest = -84.54392;
Ca_rest = 1.783825e-4;
c_m = 1.0;
g_na = 23.0;
g_ca = 0.09;
g_kp = 0.0183;
g_b = 0.03921;
PR_nak = 0.01833;
K_o = 5.4;
K_i = 145.0;
Na_o = 140.0;
Na_i = 18.0;
Ca_o = 1.8;
Ca_i=0;
g_k = 0.282*(sqrt(K_o/5.4));
g_kr=.02614*sqrt(K_o/5.4);
g_k1 = 0.6047*sqrt(K_o/5.4);
g_cat=0.05;

init_time_step = 0.01;

%% STIMULATION (you may change as desired)
s1_current = 100.0; %% [uA/cm2]
s1_ini = 100.0; %% [msec], offset from time 0 (time of delicvery of first pulse)
s1_dur = 0.5; %% [msec], duration of stimulus

s2_current = 100.0; %% [uA/cm2]
s2_ini = 600.0; %% [msec], offset from time 0 (time of delivery of second pulse)
s2_dur = 0.5; %% [msec], duration of stimulus

s3_current = 100.0; %% [uA/cm2]
s3_ini = 1200.0; %% [msec], offset from time 0 (time of delivery of second pulse)
s3_dur = 0.5; %% [msec], duration of stimulus
%% Assign variables initial values
v = Vm_rest;
ct = 0.0;
dt = init_time_step;

%% Calculate Nernst potentials based on ion concentrations
E_na = nernst*log(Na_o/Na_i);
E_k = nernst*log((K_o + PR_nak*Na_o)/(K_i + PR_nak*Na_i));
E_kr=nernst*log(K_o/K_i);
E_ks=nernst*log((K_o+PR_nak*Na_o)/(K_i+PR_nak*Na_i));
E_k1 = nernst*log(K_o/K_i);
E_kp = E_k1;

t(final_time/dt) = 0;
vv(final_time/dt) = 0;

m = 0;
h = 0;
jj = 0;
d = 0;
f = 0;
xr = 0;
xs = 0;
xss=0;
Ca_i = 0;

i_si = 0;

%% ***********Step Two: Iterate in time in a loop ***********************
for j=1:(final_time/dt);
%% Stimulate as prescribed above

if((ct>=s1_ini) && (ct<=(s1_ini+s1_dur)))
istim = s1_current;
elseif((ct>=s2_ini) && (ct<=(s2_ini+s2_dur)))
istim = s2_current;
elseif((ct>=s3_ini) && (ct<=(s3_ini+s3_dur)))
istim = s3_current;
else
istim = 0.0;
end

%% Gate variables get calculated for the current voltage values

%%*************** i_na ******************************/
if (v>=-40.0)
alpha_h = 0.0;
alpha_j = 0.0;
beta_h = 1/(0.13*(1.0 + exp((v + 10.66)/(-11.1))));
beta_j = 0.3*exp(-2.535e-7*v)/(1.0 + exp(-0.1*(v + 32)));

else
alpha_h = 0.135*exp((80.0 + v)/(-6.8));
beta_h = (3.56*exp(0.079*v)) + (3.1e5*exp(0.35*v));
alpha_j = ((-1.2714e5*exp(0.2444*v))-(3.474e-5*exp(-0.04391*v)))*((v + 37.78)/(1.0 + exp(0.311*(v+79.23))));
beta_j = 0.1212*exp(-0.01052*v)/(1.0 + exp(-0.1378*(v + 40.14)));

end

if(v==-47.13)
alpha_m = 3.2;
else
alpha_m = 0.32*(v + 47.13)/(1.0 – exp(-0.1*(v + 47.13)));
end
beta_m = 0.08*exp(-1.0*v/11.0);

m_ss = alpha_m/(alpha_m+beta_m);
h_ss = alpha_h/(alpha_h+beta_h);
j_ss = alpha_j/(alpha_j+beta_j);

tau_m = 1.0/(alpha_m+beta_m);
tau_h = 1.0/(alpha_h+beta_h);
tau_j = 1.0/(alpha_j+beta_j);

%%*************** calcium ******************************/
alpha_d = 0.095*exp(-0.01*(v – 5.0))/(1.0 + exp(-0.072*(v – 5.0)));
beta_d = 0.07*exp(-0.017*(v + 44.0))/(1.0 + exp(0.05*(v + 44.0)));
alpha_f = 0.012*exp(-0.008*(v + 28.0))/(1.0 + exp(0.15*(v + 28.0)));
beta_f = 0.0065*exp(-0.02*(v + 30.0))/(1.0 + exp(-0.2*(v + 30.0)));

d_ss = alpha_d/(alpha_d+beta_d);
f_ss = alpha_f/(alpha_f+beta_f);

tau_d = 1.0/(alpha_d+beta_d);
tau_f = 1.0/(alpha_f+beta_f);
%%*************** i_kr ******************************/
alpha_xr = 0.00138*(v+14.2)/(1-exp(-0.123*(v+14.2)));
beta_xr = 0.00061*(v+38.9)/(exp(0.145*(v+38.9))-1);

xr_ss = alpha_xr/(alpha_xr+beta_xr);
tau_xr = 1.0/(alpha_xr+beta_xr);

%%*************** i_ks ******************************/
alpha_xs = (7.19*(10^-5)*(v+shift+30))/(1-exp(-0.148*(v+shift+30)));
beta_xs = (1.31*(10^-4)*(v+shift+30))/(exp(0.0687*(v+shift+30))-1);

xs_ss = alpha_xs/(alpha_xs+beta_xs);
xss_ss=xs_ss;
tau_xs = 1.0/(alpha_xs+beta_xs);
tau_xss = 4*tau_xs;

%%*************** i_k1 ******************************/
alpha_k1 = 1.02/(1.0 + exp(0.2385*(v – E_k1 – 59.215)));
temp1 = 0.49124*exp(0.08032*(v – E_k1 + 5.476));
temp2 = exp(0.06175*(v – E_k1 – 594.31));
temp3 = (1.0 + exp(-0.5143*(v – E_k1 + 4.753)));
beta_k1 = (temp1 + temp2)/temp3;

k1_ss = alpha_k1/(alpha_k1 + beta_k1);

%%*************** i_kp ******************************/
kp = 1.0/(1.0 + exp((7.488 – v)/5.98));

%% Integrate all gates
m = m_ss-(m_ss-m)*exp(-dt/tau_m);
h = h_ss-(h_ss-h)*exp(-dt/tau_h);
jj = j_ss-(j_ss-jj)*exp(-dt/tau_j);
d = d_ss-(d_ss-d)*exp(-dt/tau_d);
f = f_ss-(f_ss-f)*exp(-dt/tau_f);
xr = xr_ss-(xr_ss-xr)*exp(-dt/tau_xr);
xs = xs_ss-(xs_ss-xs)*exp(-dt/tau_xs);
xss = xss_ss-(xss_ss-xss)*exp(-dt/tau_xss);
r = 1/(1+exp((v+9)/22.4));

%% Update intracellular calcium concentration (it’s changing in time)

Ca_i = Ca_i + dt*(-1.0e-4*i_si + 0.07*(1.0e-4 – Ca_i));
pca=3-log10(Ca_i);
g_ks=(q/100)*(.057+.19/(1+exp((-7.2+pca)/.6)));

%% Calculate all ion currents

%%*************** i_na (sodium) *********************/
i_na = g_na*m*m*m*h*jj*(v – E_na);

%%*************** i_si (slow inward calcium ) ************************/
E_si = 7.7 – 13.0287*log(Ca_i);
fff=1/(1+Ca_i/.6);
i_si = g_ca*d*f*fff*(v – E_si);

%Calcium T channel
E_ca = nernst*log(Ca_o/Ca_i);
b=1/(1+exp(-(v+14)/10.8));
g=1/(1+exp((v+60)/5.6));
i_cat=g_cat*b^2*g*(v-E_ca);

%%*************** i_k (delayed rectifier) ************************/
i_kr = g_kr*xr*r*(v – E_kr);
i_ks = g_ks*xs*xss*(v +shift – E_ks);

%%*************** i_k1 (time-independent potassium (inward rectifier) ************/
i_k1 = g_k1*k1_ss*(v – E_k1);

%%*************** i_kp (plateau potassium) ********************/
i_kp = g_kp*kp*(v – E_kp);

%%*************** i_b (background) ************************/
i_b = g_b*(v + 59.87);

%%*************** i (total) ******************************/
itot = i_na + i_si + i_cat + i_kr + i_ks + i_k1 + i_kp + i_b;

%% Update time and voltage for the next time step
ct = ct + dt;
v = v + dt*(1.0/c_m)*(istim-itot);
t(j) = ct;
vv(j) = v;
end

%% ***********Step Three: Plot results ***********************
%% plot voltage vs. time
plot(t, vv);
grid on;