% izhikevich simple model neuron network with random connectivity
% with integrate and fire postsynaptic cell
% eric zilli - june 27, 2009 using simple model MATLAB code from
% Izhikevich's book modified to be a network
%
% release version 1.0. check modeldb for updates.
%
% this source code is released into the public domain
%
% This script runs a network of optionally-noisy cells while the injected
% current is held constant. It allows the number of cells, the coupling
% type and strength, the connectivity probability, amount of noise, etc. to
% be changed.
%
% 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
warning off
% type=1 - class 1 excitable, n=100, synaptic, varying g
% type=2 - class 1 excitable, n=100, gap-junction, varying g
% type=3 - (Fig S4A-C) class 2 excitable, n=250, synaptic, varying g
% type=4 - (Fig S4D-F) class 2 excitable, n=250, gap-junction, varying g
% type=5 - (Fig 5A,D) class 2 excitable, n = 1, ISI histograms
% type=6 - (Fig 5B,E) class 2 excitable, n = 250, gap-junction, ISI histograms
% type=7 - (Fig 5C,F) class 2 excitable, n = 250, synaptic, ISI histograms
% type=8 - (Fig S8) class 2 excitable, varying n and noise, gap-junction g=0.1
% type=9 - (Discussion) class 2 excitable, n=250, gap-junctions, 100% indep. noise, 0 corr. noise
% type=10 - (Discussion) class 2 excitable, n=250, gap-junctions, 95% indep. noise var, 5% corr. noise var
% type=11 - (Discussion) class 2 excitable, n=1, 5% realistic noise var
% type=12 - (Discussion) class 2 excitable, n=2500, gap-junctions, 5% realistic corr. noise var
% type=13 - class 2 excitable, n=250, inhibitory synaptic, varying g
% type=14 - (Fig 8D,H) class 2 excitable, n = 5000, p=0.01, synaptic, ISI histograms
for type=[1]
% Note: some parameters may be overwritten in the 'if type' block:
simdur = 8000; % (ms)
dt = .1; % ms
nruns = 1;
% for statistics about each run
fileOut = 0;
% for ISI histograms:
saveFigures = 0;
numISIHists = 0;
% these allow initial transients and such to be ignored:
startAnalysisTime = dt; % (ms)
endAnalysisTime = simdur; % (ms)
% sets whether or not noise is used
useNoise = 1;
% approximate probability any two cells are connected (autoconnections will
% be removed)
% set to 0 (or set gs=0) to uncouple cells to, e.g., find noise level and
% input level such that the median matches the biological data
pcons = 1;
weightmult = 1;
% std. dev. of per-step noise term correlated among all network cells:
commonNoiseSTD = 0*useNoise;
commonNoise = 0;
if type==1
excitationClass = 1;
useGapJunctions = 0;
uniqueNoiseSTD = 130*useNoise;
I = 62;
% rough range of good gs = 150:1180
gs = 150:50:1300;
allncells = 100;
nruns = 10;
elseif type==2
excitationClass = 1;
useGapJunctions = 1;
uniqueNoiseSTD = 130*useNoise;
I = 62;
% rough range of good gs = 0.022:16
gs = 2.^[-5.5:0.5:4];
allncells = 100;
nruns = 10;
elseif type==3
excitationClass = 2;
useGapJunctions = 0;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 0:25:450;
allncells = 250;
nruns = 10;
elseif type==4
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 2.^[-6:0.5:4];
allncells = 250;
nruns = 10;
elseif type==5
excitationClass = 2;
useGapJunctions = 0;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 0;
allncells = 1;
nruns = 1;
simdur = 60000; % ms
numISIHists = 1;
elseif type==6
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 0.1;
allncells = 250;
nruns = 1;
simdur = 60000; % ms
numISIHists = 1;
elseif type==7
excitationClass = 2;
useGapJunctions = 0;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 256;
allncells = 250;
nruns = 1;
simdur = 60000; % ms
numISIHists = 1;
elseif type==8
excitationClass = 2;
useGapJunctions = 1;
%% run one at a time:
uniqueNoiseSTD = 400*useNoise;
% uniqueNoiseSTD = 200*useNoise;
% uniqueNoiseSTD = 100*useNoise; % original level
% uniqueNoiseSTD = 50*useNoise;
% uniqueNoiseSTD = 25*useNoise;
% uniqueNoiseSTD = 12.5*useNoise;
I = 125;
allncells = [16 32 64 128 256 512];
gs = 0.1;
elseif type==9
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 100*useNoise;
commonNoiseSTD = 0*useNoise;
I = 95.8;
gs = 0.1;
allncells = 250;
elseif type==10
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = sqrt(0.95)*100*useNoise;
commonNoiseSTD = sqrt(0.05)*100*useNoise;
I = 95.8;
gs = 0.1;
allncells = 250;
elseif type==11
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 0*useNoise;
commonNoiseSTD = sqrt(0.05)*100*useNoise;
I = 95.8;
gs = 0;
allncells = 1;
elseif type==12
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 0*useNoise;
commonNoiseSTD = sqrt(0.05)*100*useNoise;
I = 95.8;
gs = 0.1;
allncells = 2500;
elseif type==13
excitationClass = 2;
useGapJunctions = 0;
useNoise = 1;
uniqueNoiseSTD = 100*useNoise;
I = 110;
gs = -(2.^[9:15]); % inhibitory synapses
allncells = 1000;
nruns = 1;
simdur = 15000;
weightmult = 8;
endAnalysisTime=simdur;
elseif type==14
excitationClass = 2;
useGapJunctions = 0;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 256;
allncells = 5000;
nruns = 1;
simdur = 6000;0; % ms
pcon = 0.01;
numISIHists = 1;
end
% file to save statistics for every individual cell in network, don't
% save if indivname = ''
indivname = '';
% indivname = [strtok(outfilename,'.') '_individualcells.txt'];
Cf=100; vr=-60; vt=-40; k=0.7; % parameters used for RS
a=0.03; c=-50; d=100; % neocortical pyramidal neurons
if excitationClass==1
b=-2; % integrator
elseif excitationClass==2
b = 2; % resonator
end
vpeak=35; % spike cutoff
% store min/med/max stability time of individual cells in each simulation
mins = zeros(length(allncells),length(pcons),length(gs));
medians = zeros(length(allncells),length(pcons),length(gs));
maxs = zeros(length(allncells),length(pcons),length(gs));
for ncellsi=1:length(allncells)
ncells = allncells(ncellsi);
typeString = '';
if ncells>1
if useGapJunctions
typeString = [typeString 'g'];
else
typeString = [typeString 's'];
end
end
if useNoise
typeString = [typeString 'n'];
end
simprefix = sprintf('simple_model_RS%s_n%d_%s',typeString,ncells,datestr(now,'mmmdd'))
if ~fileOut
disp('no fileout')
end
outfilename = ['simple_model_RS' sprintf('%d',excitationClass) typeString];
outfilename = [outfilename sprintf('_LIF_%s_vary_gandnoise_%gnoise.txt',datestr(now,'mmmdd'),useNoise)];
for gi=1:length(gs)
g = gs(gi);
for pind=1:length(pcons)
pcon = pcons(pind);
stabilities = zeros(1,nruns);
ISIstabilities = zeros(1,nruns);
for run=1:nruns
% activity of postsynaptic I&F
post = zeros(1, round(simdur/dt));
% spike times of postsynaptic I&F
postspikes = spalloc(1, round(simdur/dt), 6*simdur); % expect firing at 6 Hz
% LIF time constant
tau = 4.5; % (ms)
membraneDecay = exp(-dt/tau);
postWeight=weightmult*3/ncells;
% connectivity matrix
if pcon>0
C = (rand(ncells)<pcon)>0;
C = sparse(C.*(1-eye(ncells))); % remove autoconnections
else
C = spalloc(ncells,ncells,1);
end
vhist = zeros(1,round(simdur/dt));
t = 0; % current time in simulation
v = vr*ones(ncells,1); % vr*rand(ncells,1);
u = 0*v; % initial values
% save state of one cell over simulation:
state = zeros(2,simdur/dt);
% binary array indicating which cells fire on each time steps
spikes = sparse(ncells, simdur/dt,ncells*simdur); % sparse to save memory in big simulations
tic
while t<simdur
t = t+dt; % advance to next time step
tind = 1+round(t/dt);
if commonNoiseSTD
commonNoise = commonNoiseSTD*randn;
end
% if our LIF neuron is spiking constantly, let's just skip this
% one (there is no apparent synchrony, which may either reflect the network
% itself or the LIF parameters)
if tind==500 && sum(postspikes)>(497)
% postspikes = NaN;
break
end
oldv = v;
oldu = u;
if useGapJunctions
if pcon==1
% if all-to-all connected, the gap-junctions are equivalent
% to sum(v)-ncells*v
v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(v)-ncells*v))/Cf;
else
% make matrix of differences in membrane potential
vdiffs = C.*(repmat(oldv',ncells,1) - repmat(oldv,1,ncells));
v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(vdiffs,2))/Cf;
%% slowest way:
% v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*sum(C.*vdiffs,2))/Cf;
end
else
if pcon==1
v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(sum(spikes(:,tind-1),1)-spikes(:,tind-1)))/Cf;
else
v = oldv + dt*(k*(oldv-vr).*(oldv-vt) - oldu + I + commonNoise + uniqueNoiseSTD*randn(ncells,1) + g*(C*spikes(:,tind-1)))/Cf;
end
end
u = oldu + dt*a*(b*(oldv-vr)-oldu);
vhist(tind) = v(1);
% save and reset spikes when v>=vpeak
spikes(:,tind) = v>=vpeak;
u(v>=vpeak) = u(v>=vpeak)+d;
oldv(v>=vpeak) = vpeak;
v(v>=vpeak) = c;
state(:,tind) = [oldv(1); oldu(1)];
post(:,tind) = post(:,tind-1)*membraneDecay + tau*postWeight*sum(spikes(:,tind));
postspikes(:,tind) = post(:,tind)>1;
post(postspikes(:,tind)>1,tind) = 0; % reset to 0
end
toc
clear C
% calculate mean and variance of ISIs for each postsynaptic cell
npost = 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);
for ce=1:ncells % only run this once because we only save the last value in stabilities
spikeTimes = find(spikes(ce,(startAnalysisTime/dt):(endAnalysisTime/dt))>0);
ISIs = dt*diff(spikeTimes)/1000; % convert to seconds
% count ISIs as occuring between initial
% spikes of bursts (and additional spikes as occuring in a given
% burst if they are < 50 ms apart);
fISIs = [];
csum = 0;
for is=1:length(ISIs)
csum = csum + ISIs(is);
if ISIs(is)>.050
fISIs = [fISIs csum];
csum=0;
end;
end
ISIs = fISIs;
if ce<=numISIHists
figure; hist(ISIs,linspace(.1,.3,200))
title(['VCO Cell #' num2str(ce) ' ISIs'],'fontsize',12)
xlabel('ISI duration (s)','fontsize',12)
set(gca,'fontsize',12);
set(gca,'box','off');
set(gcf,'position',[82 404 257 185]);
fname = [simprefix '_vco_isihist_' num2str(ce)];
if saveFigures && exist([fname '.eps'])~=2
print_eps([fname '.eps'])
saveas(gcf,[fname '.fig'])
end
end
% if we have 3 or fewer ISIs, let's just trash this cell:
if length(ISIs)<=3
indivmeans(ce) = NaN;
indivstds(ce) = NaN;
indivstabilities(ce) = NaN;
else
indivmeans(ce) = mean(ISIs);
indivstds(ce) = std(ISIs);
indivstabilities(ce) = 5*indivmeans(ce)^3/16/pi/pi/(eps+indivstds(ce))/(eps+indivstds(ce));
end
end
% drop NaNs and sort:
indivmeans(isnan(indivmeans)) = [];
indivstds(isnan(indivstds)) = [];
indivstabilities(isnan(indivstabilities)) = [];
[indivstabilities,ix] = sort(indivstabilities);
indivmeans = indivmeans(ix);
indivstds = indivstds(ix);
%% calculate ISI statistics for the postsynaptic LIF neurons
for ce=1:npost % only run this once because we only save the last value in stabilities
spikeTimes = find(postspikes(ce,(startAnalysisTime/dt):(endAnalysisTime/dt))>0);
ISIs = dt*diff(spikeTimes)/1000; % convert to seconds
% count ISIs as occuring between initial
% spikes of bursts (and additional spikes as occuring in a given
% burst if they are < 50 ms apart);
fISIs = [];
csum = 0;
for is=1:length(ISIs)
csum = csum + ISIs(is);
if ISIs(is)>.050
fISIs = [fISIs csum];
csum=0;
end;
end
ISIs = fISIs;
if ce==1 && numISIHists
figure; hist(ISIs,linspace(.1,.3,200))
title('Postsynaptic LIF ISIs','fontsize',12)
xlabel('ISI duration (s)','fontsize',12)
set(gca,'fontsize',12)
set(gca,'box','off')
set(gcf,'position',[82 404 257 185]);
fname = [simprefix '_post_isihist'];
if saveFigures && exist([fname '.eps'])~=2
print_eps([fname '.eps'])
saveas(gcf,[fname '.fig'])
end
end
ISImeans(ce) = mean(ISIs);
ISIstds(ce) = std(ISIs);
ISIstabilities(ce) = 5*ISImeans(ce)^3/16/pi^2/(eps+ISIstds(ce))^2;
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%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,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
Cf,vr,vt,k,a,b,c,d,vpeak,...
1/ISImeans(1),ISImeans(1),ISIstds(1),ISIstabilities(1),length(ISIs),...
min(indivmeans),median(indivmeans),max(indivmeans),...
min(indivstds),median(indivstds),max(indivstds),...
min(indivstabilities),median(indivstabilities),max(indivstabilities),...
length(find(postspikes>0)));
fclose(fid);
end
% print to screen:
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%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,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
Cf,vr,vt,k,a,b,c,d,vpeak,...
1/ISImeans(1),ISImeans(1),ISIstds(1),ISIstabilities(1),length(ISIs), ...
min(indivmeans),median(indivmeans),max(indivmeans),...
min(indivstds),median(indivstds),max(indivstds),...
min(indivstabilities),median(indivstabilities),max(indivstabilities),...
length(find(postspikes>0)));
% optinally save individual cell statistics to a file:
if ~isempty(indivname) && fileOut
fid = fopen(indivname,'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%d', ...
ncells,pcon,I,gs(gi),simdur/1000, useNoise*uniqueNoiseSTD,tau,postWeight,...
1/ISImeans(1),ISImeans(1),ISIstds(1),ISIstabilities(1),length(ISIs),...
length(find(postspikes>0)));
for ce=1:length(indivmeans)
fprintf(fid,'\t%g',indivmeans(ce));
end
for ce=1:ncells-length(indivmeans) % for alignment in ouput, put back the NaNs we removed
fprintf(fid,'\tNaN');
end
fprintf(fid,'\t');
for ce=1:length(indivmeans)
fprintf(fid,'\t%g',indivstds(ce));
end
for ce=1:ncells-length(indivmeans)
fprintf(fid,'\tNaN');
end
fprintf(fid,'\t');
for ce=1:length(indivmeans)
fprintf(fid,'\t%g',indivstabilities(ce));
end
for ce=1:ncells-length(indivmeans)
fprintf(fid,'\tNaN');
end
fprintf(fid,'\n');
fclose(fid);
end
end
medians(ncellsi,pind,gi) = median(stabilities);
maxs(ncellsi,pind,gi) = max(stabilities);
mins(ncellsi,pind,gi) = min(stabilities);
end
end
end % end ncells loop
end % end run loop
warning on
return
%% Figure 1
data = textread('simple_model_RS2sn_vary_gandncells_1noise.txt');
figure; plot(data(:,4),data(:,18),'.');
set(gca,'xlim',[0 500])
set(gca,'ylim',[0 13])
set(gca,'fontsize',12,'box','off')
xlabel('Synaptic conductance','fontsize',12); ylabel('Network frequency (Hz)','fontsize',12)
set(gcf,'position',[82 404 257 302])
if exist('fig_freq_varyg.eps')~=2
print_eps('fig_freq_varyg.eps')
saveas(gcf,'fig_freq_varyg.fig')
end
figure; plot(data(:,4),data(:,21),'.'); xlabel('Synaptic conductance','fontsize',12); ylabel('Estimated stability time (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[0 500])
set(gca,'ylim',[0 350])
set(gca,'fontsize',12,'box','off')
if exist('fig_stability_varyg.eps')~=2
print_eps('fig_stability_varyg.eps')
saveas(gcf,'fig_stability_varyg.fig')
end
figure; plot(data(:,4),data(:,20),'.'); xlabel('Synaptic conductance','fontsize',12); ylabel('Network period std. dev. (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[0 500])
set(gca,'yscale','log')
set(gca,'fontsize',12,'box','off')
if exist('fig_stdev_varyg.eps')~=2
print_eps('fig_stdev_varyg.eps')
saveas(gcf,'fig_stdev_varyg.fig')
end
data = textread('simple_model_RS2gn_vary_gandncells_1noise.txt');
figure; plot(data(:,4),data(:,18),'.');
set(gca,'xlim',[1e-2 10])
set(gca,'xtick',[.01 .1 1 10])
set(gca,'xscale','log')
set(gca,'fontsize',12,'box','off')
% set(gca,'ylim',[0 9.1])
xlabel('Gap-junction conductance','fontsize',12); ylabel('Network frequency (Hz)','fontsize',12)
set(gcf,'position',[82 404 257 302])
if exist('fig_freq_varyg.eps')~=2
print_eps('fig_freq_varyg.eps')
saveas(gcf,'fig_freq_varyg.fig')
end
figure; plot(data(:,4),data(:,21),'.'); xlabel('Gap-junction conductance','fontsize',12); ylabel('Estimated stability time (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[1e-2 20])
set(gca,'xtick',[.01 .1 1 10])
set(gca,'xscale','log')
set(gca,'fontsize',12,'box','off')
if exist('fig_stability_varyg.eps')~=2
print_eps('fig_stability_varyg.eps')
saveas(gcf,'fig_stability_varyg.fig')
end
figure; plot(data(:,4),data(:,20),'.'); xlabel('Gap-junction conductance','fontsize',12); ylabel('Network period std. dev. (s)','fontsize',12)
set(gcf,'position',[82 404 257 302])
set(gca,'xlim',[1e-2 20])
set(gca,'xtick',[.01 .1 1 10])
set(gca,'xscale','log')
set(gca,'yscale','log')
set(gca,'fontsize',12,'box','off')
if exist('fig_stdev_varyg.eps')~=2
print_eps('fig_stdev_varyg.eps')
saveas(gcf,'fig_stdev_varyg.fig')
end
%% Figure 2
% histogram plots are buried in the script above, search for numISIHists
%% Figure 3
% first run type=8 for each standard deviation of the noise, saving output
% we expect all files have the same number of rows
nfiles = 6; % i.e. number of noise levels
startrow = 0;
endrow = 6;
ISIfreqs = zeros(endrow-startrow,nfiles);
ISImeans = zeros(endrow-startrow,nfiles);
ISIstds = zeros(endrow-startrow,nfiles);
ISIstability = zeros(endrow-startrow,nfiles);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_4noise.txt');
ISIfreqs(:,1) = data((startrow+1):endrow,18);
ISImeans(:,1) = data((startrow+1):endrow,19);
ISIstds(:,1) = data((startrow+1):endrow,20);
ISIstability(:,1) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_2noise.txt');
ISIfreqs(:,2) = data((startrow+1):endrow,18);
ISImeans(:,2) = data((startrow+1):endrow,19);
ISIstds(:,2) = data((startrow+1):endrow,20);
ISIstability(:,2) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_1noise.txt');
ISIfreqs(:,3) = data((startrow+1):endrow,18);
ISImeans(:,3) = data((startrow+1):endrow,19);
ISIstds(:,3) = data((startrow+1):endrow,20);
ISIstability(:,3) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_0.5noise.txt');
ISIfreqs(:,4) = data((startrow+1):endrow,18);
ISImeans(:,4) = data((startrow+1):endrow,19);
ISIstds(:,4) = data((startrow+1):endrow,20);
ISIstability(:,4) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_0.25noise.txt');
ISIfreqs(:,5) = data((startrow+1):endrow,18);
ISImeans(:,5) = data((startrow+1):endrow,19);
ISIstds(:,5) = data((startrow+1):endrow,20);
ISIstability(:,5) = data((startrow+1):endrow,21);
data = textread('simple_model_RS2gn_LIF_vary_gandnoise_0.125noise.txt');
ISIfreqs(:,6) = data((startrow+1):endrow,18);
ISImeans(:,6) = data((startrow+1):endrow,19);
ISIstds(:,6) = data((startrow+1):endrow,20);
ISIstability(:,6) = data((startrow+1):endrow,21);
ncellsvect = data((startrow+1):endrow,1);
figure; imagesc(log10(ISIstability(1:end,:))');
set(gca,'yticklabel',[400 200 100 50 25 12.5],'xtick',[1:6],'xticklabel',[16 32 64 128 256 512]);
xlabel('Number of cells','fontsize',12); ylabel('Noise standard deviation','fontsize',12);
title('Estimated stability time (log s)','fontsize',12)
set(gca,'fontsize',12)
h=colorbar; set(h,'fontsize',12);
colormap(flipud(gray));
set(gcf,'position',[82 404 322 275])
if exist('fig_RS2gn_varynnoise.eps')~=2
print_eps('fig_RS2gn_varyn.eps')
saveas(gcf,'fig_RS2gn_varyn.fig')
end
