clear all
rng(5)
T=500; dt=.5; T=T/dt;
n=2.2; %lin/nonlin. 1=lin, >1=n for nonlin
maxopt=0; %0 for nonlin
k=.01; xn=2;
Ne=400*xn; Npv=50*xn; Nvip=25*xn; Nsom=25*xn;
Ntot=Ne+Npv+Nvip+Nsom; % # of cell groups
taus=[20 10 10 10]; %time constants, e i
conprobs=[.15, .3, .3, .3, 1;
.3, .3, 0.3, .3, 1;
.3, .3, 0.01, .3, 0;
.3, 0.01, .3, .01, 0];
W4=[0.0171 -0.9559 -0.0450 -0.5121;
0.8535 -0.9900 -0.0900 -0.3073;
1.9800 -0.1843 0 -0.7339;
1.2848 0 -0.1399 0]; %2.2379 1.3276 (L=1.2076)
NCmat=zeros(size(W4));
NCmat(:,1)=Ne; NCmat(:,2)=Npv; NCmat(:,3)=Nvip; NCmat(:,4)=Nsom;
gsyns=W4./(NCmat.*conprobs(:,1:4));
NeI=[1:Ne]; NpvI=[Ne+1:Ne+Npv]; NvipI=[Ne+Npv+1:Ne+Npv+Nvip]; NsomI=[Ne+Nvip+Npv+1:Ne+Npv+Nvip+Nsom];
NIs{1}=NeI; NIs{2}=NpvI; NIs{3}=NvipI; NIs{4}=NsomI; NIs{5}=Ntot+1;
for ni=1:5 %over networks
ni
W=zeros(Ntot, Ntot);
dvwn=9; %variance of weight equals 1/dvwn of mean
for prei=1:4 %make weight matrix
for posti=1:4
if gsyns(posti,prei)~=0
MU = log(abs(gsyns(posti,prei))^2 / sqrt((abs(gsyns(posti,prei))/dvwn)+abs(gsyns(posti,prei))^2));
SIGMA = sqrt(log((abs(gsyns(posti,prei))/dvwn)/(abs(gsyns(posti,prei)))^2 + 1));
rweights=lognrnd(MU, SIGMA, length(NIs{posti}), length(NIs{prei})).*sign(gsyns(posti,prei));
if prei==1
rweights(rweights<0)=gsyns(posti,prei)/10; %clipping
else
rweights(rweights>0)=gsyns(posti,prei)/10; %clipping
end
W(NIs{posti}, NIs{prei}) = (rand(length(NIs{posti}), length(NIs{prei}))<conprobs(posti, prei)).* rweights;
end
end
end
SpontRates=[12.8 24 8 8.9];
InpRates=[93 74 0 0];
BinpE=ones(1,T).*SpontRates(1); BinpPV=ones(1, T).*SpontRates(2); BinpVIP=ones(1, T).*SpontRates(3); BinpSOM=ones(1, T).*SpontRates(4); %no reason for this to be over time....
endT=700; stT=500;
InpEk=(rand(Ne,1)>conprobs(1,end)); InpPVk=(rand(Npv,1)>conprobs(2,end)); InpVIPk=(rand(Nvip,1)>conprobs(3,end)); InpSOMk=(rand(Nsom,1)>conprobs(4,end)); %which cells dont get inp
%run it
DV=.5;
Drates=[0 1 1 1]*5;
vn=sqrt((Drates./DV)./.083);
DepolsE=randn(Ne,1).*sqrt(Drates(1)/DV)+Drates(1); DepolsPV=randn(Npv,1).*sqrt(Drates(2)/DV)+Drates(2); DepolsVIP=randn(Nvip,1).*sqrt(Drates(3)/DV)+Drates(3); DepolsSOM=randn(Nsom,1).*sqrt(Drates(4)/DV)+Drates(4); DepolsE(DepolsE<0)=0; DepolsPV(DepolsPV<0)=0; DepolsVIP(DepolsVIP<0)=0; DepolsSOM(DepolsSOM<0)=0;
for ti=1:10 %trials
for aii=1:2 %ach or not
InpE=zeros(Ne,T); InpPV=zeros(Npv,T); InpVIP=zeros(Nvip,T); InpSOM=zeros(Nsom,T);
InpE(:,stT:endT)=InpE(:,stT:endT)+InpRates(1)+repmat(randn(Ne,1).*sqrt(2),1,endT-stT+1); %+repmat(randn(Ne,1).*InpRates(1)/40,1,endT-500+1);
InpPV(:,stT:endT)=InpPV(:,stT:endT)+InpRates(2)+repmat(randn(Npv,1).*sqrt(2),1,endT-stT+1); %+repmat(randn(Npv,1).*InpRates(2)/40,1,endT-500+1);
InpVIP(:,stT:endT)=InpVIP(:,stT:endT)+InpRates(3); %+repmat(randn(Nvip,1).*sqrt(InpRates(3)/40),1,endT-500+1);
InpSOM(:,stT:endT)=InpSOM(:,stT:endT)+InpRates(4); %+repmat(randn(Nsom,1).*sqrt(InpRates(4)/40),1,endT-500+1);
InpE(InpEk,:)=0; InpPV(InpPVk,:)=0; InpVIP(InpVIPk,:)=0; InpSOM(InpSOMk,:)=0;
Re(:,1)=zeros(Ne,1); Rpv(:,1)=zeros(Npv,1); Rvip(:,1)=zeros(Nvip,1); Rsom(:,1)=zeros(Nsom,1);
for t=2:T %over time
dRe=(-Re(:,t-1)+k.*(max((InpE(:,t-1)+BinpE(1,t-1)+ DepolsE*(aii==1)+ W(NeI,:)*[Re(:,t-1); Rpv(:,t-1); Rvip(:,t-1); Rsom(:,t-1)]),maxopt)).^n)./taus(1);
dRpv=(-Rpv(:,t-1)+k.*(max((InpPV(:,t-1)+BinpPV(1,t-1)+ DepolsPV*(aii==1)+W(NpvI,:)*[Re(:,t-1); Rpv(:,t-1); Rvip(:,t-1); Rsom(:,t-1)]),maxopt)).^n)./taus(2);
dRvip=(-Rvip(:,t-1)+k.*(max((InpVIP(:,t-1)+BinpVIP(1,t-1)+ DepolsVIP*(aii==1)+ W(NvipI,:)*[Re(:,t-1); Rpv(:,t-1); Rvip(:,t-1); Rsom(:,t-1)]),maxopt)).^n)./taus(3);
dRsom=(-Rsom(:,t-1)+k.*(max((InpSOM(:,t-1)+BinpSOM(1,t-1)+ DepolsSOM*(aii==1)+ W(NsomI,:)*[Re(:,t-1); Rpv(:,t-1); Rvip(:,t-1); Rsom(:,t-1)]),maxopt)).^n)./taus(4);
Re(:,t)= Re(:,t-1)+dRe*dt;
Rpv(:,t)= Rpv(:,t-1)+dRpv*dt;
Rvip(:,t)= Rvip(:,t-1)+dRvip*dt;
Rsom(:,t)= Rsom(:,t-1)+dRsom*dt;
end
allOTe(aii,:,:)=[Re]; %; Rpv; Rvip; Rsom];
allOTpv(aii,:,:)=[Rpv];
allOTvip(aii,:,:)=[Rvip];
allOTsom(aii,:,:)=[Rsom];
end
%peak:
OTeP(ti,:,:)=squeeze(max(allOTe(:,:,endT-stT:endT),[],3));
OTpvP(ti,:,:)=squeeze(max(allOTpv(:,:,endT-stT:endT),[],3));
OTvipP(ti,:,:)=squeeze(max(allOTvip(:,:,endT-stT:endT),[],3));
OTsomP(ti,:,:)=squeeze(max(allOTsom(:,:,endT-stT:endT),[],3));
%mean transient:
OTeM(ti,:,:)=squeeze(mean(allOTe(:,:,stT:stT+100),3));
OTpvM(ti,:,:)=squeeze(mean(allOTpv(:,:,stT:stT+100),3));
OTvipM(ti,:,:)=squeeze(mean(allOTvip(:,:,stT:stT+100),3));
OTsomM(ti,:,:)=squeeze(mean(allOTsom(:,:,stT:stT+100),3));
%steady state:
OTeT(ti,:,:)=squeeze(mean(allOTe(:,:,endT-10:endT),3));
OTpvT(ti,:,:)=squeeze(mean(allOTpv(:,:,endT-10:endT),3));
OTvipT(ti,:,:)=squeeze(mean(allOTvip(:,:,endT-10:endT),3));
OTsomT(ti,:,:)=squeeze(mean(allOTsom(:,:,endT-10:endT),3));
%baseline:
OTeTs(ti,:,:)=squeeze(mean(allOTe(:,:,200:499),3));
OTpvTs(ti,:,:)=squeeze(mean(allOTpv(:,:,200:499),3));
OTvipTs(ti,:,:)=squeeze(mean(allOTvip(:,:,200:499),3));
OTsomTs(ti,:,:)=squeeze(mean(allOTsom(:,:,200:499),3));
end
mE=squeeze(mean(OTeT,1)); mP=squeeze(mean(OTpvT,1)); mV=squeeze(mean(OTvipT,1)); mS=squeeze(mean(OTsomT,1));
mEs=squeeze(mean(OTeTs,1)); mPs=squeeze(mean(OTpvTs,1)); mVs=squeeze(mean(OTvipTs,1)); mSs=squeeze(mean(OTsomTs,1));
aveEv(ni,:)=[mean(mean(OTeT(:,2,:),3),1), mean(mean(OTpvT(:,2,:),3),1), mean(mean(OTvipT(:,2,:),3),1), mean(mean(OTsomT(:,2,:),3),1)];
aveSp(ni,:)=[mean(mean(OTeTs(:,2,:),3),1), mean(mean(OTpvTs(:,2,:),3),1), mean(mean(OTvipTs(:,2,:),3),1), mean(mean(OTsomTs(:,2,:),3),1)];
Be=[mE(1,:)-mE(2,:) ];
Bes=[mEs(1,:)-mEs(2,:) ];
Bp=[mP(1,:)-mP(2,:) ];
Bps=[mPs(1,:)-mPs(2,:) ];
Bv=[mV(1,:)-mV(2,:) ];
Bvs=[mVs(1,:)-mVs(2,:) ];
Bs=[mS(1,:)-mS(2,:) ];
Bss=[mSs(1,:)-mSs(2,:) ];
avedel(ni,:)=[sum(Be)/length(Be), sum(Bp)/length(Bp), sum(Bv)/length(Bv), sum(Bs)/length(Bs)];
avedels(ni,:)=[sum(Bes)/length(Bes), sum(Bps)/length(Bps), sum(Bvs)/length(Bvs), sum(Bss)/length(Bss)];
aveperc(ni,:)=[sum(Be<0)/length(Be), sum(Bp>0)/length(Bp), sum(Bv<0)/length(Bv), sum(Bs>0)/length(Bs)];
avepercs(ni,:)=[sum(Bes>0)/length(Bes), sum(Bps>0)/length(Bps), sum(Bvs>0)/length(Bvs), sum(Bss>0)/length(Bss)];
end
f3=figure
subplot(2,3,1)
hold on
bar(mean(aveSp,1))
errorbar(mean(aveSp,1),std(aveSp,[],1),'.k')
SpontRate(1,:)=mean(aveSp,1); SpontRate(2,:)=std(aveSp,[],1);
ax=gca;
set(ax,'XTick',[1 2 3 4])
set(ax,'XTickLabel',{'E','P','V','S'})
title('Spontaneous Firing')
subplot(2,3,4)
hold on
bar(mean(aveEv,1))
EvRatesSS(1,:)=mean(aveEv,1); EvRatesSS(2,:)=std(aveEv,[],1);
errorbar(mean(aveEv,1),std(aveEv,[],1),'.k')
ax=gca;
set(ax,'XTick',[1 2 3 4])
set(ax,'XTickLabel',{'E','P','V','S'})
title('Evoked Firing')
subplot(2,3,2)
hold on
bar(mean(avedels,1))
spDifs(1,:)=mean(avedels,1); spDifs(2,:)=std(avedels,[],1);
errorbar(mean(avedels,1),std(avedels,[],1),'.k')
ax=gca;
set(ax,'XTick',[1 2 3 4])
set(ax,'XTickLabel',{'E','P','V','S'})
title('Ach Rate - No Ach')
subplot(2,3,3)
hold on
bar(mean(avepercs,1))
spPerc(1,:)=mean(avepercs,1) ; spPerc(2,:)= std(avepercs,[],1);
errorbar(mean(avepercs,1),std(avepercs,[],1),'.k')
%scatter([1 2 3 4],[.17,.23, .28, .1]+.5, 'r')
plot([1 4], [.5 .5])
ax=gca;
set(ax,'XTick',[1 2 3 4])
set(ax,'XTickLabel',{'E','P','V','S'})
title('% cells increased')
subplot(2,3,5)
hold on
bar(mean(avedel,1))
DifsSS(1,:)=mean(avedel,1); DifsSS(2,:)=std(avedel,[],1);
errorbar(mean(avedel,1),std(avedel,[],1),'.k')
ax=gca;
set(ax,'XTick',[1 2 3 4])
set(ax,'XTickLabel',{'E','P','V','S'})
title('Ach Rate - No Ach')
subplot(2,3,6)
hold on
bar(mean(aveperc,1))
PercSS(1,:)=mean(aveperc,1); PercSS(2,:)=std(aveperc,[],1);
errorbar(mean(aveperc,1),std(aveperc,[],1),'.k')
scatter([1 2 3 4],[.17,.11, .28, .1]+.5, 'r')
plot([1 4], [.5 .5])
ylim([0 1])
ax=gca;
set(ax,'XTick',[1 2 3 4])
set(ax,'XTickLabel',{'E','P','V','S'})
title('% cells changed in the right direction')