// bulb.hoc
// Olfactory bulb network model: network specification file
// Andrew Davison, The Babraham Institute, 2000.

// ¬¬ Initialisation ¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬

objref cvode, random
cvode = new CVode(0)    // start with CVode inactive
random = new Random(seed)
random.MCellRan4(seed+1)
objref mit[nmitx][nmity], gran[ngranx][ngrany]
objref m2gAMPAvec, m2gNMDAvec, g2mvec
//m2gAMPAlist = new List()
//m2gNMDAlist = new List()
//g2mlist = new List()
m2gAMPAvec = new Vector() // indices into pnm.nclist, most often -1
m2gNMDAvec = new Vector()
g2mvec = new Vector()
objref input[nmitx][nmity]
objref outfile
outfile = new File()
strdef filename


// ¬¬  Procedures ¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬

// Create cells --------------------------------------------------------

proc create_cells() { local i,j,p
  netpar_mkpnm()
  quitmodel = 0
if (pnm.myid == 0)  print "Creating cells. This may take some time."
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      mit[i][j] = netpar_create("new Mit()", mitgid, i, j)
    }
/*
    p = 100*(i*nmity+j)/(nmitx*nmity+ngranx*ngrany)
    printf("\r=== %3d\% ===",p)
    flushf()
    doNotify() // Although this slows down cell creation, the
               // process is so long that we have to allow the
               // user to quit during the creation.
*/
  }
  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
      gran[i][j] = netpar_create("new Gran()", grangid, i, j)
    }
/*
    p = 100*(nmitx*nmity + i*ngrany+j)/(nmitx*nmity+ngranx*ngrany)
    printf("\r=== %3d\% ===",p)
    flushf()
    doNotify()
*/
  }
if (pnm.myid == 0)  printf("\n")
//  access mit[0][0].soma // because not on every machine
}

// Connect cells, set synaptic parameters ------------------------------

func wrap() {
  if ($1 < 0) {
    return $2+$1
  } else {
    if ($1 < $2) {
      return $1
    } else {
      return $1-$2
    }
  }
}

proc connect_cells() { local i,j,phi,r,ii,jj,dg,edel // 2 args - dg, fileroot
  dg = $1   // "different glomeruli"
if (pnm.myid == 0) print "Connecting cells"
  m2gAMPAvec.resize(0)
  m2gNMDAvec.resize(0)
  g2mvec.resize(0)
  sprint(filename,"%s.connect",$s2)
if (pnm.myid == 0)  outfile.wopen(filename)
  // Note: here it is possible for a mitral cell to have more than one
  // synaptic contact with any particular granule cell.
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      for k = 1, synpermit {
        phi = random.uniform(0,2*PI)
        r = random.uniform(0,rmax)
        x = dg*i*g2m + r*sin(phi)
        y = dg*j*g2m + r*cos(phi)
        ii = wrap( nint(x),ngranx )
        jj = wrap( nint(y),ngrany )
 if (pnm.myid == 0)       outfile.printf("%d %d\n%5.1f %5.1f %d %d\n\n",dg*i*g2m,dg*j*g2m,x,y,ii,jj)
        //print "Mitral cell [",i,",",j,"] connected to granule cell [",ii,",",jj,"]. "
        edel = edelay + r/rmax*conducdel
//        mit[i][j].dend m2gAMPAlist.append( new NetCon(&v(0.5),gran[ii][jj].AMPAr,thresh,edel,AMPAweight) )
//        mit[i][j].dend m2gNMDAlist.append( new NetCon(&v(0.5),gran[ii][jj].NMDAr,thresh,edel,NMDAweight) )
//        gran[ii][jj].periph g2mlist.append( new NetCon(&v(0.5),mit[i][j].GABAA,thresh,idelay,iweight) )

	m2gAMPAvec.append(netpar_connect(mitgid, i, j, grangid, ii, jj, 0, edel, AMPAweight))
	m2gNMDAvec.append(netpar_connect(mitgid, i, j, grangid, ii, jj, 1, edel, NMDAweight))
	g2mvec.append(netpar_connect(grangid, ii, jj, mitgid, i, j, 1, idelay, iweight))
      }
    }
  }
if (pnm.myid == 0)  outfile.close()
}

