%%Written by Erin Boykin, 2011
%%
%%FUNCTION runMorrisLecar
%% simulates a Morris Lecar spiking neuron
%%
%%TO USE:
%% input: numNeurons – the number of neurons in the model
%% synapseType – type of synapse, either ‘linear’, ‘ftm’, or
%% ‘nonlinear’
%% totalSimTime – the total simulation time
%% stepSize – Euler integration step size
%% I – input currents to each neuron, must be at least numNeurons
%% in length
%% k – matrix of coupling coefficients between neurons. Coupling
%% goes from row to column in the matrix. Must be at least
%% numNeurons x numNeurons in size
%% D – noise strength to each neuron. must be at least numNeurons
%% in length
%%
%%
%% output: tout – time points of output
%% yout – neuron output at each time point. Note that odd indicies
%% are outputs of the neurons, while even indicies are internal
%% variables in the neuron model.
%%
%%EXAMPLE: [tout,yout] = runMorrisLecar(2,’linear’,2000,0.01,[100,80],[0,1;0,0],[1,1]);
%%
function [tout,yout] = runMorrisLecar(numNeurons,synapseType,totalSimTime,stepSize,I,k,D)
%%Do some intial parameter checking
if numNeurons <= 0 || totalSimTime 0 || length(D) ~= numNeurons
error(‘I,k, or D vector not of correct size’)
end
%%Set step size for Euler integration
h=stepSize;
t=0:h:totalSimTime;
%%Set up variables
clear y dy;
sizeY = 2*numNeurons;
y = zeros(length(t),sizeY);
y(1,:) = rand(1,sizeY);
dy = zeros(1,sizeY);
if strcmp(synapseType, ‘nonlinear’) || strcmp(synapseType, ‘Nonlinear’)
nonlinearY = zeros(length(t),numNeurons);
nonlinearY(1,:) = rand(1,numNeurons);
nonlinearDY = zeros(1,numNeurons);
end
%%Set standard parameters of the model
V1=-1.2;
V2=18;
V3=12;
V4=17.4;
Vl=-60;
Vk=-84;
Vca=120;
gl=2;
gk=8;
gca=4;
phi = (1/15);
C = 20;
Es = -84;
S0 = 1.012658;
Vth = 0;
Vslope = 1;
tau = 7.9;
for i = 1:length(t)-1,
for nn=1:numNeurons
V=y(i,nn*2-1);
W=y(i,nn*2);
minf=.5*(1+tanh((V-V1)/V2));
winf=.5*(1+tanh((V-V3)/V4));
tauw = (1/cosh((V-V3)/(2*V4)));
dy(nn*2-1) = (I(nn) – gca*minf*(V-Vca)-gk*W*(V-Vk)-gl*(V-Vl) + D(nn)*randn)/C;
dy(nn*2) = phi*(winf-W)/tauw;
%%Calculate coupling term
couplingSum = 0;
if strcmp(synapseType, ‘linear’) || strcmp(synapseType, ‘Linear’)
for j = 1:numNeurons
if nn~=j
couplingSum = couplingSum + k(j,nn)*flinear(y(i,nn*2-1),y(i,j*2-1));
end
end
elseif strcmp(synapseType, ‘ftm’) || strcmp(synapseType, ‘FTM’)
for j = 1:numNeurons
if nn~=j
couplingSum = couplingSum + k(j,nn)*fnonlinearInstant(y(i,j*2-1))*(y(i,nn*2-1)-Es);
end
end
elseif strcmp(synapseType, ‘nonlinear’) || strcmp(synapseType, ‘Nonlinear’)
for j = 1:numNeurons
if nn~=j
couplingSum = couplingSum + k(j,nn)*nonlinearY(i,j)*(y(i,nn*2-1)-Es);
end
end
if y(i,nn*2-1) > Vth
sinf = tanh((y(i,nn*2-1)-Vth)/Vslope);
else
sinf = 0;
end
nonlinearDY(nn)=(sinf-nonlinearY(i,nn))/(tau*(S0-sinf));
else
error(‘SynapseType is not valid. Must be either linear, ftm or nonlinear’);
end
couplingSum = couplingSum/C;
dy(nn*2-1) = dy(nn*2-1) + couplingSum;
end
y(i+1,:) = y(i,:) + h.*dy;
if exist(‘nonlinearY’,'var’)
nonlinearY(i+1,:) = nonlinearY(i,:) + h.*nonlinearDY;
end
end
%%This is optional – rescale to 1ms time intervals
tout = 0:1:totalSimTime;
yout = zeros(totalSimTime+1,4);
index=1;
for i = 1:length(t)
if mod(t(i),1) == 0
yout(index,:) = y(i,:);
index=index+1;
end
end
function out = flinear(x1,x2)
out = x2 – x1;
function out = fnonlinearInstant(x2)
out = 0.5*(1+tanh((x2-0)/1));
