% Simulation  of  Cortical  Spiking  Neurons% Source:% E M Izhikevich: Simple  Models  of  Spiking  Neurons  (2004)% Step 1: Clean  workspaceclose  all; % close  open  figuresclc; % clear  command  windowclear % clear  workspace% Step 2: Input  parametersflagS = 1; % select 1 of 20 spike  train  patternsN = 5000; % number  of grid  pointstMax = 100; % max  time [ms]% Step 3:  Define  model  parameters [a b c d I] 1x5  matrix  for  each  one of 20% patterns% a, b, c, and d simulate  the  neuron% I simulates  the  current  injection% Note:% current  injection = voltage  response% voltage  injection = current  responsecoeffM = zeros (20 ,5); % pre -allocate  20x5  matrix  of zeroscoeffM (1,:) = [0.02  0.4  -65 2 10]; % tonic  spikingcoeffM (2,:) = [0.02  0.25  -65 6 1]; % phasic  spikingcoeffM (3,:) = [0.02  0.20  -50 2 15]; % tonic  burstingcoeffM (4,:) = [0.02  0.25  -55 0.05  0.6]; % phasic  burstingcoeffM (5,:) = [0.02  0.20  -55 4 10]; % mixed  modecoeffM (6,:) = [0.01  0.20  -65 8 30]; % spike  frequency  adaptationcoeffM (7,:) = [0.02  -0.1  -55 6 0]; % Class IcoeffM (8,:) = [0.20  0.26  -65 0 0]; % Class IIcoeffM (9,:) = [0.02  0.20  -65 6 7]; % spike  latencycoeffM (10,:) = [0.05  0.26  -60 0 0]; % sub -threshold  oscillations4
coeffM (11,:) = [0.10  0.26  -60  -1 0]; % resonatorcoeffM (12,:) = [0.02  -0.1  -55 6 0]; % integratorcoeffM (13,:) = [0.03  0.25  -60 4 0]; % rebound  spikecoeffM (14,:) = [0.03  0.25  -52 0 0]; % rebound  burstcoeffM (15,:) = [0.03  0.25  -60 4 0]; % threshold  variabilitycoeffM (16,:) = [1 1.5  -60 0  -65]; % bistabilitycoeffM (17,:) = [1 0.2  -60  -21 0]; % DAPcoeffM (18,:) = [0.02 1  -55 4 0]; % accomodationcoeffM (19,:) = [-0.02  -1  -60 8 80]; % inhibition -induced  spikingcoeffM (20,:) = [ -0.026  -1  -45 0 80]; % inhibition -induced  bursting% Step 4: Use  input  parameters  to  simulate  one of  twenty  patternsa = coeffM(flagS ,1);b = coeffM(flagS ,2);c = coeffM(flagS ,3);d = coeffM(flagS ,4);I = coeffM(flagS ,5);% Step 5:  Define  all  plot  titles  (20)tm = ['1: Tonic  Spiking'; ...'2:  Phasic  Spiking'; ...'3: Tonic  Bursting'; ...'4:  Phasic  Bursting'; ...'5: Mixed  Mode'; ...'6: Spike  Frequency  Adaptation'; ...'7: Class 1'; ...'8: Class 2'; ...'9: Spike  Latency'; ...'10:  Subthreshold  Oscillations'; ...'11:  Resonator'; ...'12:  Integrator'; ...'13:  Rebound  Spike'; ...'14:  Rebound  Burst'; ...'15:  Threshold  Variability'; ...'16:  Bistability'; ...'17:  DAP'; ...'18:  Accomodation'; ...'19:  Inhibition -induced  Spiking'; ...'20:  Inhibition -induced  Bursting'; ];% Step 6: Solve  Equations% define c1, c2, c3:cc(1) = 0.04; % [mv]/[ms]cc(2) = 5; % [1]/[ ms]cc(3) = 140; % [mV]/[ms]% pre -allocationt = linspace(0,tMax , N); % time  steps [ms]v = zeros(N,1); vM = zeros(N,1); % membrane  potential [mV]u = zeros(N,1); uM = zeros(N,1); % recovery  variable [mV]Iext = zeros(N,1);Iext(N/10: end) = I; % define  size of I (initialization)dt = t(2)-t(1); % one  time  step% initialize  valuesv(1) =  -65;5
vM(1) =  -65;u(1) = b * v(1);uM(1) = u(1);% begin  loopfor n = 1:N-1v(n+1) = v(n) + dt*( cc(1)*v(n)^2 + cc(2)*v(n) + cc(3) -u(n) + Iext(n) );u(n+1) = u(n) + dt*a*( b*v(n) - u(n));if v(n+1) > 30vM(n+1) = 30;v(n+1) = c;u(n+1) = u(n+1) + d;elsevM(n+1) = v(n+1);uM(n+1) = u(n+1);endend% plot  figurefigure (2); % 2 plots  in one  figureset(gcf ,'units','normalized','Position', [0.3  0.32  0.32  0.38]);% gcf is  current  plot  open (grab  current  figure)% [x1 x2 y1 y2] provides  range of x- and y- values  to zoom in onset(gca ,'fontsize', 12);% gca is  current  plot  open (grab  current  axes)hSubplot = subplot (2,1,1);% subplot  is used to put  >1 plot in one  figure  that is 2x1% h = handle  in  hSubplot (important  in  graphical  user  interfaces (GUI) )set(hSubplot ,'position', [0.13  0.46  0.78  0.45]);xP = t; % x-axis is timeyP = vM; % y-axis is  voltageplot(xP, yP ,'b','Linewidth',2);ylabel('v_M [mV]');title(tm(flagS ,:)); % choose  title  from  matrix  of 20  entries% plot  voltage  change  based  on  current  injectiongrid onset(gca ,'fontsize',12);hSubplot = subplot (2,1,2);set(hSubplot ,'position', [0.13  0.15  0.78  0.20]);xP = t;yP = Iext; % external  injectionplot(xP,yP,'r','Linewidth',2);xlabel('time t [ms]');ylabel('c_5 I_{ext} [mV]');axis ([0  tMax 0 1.2* max(Iext)]);grid onset(gca ,'fontsize', 12);