SCZ-associated variant effects on L5 pyr cell NN activity and delta osc. (Maki-Marttunen et al 2018)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:237469
" … Here, using computational modeling, we show that a common biomarker of schizophrenia, namely, an increase in delta-oscillation power, may be a direct consequence of altered expression or kinetics of voltage-gated ion channels or calcium transporters. Our model of a circuit of layer V pyramidal cells highlights multiple types of schizophrenia-related variants that contribute to altered dynamics in the delta frequency band. Moreover, our model predicts that the same membrane mechanisms that increase the layer V pyramidal cell network gain and response to delta-frequency oscillations may also cause a decit in a single-cell correlate of the prepulse inhibition, which is a behavioral biomarker highly associated with schizophrenia."
Reference:
1 . Mäki-Marttunen T, Krull F, Bettella F, Hagen E, Næss S, Ness TV, Moberget T, Elvsåshagen T, Metzner C, Devor A, Edwards AG, Fyhn M, Djurovic S, Dale AM, Andreassen OA, Einevoll GT (2019) Alterations in Schizophrenia-Associated Genes Can Lead to Increased Power in Delta Oscillations. Cereb Cortex 29:875-891 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex L5/6 pyramidal GLU cell;
Channel(s): Ca pump; I A, slow; I h; I K; I K,Ca; I K,leak; I L high threshold; I M; I Na,p; I Na,t; I T low threshold;
Gap Junctions: Gap junctions;
Receptor(s): AMPA; NMDA; Gaba;
Gene(s): Cav1.2 CACNA1C; Cav1.3 CACNA1D; Cav3.3 CACNA1I; HCN1; Kv2.1 KCNB1; Nav1.1 SCN1A; PMCA ATP2B2;
Transmitter(s): Glutamate; Gaba;
Simulation Environment: NEURON; Python; LFPy;
Model Concept(s): Schizophrenia; Oscillations;
Implementer(s): Maki-Marttunen, Tuomo [tuomo.maki-marttunen at tut.fi];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; AMPA; NMDA; Gaba; I Na,p; I Na,t; I L high threshold; I T low threshold; I K; I K,leak; I M; I h; I K,Ca; I A, slow; Ca pump; Gaba; Glutamate;
/
scznet
approxhaynet
models
README.html
Ca_HVA.mod *
Ca_LVAst.mod *
CaDynamics_E2.mod *
Ih.mod *
Im.mod *
K_Pst.mod *
K_Tst.mod *
Nap_Et2.mod *
NaTa_t.mod *
ProbAMPANMDA2.mod
ProbAMPANMDA2group.mod
ProbAMPANMDA2groupdet.mod
ProbUDFsyn2.mod
ProbUDFsyn2group.mod
ProbUDFsyn2groupdet.mod
SK_E2.mod *
SKv3_1.mod *
approxhaynetstuff.py
calcEEG.py
calcEEG_uncombined.py
calcmutgains.py
calcmutgains_comb.py
calcmutoscs.py
calcmutoscs_comb.py
calcspectra.py
calcspectra_comb.py
combinemutgains_parallel.py
combinemutgains_parallel_comb.py
combinemutgains_parallel_withLFP.py
combinemutoscs_parallel.py
combinemutoscs_parallel_comb.py
drawfig1ab.py
drawfig1c.py
drawfig2ab.py
drawfig2c.py
drawstationary_EEG.py
drawstationary_EEG_pop.py
mutation_stuff.py *
mytools.py
pars_withmids_combfs_final.sav *
runmanycellsLFP.py
runsinglecellLFP.py
saveisidists.py
savespikesshufflephases.py
scalings_cs.sav
simosc_parallel.py
simosc_parallel_comb_varconn.py
simseedburst_func.py
simseedburst_func_comb_varconn.py
simseedburst_func_withLFP.py
simseedburst_func_withLFP.pyc
                            
from mpi4py import MPI
from neuron import h
import matplotlib
matplotlib.use('Agg')
import numpy
from pylab import *
import time
import scipy.io
import pickle
import sys
import mutation_stuff
import approxhaynetstuff
import mytools
import resource

