% Acker et al 2003 model variant neuron network with random connectivity
% with time varying input with LIF-measured stability
% eric zilli - oct ??, 2009
%
% release version 1.0. check modeldb for updates.
%
% this source code is released into the public domain
%
% The Acker et al model is a biophysical model of an entorhinal cortex
% layer II stellate. The model uses voltage-shifted Hodghkin-Huxley loligo
% Na and K channels for spiking, plus a non-inactivating persistent-sodium
% current and a two-component (fast and slow) H-current which are modified
% by Acker et al. from original current models by Erik Fransen. We modify the model by
% scaling down all the conductance densities by increasing the membrane
% capacitance. This slows the minimum firing frequency of the model down so
% that the noise level can be scaled to match our biological data (ISI mean
% 0.2 s, ISI std. dev. 0.036 s).
%
% Interesting questions to continue this research: Setting slowslow=1
% decreases the adaptation's initial magnitude but the effect
% has an effect lasting many ISIs long. Is
% the corresponding time constant model realistic? Experiments looking at H
% current time constants at spiking voltages do not appear to exist, but
% gating variable "friction" in the extended HH model accounts for a
% possible non-zero base for time constants which might play a major role
% here. In-vitro experiments are needed.
% Under what conditions does the ISI period or frequency overshoot and when does it undershoot?
% What is the best way to quantify the history dependence on firing rate
% effect?
% Is it possible for the intrinsic properties of one neuron A to be exactly
% the inverse of B in the sense of this adaptation such that an
% input signal I to A would be distorted by the adaptation into an
% output I' which would then be un-distorted by B to recover a constant function
% of I? (E.g. I goes from driving at 7 Hz to 11 Hz, A's firing rate adapts from
% 11.5 Hz down to 11 Hz. Can there be a B receiving A's spikes that fires at 11 Hz
% given input at 11.5 Hz and which adapts at just the right rate that B continues to
% fire at 11 Hz as the input from A slows down slightly?)
%
% This script starts a network with a specified current level (or
% frequency) and after a specified number of spikes allows the current
% level to be changed. This can be done with a single starting and ending
% current level (Figure S2A-D), a single starting and multiple ending levels
% (Figure S2E-F), multiple starting and single ending level, or a grid of injected
% current levels (Figure S3). These correspond to three types of figures that
% might be of interest.
%
% If the FI curve of the cell or network is known, it can be used to find
% the appropriate current levels. Otherwise they must be set by hand. The
% FI curve needn't be terribly high in resolution, so generating them as
% needed should be quick and easy.
%
% Note: this is a cleaned up version of the script used to generate the
% paper figures--it is possible errors were introduced during the cleaning!
% Feel free to contact me if there is any difficulty reproducing any
% results from the manuscript.
%
clear all;
% type=1 - (Figure S2B,D) class 2 excitable, n=1, noise-free, 1 start, 1 end
% type=2 - (Figure S2F) class 2 excitable, n=1, noise-free, 1 start, 5 end, low-to-high
% type=3 - class 2 excitable, n=1, noise-free, 7:0.5:11 Hz start, 7:0.5:11 Hz end
for type=[3]
simdur = 10000; % (ms)
dt = .01; % ms
nruns =1;
saveFigures = 0;
useNoise = 0;
spikeThreshold = 20; % mV
slowslow = 0;
numDropISIs = 10;
numPreISIs = 10;
numPostISIs = 10;
if type==1
excitationClass = 2;
useGapJunctions = 1;
g = 0;
uniqueNoiseSTD = 0*useNoise;
ncells = 1;
pcon = 0;
load Acker_FI_n1.mat;
inputMags = interp1(freqs,currents,7,'cubic');
inputMags2 = interp1(freqs,currents,11,'cubic');
elseif type==2
excitationClass = 2;
useGapJunctions = 1;
g = 0;
uniqueNoiseSTD = 0*useNoise;
ncells = 1;
pcon = 0;
load Acker_FI_n1.mat;
desiredFreqs = 7:11;
inputMags = interp1(freqs,currents,7,'cubic');
inputMags2 = interp1(freqs,currents,desiredFreqs,'cubic');
elseif type==3
excitationClass = 2;
useGapJunctions = 1;
g = 0;
uniqueNoiseSTD = 0*useNoise;
ncells = 1;
pcon = 0;
load Acker_FI_n1.mat;
desiredFreqs = [7 11];
inputMags = interp1(freqs,currents,desiredFreqs,'cubic');
inputMags2 = interp1(freqs,currents,desiredFreqs,'cubic');
end
fileOut = 0;
typePrefix = [sprintf('%d',excitationClass)];
if ncells>1
if useGapJunctions
typePrefix = [typePrefix 'g'];
else
typePrefix = [typePrefix 's'];
end
end
if useNoise
typePrefix = [typePrefix 'n'];
end
outfilename = sprintf('acker_model_%s_history_dependence_%dn_%gnoise_%s.txt',typePrefix,ncells,useNoise,datestr(now,'mmmmdd'));
% file to save statistics for every individual cell in network, don't
% save if indivname = ''
indivname = '';
% indivname = [strtok(outfilename,'.') '_individualcells.txt'];
Cm = 5; 1.5; % uF/cm^2
gNa = 52; % mS/cm^2
gNap = 0.5; % mS/cm^2
gH = 1.5; % mS/cm^2
gK = 11; % mS/cm^2
gLeak = 0.5; % mS/cm^2
EH = -20; % (mv)
ELeak = -65; % (mv)
EK = -90; % (mv)
ENa = 55; % (mv)
slopes = zeros(length(inputMags));
if length(inputMags)==1 && length(inputMags2)>1
figISIs = zeros(length(inputMags2),numPreISIs+numPreISIs);
elseif length(inputMags)>1 && length(inputMags2)==1
figISIs = zeros(length(inputMags),numPreISIs+numPreISIs);
elseif length(inputMags)>1 && length(inputMags2)>1
figISIs = zeros(length(inputMags),length(inputMags2));
end
% connectivity matrix
if pcon==1
C = 1;
else
C = (rand(ncells)<pcon)>0;
C = sparse(C.*(1-eye(ncells))); % remove autoconnections
end
for inputi=1:length(inputMags)
for inputi2=1:length(inputMags2)
stabilities = zeros(1,nruns);
stabilities2 = zeros(1,nruns);
for run=1:nruns
% activity of postsynaptic I&F
npost = 1;
post = zeros(npost, round(simdur/dt)+1);
% spike times of postsynaptic I&F
postspikes = 0;
% I&F params
tau = 2.5; % (sec)
membraneDecay = exp(-dt/tau);
postWeight=4/ncells;
% magnitude of noise added to all cells on each step
commonNoiseSTD = 0;
numISIs = -1;
t = 0; % current time in simulation
tind = 1;
% initial values
v = -55*ones(1,ncells);
mNa = zeros(1,ncells);
hNa = zeros(1,ncells);
mNap = zeros(1,ncells);
nK = zeros(1,ncells);
mHfast = zeros(1,ncells);
mHslow = zeros(1,ncells);
% save state of one cell over simulation:
state = zeros(2,simdur/dt);
% binary array indicating which cells fire on each time steps
spikes = spalloc(ncells, simdur/dt,ncells*simdur); % sparse to save memory in big simulations
I = inputMags(inputi);
Is = zeros(1,round(simdur/dt));
tic
while t<simdur
t = t+dt; % advance to next time step
tind = 1+tind; %round(t/dt);
%% done below now:
oldv = v;
alphamNa = -0.1*(oldv+23)./(exp(-0.1*(oldv+23))-1);
betamNa = 4*exp(-(oldv+48)/18);
mNa = mNa + dt*(alphamNa.*(1-mNa) - betamNa.*mNa);
alphahNa = 0.07*exp(-(oldv+37)/20);
betahNa = 1./(exp(-0.1*(oldv+7))+1);
hNa = hNa + dt*(alphahNa.*(1-hNa) - betahNa.*hNa);
alphanK = -0.01*(oldv+27)./(exp(-0.1*(oldv+27))-1);
betanK = 0.125*exp(-(oldv+37)/80);
nK = nK + dt*(alphanK.*(1-nK) - betanK.*nK);
alphamNap = 1./(0.15*(1+exp(-(oldv+38)/6.5)));
betamNap = exp(-(oldv+38)/6.5)./(0.15*(1+exp(-(oldv+38)/6.5)));
mNap = mNap + dt*(alphamNap.*(1-mNap) - betamNap.*mNap);
minfHfast = 1./(1+exp((oldv+79.2)/9.78));
mtauHfast = 0.51./(exp((oldv-1.7)/10) + exp(-(oldv+340)/52)) + 1;
mHfast = mHfast + dt*((minfHfast-mHfast)./mtauHfast);
minfHslow = 1./(1+exp((oldv+71.3)/7.9));
if ~slowslow
mtauHslow = 5.6./(exp((oldv-1.7)/14) + exp(-(oldv+260)/43)) + 1;
else
mtauHslow = 65.5813 + 248.0469*exp(-(-79.2190-oldv).^2/33.5178^2);
end
mHslow = mHslow + dt*((minfHslow-mHslow)./mtauHslow);
vnoise = uniqueNoiseSTD*randn(1,ncells);
if useGapJunctions
if pcon==0
v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I)/Cm;
elseif pcon==1
v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*(sum(v)-ncells*v))/Cm;
else
vdiffs = repmat(v',ncells,1) - repmat(v,1,ncells);
v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*sum(C*vdiffs,2))/Cm;
end
else
if pcon==0
v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I)/Cm;
elseif pcon==1
v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*(sum(spikes(:,tind-1),1)-spikes(:,tind-1)'))/Cm;
else
v = v + dt*(vnoise - gH*(0.65*mHfast + 0.35*mHslow).*(v-EH) - (gNap*mNap+gNa*mNa.^3.*hNa).*(v-ENa) - gK*nK.^4.*(v-EK) - gLeak*(v-ELeak) + I + g*(C*spikes(:,tind-1))')/Cm;
end
end
state(:,tind) = [v(1) mHslow];
% spikes if upward crossing spikeThreshold mV
spikes(:,tind) = (oldv<spikeThreshold).*(v>=spikeThreshold);
%% do numPreISIs then change I and do numPostISIs
if any(spikes(:,tind))
numISIs = numISIs + 1; % starts equal to -1 so second spike is first ISI
end
if numISIs<(numPreISIs+numDropISIs)
I = inputMags(inputi);
else
I = inputMags2(inputi2);
end
if numISIs==(numPreISIs+numPostISIs+numDropISIs)
break
end
post(:,tind) = post(:,tind-1)*membraneDecay + tau*postWeight*sum(spikes(:,tind));
% faster(?) if we don't need postspikes:
post(post(:,tind)>1,tind) = -1e-10;
% slower way if postspikes is needed:
% postspikes(:,tind) = post(:,tind)>1;
% post(postspikes(:,tind)>1,tind) = 0; % reset to 0
end
toc
% calculate mean and variance of ISIs for each postsynaptic cell
npost = 1; %size(post,1);
ISImeans = zeros(1,npost);
ISIstds = zeros(1,npost);
ISIstabilities = zeros(1,nruns);
%% calculate ISI statistics for each individual unit in the
%% network:
indivmeans = zeros(1,ncells);
indivstds = zeros(1,ncells);
indivstabilities = zeros(1,ncells);
ISIs = dt*diff(find(spikes(1,:)>0))/1000;
if length(inputMags)==1 && length(inputMags2)>1 % add so we can take the average over runs
figISIs(inputi2,:) = figISIs(inputi2,:) + ISIs((numDropISIs+1):end);
elseif length(inputMags)>1 && length(inputMags2)==1
figISIs(inputi,:) = figISIs(inputi,:) + ISIs((numDropISIs+1):end);
elseif length(inputMags)>1 && length(inputMags2)>1
% using a complex number to hold two values in one matrix:
figISIs(inputi2,inputi) = figISIs(inputi2,inputi) + ISIs(numDropISIs+numPreISIs+1) + ISIs(end-1)*i;
end
% save to file:
if fileOut
fid = fopen(outfilename,'a');
fprintf(fid,'%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%d\n', ...
ncells,pcon,inputMags(inputi),inputMags2(inputi2),g,simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
mean(ISIs(2:end-numPostISIs-2)), ISIs(end-numPostISIs-2), ISIs(end-numPostISIs+1), ISIs(end-numPostISIs+2),ISIs(end),length(ISIs),...
(inputMags2(inputi2)-inputMags(inputi)), ISIs(end-numPostISIs+1)-ISIs(end-numPostISIs+2),ISIs(end-numPostISIs+1)-ISIs(end));
fclose(fid);
end
fprintf('%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%d\n', ...
ncells,pcon,inputMags(inputi),inputMags2(inputi2),g,simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
mean(ISIs(2:end-numPostISIs-2)), ISIs(end-numPostISIs-2), ISIs(end-numPostISIs+1), ISIs(end-numPostISIs+2),ISIs(end),length(ISIs),...
(inputMags2(inputi2)-inputMags(inputi)), ISIs(end-numPostISIs+1)-ISIs(end-numPostISIs+2),ISIs(end-numPostISIs+1)-ISIs(end));
slopes(inputi,inputi2) = (ISIs(end-numPostISIs-1)-ISIs(end-numPostISIs))/(inputMags2(inputi2)-inputMags(inputi));
end
end
% average ISIs over multiple runs if noisy network is in use:
if length(inputMags)>1 || length(inputMags2)>1
figISIs = figISIs./nruns;
end
end
end
state = state(:,1:tind);
if length(inputMags)==1 && length(inputMags2)==1
start=260000;stop=360000;
figure;
plot((start:stop)*dt,state(1,start:stop)');
set(gca,'xlim',[start*dt stop*dt]);
set(gcf,'position',[82 404 257 185]);
xlabel('Time (ms)','fontsize',12); ylabel('Voltage (mV)','fontsize',12);
set(gca,'fontsize',12,'box','off')
if slowslow
if saveFigures && exist(['fig_Ackerslow_history_dependence_vtrace_0a.eps'])~=2
print_eps(['fig_Ackerslow_history_dependence_vtrace_0a.eps'])
saveas(gcf,['fig_Ackerslow_history_dependence_vtrace_0a.fig'])
end
else
if saveFigures && exist(['fig_Acker_history_dependence_vtrace_0a.eps'])~=2
print_eps(['fig_Acker_history_dependence_vtrace_0a.eps'])
saveas(gcf,['fig_Acker_history_dependence_vtrace_0a.fig'])
end
end
figure;
plot((start:stop)*dt,state(2,start:stop)');
set(gcf,'position',[82 404 257 227]);
set(gca,'ygrid','on','xlim',[start*dt stop*dt]);
xlabel('Time (ms)','fontsize',12); ylabel('Slow H','fontsize',12);
set(gca,'fontsize',12,'box','off')
if slowslow
if saveFigures && exist(['fig_Ackerslow_history_dependence_htrace_0a.eps'])~=2
print_eps(['fig_Ackerslow_history_dependence_htrace_0a.eps'])
saveas(gcf,['fig_Ackerslow_history_dependence_htrace_0a.fig'])
end
else
if saveFigures && exist(['fig_Acker_history_dependence_htrace_0a.eps'])~=2
print_eps(['fig_Acker_history_dependence_htrace_0a.eps'])
saveas(gcf,['fig_Acker_history_dependence_htrace_0a.fig'])
end
end
elseif (length(inputMags)==1 && length(inputMags2)>1) || (length(inputMags)>1 && length(inputMags2)==1)
figure; plot(1./figISIs','b.-','Markersize',16);
xlabel('Interspike interval #','fontsize',12); ylabel('1/ISI duration (Hz)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'ylim',[min(desiredFreqs)-0.1 max(desiredFreqs)+0.5]);
set(gca,'fontsize',12,'box','off')
if slowslow
if saveFigures && exist(['fig_Ackerslow_history_dependence_multi_0a.eps'])~=2
print_eps(['fig_Ackerslow_history_dependence_multi_0a.eps'])
saveas(gcf,['fig_Ackerslow_history_dependence_multi_0a.fig'])
end
else
if saveFigures && exist(['fig_Acker_history_dependence_multi_0a.eps'])~=2
print_eps(['fig_Acker_history_dependence_multi_0a.eps'])
saveas(gcf,['fig_Acker_history_dependence_multi_0a.fig'])
end
end
elseif length(inputMags)>1 && length(inputMags2)>1
figure; imagesc(desiredFreqs,desiredFreqs,abs(real(figISIs)-imag(figISIs)));
title('ISI overshoot magnitude (ms)','fontsize',12)
xlabel('Starting frequency (Hz)','fontsize',12); ylabel('Ending frequency (Hz)','fontsize',12)
set(gcf,'position',[82 404 322 302])
set(gca,'fontsize',12,'box','off')
if saveFigures && exist(['fig_RS' typePrefix '_histdep_ISI_colorplot.eps'])~=2
print_eps(['fig_RS' typePrefix '_histdep_ISI_colorplot.eps'])
saveas(gcf,['fig_RS' typePrefix '_histdep_ISI_colorplot.fig'])
end
figure; imagesc(desiredFreqs,desiredFreqs,abs(1./real(figISIs)-1./imag(figISIs)));
title('Frequency overshoot magnitude (Hz)','fontsize',12)
xlabel('Starting frequency (Hz)','fontsize',12); ylabel('Ending frequency (Hz)','fontsize',12)
set(gcf,'position',[82 404 322 302])
set(gca,'fontsize',12,'box','off')
if saveFigures && exist(['fig_RS' typePrefix '_histdep_freq_colorplot.eps'])~=2
print_eps(['fig_RS' typePrefix '_histdep_freq_colorplot.eps'])
saveas(gcf,['fig_RS' typePrefix '_histdep_freq_colorplot.fig'])
end
end
