Multitarget pharmacology for Dystonia in M1 (Neymotin et al 2016)

 Download zip file 
Help downloading and running models
Accession:189154
" ... We developed a multiscale model of primary motor cortex, ranging from molecular, up to cellular, and network levels, containing 1715 compartmental model neurons with multiple ion channels and intracellular molecular dynamics. We wired the model based on electrophysiological data obtained from mouse motor cortex circuit mapping experiments. We used the model to reproduce patterns of heightened activity seen in dystonia by applying independent random variations in parameters to identify pathological parameter sets. ..."
Reference:
1 . Neymotin SA, Dura-Bernal S, Lakatos P, Sanger TD, Lytton WW (2016) Multitarget Multiscale Simulation for Pharmacological Treatment of Dystonia in Motor Cortex. Front Pharmacol 7:157 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Molecular Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex V1 L6 pyramidal corticothalamic cell; Neocortex U1 L2/6 pyramidal intratelencephalic cell; Neocortex V1 interneuron basket PV cell; Neocortex fast spiking (FS) interneuron; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron; Neocortex layer 4 neuron; Neocortex layer 2-3 interneuron; Neocortex layer 4 interneuron; Neocortex layer 5 interneuron; Neocortex layer 6a interneuron;
Channel(s): I A; I h; I_SERCA; Ca pump; I K,Ca; I Calcium; I L high threshold; I T low threshold; I N; I_KD; I M; I Na,t;
Gap Junctions:
Receptor(s): GabaA; GabaB; AMPA; mGluR;
Gene(s): HCN1;
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON; Python;
Model Concept(s): Oscillations; Activity Patterns; Beta oscillations; Reaction-diffusion; Calcium dynamics; Pathophysiology; Multiscale;
Implementer(s): Neymotin, Sam [samn at neurosim.downstate.edu]; Dura-Bernal, Salvador [salvadordura at gmail.com];
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic cell; Neocortex V1 interneuron basket PV cell; Neocortex U1 L2/6 pyramidal intratelencephalic cell; GabaA; GabaB; AMPA; mGluR; I Na,t; I L high threshold; I N; I T low threshold; I A; I M; I h; I K,Ca; I Calcium; I_SERCA; I_KD; Ca pump; Gaba; Glutamate;
/
dystdemo
readme.txt
cagk.mod
cal.mod *
calts.mod *
can.mod *
cat.mod *
gabab.mod
h_winograd.mod
HCN1.mod
IC.mod *
icalts.mod *
ihlts.mod *
kap.mod
kcalts.mod *
kdmc.mod
kdr.mod
km.mod *
mglur.mod *
misc.mod *
MyExp2SynBB.mod *
MyExp2SynNMDABB.mod
nax.mod
stats.mod *
vecst.mod *
aux_fun.inc *
conf.py
declist.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
geom.py
ghk.inc *
grvec.hoc
init.hoc
labels.hoc
labels.py *
local.hoc *
misc.h
mpisim.py
netcfg.cfg
nqs.hoc *
nqs.py
nrnoc.hoc *
pyinit.py *
python.hoc *
pywrap.hoc *
simctrl.hoc *
simdat.py
syn.py
syncode.hoc *
vector.py *
xgetargs.hoc *
                            
import sys
print sys.path
import os
import string
from neuron import *
from datetime import datetime
h("strdef simname, allfiles, simfiles, output_file, datestr, uname, osname, comment")
h.simname=simname = "dystdemo"
h.allfiles=allfiles = "pyinit.py geom.py mpisim.py"
h.simfiles=simfiles = "pyinit.py geom.py mpisim.py"
h("runnum=1")
runnum = 1.0
h.datestr=datestr = "16jun21"
h.output_file=output_file = "data/16jun21.0"
h.uname=uname = "x86_64"
h.osname=osname="linux"
h("templates_loaded=0")
templates_loaded=0
h("xwindows=1.0")
xwindows = 1.0
h.xopen("nrnoc.hoc")
h.xopen("init.hoc")
CTYPi = 60.0; STYPi = 18.0
from pyinit import *
from neuron import h, gui
from vector import *
from nqs import *
from labels import *
import random
from pylab import *
import ConfigParser

tl = tight_layout
ion()

rcParams['lines.markersize'] = 15
rcParams['lines.linewidth'] = 4
rcParams['font.size'] = 25