def simseedburst_func(Nmc=1, tstop=10200,mutID=0,rdSeed=1,Econ=0.0004,Icon=0.001,nseg=5,rateCoeff=1.0,gNoiseCoeff=1.0,gSynCoeff=1.0,Ncells2save=1,sparsedt=1.0,Nsyns2save=1,connM=[],gToBlock=[],blockEfficiency=0.0):
  rank = MPI.COMM_WORLD.Get_rank()
  dt = 0.025

  coeffCoeffs = [[0.25,0],[0.125,0],[0.5,0],[0.5,1.0/3],[0.5,2.0/3],[0.5,1.0],[-0.25,0],[-0.125,0],[-0.5,0]]
  MT = mutation_stuff.getMT(False)
  defVals = mutation_stuff.getdefvals()
  keyList = defVals.keys()
  for idefval in range(0,len(keyList)):
    if type(defVals[keyList[idefval]]) is not list:
      defVals[keyList[idefval]] = [defVals[keyList[idefval]], defVals[keyList[idefval]]] #make the dictionary values [somatic, apical]
  updatedVars = ['somatic','apical','basal'] # the possible classes of segments that defVals may apply to
  whichDefVal = [0,1,0]                      # use the defVal[0] for somatic and basal segments and defVal[1] for apical segments
  unpicklefile = open('scalings_cs.sav', 'r')
  unpickledlist = pickle.load(unpicklefile)
  unpicklefile.close()
  theseCoeffsAllAll = unpickledlist[0]

  filename = 'pars_withmids_combfs_final'
  unpicklefile = open(filename+".sav", 'r')
  unpickledlist = pickle.load(unpicklefile)
  unpicklefile.close()
  par_names = unpickledlist[0]
  par_values = unpickledlist[1]
  paramdict = {}
  for i in range(0,len(par_names)):
    paramdict[par_names[i]] = par_values[i]

  h("""
{load_file("stdlib.hoc")}
{load_file("stdrun.hoc")}

initialization_tstart = startsw()

strdef fileName
objref fileObj

fileObj = new File()

rdSeed = """+str(rdSeed)+"""
Nmc = """+str(Nmc)+"""
connectivity = 1
noisyst = 0

tstop = """+str(tstop)+"""
rcpWeightFactor = 1.5 // the factor by which reciprocal weights are stronger than unidirectional weights
pT2Tr = 0.06 //probability of reciprocating an existing connection to another L5bPC
pT2T = 0.13 //probability of a L5bPC being connected to another L5bPC
Econ = """+str(Econ)+""" //excitatory synaptic conductance
Icon = """+str(Icon)+""" //inhibitory synaptic conductance
NcontE = 5 // number of excitatory synaptic contacts per connection
NsynE = 10000 // number of excitatory synapses
NsynI = 2500 // number of inhibitory synapses
gNoiseCoeff = """+str(gNoiseCoeff)+""" // scaling of background synaptic conductances
rateE = """+str(0.72*rateCoeff)+""" // average rate of presynaptic excitatory cells
rateI = """+str(7.0*rateCoeff)+""" // average rate of presynaptic inhibitory cells
mainBifurcation = 650

{Ncells2save = """+str(Ncells2save)+"""}
sparsedt = """+str(sparsedt)+""" // recordings of [Ca], I_SK and vApical are done with low temporal resolution to save memory
{gSynCoeff = """+str(gSynCoeff)+"""}

objref tempvec
tempvec = new Vector()
{tempvec.append(Nmc)}
{Ncells2save = """+str(Ncells2save)+"""}
""")

  h("""
{load_file("models/TTC.hoc")}

objref MC_TTC
objref sl //synaptic locations list

objref rds1,rds2,rds3
{rds1 = new Random(1000*rdSeed)}
{rds1.uniform(0,1)} //random for microcircuit connectivity and noisyst

objref conMat
conMat = new Matrix(Nmc,Nmc)

for(i=0;i<Nmc;i+=1){
        conMat.x[i][i]=0
}
""")
  h('objref cvode')                                #Define an object cvode
  h('cvode = new CVode()')                         #Make cvode a time step integrator object
  h('cvode.active(0)')                             #Set the variable time step integration on/off
  h('cvode.use_fast_imem(1)')                      #Make i_membrane_ accessible

  if len(connM) == 0:
    h("""
for(i=0;i<(Nmc-2);i+=1){
        for(j=(i+1);j<Nmc;j+=1){
                if (connectivity){
                        pcon = rds1.repick()
                        if (pcon<pT2Tr){
                                conMat.x[i][j]=rcpWeightFactor*gSynCoeff
                                conMat.x[j][i]=rcpWeightFactor*gSynCoeff
                        } else {
                                if (pcon<(pT2Tr + 0.5*pT2T)){
                                        conMat.x[i][j]=gSynCoeff
                                        conMat.x[j][i]=0
                                } else {
                                        if (pcon<(pT2Tr + pT2T)){
                                                conMat.x[i][j]=0
                                                conMat.x[j][i]=gSynCoeff
                                        } else {
                                                conMat.x[i][j]=0
                                                conMat.x[j][i]=0
                                        }
                                }
                        }
                } else {
                        conMat.x[i][j]=0
                        conMat.x[j][i]=0
                }
        }
}
""")
  else:
    for i in range(0,Nmc):
      for j in range(0,Nmc):
        if connM[i][j]:
          h("conMat.x["+str(i)+"]["+str(j)+"]="+str(gSynCoeff*connM[i][j])) # Remember that conMat.x[i][j] is connection FROM j TO i  print "Connectivity OK!"

  h("""
{load_file("netparmpi.hoc")}
objref epnm

epnm = new ParallelNetManager(Nmc)
{epnm.round_robin()}

strdef treename
objref NsynsE, NsynsI, preTrainList
""")

  for i in range(0,Nmc):
    h("""
  i = """+str(i)+"""
  if (epnm.gid_exists(i)) {
        print \"rank = """+str(rank)+""", gid \", i, \" exists\\n\"
        MC_TTC = new TTC()
        epnm.register_cell(i,MC_TTC)
        epnm.pc.gid2cell(i).initRand(1000*rdSeed+i)
        epnm.pc.gid2cell(i).setnetworkparameters(rcpWeightFactor,Econ,Icon,NsynE,NsynI,NcontE,1.0,1.0,1.0,gNoiseCoeff)
  }""")

  approxhaynetstuff.setparams(paramdict,Nmc)

  for i in range(0,Nmc):
    h("""
  i = """+str(i)+"""
  if (epnm.gid_exists(i)) {
    lengthA = epnm.pc.gid2cell(i).apic[0].L + epnm.pc.gid2cell(i).apic[1].L
    lengthB = epnm.pc.gid2cell(i).dend.L
    pA = lengthA/(lengthA + lengthB)
    {NsynsE = new List()}
    {NsynsI = new List()}
    for i1 = 0, 2 {
      {NsynsE.append(new Vector("""+str(nseg)+"""))}
      {NsynsI.append(new Vector("""+str(nseg)+"""))}
    }
    for(i1=0;i1<(NsynE+NsynI);i1+=1){
      if (epnm.pc.gid2cell(i).rd1.repick()<pA){
        treename = "apic"
        compInd = 1
      } else {
        treename = "dend"
        compInd = 0
      }
      sl = epnm.pc.gid2cell(i).locateSites(treename,epnm.pc.gid2cell(i).rd1.repick()*epnm.pc.gid2cell(i).getLongestBranch(treename))
      sitenum = int((sl.count()-1)*epnm.pc.gid2cell(i).rd1.repick())
      compInd = compInd + sl.o[sitenum].x[0] // if we are at apical, and sl.o[sitenum].x[0]=1, then compInd = 2, otherwise 1 at apical, and 0 at basal
      segInd = int(sl.o[sitenum].x[1]*"""+str(nseg)+""")
      if (i1<NsynE) {
        NsynsE.o[compInd].x[segInd] = NsynsE.o[compInd].x[segInd] + 1
      } else {
        NsynsI.o[compInd].x[segInd] = NsynsI.o[compInd].x[segInd] + 1
      }
    }
    {epnm.pc.gid2cell(i).distributeSyn(NsynsE,NsynsI)}

    {preTrainList = new List()}
    {rds2 = new Random(1000*rdSeed+i)}//random for presynaptic trains
    {rds3 = new Random(1000*rdSeed+i)}//random for presynaptic trains
    {rds2.negexp(1/rateE)}
    {rds3.negexp(1/rateI)}

    for(compInd=0;compInd<3;compInd+=1){
      for(segInd=0;segInd<"""+str(nseg)+""";segInd+=1){
        {preTrainList.append(new Vector())}
        pst=0 //presynaptic spike time
        if(NsynsE.o[compInd].x[segInd]==0) {
          print "Warning: NsynsE.o[",compInd,"].x[",segInd,"] = 0!!!!"
          pst = 1e6
        }
        while(pst < tstop){
          pst+= 1000*rds2.repick()/NsynsE.o[compInd].x[segInd]
          {preTrainList.o[preTrainList.count()-1].append(pst)}
        }
      }
    }
    for(compInd=0;compInd<3;compInd+=1){
      for(segInd=0;segInd<"""+str(nseg)+""";segInd+=1){
        {preTrainList.append(new Vector())}
        pst=0 //presynaptic spike time
        if(NsynsI.o[compInd].x[segInd]==0) {
          print "Warning: NsynsI.o[",compInd,"].x[",segInd,"] = 0!!!!"
          pst = 1e6
        }
        while(pst < tstop){
          pst+= 1000*rds3.repick()/NsynsI.o[compInd].x[segInd]
          {preTrainList.o[preTrainList.count()-1].append(pst)}
        }
      }
    }
    {epnm.pc.gid2cell(i).setpretrains(preTrainList)}
    {epnm.pc.gid2cell(i).queuePreTrains()}

    //print \"distributeSyn(), setpretrains(), queuePreTrains() OK! i=\", i, \".\"
  }
""")

  #LOAD MUTATIONS
  counter = -1
  found = 0
  for igene in range(0, len(MT)):
    for imut in range(0, len(MT[igene])):
      theseCoeffs = theseCoeffsAllAll[0][igene][imut]

      nVals = len(MT[igene][imut])*[0]
      thesemutvars = []
      for imutvar in range(0,len(MT[igene][imut])):
        thesemutvars.append(MT[igene][imut][imutvar][0])
        if type(MT[igene][imut][imutvar][1]) is int or type(MT[igene][imut][imutvar][1]) is float:
          MT[igene][imut][imutvar][1] = [MT[igene][imut][imutvar][1]]
        nVals[imutvar] = len(MT[igene][imut][imutvar][1])
      cumprodnVals = cumprod(nVals)
      allmutvars = cumprodnVals[len(MT[igene][imut])-1]*[thesemutvars]
      allmutvals = []
      for iallmutval in range(0,cumprodnVals[len(MT[igene][imut])-1]):
        allmutvals.append([0]*len(thesemutvars))
      for iallmutval in range(0,cumprodnVals[len(MT[igene][imut])-1]):
        for imutvar in range(0,len(MT[igene][imut])):
          if imutvar==0:
            allmutvals[iallmutval][imutvar] = MT[igene][imut][imutvar][1][iallmutval%nVals[imutvar]]
          else:
            allmutvals[iallmutval][imutvar] = MT[igene][imut][imutvar][1][(iallmutval/cumprodnVals[imutvar-1])%nVals[imutvar]]

      for iallmutval in range(0, len(theseCoeffs)):
        if igene == 0 and imut == 0 and iallmutval == 0:
          iters = [-1,0,2,6,8]
        else:
          iters = [0,2,6,8]
        for iiter in range(0,len(iters)):
          iter = iters[iiter]
          counter = counter + 1
          if counter == mutID:
            found = 1
            break
        if found:
          break
      if found:
        break
    if found:
      break
  if not found:
    print "Not found corresponding mutID, mutID = " + str(mutID)
    sys.exit()

  if iter >= 0:
    thisCoeff = coeffCoeffs[iter][0]*theseCoeffs[iallmutval] + coeffCoeffs[iter][1]*(1.0 - 0.5*theseCoeffs[iallmutval])
  else:
    thisCoeff = 0

  mutText = ""
  for imutvar in range(0,len(MT[igene][imut])):
    if imutvar > 0 and imutvar%2==0:
      mutText = mutText+"\n"
    mutvars = allmutvars[iallmutval][imutvar]
    mutvals = allmutvals[iallmutval][imutvar]
    if type(mutvars) is str:
      mutvars = [mutvars]
    mutText = mutText + str(mutvars) + ": "
    for kmutvar in range(0,len(mutvars)):
      mutvar = mutvars[kmutvar]
      if mutvar.find('offm') > -1 or mutvar.find('offh') > -1 or mutvar.find('ehcn') > -1:
        newVal =  [x+mutvals*thisCoeff for x in defVals[mutvar]]
        if mutvals >= 0 and kmutvar==0:
          mutText = mutText + "+" + str(mutvals) +" mV"
        elif kmutvar==0:
          mutText = mutText  + str(mutvals) +" mV"
      else:
        newVal = [x*(mutvals**thisCoeff) for x in defVals[mutvar]]
        if kmutvar==0:
          mutText = mutText + "*" + str(mutvals)
      if kmutvar < len(mutvars)-1:
        mutText = mutText + ", "
      if mutvar.find('_Ih') > -1:
        updateThese = [1,1,1]
      elif mutvar.find('_Ca_HVA') > -1 or mutvar.find('_Ca_LVAst') > -1 or mutvar.find('_SKv3.1') > -1 or mutvar.find('_Ca_HVA') > -1 or mutvar.find('_SK_E2') > -1 or mutvar.find('_NaTa_t') > -1 or mutvar.find('_CaDynamics_E2') > -1:
        updateThese = [1,1,0]
      elif mutvar.find('_K_Pst') > -1 or mutvar.find('_K_Tst') > -1 or mutvar.find('_Nap_Et2') > -1:
        updateThese = [1,0,0]
      elif mutvar.find('_Im') > -1:
        updateThese = [0,1,0]
      else:
        print "Error: str=" + str(mutvar)
        updatedThese = [0,0,0]
      for iupdated in range(0,3):
        if updateThese[iupdated]:
          print """forsec epnm.pc.gid2cell(i)."""+str(updatedVars[iupdated])+""" {
"""+mutvar+""" = """+str(newVal[whichDefVal[iupdated]])+"""
}""" 
          for i in range(0,Nmc):
            h("""
i = """+str(i)+"""
if (epnm.gid_exists(i)) {
  forsec epnm.pc.gid2cell(i)."""+str(updatedVars[iupdated])+""" {
    """+mutvar+""" = """+str(newVal[whichDefVal[iupdated]])+"""
  }
}""")

  print mutText
  h("""
thisCa = 0.0001
for(i=0;i<Nmc;i+=1){
  if (epnm.gid_exists(i)) {
    thisCa = epnm.pc.gid2cell(i).soma.minCai_CaDynamics_E2
  }
}
""")
  thisCa = h.thisCa

  myMechs = ['Ca_HVA','Ca_LVAst','Ih','Im','K_Pst','K_Tst','NaTa_t','Nap_Et2','SK_E2','SKv3_1','']
  myMechToAdd = ""
  for iblock in range(0,len(gToBlock)):
    for iMech in range(0,len(myMechs)):
      if gToBlock[iblock] in myMechs[iMech]:
        break
    if iMech <= 9:
      if type(blockEfficiency) is list:
        for i in range(0,Nmc):
          h("""
i = """+str(i)+"""
if (epnm.gid_exists(i)) {
  epnm.pc.gid2cell(i).soma g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" = g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" * """+str(blockEfficiency[0])+"""
  forsec epnm.pc.gid2cell(i).apical { g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" = g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" * """+str(blockEfficiency[1])+""" }
}""")
        print("""
i = """+str(i)+"""
if (epnm.gid_exists(i)) {
  epnm.pc.gid2cell(i).soma g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" = g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" * """+str(blockEfficiency[0])+"""
  forsec epnm.pc.gid2cell(i).apical { g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" = g"""+str(myMechs[iMech])+"""bar_"""+str(myMechs[iMech])+""" * """+str(blockEfficiency[1])+""" }
}""")
        myMechToAdd = myMechToAdd+myMechs[iMech]+'x'+str(blockEfficiency[0])+'-'+str(blockEfficiency[1])+'_'        
      else:
        print("forall if(ismembrane(\""+str(myMechs[iMech])+"\")) { g"+str(myMechs[iMech])+"bar_"+str(myMechs[iMech])+" = g"+str(myMechs[iMech])+"bar_"+str(myMechs[iMech])+" * "+str(blockEfficiency)+" }")
        h("forall if(ismembrane(\""+str(myMechs[iMech])+"\")) { g"+str(myMechs[iMech])+"bar_"+str(myMechs[iMech])+" = g"+str(myMechs[iMech])+"bar_"+str(myMechs[iMech])+" * "+str(blockEfficiency)+" }")
        myMechToAdd = myMechToAdd+myMechs[iMech]+'x'+str(blockEfficiency)+'_'
    else:
      print "Error: No mechanism recognized"


  h("""
v_init = -80
cai0_ca_ion = thisCa
dt = """+str(dt)+"""
objref syninds, conMatRows

for(i=0;i<Nmc;i+=1){
  if (epnm.gid_exists(i)) {
    //epnm.pc.gid2cell(i).geom_nseg()
    epnm.pc.gid2cell(i).insertMCcons(conMat.getcol(i))
  }
}

{syninds = new Vector()}
{conMatRows = new List()}
for(i=0;i<Nmc;i+=1){
        syninds.append(2*3*"""+str(nseg)+""")
}

// appending the microcircuit connections
for(i=0;i<Nmc;i+=1){
        conMatRows.append(new Vector())
        for(j=0;j<Nmc;j+=1){
                conMatRows.o[i].insrt(j,conMat.x[j][i])
                if (conMat.x[j][i] != 0){
                        for(jj=0;jj<NcontE;jj+=1){
                                epnm.nc_append(j,i,syninds.x[i],1,0.5)
                                syninds.x[i] +=1
                        }
                }
        }
}
""")
  h("forall nseg="+str(nseg))

  h("""
objref vSomaList, tvecList, caSomaList, skSomaList, cahvaSomaList, calvaSomaList
objref natSomaList, napSomaList, ihSomaList, kv31SomaList, ktSomaList, kpSomaList, IList
objref apcvecList, apcList, netcon, nil, spikes, spikedCells
{spikes = new Vector()}
{spikedCells = new Vector()}

{apcvecList = new List()}
{apcList = new List()}
{vSomaList = new List()}
{caSomaList = new List()}
{skSomaList = new List()}
{cahvaSomaList = new List()}
{calvaSomaList = new List()}
{natSomaList = new List()}
{napSomaList = new List()}
{ihSomaList = new List()}
{kv31SomaList = new List()}
{ktSomaList = new List()}
{kpSomaList = new List()}
{IList = new List()}
{tvecList = new List()}
""")

  Nsyns = numpy.array(h.syninds)-2*3*nseg
  cumpNsyns = numpy.cumsum(Nsyns)/numpy.sum(Nsyns)
  randVec = [rand() for x in range(0,Nsyns2save)]
  if sum(Nsyns) > 0:
    cellsSynRecorded = [next(i for i,x in enumerate(cumpNsyns) if x > randVec[j]) for j in range(0,Nsyns2save)]
    synIndsRecorded = [int(2*3*nseg+rand()*Nsyns[i]) for i in cellsSynRecorded]
  else:
    cellsSynRecorded = []
    synIndsRecorded = []

  for i in range(0,int(h.Ncells2save)):
    h("""
  i = """+str(i)+"""
  if (epnm.gid_exists(i)) {
    {vSomaList.append(new Vector())}
    {caSomaList.append(new Vector())}
    {skSomaList.append(new Vector())}
    {cahvaSomaList.append(new Vector())}
    {calvaSomaList.append(new Vector())}
    {natSomaList.append(new Vector())}
    {napSomaList.append(new Vector())}
    {ihSomaList.append(new Vector())}
    {kv31SomaList.append(new Vector())}
    {ktSomaList.append(new Vector())}
    {kpSomaList.append(new Vector())}
    {tvecList.append(new Vector())}

    access epnm.pc.gid2cell(i).soma
    {vSomaList.o[vSomaList.count()-1].record(&v(0.5),dt)}
    {caSomaList.o[caSomaList.count()-1].record(&cai(0.5),sparsedt)}
    {skSomaList.o[skSomaList.count()-1].record(&ik_SK_E2(0.5),sparsedt)}
    {cahvaSomaList.o[skSomaList.count()-1].record(&ica_Ca_HVA(0.5),sparsedt)}
    {calvaSomaList.o[skSomaList.count()-1].record(&ica_Ca_LVAst(0.5),sparsedt)}
    {natSomaList.o[skSomaList.count()-1].record(&ina_NaTa_t(0.5),sparsedt)}
    {napSomaList.o[skSomaList.count()-1].record(&ina_Nap_Et2(0.5),sparsedt)}
    {ihSomaList.o[skSomaList.count()-1].record(&ihcn_Ih(0.5),sparsedt)}
    {kv31SomaList.o[skSomaList.count()-1].record(&ik_SKv3_1(0.5),sparsedt)}
    {ktSomaList.o[skSomaList.count()-1].record(&ik_K_Tst(0.5),sparsedt)}
    {kpSomaList.o[skSomaList.count()-1].record(&ik_K_Pst(0.5),sparsedt)}
  }
""")
    indSynIndsRecorded = [ix for ix,x in enumerate(cellsSynRecorded) if x==i]
    for isyn in range(0,len(indSynIndsRecorded)):
      h("""
  if (epnm.gid_exists(i)) {
    {IList.append(new Vector())}
    {IList.o[IList.count()-1].record(&epnm.pc.gid2cell(i).synlist.o["""+str(synIndsRecorded[indSynIndsRecorded[isyn]])+"""].i, sparsedt)}
  }
""")
    if i < Ncells2save:
      h("""
  if (epnm.gid_exists(i)) {
    {vSomaList.append(new Vector())}
    {vSomaList.o[vSomaList.count()-1].record(&v(0.5),dt)}
  }
""")

  for i in range(0,Nmc):
    h("""
  i = """+str(i)+"""
  if (epnm.gid_exists(i)) {
          access epnm.pc.gid2cell(i).soma
          {apcList.append(new APCount(0.5))}
          {apcvecList.append(new Vector())}
          apcList.o[apcList.count()-1].thresh= -40
          {apcList.o[apcList.count()-1].record(apcvecList.o[apcList.count()-1])}
          {netcon = new NetCon(&v(0.5), nil)}
          netcon.threshold = -20
          {netcon.record(spikes, spikedCells)  }
}""")

  h("""
{epnm.set_maxstep(100)}

stdinit()

if (epnm.gid_exists(0)) {
        print \"\\n\"
        sim_tstart = startsw()
        initializationtime = (sim_tstart-initialization_tstart)/3600
        print \"Initialization completed. Initialization took \", initializationtime, \" hours\\n\"
        print \"Starting simulation\\n\"
        print \"\\n\"
}
""")
  #h("""
  #{epnm.psolve(tstop)}
  #
  #if (epnm.gid_exists(0)) {
  #        simruntime = (startsw() - sim_tstart)/3600
  #        print \"Simulation took \", simruntime, \" hours\\n\"
  #}
  #""")
  tstep = 0
  use_ipas = False
  use_icap = False
  use_isyn = False
  rec_current_dipole_moment = True
  h("objref allseclist,memirec,memireclist")
  h("allseclist = new SectionList()")
  h("memireclist = new SectionList()")
  h("memirec = new Vector()")

  dlring = 56.0; dlringI = 110.0
  pos = []; posI = []
  for iring in range(0,40):
    r = iring * dlring
    rI = iring * dlringI
    if iring==0:
      NcellsThisRing=1
      NcellsThisRingI=1
    else:
      NcellsThisRing = int(3.1416*2*r/dlring+0.5)
      NcellsThisRingI = int(3.1416*2*rI/dlringI+0.5)
    for icell in range(0,NcellsThisRing):
      pos.append([r*cos(2*3.1416*(icell+(iring%2)*0.5)/NcellsThisRing), r*sin(2*3.1416*(icell+(iring%2)*0.5)/NcellsThisRing)])
    for icell in range(0,NcellsThisRingI):
      posI.append([rI*cos(2*3.1416*(icell+(iring%2)*0.5)/NcellsThisRingI), rI*sin(2*3.1416*(icell+(iring%2)*0.5)/NcellsThisRingI)])
    
  midpointsSecs = []
  midpoints = []
  for i in range(0,Nmc):
    h("""
  i = """+str(i)+"""
  found = 0
  if (epnm.gid_exists(i)) {
    epnm.pc.gid2cell(i).soma allseclist.append()
    epnm.pc.gid2cell(i).apic[0] allseclist.append()    
    epnm.pc.gid2cell(i).apic[1] allseclist.append()    
    epnm.pc.gid2cell(i).dend allseclist.append()    
    found = 1
  }
  """)
    if h.found > 0:
      midpointsSecs.append([pos[i][0],0,pos[i][1]])    
      midpointsSecs.append([pos[i][0],paramdict['L_soma']*0.5+paramdict['L_apic[0]']*0.5,pos[i][1]])                        #apic0 midpoint
      midpointsSecs.append([pos[i][0],paramdict['L_soma']*0.5+paramdict['L_apic[0]']+paramdict['L_apic[1]']*0.5,pos[i][1]]) #apic1 midpoint
      midpointsSecs.append([pos[i][0],-paramdict['L_soma']*0.5-paramdict['L_dend']*0.5,pos[i][1]]) #apic1 midpoint

  iseg = 0
  isec = 0
  for sec in h.allseclist:
    for seg in sec:
      if iseg < nseg*4*Nmc:
        if sec.name().find('soma') > -1 or sec.name().find('apic') > -1:
          midpoints.append([midpointsSecs[isec][0], midpointsSecs[isec][1]-sec.L/2+sec.L*seg.x, midpointsSecs[isec][2]])
        else:
          midpoints.append([midpointsSecs[isec][0], midpointsSecs[isec][1]+sec.L/2-sec.L*seg.x, midpointsSecs[isec][2]])
      else:
        midpoints.append([midpointsSecs[isec][0], midpointsSecs[isec][1], midpointsSecs[isec][2]])
      #h.memirec.record(seg._ref_i_membrane_, dt) #(not needed in NEURON version 7.5!)
      iseg = iseg+1
    isec = isec+1
  NsegTot = iseg
  
  imem = {'imem': numpy.zeros([NsegTot,])}
  DIPOLE_MOMENT = []
  dtPrint = 500
  nextPrint = dtPrint
  timeNow = time.time()
  while h.t < tstop:
    if h.t >= 0:
      i = 0
      totnsegs = 0
      if use_isyn:
        imem['isyn_e'] = 0. # need to reset these for every iteration                                                                                                                                                              
        imem['isyn_i'] = 0. # because we sum over synapses                                                                                                                                                                         
      for sec in h.allseclist:
        for seg in sec:
          imem['imem'][i] = seg.i_membrane_
          if use_ipas:
            imem['ipas'][i] = seg.i_pas
          if use_icap:
            imem['icap'][i] = seg.i_cap
          i += 1

      if rec_current_dipole_moment:
        DIPOLE_MOMENT.append(numpy.dot(imem['imem'], midpoints).tolist())
    h.fadvance()
    if h.t > nextPrint:
      print "time="+str(h.t)+", used "+str(time.time()-timeNow)+" seconds"
      nextPrint = nextPrint + dtPrint
  print "Simulation OK!"
  print h.tstop

  #times = numpy.array(h.tvecList)
  vSoma = numpy.array(h.vSomaList)
  spikes = numpy.array(h.spikes)
  spikedCells = numpy.array(h.spikedCells)

  picklelist = [spikes,spikedCells,vSoma,DIPOLE_MOMENT]
  file = open('spikes_parallel_'+myMechToAdd+str(nseg)+'_'+str(tstop)+'_'+str(mutID)+'_'+str(Econ)+'_'+str(Icon)+'_'+str(rateCoeff)+'_'+str(gNoiseCoeff)+'_'+str(gSynCoeff)+'_'+str(rdSeed)+'_'+str(rank)+'_of_'+str(Nmc)+'_withLFP.sav', 'w')
  pickle.dump(picklelist,file)
  file.close()
  print 'Saved to spikes_parallel_'+myMechToAdd+str(nseg)+'_'+str(tstop)+'_'+str(mutID)+'_'+str(Econ)+'_'+str(Icon)+'_'+str(rateCoeff)+'_'+str(gNoiseCoeff)+'_'+str(gSynCoeff)+'_'+str(rdSeed)+'_'+str(rank)+'_of_'+str(Nmc)+'_withLFP.sav'
  caSoma = numpy.array(h.caSomaList)
  skSoma = numpy.array(h.skSomaList)
  cahvaSoma = numpy.array(h.cahvaSomaList)
  calvaSoma = numpy.array(h.calvaSomaList)
  natSoma = numpy.array(h.natSomaList)
  napSoma = numpy.array(h.napSomaList)
  ihSoma = numpy.array(h.ihSomaList)
  kv31Soma = numpy.array(h.kv31SomaList)
  ktSoma = numpy.array(h.ktSomaList)
  kpSoma = numpy.array(h.kpSomaList)
  Is = numpy.array(h.IList)
  if len(caSoma) > 0:
    sparseTimes = [sparsedt*x for x in range(0,len(caSoma[0]))]
  else:
    sparseTimes = []


  if len(vSoma) > 0:
    times = [dt*x for x in range(0,len(vSoma[0]))]
    if len(vSoma)==1:
      spikes = mytools.spike_times(times,vSoma.T)
    else:
      spikeTimes = []
      spikeNeurs = []
      for i in range(0,len(vSoma)):
        theseSpikes = mytools.spike_times(times,vSoma[i])
        spikeTimes = spikeTimes + theseSpikes
        spikeNeurs = spikeNeurs + [i for x in theseSpikes]
      spikes = [spikeTimes, spikeNeurs]
  else:
    times = []


  h("""
{epnm.pc.runworker()}
{epnm.pc.done()}
""")

  return [times,vSoma,spikes]



Loading data, please wait...