Computational Neuroscience Matlab Plotting Basics Notes

Basic Settings

dt = 0.01; 			% Time-step for plot
tmax = 10;			% Maximum time for plot
t = 0:dt:tmax;		% Vector of time-points

Parameters

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

Membrane potential exponential decay

y = E_L + (starting_voltage_in_mV - E_L)*exp(-t/tau);

Simple Simulation

y_sim = zeros(size(t));
for i = 2:length(t)
	y_sim(i) = y_sim(i-1) + dt;			% Updating y value
end

Simulating Leaky Integrate-and-Fire Model (LIF Model)

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

Simulating Exponential Leaky Integrate-and-Fire Model (ELIF Model)

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;

Simulating Adaptive Exponential Leaky Integrate-and-Fire Model (AELIF Model)

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

Setting the Plotting Parameters

set(0,'DefaultLineLineWidth',2,...
    'DefaultLineMarkerSize',8, ...
    'DefaultAxesLineWidth',2, ...
    'DefaultAxesFontSize',14,...
    'DefaultAxesFontWeight','Bold');

Basic Plotting

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')