defbinsz=10
stimdel = 500 # delay from stim to measured firing rate peak response

defnCPU=8
print 'Using ' , defnCPU, ' cpus by default.'

config = ConfigParser.ConfigParser()

# determine config file name
def setfcfg ():
  fcfg = "physiol.cfg" # default config file name
  for i in xrange(len(sys.argv)):
    if sys.argv[i].endswith(".cfg") and os.path.exists(sys.argv[i]):
      fcfg = sys.argv[i]
  print "config file is " , fcfg
  return fcfg

fcfg=setfcfg() # config file name

config.read(fcfg)

def conffloat (base,var): return float(config.get(base,var))
def confint (base,var): return int(config.get(base,var))
def confstr (base,var): return config.get(base,var)

tstart = 0
tstop = conffloat('run','tstop') + tstart
binsz = conffloat("run","binsz") # bin size (in milliseconds) for MUA
recdt = conffloat('run','recdt')
recvdt = conffloat('run','recvdt')
vtslow=Vector(); vtslow.indgen(tstart,tstop-recdt,recdt); vtslow=numpy.array(vtslow.to_python())
vtfast=Vector(); vtfast.indgen(tstart,tstop-recvdt,recvdt); vtfast=numpy.array(vtfast.to_python())
simstr = config.get('run','simstr')

def makefname (simstr,pcid): return 'data/' + simstr + '_pc_' + str(pcid) + '.npz'

ix,ixe=[1e9 for i in xrange(CTYPi)],[-1e9 for i in xrange(CTYPi)]
numc=[0 for i in xrange(CTYPi)]
allcells,icells,ecells=0,0,0

# reload the variables from config file
def reloadvars ():
  global tstart,loadstate,tstop,binsz,recdt,recvdt,vtslow,vtfast,simstr,numc,allcells,ecells,icells,startt
  tstart = conffloat('run','loadtstop')
  loadstate = confint('run','loadstate')
  if loadstate == 0: tstart = 0 # only use previous end time if loading state
  tstop = conffloat('run','tstop') + tstart
  binsz = conffloat("run","binsz") # bin size (in milliseconds) for MUA
  recdt = conffloat('run','recdt')
  recvdt = conffloat('run','recvdt')
  vtslow=Vector(); vtslow.indgen(tstart,tstop-recdt,recdt); vtslow=numpy.array(vtslow.to_python())
  vtfast=Vector(); vtfast.indgen(tstart,tstop-recvdt,recvdt); vtfast=numpy.array(vtfast.to_python())
  simstr = config.get('run','simstr')
  numc=[0 for i in xrange(CTYPi)]
  allcells,ecells,icells=0,0,0

# setup start/end indices for cell types
def makeix (lctyID):
  global ix,ixe,allcells,ecells,icells
  allcells,ecells,icells=0,0,0
  for i in xrange(CTYPi):
    ix[i]=1e9
    ixe[i]=-1e9
    numc[i]=0
  for i in xrange(len(lctyID)):
    ty = lctyID[i]
    numc[ty]+=1
    allcells+=1
    ix[ty] = min(ix[ty],i)
    ixe[ty] = max(ixe[ty],i)
    if h.ice(ty): icells+=1
    else: ecells+=1

lfastvar = ['volt', 'iAM', 'iNM', 'iGB', 'iGA','ina','ik','ica','ih']

# reads data from files saved by mpisim.py - must get nhost correct
def readdat (simstr,nhost):
  ld = {} # dict for data from a single host
  for pcid in xrange(nhost):
    fn = makefname(simstr,pcid)
    ld[pcid]=numpy.load(fn)
  return ld
  
