// 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 m2gAMPAlist, m2gNMDAlist, g2mlist
m2gAMPAlist = new List()
m2gNMDAlist = new List()
g2mlist = new List()
objref input[nmitx][nmity]
objref outfile
outfile = new File()
strdef filename


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

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

proc create_cells() { local i,j,p
  quitmodel = 0
  print "Creating cells. This may take some time."
  for i = 0, nmitx-1 {
    for j = 0, nmity-1 {
      mit[i][j] = new Mit()
    }
    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] = new Gran()
    }
    p = 100*(nmitx*nmity + i*ngrany+j)/(nmitx*nmity+ngranx*ngrany)
    printf("\r=== %3d\% ===",p)
    flushf()
    doNotify()
  }
  printf("\n")
  access mit[0][0].soma
}

// 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"
  print "Connecting cells"
  m2gAMPAlist.remove_all()
  m2gNMDAlist.remove_all()
  g2mlist.remove_all()
  sprint(filename,"%s.connect",$s2)
  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 )
        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) )
      }
    }
  }
  outfile.close()
}

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

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

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

proc randomise_NMDA() { local i // 2 args - mean weight, variance
  m2gNMDAlist.object(0).weight = random.normal($1,$2)
  //m2gNMDAlist.object(0).weight = random.poisson($1)
  for i = 1,m2gNMDAlist.count()-1 {
    m2gNMDAlist.object(i).weight = 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 {
      gran[i][j].NMDAr.Alpha = alpha
      gran[i][j].NMDAr.Beta = beta
    }
  }
}

proc set_GABAA_time_constant() { local i,j, tau
				           // 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 {
      mit[i][j].GABAA.tau = tau
    }
  }
}

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

proc insert_iclamps() { local i,j // 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 {
      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 {
         input[i][j].del = random.repick()
      }
    }
  }
}

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

proc random_init() { local i,j
  random.normal(-65,25)
  for i = 0,nmitx-1 {
    for j = 0, nmity-1 {
      mit[i][j].soma.v(0.5) = random.repick()
      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 {
      gran[i][j].soma.v(0.5) = random.repick()
      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)