function [Synch, Vm, S, D, instISI, spklist, LamRatio, Lambda2, LambdaVar] = ML_Net_Sweep(T, N, k, rho, I, SynStren, coefs, Syns);
rand('seed',sum(100*clock));
randn('seed',sum(100*clock));
load('ber_3000_0.010_4.0_1.4_1.3_1.2.dat');
Syns=cell(3000,1);
for i=1:1:3000
Syns{i}=find(ber_3000_0_010_4_0_1_4_1_3_1_2(i,:));
end
if(nargin<1)
T = 1000;
end;
if(nargin<2)
N=3000;
end;
if(nargin<3)
k=30;
end;
if(nargin<4)
rho=0.01;
end;
if(nargin<5)
I=0;
end;
if(nargin<6)
SynStren=0.01;
end;
dt = 0.02;
NumSamples = T/dt;
NoiseScalar = (sqrt(dt)*0.2);
lambda2=0;
LambdaN=0;
LambdaVar=0;
if(0)
if(iscell(Syns))
K=length(Syns{1});
CellNum = [];
PostSyns=[];
for i=1:N
K=length(Syns{i});
CellNum=[CellNum; i*ones(K,1)];
PostSyns = [PostSyns; Syns{i}'];
end;
SparseSyns=sparse(CellNum, PostSyns,ones(length(CellNum),1));
else
CellNum=repmat(1:N,1,k);
SparseSyns=sparse([1:N, CellNum], [1:N, Syns(:)'],[-30*ones(N,1); ones(k*N,1)], N, N);
SparseSyns=sparse(CellNum, Syns(:), ones(k*N,1), N, N);
SparseSyns(find(SparseSyns>1))=1;
end;
SparseSyns=SparseSyns-sparse(diag(sum(SparseSyns,1)));
D = eig(full(SparseSyns));
figure
plot(D,'.');
xlabel('Real(\lambda)');
ylabel('Imag(\lambda)');
LambdaVar = sum(abs(D(2:end)-mean(D(2:end))).^2)/(mean(D(2:end))^2*(N-1));
LamRatio=min(real(D(abs(D)>1e-12)))/max(real(D(abs(D)>1e-12)))
Lambda2=max(real(D(abs(D)>1e-12)));
end;
if(nargin<7)
coefs.C=(1);
coefs.gL=(8);
coefs.EL=(-53.23878);
coefs.gNa= (18.22315);
coefs.ENa=(60);
coefs.gK=(4);
coefs.EK=(-95.52116);
coefs.Vhalfm=(-7.37257);
coefs.km=(11.97239);
coefs.Vhalfn=(-16.34877);
coefs.kn=(4.21133);
coefs.tau=(1);
coefs.E_EPSP = (0);
coefs.tauEPSPr = (0.25);
coefs.tauEPSPf = (.5);
coefs.gEPSP=SynStren/( coefs.tauEPSPf - coefs.tauEPSPr);
end;
[t,vinit] = ode23(@ML_derivs,[0:0.01:100],[-60, 0, 0, 0],[],I, coefs);
[ind, spkind] = findpeaks(vinit(:,1),'minpeakheight',-20);
PredNumSpikes = length(ind)*N*1.2*T/100;
ISI = 100/length(ind);
ExpSynCurrent = 1/ISI*k*SynStren*(-mean(vinit(:,1)));
t = t(spkind(end-1):spkind(end));
t=t-t(1);
vinit=vinit(spkind(end-1):spkind(end),:);
ind = ceil((length(t))*rand(N,1));
V=(vinit(ind,1));
n=(vinit(ind,2));
lastSpkTime = -t(ind);
Ser = zeros(N,1);
Sef = zeros(N,1);
if(nargout>5)
LengthSpkList = PredNumSpikes;
spklist = zeros(LengthSpkList,2);
end
Vm = zeros(NumSamples,1);
S = zeros(NumSamples,1);
Synch = zeros(NumSamples,1);
spknum=0;
t=0;
dVdt = zeros(N,1);
dndt = zeros(N,1);
dssdt = zeros(N,1);
dsfdt = zeros(N,1);
lastV = zeros(N,1);
ISI = ISI*ones(N,1);
D = zeros(T/dt,1);
instISI = zeros(T/dt,1);
thresh = -20;
tempV=linspace(coefs.EK, coefs.ENa, 1024);
tempV = tempV';
NInfLUT = 1.0./(1.0+exp((coefs.Vhalfn - tempV)/coefs.kn ));
NTauLUT = exp(-0.07*tempV-3);
MInfLUT = (1.0./(1.0 + exp((coefs.Vhalfm - tempV)/coefs.km )));
LUTLength = length(MInfLUT);
PhaseLUT = exp(1i*linspace(0,2*pi,1024));
PhaseLUTLen=length(PhaseLUT);
SynStrenSweep=linspace(0.01,0.01,NumSamples);
Tdrive=ones(1,ceil(NumSamples/3))*8;
I=linspace(8,-4,floor(2*NumSamples/3));
I=horzcat(Tdrive,I);
for i=1:NumSamples;
coefs.gEPSP=SynStren/( coefs.tauEPSPf - coefs.tauEPSPr);
Vind = ceil(LUTLength*(V-coefs.EK)/(coefs.ENa-coefs.EK));
dVdt = 1.0./coefs.C*(I(i)-ExpSynCurrent-coefs.gL.*(V-coefs.EL) - ...
coefs.gNa.*MInfLUT(Vind).*(V-coefs.ENa) - ...
coefs.gK.*n.*(V-coefs.EK) + ...
coefs.gEPSP.*(Sef-Ser).*(coefs.E_EPSP-V));
dndt = (NInfLUT(Vind)-n)./NTauLUT(Vind);
dssdt = -Ser/coefs.tauEPSPr;
dsfdt = -Sef/coefs.tauEPSPf;
lastV = V;
V = V + dVdt*dt + (randn(N,1))*NoiseScalar;
n = n + dndt*dt;
Ser = Ser + dssdt*dt;
Sef = Sef + dsfdt*dt;
spks = find((V>thresh).*(lastV<thresh));
if(~isempty(spks))
ISI(spks)=t-lastSpkTime(spks);
lastSpkTime(spks) = t;
for spk=1:length(spks);
if(iscell(Syns))
Ser(Syns{spks(spk),:})=Ser(Syns{spks(spk),:})+1;
Sef(Syns{spks(spk),:})=Sef(Syns{spks(spk),:})+1;
else
Ser(Syns(spks(spk),:))=Ser(Syns(spks(spk),:))+1;
Sef(Syns(spks(spk),:))=Sef(Syns(spks(spk),:))+1;
end;
end;
if(nargout>5)
spklist(spknum+1:spknum+1+length(spks),1)=t;
spklist(spknum+1:spknum+length(spks),2)=spks;
spknum=spknum+length(spks);
end;
end;
Vm(i)=V(1);
S(i)=(coefs.gEPSP.*(Sef(1)-Ser(1)).*(coefs.E_EPSP-V(1)));
Phase = mod(ceil((PhaseLUTLen-1)*(t-lastSpkTime)/mean(ISI)),PhaseLUTLen)+1;
Synch(i) = abs(mean(PhaseLUT(Phase)));
D(i)=length(spks);
instISI(i)=mean(ISI);
t=t+dt;
end;
if(nargout>5)
if spknum<length(spklist);
spklist=spklist(1:spknum,:);
end;
end;
figure;
subplot 311
plot(dt:dt:T,D,'k');
ylabel('Firing Density');
subplot 312
plot(dt:dt:T,instISI,'k');
ylabel('|ISI|');
subplot 313
plot(dt:dt:T, Synch,'k');
ylabel('Synch');