proc set_GABAA_weights() { local i // 1 arg - weight
  for i = 0,g2mvec.size()-1 {
//    g2mlist.object(i).weight = $1
	netpar_weight(g2mvec, i, $1)
  }
}

proc set_AMPA_weights() { local i // 1 arg - weight
  for i = 0,m2gAMPAvec.size()-1 {
//    m2gAMPAlist.object(i).weight = $1
	netpar_weight(m2gAMPAvec, i, $1)
  }
}

proc set_NMDA_weights() { local i // 1 arg - weight
  for i = 0,m2gNMDAvec.size()-1 {
//    m2gNMDAlist.object(i).weight = $1
	netpar_weight(m2gNMDAvec, i, $1)
  }
}

proc randomise_NMDA() { local i // 2 args - mean weight, variance
//  m2gNMDAlist.object(0).weight = random.normal($1,$2)
  netpar_weight(m2gNMDAvec, 0, random.normal($1,$2))
  //m2gNMDAlist.object(0).weight = random.poisson($1)
  for i = 1,m2gNMDAvec.size()-1 {
//    m2gNMDAlist.object(i).weight = random.repick()
	netpar_weight(m2gNMDAvec, i, random.repick())
  }
}

proc set_NMDA_time_constants() { local i,j, alpha, beta
				           // 2 args - rise time constant (ms), 
				           //     decay time constant (ms)
  if ($1 < 1e-9) { $1 = 1e-9 }
  if ($2 < 1e-9) { $2 = 1e-9 }
  beta = 1/$2
  alpha = 1/$1 - beta
   

  for i = 0, ngranx-1 {
    for j = 0, ngrany-1 {
	if (pnm.gid_exists(grangid.x[i][j])) {
      gran[i][j].NMDAr.Alpha = alpha
      gran[i][j].NMDAr.Beta = beta
	}
    }
  }
}

proc set_GABAA_time_constant() { local i,j, tau, gid
				           // 1 args - decay time constant (ms)
  if ($1 < 1e-9) { $1 = 1e-9 }
  tau = $1

  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
	if (pnm.gid_exists(mitgid.x[i][j])) {
      mit[i][j].GABAA.tau = tau
	}
    }
  }
}

// Add input currents --------------------------------------------------

proc insert_iclamps() { local i,j, ran // 2 args - del dur
  // if $1 is negative, delay is randomly chosen in the uniform interval 0,$1
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
	if (pnm.gid_exists(mitgid.x[i][j])) {
      mit[i][j].glom input[i][j] = new IClamp(0.5)
      input[i][j].dur = $2
      input[i][j].del = abs($1)
	}
    }
  }
  random.uniform(0,abs($1))
  if ($1 < 0) {
    for i = 0, nmitx-1 {
      for j = 0, nmity-1 {
         ran = random.repick()
	if (pnm.gid_exists(mitgid.x[i][j])) {
         input[i][j].del = ran
	}
      }
    }
  }
}

// Randomise initial conditions ----------------------------------------

proc random_init() { local i,j, ran
  random.normal(-65,25)
  for i = 0,nmitx-1 {
    for j = 0, nmity-1 {
      ran = random.repick()
	if (pnm.gid_exists(mitgid.x[i][j])) {
      mit[i][j].soma.v(0.5) = ran
      mit[i][j].dend.v(0.5) = mit[i][j].soma.v(0.5)
      mit[i][j].prim.v(0.5) = mit[i][j].soma.v(0.5)
      mit[i][j].glom.v(0.5) = mit[i][j].soma.v(0.5)
	}
    }
  }
  for i = 0,ngranx-1 {
    for j = 0, ngrany-1 {
      ran = random.repick()
	if (pnm.gid_exists(grangid.x[i][j])) {
      gran[i][j].soma.v(0.5) = ran
      gran[i][j].deep.v(0.5) = gran[i][j].soma.v(0.5)
      gran[i][j].periph.v(0.5) = gran[i][j].soma.v(0.5)
	}
    }
  }
}

// ¬¬ Create the model ¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬

create_cells()
connect_cells(diffglom,fileroot)
insert_iclamps(-200,tstop)

// set synaptic properties
Cdur_NMDA = NMDArisetime
mg_NMDA = mgconc	
set_AMPA_weights(AMPAweight)
set_NMDA_weights(NMDAweight)
set_GABAA_weights(iweight)
set_NMDA_time_constants(NMDArise,NMDAdecay)