function [t, vv] = LR1(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;
g_k = 0.282*(sqrt(K_o/5.4));
g_k1 = 0.6047*sqrt(K_o/5.4);
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
%% 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_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;
x = 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;
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_k ******************************/
alpha_x = 0.0005*exp(0.083*(v + 50.0))/(1.0 + exp(0.057*(v + 50.0)));
beta_x = 0.0013*exp(-0.06*(v + 20.0))/(1.0 + exp(-0.04*(v + 20.0)));
x_ss = alpha_x/(alpha_x+beta_x);
tau_x = 1.0/(alpha_x+beta_x);
if(v>-100.0)
x_i = 2.387*(exp(0.04*(v + 77.0)) - 1.0)/((v+77.0)*exp(0.04*(v + 35.0)));
else
x_i = 1.0;
end
%%*************** 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);
x = x_ss-(x_ss-x)*exp(-dt/tau_x);
%% 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));
%% 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);
i_si = g_ca*d*f*(v - E_si);
%%*************** i_k (delayed rectifier) ************************/
i_k = g_k*x*x_i*(v - E_k);
%%*************** 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_k + 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;
Modifications:
First line:
function [t, vv] = LR1(final_time)
was added an extra input parameter:
function [t, vv] = LR1(final_time, Q)
In Step One: Initailize,
g_k1 = 0.6047*sqrt(K_o/5.4);
the parameter was applied to the equation as a scaling factor:
g_k1 = Q*(0.6047*sqrt(K_o/5.4));
The S2 stimulation can be deleted or delayed.