run_pDpsim.m
tit=1;
twoorthree=input('network architecture: one inhibitory population (1) or two inhibitory populations, one feedback and one feedworward (2) ');
randortuned=input('network type: random (1), scaled (2), tuned E+I (3) or tuned I (4) ');
analys=input('analysis: yes (1) or no (2) ');
if twoorthree==1
pagexls=1;
if randortuned==1
ext='';
elseif randortuned==2
ext='';
ext1='_tI';
fact=1.4;
elseif randortuned==3
ext='_tEI';
elseif randortuned==4
ext='_tI';
end
ods=input('odor set: learned (1), novel (2), subspace fig.4 (3) or subspace fig. 6 (4) ');
if ods==1
num_odors=10;
load('odorset_10learned.mat')
elseif ods==2
num_odors=10;
load('odorset_10novel.mat')
elseif ods==3
load('odorspace_fig4.mat')
elseif ods==4
load('odorspace_fig6.mat')
else
disp('error')
end
PIN=0;
else
pagexls=2;
load('odorset_20odors.mat');num_odors=20;ods=1;
PIN=input('inhibiting I neurons (vPIN): no (0), FFI (1) or FBI (7)');
if randortuned==1
ext='_r';
elseif randortuned==2
ext='_r';
ext1='_s';
fact=input('scaling factor I to E strength?');
elseif randortuned==4
ext='_s';
elseif randortuned==3
disp('network not available')
end
end
[param,connec,tot]=xlsread('parameters.xlsx',pagexls);
if twoorthree==2 && randortuned<3
tot_mat=10;
else
tot_mat=length(param);
end
r_ob=int8(r_olfbs);clear r_olfbs
for mat=1:tot_mat %can be changed if simulation of only subsets of networks
dt=0.1;
Ninh=1000;
if twoorthree==2
Istim=ones(length(r_ob),Ninh);
PA1=0;
if PIN>0
[inh_neur]=xlsread('parameters.xlsx',3); %neurons identity that are silenced can be changed: these are the ones used in the paper
for oo=1:num_odors
times=[500/dt+3000/dt*(oo-1);3000/dt+3000/dt*(oo-1)];
Istim(times(1):times(2),[inh_neur(oo,PIN):inh_neur(oo,PIN+1),inh_neur(mat,PIN+2):inh_neur(mat,PIN+3)])=PA1;
clear times
end
end
end
struct=1;
if randortuned==2
load(strcat(connec{mat},ext1))
if twoorthree==1
load(strcat(connec{mat},ext),'w_ie');
else
conn=connec{mat};
load(strcat(conn(1:9),num2str(1),ext),'w_ie');
end
w_ie=fact*w_ie;
else
load(strcat(connec{mat},ext))
end
Nexc=4000;
Ninh=length(w_ii);
Nob=1500;
y=param(mat,:);
%neuronal parameters
tau_exc=30;
tau_inh=10;
E_e=0;
E_i=-70;
g_rest_e=1.35;
g_rest_i=0.9;
tau_e=85;
tau_i=50;
tau_ref_e=8/dt;
tau_ref_i=8/dt;
tau_w=40;
ada=1;
E_th=-38;
E_th_i=-45;
if twoorthree==1
E_rest=-60; %
E_rest_ii=-65;
beka=10;
elseif twoorthree==2
E_rest=-60;
E_rest_ii=-60;
beka=5;
end
%activity neuron
if ods>2
ds=0.01;
else
ds=0.1;
end
downs=floor(length(r_ob)*ds);
perc=1;
NumNe=[1:Nexc/perc];
NumNi=[1:Ninh/perc];
V_id=E_rest_ii.*ones(downs,Ninh/perc);
V_ed=E_rest.*ones(downs,Nexc/perc);
g_oed=zeros(downs,Nexc/perc);
g_eed=zeros(downs,Nexc/perc);
g_ied=zeros(downs,Nexc/perc);
g_ief=zeros(downs,Nexc/perc);
g_ieb=zeros(downs,Nexc/perc);
%% solve
t=1;
ts=1;
tsc=1;
tsci=1;
tsci_b=1;
VI=E_rest_ii.*ones(1,Ninh);
VE=E_rest.*ones(1,Nexc);
GEE=zeros(1,Nexc);
GIE=zeros(1,Nexc);
GIEB=zeros(1,Nexc);
GIEF=zeros(1,Nexc);
GEI=zeros(1,Ninh);
GII=zeros(1,Ninh);
GOE=zeros(1,Nexc);
GOI=zeros(1,Ninh);
refrac_E=single(zeros(1,Nexc));
refrac_I=single(zeros(1,Ninh));
WE=zeros(1,Nexc);
spikecount_E=zeros(1000000,2);
spikecount_I=zeros(3000000,2);
while t<=length(r_ob)
V_e=VE;
V_i=VI;
g_ee=GEE;
g_ibe=GIEB;
g_ife=GIEF;
if twoorthree==1
g_ie=GIE;
else
g_ie=GIEB+GIEF;
end
g_ei=GEI;
g_ii=GII;
g_oe=GOE;
g_oi=GOI;
w_e=WE;
VE=V_e+dt/tau_e.*((E_rest-V_e)+((g_oe+g_ee).*(E_e-V_e)+g_ie.*(E_i-V_e)-w_e)/g_rest_e);
VI=V_i+dt/tau_i.*((E_rest_ii-V_i)+((g_oi+g_ei).*(E_e-V_i)+g_ii.*(E_i-V_i))/g_rest_i);
WE=w_e+dt/tau_w*(ada*(V_e-E_rest)-w_e);
refre=find(refrac_E(1,:)>0);
refri=find(refrac_I(1,:)>0);
VE(refre)=E_rest;
VI(refri)=E_rest_ii;
neufire=find(VE>E_th);
neufiri=find(VI>E_th_i);
spikecount_E(tsc:tsc+length(neufire)-1,1)=t;
spikecount_E(tsc:tsc+length(neufire)-1,2)=neufire;
tsc=tsc+length(neufire);
spikecount_I(tsci:tsci+length(neufiri)-1,1)=t;
spikecount_I(tsci:tsci+length(neufiri)-1,2)=neufiri;
tsci=tsci+length(neufiri);
obfire=find(r_ob(t,:)>0);
refrac_E=refrac_E-1;
refrac_I=refrac_I-1;
refrac_E(neufire)=tau_ref_e;
refrac_I(neufiri)=tau_ref_i;
T=2;
WE(1,neufire)=WE(1,neufire)+beka;
GEE=g_ee.*exp(-dt/tau_exc)+y(1).*sum(w_ee(neufire,:),1);
GEI=g_ei.*exp(-dt/tau_exc)+y(4).*sum(w_ei(neufire,:),1);
if twoorthree==1
GIE=g_ie.*exp(-dt/tau_inh)+y(2).*sum(w_ie(neufiri,:),1);
GII=g_ii.*exp(-dt/tau_inh)+y(6).*sum(w_ii(neufiri,:),1);
else
GIEB=g_ibe.*exp(-dt/tau_inh);
GIEF=g_ife.*exp(-dt/tau_inh);
for in=1:length(neufiri)
if neufiri(in)<Ninh/2+1
GIEB=GIEB+y(7).*squeeze(Istim(t,neufiri(in)))*w_ie(neufiri(in),:);
else
GIEF=GIEF+y(2).*squeeze(Istim(t,neufiri(in)))*w_ie(neufiri(in),:);
end
end
GIE=GIEB+GIEF;
GII=g_ii.*exp(-dt/tau_inh)+y(6).*sum(squeeze(Istim(t,neufiri))*w_ii(neufiri,:),1);
end
GOE=g_oe.*exp(-dt/tau_exc)+y(3).*sum(w_oe(obfire,:),1);
GOI=g_oi.*exp(-dt/tau_exc)+y(5).*sum(w_oi(obfire,:),1);
if mod(t-1,1/ds)==0 || t==1
V_id(ts,:)=VI(1,NumNi);
V_ed(ts,:)=VE(1,NumNe);
g_oed(ts,:)=GOE(1,NumNe);
g_eed(ts,:)=GEE(1,NumNe);
g_ied(ts,:)=GIE(1,NumNe);
if twoorthree==2
g_ieb(ts,:)=GIEB(1,NumNe);
g_ief(ts,:)=GIEF(1,NumNe);
end
ts=ts+1;
end
t=t+1;
clear refre refri neufire neufiri obfire
end
clear refrac_E refrac_I
spikecount_E=spikecount_E(1:tsc-1,:);
spikecount_I=spikecount_I(1:tsci-1,:);
g_ed=g_oed+g_eed;
%% ANALYSIS
ACT{tit}=spikecount_E;
ACT_I{tit}=spikecount_I;
if analys==1
[obs loss] = observables(num_odors,Nexc,Ninh,dt,spikecount_E,spikecount_I,g_oed,g_ed,g_ied);
obs_mat(mat,:)=obs
loss_mat(mat,:)=loss;
[res_cot] = cotuning(num_odors,Nexc,dt,g_oed,g_ed,g_ied,EnsE);
cot_mat(mat,:)=res_cot;
end
%figure,plot(spikecount_E(:,1),spikecount_E(:,2),'.k')
tit=tit+1
clearvars -except mat twoorthree tit randortuned analys pagexls ext ext1 fact ods num_odors r_ob tot_mat dt PIN ACT ACT_I param connec tot obs_mat cot_mat loss_mat
%save()
end