# load/concat data from mpisim.py (which saves data from each host separately) - must get nhost correct
def loaddat (simstr,nhost,quiet=False):
  ld = readdat(simstr,nhost) # dict for data from a single host
  lids,lspks=array([]),array([]) # stitch together spike times
  for pcid in xrange(nhost):
    lids=concatenate((lids,ld[pcid]['lids']))
    lspks=concatenate((lspks,ld[pcid]['lspks']))
  ldout={} # concatenated output from diff hosts
  ldout['lids']=lids; ldout['lspks']=lspks;
  for k in ['lctyID','vt','lX', 'lY', 'lZ']: ldout[k]=ld[0][k][:]
  makeix(ldout['lctyID'])
  ncells = len(ld[0]['lctyID'])
  vt = ld[0]['vt'][:]
  locvarlist = []
  for loc in ['soma','Adend3']:
    for var in ['volt']:
      locvarlist.append(loc + '_' + var)
  locvarlist.append('Bdend_volt')
  for k in locvarlist:
    if not k in ld[0]: continue # skip?
    if not quiet: print k
    if k.endswith('volt') or lfastvar.count(k.split('_')[1])>0: sz = len(vtfast)
    else: sz = len(vtslow)
    ldout[k] = zeros((ncells,sz))
    ldk = ldout[k]
    for pcid in xrange(nhost):
      d = ld[pcid][k]
      lE = ld[pcid]['lE']; 
      for row,cdx in enumerate(lE): ldk[cdx,:] = d[row,0:sz]
      if k == 'soma_volt': # special case for soma voltage of I cells
        lI=ld[pcid]['lI']; dI=ld[pcid]['soma_voltI']; row=0
        for row,cdx in enumerate(lI): ldk[cdx,:] = dI[row,0:sz]
  for pcid in xrange(nhost):
    del ld[pcid].f
    ld[pcid].close()
    del ld[pcid]
  return ldout

# colors for spikes for raster -- based on cell type
def getcolors (lspks,lids,lctyID):
  lclrs = []
  for i in xrange(len(lids)):
    ty = lctyID[int(lids[i])]
    if IsLTS(ty):
      lclrs.append('b')
    elif ice(ty):
      lclrs.append('g')
    else:
      lclrs.append('r')
  return lclrs
			
#
def drawraster (ld,sz=2):
  lspks,lids,lctyID=ld['lspks']/1e3,ld['lids'],ld['lctyID']
  try:
    lclrs = ld['lclrs']
  except:    
    lclrs = getcolors(lspks,lids,lctyID)
    ld['lclrs'] = lclrs
  scatter(lspks,lids,s=sz**2,c=lclrs,marker='+')
  allcells = len(lctyID)
  xlim((tstart/1e3,tstop/1e3)); ylim((0,allcells)); tight_layout(); xlabel('Time (s)',fontsize=30); 
  ax = gca()
  ax.set_yticks([])

# get cells in vid if they're of specified type (in lty)
def getall (ix,ixe,lty=[E2,E4,E5R,E5B,E5P,E6,E6C]):
  vact = Vector()
  for ct in lty:
    vtmp = Vector(); vtmp.indgen(ix[ct],ixe[ct],1); vact.append(vtmp)
  return vact

# get an array of times of interest
def gettimes (tstart,tstop,binsz):
  vt = h.Vector()
  vt.indgen(tstart,tstop,binsz)
  return vt

# get an NQS with spikes (snq)
def getsnq (ld):
  lspks,lids,lctyID=ld['lspks'],ld['lids'],ld['lctyID']
  snq = NQS('id','t','ty','ice')
  snq.v[0].from_python(lids)
  snq.v[1].from_python(lspks)
  for i in xrange(len(lids)):
    snq.v[2].append(lctyID[int(lids[i])])
    snq.v[3].append(h.ice(lctyID[int(lids[i])]))
  ld['snq']=snq
  return snq

# print spike rates in each population in the periods of interest (baseline, signal, recall)
def prspkcount (snq,tstart,tstop,binsz,times=None,quiet=False):
  snq.verbose=0
  if times is None:
    times = gettimes(tstart,tstop,binsz) # periods of interest
  lfa=[]
  vact = getall(ix,ixe,lty=[E2,E4,E5R,E5B,E5P,E6,E6C])
  for startt in times:
    endt = startt + binsz
    if endt > tstop: break
    fa = round(1e3*snq.select("id","EQW",vact,"t","[]",startt,endt) / ((endt-startt)*vact.size()),2)
    lfa.append(fa)
    if not quiet:
      stro = "t=" + str(startt) + "-" + str(endt) + ". E:" + str(fa) + " Hz."
      for ct in [I2, I2L, I5, I5L, I6, I6L]:
        fb = round(1e3*snq.select("id","[]",ix[ct],ixe[ct],"t","[]",startt,endt) / ((endt-startt)*numc[ct]),2)
        stro += ' ' + CTYP[ct] + ':' + str(fb)
      fi = round(1e3*snq.select("ice",1,"t","[]",startt,endt) / ((endt-startt)*icells),2)
      stro += ' I:' + str(fi)
      for ct in [E2, E4, E5R, E5B, E5P, E6, E6C]:
        fb = round(1e3*snq.select("id","[]",ix[ct],ixe[ct],"t","[]",startt,endt) / ((endt-startt)*numc[ct]),2)
        stro += ' ' + CTYP[ct] + ':' + str(fb)
      fe = round(1e3*snq.select("ice",0,"t","[]",startt,endt) / ((endt-startt)*ecells),2)
      stro += ' E:' + str(fe)
      print stro
  snq.verbose=1
  return lfa

