(Click to enlarge.)
Read More
Superimposed Action Potential Graph
Sensitivity Analysis
Resting Membrane Potential Shift
Archive for May 2011
Codes for GUI:
function varargout = BME301GUI(varargin)
% BME301GUI MATLAB code for BME301GUI.mov
% BME301GUI, by itself, creates a new BME301GUI or raises the existing
% singleton*.
%
% H = BME301GUI returns the handle to a new BME301GUI or the handle to
% the existing singleton*.
%
% BME301GUI('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in BME301GUI.M with the given input arguments.
%
% BME301GUI('Property','Value',...) creates a new BME301GUI or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before BME301GUI_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to BME301GUI_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help BME301GUI
% Last Modified by GUIDE v2.5 11-May-2011 11:47:52
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @BME301GUI_OpeningFcn, ...
'gui_OutputFcn', @BME301GUI_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before BME301GUI is made visible.
function BME301GUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to BME301GUI (see VARARGIN)
% Choose default command line output for BME301GUI
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes BME301GUI wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = BME301GUI_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = handles.output;
% --- Executes on button press in Play.
function Play_Callback(hObject, eventdata, handles)
% hObject handle to Play (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
Pc = handles.InputB;
g = [1:-0.05:Pc];
z = zeros(100000,length(g));
i = 1;
while i <= length(g)
[ti, vt] = LR1M(1000, g(i));
z(:, i) = vt;
i = i + 1;
end
for n = 1:(length(g))
plot(ti, z(:,n)); axis([0 1000 -100 50]);
xlabel('Time (ms)'); ylabel('Transmembrane Voltage (mV)');
M(:,n) = getframe;
end
moviereps = -5; % number of repetitions of the movie
speedratio = 10; % ratio of movie speed to actual speed
handles.Fig = M;
movie(handles.Fig,moviereps,speedratio)
guidata(hObject, handles);
function InputB_Callback(hObject, eventdata, handles)
% hObject handle to InputB (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of InputB as text
% str2double(get(hObject,'String')) returns contents of InputB as a double
handles.InputB = str2num(get(hObject,'String'))*0.01;
if(isempty(handles.InputB))
set(hObject,'String','0')
end
guidata(hObject,handles);
% --- Executes during object creation, after setting all properties.
function InputB_CreateFcn(hObject, eventdata, handles)
% hObject handle to InputB (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
% --- Executes on button press in plot.
function plot_Callback(hObject, eventdata, handles)
% hObject handle to plot (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
[time, voltage] = LR1M(1000, handles.InputB);
handles.Fig = plot(time, voltage);
axis([0 1000 -100 50]);
xlabel('Time (ms)')
ylabel('Transmembrane Voltage (mV)')
guidata(hObject,handles);
% --- Executes on button press in Save.
function Save_Callback(hObject, eventdata, handles)
% hObject handle to Save (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
[file_name, path_name] = uiputfile({ '*.png','PNG File (*.png)';...
'*.bmp','Bitmap (*.bmp)'; '*.fig','Figure (*.fig)'}, ...
'Save picture as','default');
if isequal(file_name,0) || isequal(path_name,0)
return
end
% --- Executes on button press in Exit.
function Exit_Callback(hObject, eventdata, handles)
% hObject handle to Exit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
close
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.