Biophysically realistic neural modeling of the MEG mu rhythm (Jones et al. 2009)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:136803
"Variations in cortical oscillations in the alpha (7–14 Hz) and beta (15–29 Hz) range have been correlated with attention, working memory, and stimulus detection. The mu rhythm recorded with magnetoencephalography (MEG) is a prominent oscillation generated by Rolandic cortex containing alpha and beta bands. Despite its prominence, the neural mechanisms regulating mu are unknown. We characterized the ongoing MEG mu rhythm from a localized source in the finger representation of primary somatosensory (SI) cortex. Subjects showed variation in the relative expression of mu-alpha or mu-beta, which were nonoverlapping for roughly 50% of their respective durations on single trials. To delineate the origins of this rhythm, a biophysically principled computational neural model of SI was developed, with distinct laminae, inhibitory and excitatory neurons, and feedforward (FF, representative of lemniscal thalamic drive) and feedback (FB, representative of higher-order cortical drive or input from nonlemniscal thalamic nuclei) inputs defined by the laminar location of their postsynaptic effects. ..."
Reference:
1 . Jones SR, Pritchett DL, Sikora MA, Stufflebeam SM, Hämäläinen M, Moore CI (2009) Quantitative analysis and biophysically realistic neural modeling of the MEG mu rhythm: rhythmogenesis and modulation of sensory-evoked responses. J Neurophysiol 102:3554-72 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex V1 L6 pyramidal corticothalamic cell; Neocortex V1 L2/6 pyramidal intratelencephalic cell;
Channel(s): I Na,t; I T low threshold; I K; I h;
Gap Junctions:
Receptor(s): GabaA; GabaB; AMPA; NMDA;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Touch;
Implementer(s):
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic cell; Neocortex V1 L2/6 pyramidal intratelencephalic cell; GabaA; GabaB; AMPA; NMDA; I Na,t; I T low threshold; I K; I h;
begintime = startsw()
// 2D init file
X_DIM = 10 // number of Pyramidal cells along the X dimension 
Y_DIM = 10 // number of Pyramidal cells along the Y dimension

N_TYPE = 4 // number of cell types
{load_file("parlib.hoc")} // create pc and all the parallel specific functions, etc.
declare_gids()
stim_delay_offset = 1

Tfactor=0.0002
Ttau=0.0
Hfactor=0.000001
Htau=0.003

FSx=4 // Coordinates for *all* feeds 
FSy=4

{load_file("sj10-cortex.hoc")}
{load_file("wiring_proc_2Dv2.hoc")}

/*
load_file("noise2D_V2.hoc") 
UnoiseV(-0.3,0.3)
UnoiseII(-0.3,0.3)
UnoiseIPL5(-0.3,0.3)
UnoiseIPL2(-0.3,0.3)
*/

// Intra-cortical wiring
{load_file("wiring-SmlFeed-3_7.hoc")}
//----------------------------------------
// Mu timing
{load_file("MuBurst_10.hoc")}
//----------------------------------------
// Mu weights
{load_file("E-FFFBx_fixed_10.hoc")}
//----------------------------------------
// Evoked signal below



// Redefinition of init()
// Read files containing segment rest voltages and state variables
// at steady-state. **PL2 not comprhensive**
objref f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16
objref PL5rest,PL2rest,MAR,HCAT,MCAT,HCA,MCA,PL5CAI
objref HHH,MHH,NHH,NKCA,NKM,PL5_DPL,PL2_DPL,PL5_Q
{
f1 = new File()
{f1.ropen("./STATES/PL5_Rest.dat")}
PL5rest = new Vector()
NumF1=PL5rest.scanf(f1,1,1)
f1.close()
f2 = new File()
{f2.ropen("./STATES/PL2_Rest.dat")}
PL2rest = new Vector()
NumF2=PL2rest.scanf(f2,1,1)
f2.close()
f3 = new File()
{f3.ropen("./STATES/M_AR.dat")}
MAR = new Vector()
NumF3=MAR.scanf(f3,1,1)
f3.close()
f4 = new File()
{f4.ropen("./STATES/H_CAT.dat")}
HCAT = new Vector()
NumF4=HCAT.scanf(f4,1,1)
f4.close()
f5 = new File()
{f5.ropen("./STATES/M_CAT.dat")}
MCAT = new Vector()
NumF5=MCAT.scanf(f5,1,1)
f5.close()
f6 = new File()
{f6.ropen("./STATES/H_CA.dat")}
HCA = new Vector()
NumF6=HCA.scanf(f6,1,1)
f6.close()
f7 = new File()
{f7.ropen("./STATES/M_CA.dat")}
MCA = new Vector()
NumF7=MCA.scanf(f7,1,1)
f7.close()
f8 = new File()
{f8.ropen("./STATES/H_HH.dat")}
HHH = new Vector()
NumF8=HHH.scanf(f8,1,1)
f8.close()
f9 = new File()
{f9.ropen("./STATES/M_HH.dat")}
MHH = new Vector()
NumF9=MHH.scanf(f9,1,1)
f9.close()
f10 = new File()
{f10.ropen("./STATES/N_HH.dat")}
NHH = new Vector()
NumF10=NHH.scanf(f10,1,1)
f10.close()
f11 = new File()
{f11.ropen("./STATES/N_KCA.dat")}
NKCA = new Vector()
NumF11=NKCA.scanf(f11,1,1)
f11.close()
f12 = new File()
{f12.ropen("./STATES/N_KM.dat")}
NKM = new Vector()
NumF12=NKM.scanf(f12,1,1)
f12.close()
f13 = new File()
{f13.ropen("./STATES/PL5_Dipole.dat")}
PL5_DPL = new Vector()
NumF13=PL5_DPL.scanf(f13,1,1)
f13.close()
f14 = new File()
{f14.ropen("./STATES/PL5_CAI.dat")}
PL5CAI = new Vector()
NumF14=PL5CAI.scanf(f14,1,1)
f14.close()
f15 = new File()
{f15.ropen("./STATES/PL2_Dipole.dat")}
PL2_DPL = new Vector()
NumF15=PL2_DPL.scanf(f15,1,1)
f15.close()
//
f16 = new File()
{f16.ropen("./STATES/PL5_Q.dat")}
PL5_Q = new Vector()
NumF16=PL5_Q.scanf(f16,1,1)
f16.close()
}
// Redefine init() to use those voltages & state variables