# getfnq - make an NQS with ids, firing rates, types
def getfnq(snq,lctyID,skipms=200):
  snq.verbose=0; snq.tog("DB"); 
  fnq = h.NQS("id","freq","ty")
  tf = tstop - skipms # duration we're considering for frequency calc
  for i in xrange(allcells):
    n = float( snq.select("t",">",skipms,"id",i) ) # number of spikes
    fnq.append(i, n*1e3/tf, lctyID[i])
  snq.verbose=1
  return fnq

# pravgrates - print average firing rates using fnq
def pravgrates(snq,skipms=500):
  snq.verbose=0
  for ty in xrange(CTYPi):
    if numc[ty] < 1: continue
    print CTYP[ty], ' avg rate = ', getrate(snq,ty,500.0)[2]
  snq.verbose=1

###########################################################3

# reload variables, reset config file name
def myreload (Simstr,ncpu=defnCPU, quiet=False):
  global fcfg,simstr
  if not quiet: print 'simstr was ' , simstr
  simstr = Simstr
  print 'simstr: ' , simstr
  fcfg = 'backupcfg/' + simstr + '.cfg'
  if len(config.read(fcfg))==0:
    print 'myreload ERRA: could not read ', fcfg, '!!!'
    return None
  reloadvars()
  if not quiet: print fcfg
  try:
    ld=loaddat(simstr,ncpu,quiet);
    return ld
  except:
    print 'could not reload data from' , Simstr
    return None

ld=None
  
#
def mydrawrast (lld):
  for i,ld in enumerate(lld):
    figure(); 
    drawraster(ld,sz=4); title(lsimstr[i])
    xlabel('')
    xlim((0.5,2))  
    fsz= 20
    ax = gca(); 
    ax.set_yticks([]); ax.set_xticks([]); ax.set_xlabel('');
    tx=-0.03
    lty = [(ix[ct]+(ixe[ct]-ix[ct])*0.3)/allcells for ct in [E2,E4,E5R,E5B,E5P,E6C,E6]]
    for ty,lbl in zip(lty,['E2','E4','E5a','E5b','E5P','E6C','E6']): addtext(1,1,[1],[lbl],tx=tx,ty=ty,c='r',fsz=fsz+5)
    lty = [(ix[ct]+(ixe[ct]-ix[ct])*0.075)/allcells for ct in [I2,I5,I6]]
    for ty,lbl in zip(lty,['I2','I5','I6']): addtext(1,1,[1],[lbl],tx=tx,ty=ty,c='g',fsz=fsz)
    lty = [(ix[ct]+(ixe[ct]-ix[ct])*0.25)/allcells for ct in [I2L,I5L,I6L]]
    for ty,lbl in zip(lty,['I2L','I5L','I6L']): addtext(1,1,[1],[lbl],tx=tx,ty=ty,c='b',fsz=fsz)

#
def addtext (row,col,lgn,ltxt,tx=-0.025,ty=1.03,c='k',fsz=30):
  for gn,txt in zip(lgn,ltxt):
    if row == col and col == 1:
      ax = gca()
    else:
      ax = subplot(row,col,gn)
    text(tx,ty,txt,fontweight='bold',transform=ax.transAxes,fontsize=fsz,color=c);

#
def Gaddtext (G,row,col,txt,tx=-0.025,ty=1.03,c='k',fsz=30):
  ax = subplot(G[row,col])
  text(tx,ty,txt,fontweight='bold',transform=ax.transAxes,fontsize=fsz,color=c);

def naxbin (ax,nb): ax.locator_params(nbins=nb);

lsimstr = ['d2', 'dystonia', 'latchup']
lld,lsnq = [],[]
for simstr in lsimstr:
  ld = loaddat(simstr,defnCPU)
  lld.append(ld)
  snq = getsnq(ld)
  lsnq.append(snq)

mydrawrast(lld)


Loading data, please wait...