// $Id: network.hoc,v 1.352 2012/03/08 17:56:22 samn Exp $ //* Numbers and connectivity params DOPE_INTF6=1 EDOPE_INTF6=1 ESTDP_INTF6 = ISTDP_INTF6 = 0 declare("scale",2) // Sam had this on 4. declare("colW",1,"colH",1) declare("numcols",colW*colH) declare("dbgcols",0) // whether to debug columns by making them have the same wiring and inputs declare("colr",2) // maximal trans-column projection distance; 0 within col; 1 next col etc declare("colnq","o[5]","lcol",new List()) declare("autotune",0) // declare("useSTDP",1) // {sprint(tstr,"o[%d]",numcols) declare("col",tstr)} {sprint(tstr,"o[%d][%d]",colH,colW) declare("gcol",tstr)} // 2D column grid double div[numcols][numcols][CTYPi][CTYPi]//div[i][j]==# of outputs from type i->j double wmat[numcols][numcols][CTYPi][CTYPi][STYPi] // wmat[i][j][k]==weight from type i->j for synapse k double delm[numcols][numcols][CTYPi][CTYPi]//avg. delay from type i->j double deld[numcols][numcols][CTYPi][CTYPi]//delay variance from type i->j double conv[numcols][numcols][CTYPi][CTYPi] dosetpmat=name_declared("pmat")==0 {sprint(tstr,"d[%d][%d][%d][%d]",numcols,numcols,CTYPi,CTYPi) declare("pmat",tstr)} double synloc[CTYPi][CTYPi]//location of synapses declare("cedp",new List()) // list for DP cells //* swire & related column variables declare("colside",30) // column diameter in micrometers declare("swire",0) // whether to use 'spatial' wiring, 2 means use swirecut, 3 means use swirecutfl (Sam had set to 3) declare("checkers",0) // whether to arrange cells in checkerboard pattern declare("cxinc",3,"cyinc",3,"crad",1) declare("EEsq",(colside/2)^2,"EIsq",(colside/2)^2,"IEsq",(colside/1)^2,"IIsq",(colside/2)^2) // params for swirecut double dels[CTYPi][CTYPi]//[STYPi] //stdev of delays double delv[CTYPi][CTYPi]//variance of delays double vcond[CTYPi][CTYPi]//[STYPi] //conduction velocities declare("layerzvar",25) // range of z location of cells in a layer (in micrometers) declare("colxydist",50) // distance between adjacent columns (x,y) declare("vlayerz",new Vector(CTYPi)) declare("pseed",4321) declare("dvcond","d[2]") //intra layer conduction delays declare("maxsfall",0.001) // max fall-off in prob of connection @ opposite side of column (used by swire) declare("slambda",colside/3) // spatial length constant for probability of connections, used in swirecut dvcond[0] = 1 //if presynaptic cell is excit, conduction delay == 1.0 m/s dvcond[1] = 0.4 //if presynaptic cell is inhib, conduction delay == 0.4 m/s declare("lambda",30,"hval",0) //* other variables declare("mknetnqss",1) // whether COLUMN should make connsnq,cellsnq declare("fc",0.55) declare("EEGain",fc*4*15/11,"EIGain",fc*15,"IEGain",fc*4*15/11,"IIGain",fc*4*15/11) declare("NMAMR",0.1,"EENMGain",1) declare("EIGainInterC0",0.125,"EIGainInterC1",1,"EEGainInterC0",1,"EEGainInterC1",0.25) double EIGainInterC[2],EEGainInterC[2] EIGainInterC[0] = EIGainInterC0 // 0->1 EIGainInterC[1] = EIGainInterC1 EEGainInterC[0] = EEGainInterC0 EEGainInterC[1] = EEGainInterC1 declare("TCNM",0) // whether TC->NM2 on declare("NMUSCLES",4) // number of muscle groups in model declare("sepDPpos",1) // whether to position DP cells in each group separately if(autotune) { seadsetting_INTF6 = 2 wsetting_INTF6 = 1 // use sywv during sim } else if(useSTDP) { seadsetting_INTF6 = 3 // use plasticity } else { seadsetting_INTF6 = 2 // fixed weights } declare("wnqstr","") // intracol wiring NQS (single column) //* prdiv() - print div proc prdiv () { local ii,jj,kk,ll for ii=0,numcols-1 for jj=0,numcols-1 for kk=0,CTYPi-1 for ll=0,CTYPi-1 if(div[ii][jj][kk][ll]) { printf("div[%d][%d][%s][%s]=%g\n",ii,jj,CTYP.o(kk).s,CTYP.o(ll).s,div[ii][jj][kk][ll]) } } // %con (con/pre) = %div (div/post) DEAD_DIV_INTF6=0 declare("jcn",1) declare("disinhib",0) //iff==1 , turn off inhibition, by setting wmat[I%d][...]==0 in inhiboff() declare("pmatscale",1) // scale for pmat - allows keeping it fixed while changing # of cells in network declare("wmatscale",1) // scale for wmat - called after setwmat batch_flag=declare("dstr",datestr,"setdviPT",NORM) declare("params","not batch","ofile",output_file) declare("dvseed",534023) // seed for wiring dosetcpercol=name_declared("cpercol")==0 // whether to set values in cpercol or use user-supplied values {sprint(tstr,"d[%d][%d]",numcols,CTYPi) declare("cpercol",tstr)} // cells of a specific type per column {sprint(tstr,"o[%d]",numcols) declare("vcpercol",tstr)} for i=0,numcols-1 vcpercol[i]=new Vector(CTYPi) declare("DPESFeedForward",0.15) // level of feedforward connectivity from DP -> ES // declare("DPESFeedForward",0.1) // level of feedforward connectivity from DP -> ES declare("EMRecur",0.0) //declare("EMRecur",0.07) // level of recurrence between EM cells declare("ESRecur",0.0)//declare("ESRecur",0.191)// level of recurrence between ES cells declare("EMESFeedBack",0.0)//declare("EMESFeedBack",0.017)// level of feedback from EM -> ES declare("ESEMFeedForward",0.08)// level of feedforward connectivity from ES -> EM sprint(tstr,"d[%d][%d]",CTYPi,CTYPi) if(!(i=name_declared("delmscale"))) { declare("delmscale",tstr) // scale delm values by this # } if(!i) { for i=0,CTYPi-1 for j=0,CTYPi-1 delmscale[i][j]=1 } //** mklayerz -- get average z location of a layer, based on frana and tononi // returns value in micrometers. 0 is at top (layer 1), max val // is at layer 6 ... actual #s don't matter, only distances between // the layers matters... proc mkvlayerz () { local i,vl vlayerz.indgen(0,CTYPi-1,1) vlayerz.apply("layer") for vtr(&vl,vlayerz,&i) { if(vl == 2 || vl == 3){ vlayerz.x(i)=1740 // average of 1540+1940 from frana } else if(vl == 4) { vlayerz.x(i) = (1740+1130) / 2 } else if(vl == 5 ){ vlayerz.x(i) = 1130 } else if(vl == 6 ){ vlayerz.x(i) = 488 } else vlayerz.x(i) = 200 } } //* setcpercol - set # of cells per column proc setcpercol () { local i,j // (/u/samn/vcsim/notebook.dol_1:24562)(notebook.dol_1:24492) if(dosetcpercol) { // if user didn't supply values (default), set # of cells of a type per column cpercol[0][ES] = int(48 * scale) // sensory layer cpercol[0][ISL] = 5 * scale cpercol[0][IS] = 11 * scale cpercol[0][EM] = int(48 * scale) // motor layer cpercol[0][IML] = 5 * scale cpercol[0][IM] = 11 * scale cpercol[0][DP] = cpercol[0][ES] // arm-cell/spindle layer for i=0,CTYPi-1 if(cpercol[0][i]) cpercol[0][i]=int(cpercol[0][i]+0.5) // round //store values in vcpercol Vectors for i=0,numcols-1 { {vcpercol[i].resize(CTYPi) vcpercol[i].fill(0)} // store the values in a vector for j=0,CTYPi-1 vcpercol[i].x(j)=cpercol[i][j] } } } //* setpmat() proc setpmat () { local pre,po,ii,jj,kk,ll if(!dosetpmat) return // if pmat setup by user (in notebook), then don't reset its values for ii=0,numcols-1 for jj=0,numcols-1 for kk=0,CTYPi-1 for ll=0,CTYPi-1 pmat[ii][jj][kk][ll]=0 // intracol connections ii=0 pmat[0][0][DP][ES]=DPESFeedForward // 0.1, DP = drawspike cells (muscle/arm layer) // sensory layer pmat[ii][ii][ES][EM]=ESEMFeedForward pmat[ii][ii][ES][ES]=ESRecur pmat[ii][ii][ES][ISL]=0.51 pmat[ii][ii][ES][IS]=0.43 pmat[ii][ii][ISL][ES]=0.35 pmat[ii][ii][ISL][ISL]=0.09 pmat[ii][ii][ISL][IS]=0.53 pmat[ii][ii][IS][ES]=0.44 pmat[ii][ii][IS][ISL]=0.34 pmat[ii][ii][IS][IS]=0.62 // motor layer pmat[ii][ii][EM][EM]=EMRecur pmat[ii][ii][EM][ES]=EMESFeedBack pmat[ii][ii][EM][IML]=0.51 pmat[ii][ii][EM][IM]=0.43 pmat[ii][ii][IML][EM]=0.35 pmat[ii][ii][IML][IML]=0.09 pmat[ii][ii][IML][IM]=0.53 pmat[ii][ii][IM][EM]=0.44 pmat[ii][ii][IM][IML]=0.34 pmat[ii][ii][IM][IM]=0.62 } //* scalepmat(fctr) - multiply values in pmat by fctr proc scalepmat () { local fctr,from,to,c1,c2 fctr=$1 for c1=0,numcols-1 for c2=0,numcols-1 { for from=0,CTYPi-1 for to=0,CTYPi-1 pmat[c1][c2][from][to] *= fctr } } //* pmat2nq - return an NQS with info in pmat obfunc pmat2nq () { local i,j,k,l localobj nqpmat nqpmat=new NQS("col1","col2","froms","tos","from","to","pij") {nqpmat.strdec("froms") nqpmat.strdec("tos")} for i=0,numcols-1 for j=0,numcols-1 for k=0,CTYPi-1 for l=0,CTYPi-1 if(pmat[i][j][k][l]) { nqpmat.append(i,j,CTYP.o(k).s,CTYP.o(l).s,k,l,pmat[i][j][k][l]) } return nqpmat } //* nq2pmat - load NQS ($o1) into pmat proc nq2pmat () { local i,j,k,l localobj nq,vf,vto,vc1,vc2,vpij {nq=$o1 nq.tog("DB") vc1=nq.getcol("col1") vc2=nq.getcol("col2")} {vf=nq.getcol("from") vto=nq.getcol("to") vpij=nq.getcol("pij")} for i=0,numcols-1 for j=0,numcols-1 for k=0,CTYPi-1 for l=0,CTYPi-1 pmat[i][j][k][l]=0 for i=0,vf.size-1 pmat[vc1.x(i)][vc2.x(i)][vf.x(i)][vto.x(i)]=vpij.x(i) print "loaded " , nq , " into pmat" } //* synapse locations DEND SOMA AXON proc setsynloc () { local from,to for from=0,CTYPi-1 for to=0,CTYPi-1 { if(ice(from)) { if(IsLTS(from)) { synloc[from][to]=DEND // distal [GA2] - from LTS } else { synloc[from][to]=SOMA // proximal [GA] - from FS } } else { synloc[from][to]=DEND // E always distal. use AM2,NM2 } } } //* setdelmats -- setup delm,deld proc setdelmats () { local from,to,ii,jj,fcm,fcd for ii=0,numcols-1 for jj=0,numcols-1 { if(ii==jj) { fcm=fcd=1 } else {fcm=3 fcd=6} for from=0,CTYPi-1 for to=0,CTYPi-1 { if(synloc[from][to]==DEND) { delm[ii][jj][from][to]=4 * delmscale[from][to] * fcm deld[ii][jj][from][to]=1 * fcd } else { delm[ii][jj][from][to]=2.0 * delmscale[from][to] * fcm deld[ii][jj][from][to]=0.2 * fcd } } } } //* weight params //** delay all 2+/-0.02 within column for now proc setwmat () { local from,to,sy,gn,c,ii,jj for ii=0,numcols-1 for jj=0,numcols-1 for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=0,STYPi-1 wmat[ii][jj][from][to][sy]=0 ii=0 // wire DP directly to ES wmat[0][0][DP][ES][AM2]= 1.75 * 5.01 / EEGain // 1.67 // sensory layer wmat[ii][ii][ES][ES][AM2]=0.66 wmat[ii][ii][ES][ISL][AM2]=0.23 wmat[ii][ii][ES][IS][AM2]=0.23 wmat[ii][ii][ES][EM][AM2]=0.88 * 2 wmat[ii][ii][IS][ES][GA]=1.5 wmat[ii][ii][IS][ISL][GA]=1.5 wmat[ii][ii][IS][IS][GA]=1.5 wmat[ii][ii][ISL][ES][GA2]=0.83 wmat[ii][ii][ISL][ISL][GA2]=1.5 wmat[ii][ii][ISL][IS][GA2]=1.5 // motor layer wmat[ii][ii][EM][EM][AM2]=0.71 wmat[ii][ii][EM][IML][AM2]=0.23 wmat[ii][ii][EM][IM][AM2]=0.23 wmat[ii][ii][EM][ES][AM2]=0.24 wmat[ii][ii][IM][EM][GA]=1.5 wmat[ii][ii][IM][IM][GA]=1.5 wmat[ii][ii][IM][IML][GA]=1.5 wmat[ii][ii][IML][EM][GA2]=0.83 wmat[ii][ii][IML][IM][GA2]=1.5 wmat[ii][ii][IML][IML][GA2]=1.5 //set NMDA weights for ii=0,numcols-1 for jj=0,numcols-1 { for from=0,CTYPi-1 for to=0,CTYPi-1 wmat[ii][jj][from][to][NM2]=NMAMR*wmat[ii][jj][from][to][AM2] } //gain control for ii=0,numcols-1 for jj=0,numcols-1 for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=AM,GA2 if(wmat[ii][jj][from][to][sy]>0) { if(ice(from)) { if(ice(to)) { gn = IIGain } else { gn = IEGain } if(IsLTS(from) && !IsLTS(to)) gn *= 0.5 } else { if(ice(to)) { gn = EIGain if(IsLTS(to)) gn *= 0.5 } else { gn = EEGain if(sy==NM || sy==NM2) gn *= EENMGain // E->E NMDA gain } } wmat[ii][jj][from][to][sy] *= gn } } //* scalewmat(fctr) - multiply values in wmat by fctr proc scalewmat () { local fctr,from,to,c1,c2,sy fctr=$1 for c1=0,numcols-1 for c2=0,numcols-1 { for from=0,CTYPi-1 { if(from==DP || IsTHAL(from)) continue for to=0,CTYPi-1 { if(IsTHAL(to)) continue for sy=0,STYPi-1 wmat[c1][c2][from][to][sy] *= fctr } } } } //* %con (con/pre) = %div (div/post) {sp = new NQS() cp = new NQS()} //** gethublims(col,hubtype,hubfactor,numhubs,mode) // get a matrix of size CTYPi X CTYPi, specifying div with mat.x(hubtype,othertype) // and conv with mat.x(othertype,hubtype) // hubtype = type of hub. hubfactor = desired ratio of hub div/conv vs non-hub div/conv // numhubs = # of hubs. col = COLUMN for which to set hubs. // mode == 0 <-- hub div(conv) is set to hubfactor * original div(conv) // mode == 1 <-- hub div(conv) is set so that final hub div = hubfactor * final non_hub div (same for conv) // formula is based on: m / ((N-H*m) / (C-H)) = F , and then solving for m // m = div for the hubs, F = desired ratio of final hub div to final non-hub div // N = # of synapses (links), C = total # of postsynaptic cells (including hubs) , H = # of hubs // similarly done for conv , but replace N with appropriate values // (/u/samn/intfcol/notebook.dol_1:21933) obfunc gethublims () { local ct,mode,from,to,lim,nc,nhubs,fctr localobj col,mat {col=$o1 ct=$2 fctr=$3 nhubs=$4 mode=$5 mat=new Matrix(CTYPi,CTYPi)} for to=0,CTYPi-1 if(col.numc[to] && col.div[ct][to]) { {nc=col.numc[to] if(ct==to)nc-=1} // deduct for self-link if(mode==0) { lim = int( 0.5 + col.div[ct][to]*fctr ) } else { lim = int( 0.5 + col.div[ct][to]*col.numc[ct]*fctr/(col.numc[ct]-nhubs+fctr*nhubs) ) } mat.x(ct,to) = MINxy(lim, nc) // at most div to all postsynaptic cells } for from=0,CTYPi-1 if(col.numc[from] && col.div[from][ct]) { {nc=col.numc[from] if(ct==from)nc-=1} // deduct for self-link if(mode==0) { lim = int( 0.5 + col.conv[from][ct]*fctr ) } else { lim = int( 0.5 + col.div[from][ct]*col.numc[from]*fctr/(col.numc[ct]-nhubs+fctr*nhubs) ) } mat.x(from,ct) = MAXxy(MINxy(lim, nc),1) // at most conv from all presynaptic cells, but at least 1 } return mat } //** addhubs(column,cell-type,numhubs,scaling factor,skipI[,seed,allowz,hubmode,verbose]) // add hubs to the network by stealing wires from other neurons // $o1 == column // $2 == cell type of hub // $3 == number of hubs to add // $4 == scaling factor (should be > 1.0) for conv,div of hub // $5 == skip div/conv of I cells // $6 == seed - optional // $7 == allowz - whether to allow pulling all links from/to another cell // $8 == hubmode - which mode to use for gethublims (see above) // $9 == verbose - optional // function returns a Vector containing the ids of the cells selected as hubs (within column ids) obfunc addhubs () { local a,ct,fctr,nhubs,idx,jdx,lseed,hubid,szorig,cursz,preid,poid,lim,skipI,to,from,vrb,changed,allowz,hmode\ localobj col,ce,vin,vout,nq,vd,vc,vdd,vdt,vddt,vpicked,vhubid,vw1,vw2,vsyn,vprob,vsynt,vtmp,vdsz,vcsz,mhlim col=$o1 ct=$2 nhubs=$3 fctr=$4 skipI=$5 if(numarg()>5) lseed=$6 else lseed=1234 if(numarg()>6) allowz=$7 else allowz=1 if(numarg()>7) hmode=$8 else hmode=0 if(numarg()>8) vrb=$9 else vrb=0 {ce=col.ce hashseed_stats(lseed,lseed,lseed)} a=allocvecs(vin,vout,vd,vc,vdd,vdt,vddt,vpicked,vw1,vw2,vsyn,vprob,vsynt,vtmp,vdsz,vcsz) vrsz(col.allcells,vin,vout,vd,vc,vdd,vdt,vddt,vpicked,vw1,vw2,vsyn,vprob,vsynt,vdsz,vcsz,vtmp) mhlim=gethublims(col,ct,fctr,nhubs,hmode) //vin,vout = input/output markers. vd,vc = div/conv. //vdd div/conv delays, vdt div/conv temp. vddt=div/conv delay temp //vpicked=which cells already picked as hubs vhubid=new Vector() {vhubid.indgen(col.ix[ct],col.ixe[ct],1) vhubid.shuffle() vhubid.resize(nhubs)} if(vrb) vlk(vhubid) for idx=0,vhubid.size-1 vpicked.x(vhubid.x(idx))=1 for idx=0,vhubid.size-1 { hubid=vhubid.x(idx) if(vrb) printf("hub%d id = %d\n",idx+1,hubid) {ce.o(hubid).getdvi(vd,vdd,vw1,vw2,vprob,vsyn) ce.o(hubid).getconv(vc)}//IDs of post/presynaptic cells {ce.o(hubid).getconv(1.2,vcsz) vdsz.resize(CTYPi) vdsz.fill(0)}//counts of post/pre types for jdx=0,vd.size-1 vdsz.x(ce.o(vd.x(jdx)).type)+=1 {vout.fill(0) vin.fill(0)} //init as 0 for jdx=0,vd.size-1 vout.x(vd.x(jdx))=1 //mark current postsynaptic cells for jdx=0,vc.size-1 vin.x(vc.x(jdx))=1 //mark current presynaptic cells for to=0,CTYPi-1 if(col.numc[to] && col.div[ct][to] && (!skipI || !ice(to))) { cursz=szorig=vdsz.x(to) // update divergence if(vrb) print "\torig div -> " , CTYP.o(to).s, " = " , szorig {lim=mhlim.x(ct,to) changed=1} while(cursz " , CTYP.o(to).s, " = " , cursz } ce.o(hubid).setdvi(vd,vdd,vsyn,1) // update hub dvi for from=0,CTYPi-1 if(col.numc[from] && col.div[from][ct] && (!skipI || !ice(from))) { cursz=szorig=vcsz.x(from) // update convergence {lim=mhlim.x(from,ct) changed=1} if(vrb) print "\torig conv <- ", CTYP.o(from).s, " = " , szorig while(cursz1) { // make sure not to remove all inputs of a type to a cell vdt.x( jdx ) = hubid // reassign input to hub ce.o(preid).setdvi(vdt,vddt,vsynt,1) // reset presynaptic cell's div vin.x( preid ) = changed = 1 // mark input cursz += 1 break } } } } } if(vrb) print "\tnew conv <- " , CTYP.o(from).s, " = " , cursz } } ce.o(0).finishdvir() {dealloc(a) return vhubid} } //* mkcolnqs - make an nqs with current pmat,wmat,delm,deld info for use by both COLUMNs(AREAS) for wiring // "dist" represents distance between columns: dist==0 for intra-COLUMN setup, dist>0 for INTER-COLUMN setup proc mkcolnqs () { local from,to,sy,i,j,d localobj froms,tos,sys for i=0,numcols-1 { {nqsdel(colnq[i]) colnq[i]=new NQS("froms","tos","sys","from","to","sy","w","pij","delm","deld","loc","dist")} colnq[i].strdec("froms","tos","sys") d=0 //intracol for from=0,CTYPi-1 { froms=CTYP.o(from) for to=0,CTYPi-1 { tos=CTYP.o(to) if(pmat[i][i][from][to]>0) for sy=0,STYPi-1 if(wmat[i][i][from][to][sy]>0) { sys=STYP.o(sy) colnq[i].append(froms.s,tos.s,sys.s,from,to,sy,wmat[i][i][from][to][sy],pmat[i][i][from][to],delm[i][i][from][to],deld[i][i][from][to],synloc[from][to],d) } } } for j=0,numcols-1 if(i!=j) { d=1 //intercol for from=0,CTYPi-1 { froms=CTYP.o(from) for to=0,CTYPi-1 { tos=CTYP.o(to) if(pmat[i][j][from][to]>0) for sy=0,STYPi-1 if(wmat[i][j][from][to][sy]>0) { sys=STYP.o(sy) colnq[i].append(froms.s,tos.s,sys.s,from,to,sy,wmat[i][j][from][to][sy],pmat[i][j][from][to],delm[i][j][from][to],deld[i][j][from][to],synloc[from][to],d) } } } } } } //** swirecut (col,EEsq,EIsq,IEsq,IIsq[,seed,slambda]) - spatial wiring: wires the column using pmat and cell positions // (wiring probability effected by distance btwn cells) // seed is random # seed // lambda specifies length-constant for spatially-dependent fall-off in wiring probability // at one lambda away, probability of connect is value in pmat proc swirecut () { local x,y,z,ii,jj,a,del,prid,poid,prty,poty,dv,lseed,h,prob,slambda,dsq,dist,slambdasq,\ EEsq,EIsq,IEsq,IIsq,ic1,ic2,pdx\ localobj col,ce,vidx,vdel,vdist,vwt1,vwt2,vtmp,opr,opo,st,rdm,vprob a=allocvecs(vidx,vdel,vtmp,vdist,vwt1,vwt2,vprob) z=0 col=$o1 ce=col.ce EEsq=$2 EIsq=$3 IEsq=$4 IIsq=$5 if (argtype(6)==0) lseed=$6 else lseed=1234 if(numarg()>6) slambda=$7 else slambda=10 if(slambda<=0){ printf("swirecut WARN: invalid slambda=%g, setting slambda to %g\n",colside/3) slambda=colside/3 } slambdasq = slambda^2 // using squared distance vrsz(1e4,vidx,vdel,vdist,vtmp) rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator rdm.uniform(0,1) vprob.resize(ce.count^2) vprob.setrand(rdm) pdx=0 // Create the connectivity NQS table. if(col.mknetnqss) {nqsdel(col.connsnq) col.connsnq=new NQS("id1","id2","del","dist","wt1","wt2")} for y=0,ce.count-1 { opr=ce.o(y) vrsz(0,vidx,vdel,vdist,vwt1,vwt2) prid=opr.id prty=opr.type ic1=ice(opr.type) for poty=0,CTYPi-1 if (col.numc[poty]!=0 && (h=col.pmat[prty][poty])>0) { for poid=col.ix[poty],col.ixe[poty] if(prid!=poid) { // go thru postsynaptic cells opo = ce.o(poid) ic2=ice(ce.o(poid).type) dx = opr.xloc - opo.xloc dy = opr.yloc - opo.yloc dsq = dx^2 + dy^2 if(ic1 && ic2 && dsq > IIsq) { // check max dists continue } else if(!ic1 && !ic2 && dsq > EEsq) { continue } else if(!ic1 && ic2 && dsq > EIsq) { continue } else if(ic1 && !ic2 && dsq > IEsq) continue if(dsq <= slambdasq) { prob=h } else { prob = h * exp(1 - sqrt(dsq)/slambda) // probability of connect // print slambda,dx,dy,dsq,h,prob } if( prob >= vprob.x(pdx) ) { // rdm.uniform(0,1) del = rdm.uniform(col.delm[prty][poty]-col.deld[prty][poty],col.delm[prty][poty]+col.deld[prty][poty]) {vidx.append(poid) vdel.append(del)} if(synloc[prty][poty]==DEND) vdist.append(1) else vdist.append(0) if(col.mknetnqss || wsetting_INTF6==1) { if(col.syty1[prty][poty]>=0) vwt1.append(col.wmat[prty][poty][col.syty1[prty][poty]]) else vwt1.append(0) if(col.syty2[prty][poty]>=0) vwt2.append(col.wmat[prty][poty][col.syty2[prty][poty]]) else vwt2.append(0) } } pdx += 1 } } if(vidx.size>0) { if(wsetting_INTF6==1) opr.setdvi(vidx,vdel,vdist,1,vwt1,vwt2) else opr.setdvi(vidx,vdel,vdist,1) if(col.mknetnqss) { for ii=0,vidx.size-1 col.connsnq.append(prid,vidx.x(ii),vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii)) } } } dealloc(a) if(verbose) printf("\n") } //* conprob func conprob () { // calculate connection probability using global slambda and hval if ($1<=slambda) return hval else return hval*exp(1-$1/slambda) // probability of connect } //** swirecutfl (col,EEsq,EIsq,IEsq,IIsq[,seed,slambda]) - spatial wiring: wires the column using pmat and cell positions // (wiring probability effected by distance btwn cells) // seed is random # seed // slambda specifies length-constant for spatially-dependent fall-off in wiring probability // at one slambda away, probability of connect is value in pmat proc swirecutfl () { local x,y,z,ii,jj,a,del,prid,poid,prty,poty,dv,lseed,h,prob,slambda,dsq,dist,slambdasq,\ EEsq,EIsq,IEsq,IIsq,ic1,ic2,pdx,n,nn,maxd,rsz\ localobj col,ce,vidx,vdel,vdist,vwt1,vwt2,vtmp,vpre,opr,opo,st,rdm,vprob,oq,ol,mmaxd oq=new NQS("id1","id2","del","dist","wt1","wt2","dist0","prob") oq.verbose=0 ol=new List() mmaxd=new Matrix(2,2) a=allocvecs(vidx,vdel,vtmp,vdist,vwt1,vwt2,vprob,vpre) z=0 {ol.append(vpre) ol.append(vidx) ol.append(vdel) ol.append(vdist) ol.append(vwt1) ol.append(vwt2)} col=$o1 ce=col.ce EEsq=$2 EIsq=$3 IEsq=$4 IIsq=$5 if (argtype(6)==0) lseed=$6 else lseed=1234 if(numarg()>6) slambda=$7 else slambda=10 if(slambda<=0){ printf("swirecut WARN: invalid slambda=%g, setting slambda to %g\n",colside/3) slambda=colside/3 } slambdasq = slambda^2 // using squared distance vrsz(1e4,vidx,vdel,vdist,vtmp,vpre) rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator rdm.uniform(0,1) rsz=ce.count^2 vprob.resize(rsz) vprob.setrand(rdm) pdx=0 {mmaxd.x(0,0)=sqrt(EEsq) mmaxd.x(0,1)=sqrt(EIsq) mmaxd.x(1,0)=sqrt(IEsq) mmaxd.x(1,1)=sqrt(IIsq)} // Create the connectivity NQS table. if(col.mknetnqss) {nqsdel(col.connsnq) col.connsnq=new NQS("id1","id2","del","dist","wt1","wt2")} for y=0,ce.count-1 { opr=ce.o(y) vrsz(0,vidx,vdel,vdist,vwt1,vwt2) prid=opr.id prty=opr.type ic1=ice(opr.type) for poty=0,CTYPi-1 if (col.numc[poty]!=0 && (h=hval=col.pmat[prty][poty])>0) { ic2 = ice(poty) oq.clear() if ((n=opr.floc(opr.xloc,opr.yloc,1e9,oq.getcol("id2"),oq.getcol("dist"),mmaxd.x(ic1,ic2),poty))==0) { if(verbose) printf("swirecutfl: No possibilities for connections found from %d\n",opr.id) } else { oq.getcol("dist0").copy(oq.getcol("dist")) oq.pad(n) oq.getcol("dist0").apply("conprob") // turn distances into probabilities if (pdx+n>=rsz) { rdm.uniform(0,1) vprob.setrand(rdm) pdx=0 } oq.getcol("prob").copy(vprob,pdx,pdx+n-1) // pick out some random probabilities to compare to pdx+=n oq.calc(".sub()") if ((nn=oq.select("dist0",">",0,"id2","!=",prid))==0) { // looks for dist_prob > prob if(verbose) printf("swirecutfl: No connections found from %d out of %d assayed\n",opr.id,n) } else { oq.cpout() //oq.fill("id1",prid) rdm.uniform(col.delm[prty][poty]-col.deld[prty][poty],col.delm[prty][poty]+col.deld[prty][poty]) oq.getcol("del").setrand(rdm) if (col.syty1[prty][poty]>=0) oq.getcol("wt1").fill(col.wmat[prty][poty][col.syty1[prty][poty]]) else oq.getcol("wt1").fill(0) if (col.syty2[prty][poty]>=0) oq.getcol("wt2").fill(col.wmat[prty][poty][col.syty2[prty][poty]]) else oq.getcol("wt2").fill(0) {vwt1.append(oq.getcol("wt1")) vwt2.append(oq.getcol("wt2"))} {vidx.append(oq.getcol("id2")) vdel.append(oq.getcol("del"))} vdist.resize(vdist.size+nn) if(synloc[prty][poty]==DEND) vdist.fill(1,vdist.size-nn,vdist.size-1) else vdist.fill(0,vdist.size-nn,vdist.size-1) } } } if(vidx.size>0) { if(wsetting_INTF6==1) opr.setdvi(vidx,vdel,vdist,1,vwt1,vwt2) else opr.setdvi(vidx,vdel,vdist,1) if(col.mknetnqss){vpre.resize(vidx.size) vpre.fill(prid) col.connsnq.append(ol)} // for ii=0,vidx.size-1 col.connsnq.append(prid,vidx.x(ii),vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii)) } } col.ce.o(0).finishdvir() dealloc(a) nqsdel(oq) if(verbose) printf("\n") } //* swireDPs(seed,slambda,) - wire the DP (muscle proprioceptive) cells - they're only in COLUMN[0] (sensory) proc swireDPs () { local x,y,z,ii,jj,a,n,del,prid,poid,prty,poty,dv,lseed,prob,dsq,dist,ic1,ic2,pdx,rsz,cdx\ localobj cl,lce,vidx,vdel,vdend,vwt1,vwt2,vprob,vtmp,opr,opo,st,nc,oq,ol,rdm,vpoty //prid,poid,vdel.x[ii],vdist.x[ii],vwt1.x[ii],vwt2.x[ii] oq=new NQS("id1","id2","del","dist","wt1","wt2","dist0","prob") oq.verbose=0 a=allocvecs(vprob,vpoty) ol=new List() for ii=0,5 ol.append(oq.v[ii]) // will use ol for the appends vidx=oq.getcol("id1") // vdel,vdend,vwt1,vwt2,vtmp,vprob cdx=0 cl=col[cdx] lce=cl.ce if (argtype(1)==0) lseed=$1 else lseed=1234 if(argtype(2)==0) slambda=$2 else slambda=10 if(slambda<=0){slambda=colside/3 printf("swireDPs WARN: invalid lambda=%g, setting lambda to %g\n",slambda)} rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator rdm.uniform(0,1) vprob.resize(rsz=lce.count^2) vprob.setrand(rdm) pdx=0 vpoty.append(ES) //** wire in the DPs prty=DP for vtr(&poty,vpoty) { hval=pmat[cdx][cdx][prty][poty] for ltr(opr,cedp,&ii) { prid=opr.id oq.clear() if ((n=col[cdx].ce.o(0).floc(opr.xloc,opr.yloc,1e9,oq.getcol("id2"),oq.getcol("dist"),4*slambda,poty))==0) { printf("swireDPs: No possibilities for connections found from %d\n",opr.id) } else { oq.getcol("dist0").copy(oq.getcol("dist")) oq.pad(n) oq.getcol("dist0").apply("conprob") // turn distances into probabilities if (pdx+n>=rsz) { vprob.setrand(rdm) pdx=0 } oq.getcol("prob").copy(vprob,pdx,pdx+n-1) // pick out some random probabilities to compare to pdx+=n oq.calc(".sub()") if ((nn=oq.select("dist0",">",0))==0) { // looks for dist_prob > prob printf("swireDPs: No connections found from %d out of %d assayed\n",opr.id,n) } else { oq.cpout() oq.fill("id1",prid) rdm.uniform(cl.delm[prty][poty]-cl.deld[prty][poty],cl.delm[prty][poty]+cl.deld[prty][poty]) oq.getcol("del").setrand(rdm) if (cl.syty1[prty][poty]>=0) oq.getcol("wt1").fill(cl.wmat[prty][poty][cl.syty1[prty][poty]]) else oq.getcol("wt1").fill(0) if (cl.syty2[prty][poty]>=0) oq.getcol("wt2").fill(cl.wmat[prty][poty][cl.syty2[prty][poty]]) else oq.getcol("wt2").fill(0) if (cl.mknetnqss) cl.connsnq.append(ol) for vtr(&poid,oq.out.getcol("id2"),&jj) { ncl.append(nc=new NetCon(opr.drspk,lce.o(poid))) nc.delay=oq.getcol("del").x[jj] nc.weight[cl.syty1[prty][poty]]=oq.getcol("wt1").x[jj] nc.weight[cl.syty2[prty][poty]]=oq.getcol("wt2").x[jj] } } } } } dealloc(a) nqsdel(oq) } //* checkerpos(col,rad,xinc,yinc,seed) - arrange cells in checkerboard pattern proc checkerpos () { local i,j,x,y,csz,xinc,yinc,seed,rad,a localobj col,ce,rdm,xo,nqc col=$o1 csz=col.cside rad=$2 xinc=$3 yinc=$4 rdm=new Random() if(numarg()>4) seed=$5 else seed=1234 rdm.ACG(seed) ce=col.ce nqc=new NQS("x","y") for(y=0;y=csz || y>=csz) { r = rdm.uniform(0,nqc.v.size-1) x = nqc.v[0].x(r) + rdm.uniform(-rad,rad) y = nqc.v[1].x(r) + rdm.uniform(-rad,rad) } {xo.xloc=x xo.yloc=y} if(col.cellsnq!=nil) {col.cellsnq.v[5].x(i)=xo.xloc col.cellsnq.v[6].x(i)=xo.yloc} } nqsdel(nqc) } //* mkcols - make the COLUMNs proc mkcols () { local id,x,y,seed,svnum localobj nq id=nextGID_INTF6=0 for y=0,colH-1 for x=0,colW-1 { if(dbgcols)seed=dvseed else seed=(id+1)*dvseed svnum=vcpercol[id].x[DP] vcpercol[id].x[DP]=0 // ==== REMOVE DPs so can create separately ==== lcol.append(gcol[y][x]=new COLUMN(id,vcpercol[id],colnq[id],seed,x,y,setdviPT,mknetnqss,1)) print "made COLUMN, wiring" // vcpercol[id].x[DP]=svnum // replace col[id]=gcol[y][x] col[id].verbose=verbose_INTF6 nextGID_INTF6 += cpercol[0][DP] // inc by # of DPs if(strlen(wnqstr)>0 && id==1) { // wire from NQS nq=new NQS(wnqstr) col[id].wirenq(nq) nqsdel(nq) } else if(swire>0) { // spatial-wiring col[id].setcellpos(vlayerz,layerzvar,pseed*(id+1),colside) if(checkers) checkerpos(col[id],crad,cxinc,cyinc,pseed*(id+1)) if(swire==1) { col[id].swire(col[id].wseed,maxsfall) } else if(swire==2) { swirecut(col[id],EEsq,EIsq,IEsq,IIsq,col[id].wseed,slambda) } else if(swire==3) { swirecutfl(col[id],EEsq,EIsq,IEsq,IIsq,col[id].wseed,slambda) // uses intf6::floc for speedup } } else col[id].wire(col[id].wseed) // random wiring (no spatial dependence) id+=1 } } //** mkDPs() create and locate the Stimulators -- these are DRSPK cells rather than INTF6 proc mkDPs () { local ii,x,y,c,ndp,suty,rw,cl,rs,cs,pad,cdx localobj cellsnq,xo,rdm,ce cdx=0 ndp=col[cdx].numc[DP]=cpercol[cdx][DP] // will be outside of col now col[cdx].ix[DP]=col[cdx].allcells col[cdx].ixe[DP]=col[cdx].allcells+ndp-1 ce=col[cdx].ce gid=0 rs=2 cs=2 // 2 rows and 2 columns for ii=0,col[cdx].numc[DP]-1 cedp.append(new DPC()) // create the DPs rdm=new Random() rdm.ACG(pseed+DP) pad=colside/cs/8 if (col[cdx].cellsnq==nil) c=-1 else { cellsnq=col[cdx].cellsnq cellsnq.tog("DB") c=cellsnq.fi("xloc") gid=cellsnq.applf(".max","gid")+1 if(gid!=col[cdx].ix[DP]) {printf("mkDPS ERRA %d\n",gid) return } } if (gid==0) gid=col[cdx].ix[DP] //** assign DP identity for ltr(xo,cedp,&ii) { xo.id=gid xo.type=DP suty=xo.subtype=int(xo.id%NMUSCLES) // which muscle it belongs to rw=int(suty/cs) cl=suty%cs // xloc,yloc // print suty,rw,cl,cl*colside/cs,(cl+1)*colside/cs,rw*colside/cs,(rw+1)*colside/rs if(sepDPpos) { // each set of DP cells is spatially separated xo.xloc=x=rdm.uniform(cl*colside/cs+pad,(cl+1)*colside/cs-pad) xo.yloc=y=rdm.uniform(rw*colside/rs+pad,(rw+1)*colside/rs-pad) } else { // dont segregate DP cells xo.xloc=x=rdm.uniform(0,colside) xo.yloc=y=rdm.uniform(0,colside) } xo.zloc=z=suty ce.o(col[cdx].ix[ES]+ii).zloc=ce.o(col[cdx].ix[EM]+ii).zloc=z//same group // add DPs to NQS; gid(0) id(1) col(2) ty(3) ice(4) xloc(5) yloc(6) zloc(7) // gid id col ty ice x,y,z if(c!=-1) col[cdx].cellsnq.append(gid,xo.id,0,xo.type,0,x,y,suty) // use zloc to store muscle # gid+=1 } } //* wirecols - setup inter-COLUMN connectivity with NetCon proc wirecols () { local i,j if(numcols==1) return for case(&i,0,1) if(EEGainInterC[i] || EIGainInterC[i]) { j = 1 - i gcol[0][i].wire2col(gcol[0][j],colnq[i],1,ncl) // unidirectional wiring } } //* intercoloff(c0,c1) - turn off all weights between COLUMNs proc intercoloff () { local i,c0,c1 localobj xo c0=$1 c1=$2 for ltr(xo,ncl) if(isojt(xo.pre,col.ce.o(0)) && isojt(xo.syn,col.ce.o(0))) { if(xo.pre.col!=c0 || xo.syn.col!=c1) continue for i=0,6 xo.weight(i)=0 } } //* intercolmul(c0,c1,from,to,sy,w) proc intercolsyw () { local c0,c1,from,to,sy,w localobj xo c0=$1 c1=$2 from=$3 to=$4 sy=$5 w=$6 for ltr(xo,ncl) if(isojt(xo.pre,col.ce.o(0)) && isojt(xo.syn,col.ce.o(0))) { if(xo.pre.col!=c0 || xo.syn.col!=c1) continue if(xo.pre.type==from && xo.syn.type==to) xo.weight(sy)=w } } //** drdp() DP locations proc drdp () { local cdx localobj nq nq=col[0].cellsnq gvmarkflag=1 ge(0) nq.marksym="o" for ii=0,5 if (nq.select("ty",DP,"zloc",ii)) nq.gr("yloc","xloc",0,ii+1,8) } //* drcelllocs - draw location of cells proc drcelllocs () { local i,gvt,ty gvt=gvmarkflag {gvmarkflag=1 col[CDX].cellsnq.marksym="O"} if(numarg()>0) { if(col[CDX].cellsnq.select("ty",$1)) col[CDX].cellsnq.gr("yloc","xloc",0,1,3) } else for i=0,1 if(col[CDX].cellsnq.select("ice",i)) col[CDX].cellsnq.gr("yloc","xloc",0,i+2,3) gvmarkflag=gvt } //* drcelldiv(cellid[,skipI,skipE]) - draw outputs from a cell proc drcelldiv () { local id1,id2,x1,y1,x2,y2,ic1,ic2,clr1,clr2,gvt,skipI,skipE localobj cnq,nq gvt=gvmarkflag gvmarkflag=1 cnq=col[CDX].connsnq nq=col[CDX].cellsnq id1=$1 nq.verbose=0 if(!nq.select(-1,"id",id1)) return if(numarg()>1) skipI=$2 else skipI=0 if(numarg()>2) skipE=$3 else skipE=0 x1=nq.v[5].x(nq.ind.x(0)) y1=nq.v[6].x(nq.ind.x(0)) ic1=nq.v[4].x(nq.ind.x(0)) if(ic1) clr1=3 else clr1=2 nq.select("id",id1) nq.gr("yloc","xloc",0,clr1,15) cnq.select(-1,"id1",id1) cnq.verbose=0 for vtr(&i,cnq.ind) { id2 = cnq.v[1].x(i) if(nq.select("id",id2)) { x2=nq.fetch("xloc") y2=nq.fetch("yloc") ic2=nq.fetch("ice") if(ic2) { clr2=3 if(skipI) continue } else { clr2=2 if(skipE) continue } nq.gr("yloc","xloc",0,clr2,8) drline(x1,y1,x2,y2,g,clr1,7) } } cnq.verbose=nq.verbose=1 gvmarkflag=gvt } //* drcellconv(cellid) - draw inputs to a cell proc drcellconv () { local id1,id2,x1,y1,x2,y2,ic1,ic2,clr1,clr2,gvt localobj cnq,nq gvt=gvmarkflag gvmarkflag=1 cnq=col[CDX].connsnq nq=col[CDX].cellsnq id2=$1 nq.verbose=0 nq.select(-1,"id",id2) x2=nq.v[5].x(nq.ind.x(0)) y2=nq.v[6].x(nq.ind.x(0)) ic2=nq.v[4].x(nq.ind.x(0)) if(ic2) clr2=3 else clr2=2 nq.select("id",id2) nq.gr("yloc","xloc",0,clr2,15) cnq.select(-1,"id2",id2) cnq.verbose=0 for vtr(&i,cnq.ind) { id1 = cnq.v[0].x(i) if(nq.select("id",id1)) { x1=nq.fetch("xloc") y1=nq.fetch("yloc") ic1=nq.fetch("ice") if(ic1) clr1=3 else clr1=2 nq.gr("yloc","xloc",0,clr1,8) drline(x1,y1,x2,y2,g,clr1,7) } } cnq.verbose=nq.verbose=1 gvmarkflag=gvt } //* mkwgnq() - return nqs with intracolumnar weight info obfunc mkwgnq () { local a,i,j,id1,id2,wg,ic,ty1,ty2,wt1,wt2,del,dist,c\ localobj mycel,vwt1,vwt2,myv,myv1,myv2,myv3,myv4,myv5,myv6,wgnq,col a=allocvecs(vwt1,vwt2,myv,myv1,myv2,myv3,myv4,myv5,myv6) {nqsdel(wgnq) wgnq=new NQS("col","id1","id2","wg","ice","ty1","ty2","wt1","wt2","del","dist")} for ltr(col,lcol,&c) for i=0,col.allcells-1 { mycel=col.ce.o(i) mycel.getdvi(myv,myv1,myv2,myv3,myv4,myv5,myv6) id1=i ic=ice(mycel.type) ty1=mycel.type vrsz(mycel.getdvi,vwt1,vwt2) if(wsetting_INTF6==1) { mycel.getsywv(vwt1,vwt2) } else { vwt1.copy(myv3) vwt2.copy(myv4) } for j=0,myv6.size-1 { id2=myv.x(j) ty2=col.ce.o(id2).type wg = myv6.x(j) wt1 = wg * vwt1.x(j) if(ic) wt2 = vwt2.x(j) else wt2 = wg * vwt2.x(j) del = myv1.x(j) dist = myv5.x(j) wgnq.append(c,id1,id2,wg,ic,ty1,ty2,wt1,wt2,del,dist) } } dealloc(a) return wgnq } //* function calls to setup network //** z location of each layer in microns mkvlayerz() //** # of cells per column setcpercol() //** setup pmat if(name_declared("nqpmat")==2) { // read pmat from NQS if available, else set to default if(nqpmat!=nil) nq2pmat(nqpmat) else setpmat() } else setpmat() if(pmatscale!=1) scalepmat(pmatscale) //** setup synapse locations,delays,wmat setsynloc() setdelmats() setwmat() // new KMJ version if(wmatscale!=1) scalewmat(wmatscale) scrsz=50*1e3 double scr[scrsz] //** make cells, columns, wire columns {mkcolnqs() printf("finished mkcolnqs\n")} {mkcols() printf("finished mkcols\n")} wirecols(1) // only 1 col at present {mkDPs() printf("finished mkDPs\n")} {swireDPs(col[0].wseed,lambda) printf("finished swireDPs\n")}