proc init(){ local i,j,Vnum,v2num,DPnum
finitialize(-64.9737) // initialize all cells to Inhibitory rest potential
//

// Initialize Layer V Pyramidal voltage and states
for i=0,XD{
for j=0,YD{
 v2num=0 Vnum=0 
	 if (cell_exist(PL5_type, i, j)) {
 PL5[i][j].soma {v=-70.197014 m_ar=0.29453859 h_cat=0.079413855 m_cat=0.091981604\
    h_hh=0.75925494 m_hh=0.028240063 n_hh=0.24191157 n_km=0.011362846 h_ca=0.63961253\
    m_ca=3.8876755e-05 n_kca=4.9998916e-05 cai=0.00010000285}
 PL5[i][j].soma distance()
 forsec PL5[i][j].dendritic for(x,0){v=PL5rest.x(Vnum) m_ar=MAR.x(Vnum)\
  h_cat=HCAT.x(Vnum) m_cat=MCAT.x(Vnum) h_ca=HCA.x(Vnum) m_ca=MCA.x(Vnum)\
  h_hh=HHH.x(Vnum) m_hh=MHH.x(Vnum) n_hh=NHH.x(Vnum) n_kca=NKCA.x(Vnum)\
  n_km=NKM.x(Vnum) Qsum_dipole=PL5_DPL.x(Vnum) Q_dipole=PL5_Q.x(Vnum) cai=PL5CAI.x(Vnum)
  Vnum=Vnum+1}
//
// dend[5] not properly set by above ???
  PL5[i][j].dend[5].v = -70.22481
	}
//
// Initialize Layer II Pyramidal 
	if (cell_exist(PL2_type, i, j)) {
 PL2[i][j].soma {v=-71.462514 n_km=0.0098896622}
   PL2[i][j].soma distance()
   forsec PL2[i][j].dendritic for(x,0){v=PL2rest.x(v2num) Qsum_dipole=PL2_DPL.x(v2num)
    v2num=v2num+1}
	}
} }
// if (cvode.active()) {cvode.re_init() }else {fcurrent()}
fcurrent()
frecord_init()
}
/////////////////////////////////////////////////////////////
// Define the evoked signal with initial start times
FFgid = 1+gidAlphaend
FBgid = 2+gidAlphaend
FF2gid = 3+gidAlphaend
objref  FF, FB, FF2
if (pc.id == 0) {
FF = new FeedX() // Thalamic input (Feed-forward) 
FF.pp.number=1
FF.pp.start=454 - stim_delay_offset
FB = new FeedX() // Feed-Back (eg. Pre-frontal input) 
FB.pp.number=1
FB.pp.start=499 - stim_delay_offset
FF2 = new FeedX() // Thalamic input (Feed-forward) 
FF2.pp.number=1
FF2.pp.start=564 - stim_delay_offset
pc.set_gid2node(FFgid, pc.id)
pc.set_gid2node(FBgid, pc.id)
pc.set_gid2node(FF2gid, pc.id)
pc.cell(FFgid, new NetCon(FF.pp, nil))
pc.cell(FBgid, new NetCon(FB.pp, nil))
pc.cell(FF2gid, new NetCon(FF2.pp, nil))
}
// and connect it
{load_file("scale_ep_thresh.hoc")}
///////////////////////////////////////////////////////////////
//Total time for one run
tstop=1500
// time increment
dt=0.05
// Number of runs
NRUN = 1 // orig 26
proc runonce(){local runtime
//
	runtime = startsw()
	if (0) {
init()
while (t<tstop){
          fadvance()

// below put a '%f' for every data type defined
          fprint( "%f %f  \n", \
 t, \
 L2_dipole() + L5_dipole() ) // end print statement
}
	}else{
		psolve(tstop)
		// end of run print the total dipole
		calc_total_dipole()
		if (pc.id == 0) {
			print_total_dipole()
		}
	}
	if (pc.id == 0) printf("runtime = %g\n", startsw() - runtime)
} //end runonce()

make_MuBursts(400)
e_fffbx(0.00004,0.00004,5,1)

proc runit(){ local i
if (pc.id == 0) wopen("Mu_output.dat")
for i=0, NRUN-1{ runonce() //increments the evoked response starttime each run
if (pc.id == 0) {
FF.pp.start = FF.pp.start + 3.85
FB.pp.start = FB.pp.start + 3.85
FF2.pp.start = FF2.pp.start + 3.85
}
if (pc.id == 0) print(i)}
if (pc.id == 0) wopen()
}

setuptime = startsw() - begintime
if (pc.id == 0) printf("setuptime =  %g\n", setuptime)
//

want_all_spikes()

runit()

spikeout()

finishtime = startsw() - begintime
if (pc.id == 0) printf("finishtime = %g\n", finishtime)

if (pc.nhost > 1) {
	pc.runworker()
	pc.done()
	quit()
}