% 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 histogramsfortype=[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;
iftype==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;
elseiftype==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;
elseiftype==3
excitationClass = 2;
useGapJunctions = 0;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 0:25:450;
allncells = 250;
nruns = 10;
elseiftype==4
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 2.^[-6:0.5:4];
allncells = 250;
nruns = 10;
elseiftype==5
excitationClass = 2;
useGapJunctions = 0;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 0;
allncells = 1;
nruns = 1;
simdur = 60000; % ms
numISIHists = 1;
elseiftype==6
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 0.1;
allncells = 250;
nruns = 1;
simdur = 60000; % ms
numISIHists = 1;
elseiftype==7
excitationClass = 2;
useGapJunctions = 0;
uniqueNoiseSTD = 100*useNoise;
I = 95.8;
gs = 256;
allncells = 250;
nruns = 1;
simdur = 60000; % ms
numISIHists = 1;
elseiftype==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 = [163264128256512];
gs = 0.1;
elseiftype==9
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 100*useNoise;
commonNoiseSTD = 0*useNoise;
I = 95.8;
gs = 0.1;
allncells = 250;
elseiftype==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;
elseiftype==11
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 0*useNoise;
commonNoiseSTD = sqrt(0.05)*100*useNoise;
I = 95.8;
gs = 0;
allncells = 1;
elseiftype==12
excitationClass = 2;
useGapJunctions = 1;
uniqueNoiseSTD = 0*useNoise;
commonNoiseSTD = sqrt(0.05)*100*useNoise;
I = 95.8;
gs = 0.1;
allncells = 2500;
elseiftype==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;
elseiftype==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 neuronsif excitationClass==1
b=-2; % integratorelseif excitationClass==2
b = 2; % resonatorend
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>1if useGapJunctions
typeString = [typeString 'g'];
else
typeString = [typeString 's'];
endendif 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 matrixif pcon>0
C = (rand(ncells)<pcon)>0;
C = sparse(C.*(1-eye(ncells))); % remove autoconnectionselse
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;breakend
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;endelseif 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;
endend
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 0end
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',[82404257185]);
fname = [simprefix '_vco_isihist_' num2str(ce)];
if saveFigures && exist([fname '.eps'])~=2
print_eps([fname '.eps'])
saveas(gcf,[fname '.fig'])
endend% if we have 3 or fewer ISIs, let's just trash this cell:iflength(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));
endend% 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 neuronsfor 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',[82404257185]);
fname = [simprefix '_post_isihist'];
if saveFigures && exist([fname '.eps'])~=2
print_eps([fname '.eps'])
saveas(gcf,[fname '.fig'])
endend
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));
endfor 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));
endfor ce=1:ncells-length(indivmeans)
fprintf(fid,'\tNaN');
end
fprintf(fid,'\t');
for ce=1:length(indivmeans)
fprintf(fid,'\t%g',indivstabilities(ce));
endfor ce=1:ncells-length(indivmeans)
fprintf(fid,'\tNaN');
end
fprintf(fid,'\n');
fclose(fid);
endend
medians(ncellsi,pind,gi) = median(stabilities);
maxs(ncellsi,pind,gi) = max(stabilities);
mins(ncellsi,pind,gi) = min(stabilities);
endendend% end ncells loopend% 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',[0500])
set(gca,'ylim',[013])
set(gca,'fontsize',12,'box','off')
xlabel('Synaptic conductance','fontsize',12); ylabel('Network frequency (Hz)','fontsize',12)
set(gcf,'position',[82404257302])
if exist('fig_freq_varyg.eps')~=2
print_eps('fig_freq_varyg.eps')
saveas(gcf,'fig_freq_varyg.fig')
endfigure; plot(data(:,4),data(:,21),'.'); xlabel('Synaptic conductance','fontsize',12); ylabel('Estimated stability time (s)','fontsize',12)
set(gcf,'position',[82404257302])
set(gca,'xlim',[0500])
set(gca,'ylim',[0350])
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')
endfigure; plot(data(:,4),data(:,20),'.'); xlabel('Synaptic conductance','fontsize',12); ylabel('Network period std. dev. (s)','fontsize',12)
set(gcf,'position',[82404257302])
set(gca,'xlim',[0500])
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-210])
set(gca,'xtick',[.01.1110])
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',[82404257302])
if exist('fig_freq_varyg.eps')~=2
print_eps('fig_freq_varyg.eps')
saveas(gcf,'fig_freq_varyg.fig')
endfigure; plot(data(:,4),data(:,21),'.'); xlabel('Gap-junction conductance','fontsize',12); ylabel('Estimated stability time (s)','fontsize',12)
set(gcf,'position',[82404257302])
set(gca,'xlim',[1e-220])
set(gca,'xtick',[.01.1110])
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')
endfigure; plot(data(:,4),data(:,20),'.'); xlabel('Gap-junction conductance','fontsize',12); ylabel('Network period std. dev. (s)','fontsize',12)
set(gcf,'position',[82404257302])
set(gca,'xlim',[1e-220])
set(gca,'xtick',[.01.1110])
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',[400200100502512.5],'xtick',[1:6],'xticklabel',[163264128256512]);
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',[82404322275])
if exist('fig_RS2gn_varynnoise.eps')~=2
print_eps('fig_RS2gn_varyn.eps')
saveas(gcf,'fig_RS2gn_varyn.fig')
end