clear global Network
global Network
disp('Preparing simulations')
tic
% --- Time integration ---
Network.dt = 0.01; % Integration step (in ms)
simdur = 2e3; % ms of simulation time
Network.nStep = round(simdur/Network.dt); % Number of time steps
Network.OutputVm = 0; % 1=Output membrane potentials to Vm.txt, 0=no such output
Network.OutputSp = 1; % 1=Output spike times to Spikes.txt, 0=no such output
Network.OutputSynState = 0; % 1=Output synapse state to SynState.txt, 0=no such output
fInit = fopen('Inits.txt','w');
fprintf(fInit,'%u\t%g\n%u\t%u\t%u',Network.nStep,Network.dt,Network.OutputVm,Network.OutputSp,Network.OutputSynState);
fclose(fInit);
% --- Populations ---
% generate Ncell cells, unconnected, that all have a constant synaptic
% input
Ncell = 400;
Iinput = 0;
gsyn = sort([linspace(-0.2,3,(Ncell-100)), linspace(-0.1,0.1,100)]);
AddPopulation(...
'Name','Ex',...
'nCell',Ncell,...
'Type','HH',...
'Iinput',Iinput ,'Cm',10.0,'g_na',100.0,'g_naL',0.0175,'g_k',40.0,'g_kL',0.05,'g_clL',0.05,'phi',3.0,'E_k',Ek,'E_na',53,'E_cl',-82,...
'gsyntest', gsyn, 'Esyntest', 50,...
'V',-50,'n',0.07,'h',0.97,...
'Position',repmat(linspace(0,10,Ncell)',1,3)...
);
% --- Output to files
WriteCells;
WriteConnectivity;
toc
% --- run simulation in C++ environment
disp('Running simulations')
tic
if(ismac)
system('VerdandiLite_v1_0/VerdandiLite_OSX');
elseif(ispc)
%system('VerdandiLite_v1_0\VerdandiLite_WIN32.exe');
system('VerdandiLite_v1_0\VerdandiLite_CPP\BJZ.exe');
else
error('No binaries available for VerdandiLite. Please compile C++ source files and add modify these Matlab lines.')
end
SkuldOut % create output for analysis in Skuld (not used)
toc
Tana = [0,Inf]; % [ms] time span for analysis can be changed
Cellnums = (1:Ncell)-1;
method = 2; %use method 2, when input is constant, else use 1
Spikes = importdata('Spikes.txt');
% use only spikes in time span to analyze
ind = find(Spikes(:,1)>Tana(1) & Spikes(:,1)<Tana(2));
S = Spikes(ind,:);
S2 = sortrows(S,[2,1]);
switch method
case 1 % average firing rate during entire simulation
Nspike = hist(S(:,2),Cellnums);
DT = (Tana(2)-Tana(1))/1e3;
f = Nspike/DT;
case 2 % spike rate at and of simulation ( 1/(last interspike interval) )
f = zeros(size(Cellnums));
Nspike = hist(S(:,2),Cellnums);
for n = find(Nspike>2)
lastspikepos = find(S2(:,2) == Cellnums(n),1,'last');
f(n) = 1e3/(S2(lastspikepos,1)-S2(lastspikepos-1,1)); % calculate difference in spike times between last and secondlast spike
end
end
%% plot result
g = sort([linspace(-0.2,3,(Ncell*10-1000)), linspace(-0.1,0.1,1000)]);
F = interp1(gsyn,f,g,'pchip');
plot(gsyn,f,'r.'); hold on; xlabel('Synaptic conductane (mS/cm^2)'), ylabel('Spike rate')
plot(g,F); hold on; xlabel('Synaptic conductane (?/cm^2)'), ylabel('Spike rate')
%% save result
%