dt = 0.01; % Time-step for plot
tmax = 10; % Maximum time for plot
t = 0:dt:tmax; % Vector of time-points
tau = 0.010; % Membrane time constant
E_L = -0.070; % Leak potential (steady state or resting potential)
E_K = -0.080; % Potassium reversal potential (for refractory conductance)
Vth = -0.050; % Threshold potential to produce spike
Vreset = -0.080; % Reset potential (post-spike)
Vmax = 0.050; % Clip spikes at this value of membrane potential
Cm = 100e-12; % Total membrane capacitance
G_L = Cm/tau; % Total membrane conductance (leak conductance)
Iapp = 180e-12; % Applied current
tauw = 200e-3; % Adaptation time constant (s)
a = 2e-9; % Adaptation recovery (S)
b = 0.02e-9; % Adaptation strength (A)
ton = 0.15; % Time to begin applied current (onset)
t0ff = 0.35; % Time to end applied current (offset)
non = round(ton/dt); % Time-point index of current onset
noff = round(toff/dt); % Time-point index of current offset
y = E_L + (starting_voltage_in_mV - E_L)*exp(-t/tau);
y_sim = zeros(size(t));
for i = 2:length(t)
y_sim(i) = y_sim(i-1) + dt; % Updating y value
end
The formula for the leaky integrate-and-fire model is as follows:
C_m*(dV_m/dt) = -G_L(V_m - E_L) + I_app
I = zeros(size(t)); % vector for current at each time-point
I(non:noff) = Iapp; % add the applied current
V = E_L*ones(size(t)); % initialize the membrane potential vector
spikes = zeros(size(t)); % initialize a vector to record spikes
for i = 2:length(t); % loop through all time points
% next line: Forward Euler method to update membrane potential
V(i) = V(i-1) + dt*(I(i) +G_L*(E_L-V(i-1)))/Cm;
if V(i) > Vth; % if potential is above threshold
spikes(i) = 1; % record the spike at that time-point
V(i) = Vreset; % reset the potential
end;
end; % end the loop & go to next time-point
The formula for the exponential leaky integrate-and-fire model is as follows:
C_m(dV_m/dt) = -G_L(V_m - E_L) + I_app + dthexp((V(i-1)-Vth)/dth))
I = zeros(size(t)); % vector for current at each time-point
I(non:noff) = Iapp; % add the applied current for the trial
V = E_L*ones(size(t)); % initialize the membrane potential vector
spikes = zeros(size(t)); % initialize a vector to record spikes
for i = 2:length(t); % loop through all time points
% next line: Forward Euler update of membrane potential
V(i) = V(i-1) + dt*(I(i) +G_L*(E_L - V(i-1) + Delta_th*exp((V(i-1)-Vth)/Delta_th)))/Cm;
if ( V(i) > Vmax ) % if potential is greater than Vmax
spikes(i) = 1; % record the spiketime
V(i) = Vreset; % Reset the membrane potential
V(i-1) = Vmax; % add a spike on prior time-point for visual purposes
end;
end;
The formula for the adaptive exponential leaky integrate-and-fire model is as follows:
C_m(dV_m/dt) = -G_L(V_m - E_L) + I_app + dTexp((v(j)-Vth/dT)) - w(j)
I = I0*ones(size(tvector)); % applied current initialized to baseline
I(non:noff) = Iapp; % add the step to the applied current vector
v = zeros(size(tvector)); % initialize membrane potential at all time-points
v(1) = E_L; % set initial value to be the leak potential
w = zeros(size(tvector)); % initialize adaptation variable at all time-points
spikes = zeros(size(tvector)); % intialize a vector to record spiketimes
for j = 1:length(tvector)-1 % simulation through all time-points
if ( v(j) > Vmax ) % if there is a spike
v(j) = V_Reset; % reset the voltage
w(j) = w(j) + b; % increase the adaptation variable by b
spikes(j) = 1; % record the spike
end
% next line integrates the membrane potential over time, using the
% Forward Euler method.
% first term within parentheses is like the LIF model
% second term is an exponential spiking term
% third term includes adaptation
v(j+1) = v(j) + dt*(G_L*(E_L-v(j) + deltaT*exp((v(j)-V_Thresh)/deltaT)) - w(j) + I(j))/C;
% next line decys the adaptation toward a steady state in between spikes
w(j+1) = w(j) + dt*(a*(v(j)-E_L) - w(j))/tauw;
end
set(0,'DefaultLineLineWidth',2,...
'DefaultLineMarkerSize',8, ...
'DefaultAxesLineWidth',2, ...
'DefaultAxesFontSize',14,...
'DefaultAxesFontWeight','Bold');
figure(1)
clf
% Creating Subplot
subplot(2,1,1)
plot(t,v,'k')
hold on
xlabel('x_label_name')
ylabel('y_label_name')
plot([0 1 2 3 4 5],[0 -10 -20 -30 -40], 'xk')
axis([0 5 -50 0])
% Creating Annotation for Subplot
annotation('textbox',[0.00 0.99 0.02 0.02],'LineStyle','none', ... 'FontSize',16,'FontWeight','Bold','String','subplot_name')