%% SFO Model
% Burst Firing:
% For values of gNa=150 & gK=100 (in mS/cm^2) we get burst firing,
% gKS ranges from 2-5mS/cm^2 depending on the burst regime (B1 vs B2).
% Tonic Firing:
% Transition burst firing regime to tonic firing by either:
% 1. Increasing gK from 100mS/cm^2 to 280mS/cm^2.
% 2. Increasing gNa from from 150mS/cm^2 to 170mS/cm^2.
clear
%% Time:
dt = 0.01; % Step Size for Eulers Method
timesec = 30; % End Time (s)
tfinal = timesec*1000; % End Time (ms)
Tstep =(0:dt:tfinal); % Time Array (ms)
%% Initial Values:
V=-60; % Membrane Potential (mV)
nA=0; % IK Activation
mA=0.35; % IKS Activation
hA=0; % IKS Inactivation
mKA=0; % IA Activation
hKA=0; % IA Inactivation
sA=0; % INa Activation
kA=0; % INa Inactivation
mNaP=0; % INaP Activation
hNaP=0; % INaP Inactivation
mNCa=0; % ICa Activation
%% Cell Properties
% Cell Morphology:
CellSize = 10; % Model Cell Diameter in microns
CellArea = 4*pi*((CellSize/2).^2); % Model Cell Area for Spherical Shape
% Cell Capacitance:
Cr = 5; % Real Cell Capacitance in pF
C = (Cr/CellArea)*100; % Model Capacitance in uF/cm^2
% Cell Resistance:
Rr = 1; % Real Cell Resistance in gigaohms
R = (Rr * CellArea)*10; % Model Cell Resistance in ohms*cm^2
% Reversal Potential (in mV):
Er = -65; % Reversal Potential for leak
ENa = 107; % Reversal Potential for Na+
EK = -88; % Reversal Potential for K+
ENSCC = -35; % Reversal Potential for NSCC
ECa = 120; % Reversal Potential for Ca2+
% Conductance (in mS/cm^2):
gLeak = (1/R)*1000; % Conductance of IL
gK = 100; % Conductance of IK
gKS = 3; % Conductance of IKS
gA = 3; % Conductance of IA
gNa = 150; % Conductance of INa
gNaP = 0.13; % Conductance of INaP
gNSCC = 0.2; % Conductance of INSCC
gCa = 0.3; % Conductance of ICa
%% Zeros Vectors
Vm=zeros(1,length(Tstep));
Iz=zeros(1,length(Tstep));
IN=zeros(1,length(Tstep));
%% Looping Code
% Solving differential equations with Eulers method:
tidx=0;
for z=Tstep
tidx=tidx+1;
%% Injected Current
% Where 1 pA = 0.3183 uA/cm^2 & 10 pA = 3.18 uA/cm^2
I=0;
%% Time Constants (in ms):
tau_nA = 7.2-(6.4/(1+exp((V+28.3)/-19.2))); % Time Constant for IK activation
tau_mA = 3000; % Time Constant for IKS activation
tau_hA = 10; % Time Constant for IKS inactivation
tau_sA = 0.1; % Time Constant for INa activation
tau_kA = 1; % Time Constant for INa inactivation
tau_mNaP = 5; % Time Constant for INaP activation
tau_hNaP = 50; % Time Constant for INaP inactivation
tau_mNCa = 5; % Time Constant for ICa activation
tau_mKA = 5; % Time Constant for IA activation
tau_hKA = 30; % Time Constant for IA inactivation
%% Current Equations:
% Potassium Currents:
% Delayed-rectifier K+ Current (IK):
nA0=1/(1+exp((V-2)/-8));
nA = nA + dt*((nA0-nA)/tau_nA);
IK = gK*nA^4*(V-EK);
% Slow-activating K+ Current (IKS):
mA0 = 1/(1+exp((V+44)/-18));
mA = mA + dt*((mA0-mA)/tau_mA);
hA0 = 1/(1+exp((V+61)/8));
hA = hA + dt*((hA0-hA)/tau_hA);
IKS = gKS*mA^3*hA*(V-EK);
% Transient K+ Current (IA):
mKA0 = 1/(1+exp((V+44)/-18));
mKA = mKA + dt*((mKA0-mKA)/tau_mKA);
hKA0 = 1/(1+exp((V+61)/8));
hKA = hKA + dt*((hKA0-hKA)/tau_hKA);
IA = gA*mKA^3*hKA*(V-EK);
% Sodium Currents:
% Transient Na+ Current (INa):
sA0 = 1/(1+exp((V+31)/-6.11));
sA = sA + dt*((sA0-sA)/tau_sA);
kA0 = 1/(1+exp((V+62)/6.15));
kA = kA + dt*((kA0-kA)/tau_kA);
INa = gNa*sA^3*kA*(V-ENa);
% Persisent Na+ Current (INaP):
mNaP0 = 1/(1+exp((V+55)/-4));
mNaP = mNaP + dt*((mNaP0-mNaP)/tau_mNaP);
hNaP0 =1/(1+exp((V+45)/6.1));
hNaP = hNaP + dt*((hNaP0-hNaP)/tau_hNaP);
INaP = gNaP*mNaP^3*hNaP*(V-ENa);
% Calcium Currents:
% N-Type Ca2+ Current (ICa):
mNCa0 = 1/(1+exp((V+14)/-5.8));
mNCa = mNCa + dt*((mNCa0-mNCa)/tau_mNCa);
ICa = gCa*mNCa^2*(V-ECa);
% Leak Current (IL):
IL = gLeak*(V-Er);
% Non-selective Cation Current (INSCC):
INSCC = gNSCC*(V-ENSCC);
% Noise (IN):
stdnoise = 9; % Standard Deviation for Noise
IN=stdnoise*randn; % Noise Current (Gaussian Distribution)
%% Voltage Calculation:
V = V + ((dt/C)*(I + IN - IL - INa - IK - INSCC - INaP - IA - ICa - IKS));
%% Voltage Array:
Vm(tidx)= V;
end
%% Membrane Potential vs Time Plot:
figure('Renderer', 'painters', 'Position', [100 100 1200 500])
sh(1)=subplot(1,1,1);
plot(Tstep,Vm,'k-')
% Figure Settings:
box off
ax = gca;
ax.FontSize = 14;
xlabel('Time (s)','FontSize',14)
ylabel('Membrane Potential (mV)','FontSize',14)
xt=arrayfun(@num2str,get(gca,'xtick')/1000,'un',0);
set(gca,'xticklabel',xt)
| 
