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

 Download zip file 
Help downloading and running models
"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. ..."
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]
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 L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell;
Channel(s): I Na,t; I T low threshold; I K; I h;
Gap Junctions:
Receptor(s): GabaA; GabaB; AMPA; NMDA;
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Touch;
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; GabaA; GabaB; AMPA; NMDA; I Na,t; I T low threshold; I K; I h;
objref pc, nil
pc = new ParallelContext()
{load_file("stdlib.hoc") } // need String

// map functions between gid and (type, iXD, jYD)
// original serial program creates
//objref PL5[X_DIM][Y_DIM]
//objref PL2[X_DIM][Y_DIM]
//objref IPL2[X_DIM][Y_DIM]
//objref IPL5[X_DIM][Y_DIM]
//Although some of the IPL are empty because of zigzag there is no reason
//to have contiguous gids so we start with a simple map.

// type ranges from 0 to 3 and refers to PL5, PL2, IPL2, IPL5
gid_end = N_TYPE*X_DIM*Y_DIM // NOT the number of cells (because of zigzag)
PL5_type = 0
PL2_type = 1
IPL2_type = 2
IPL5_type = 3

//gid = type2gid(type, ixd, jyd)
func type2gid() {local type, ixd, jyd, gid
	type = $1  ixd = $2  jyd = $3
	gid = type * X_DIM * Y_DIM  +  ixd * Y_DIM + jyd
	return gid

//type = gid2type(gid, &ixd, &jyd)
func gid2type() { local gid, x, type, ixd, jyd
	gid = $1
	jyd = gid%Y_DIM
	x = (gid - jyd)/Y_DIM
	ixd = x%X_DIM
	type = (x - ixd)/X_DIM
	$&2 = ixd
	$&3 = jyd
	return type

// only if cell exists
//1or0 = cell_exist(type, ixd, jyd) {
func cell_exist() {
	return (pc.gid_exists(type2gid($1, $2, $3)) > 1)

// round robin iterator
// for pcitr(&gid) { ... }
iterator pcitr() {local gid
	for (gid =; gid < gid_end; gid += pc.nhost) {
		$&1 = gid
		// will generally have to check if the gid is really associated
		// with a cell by checking the return value of pc.gid_exists(gid)

// gids owned by this host (because of zigzag, some the gids for
// inhibitory cells will not be associated with cells.
proc declare_gids() {local gid
	for pcitr(&gid) {

//cell = cell_create("create statement", type, iXD, jYD)
// assumes gid exists
objref tmpcell
obfunc create_cell() { local gid   localobj cell, nc, nil, s
	gid = type2gid($2, $3, $4)
	if (pc.gid_exists(gid)) {
		s = new String()
		sprint(s.s, "tmpcell = %s", $s1)	
		tmpcell.connect2target(nil, nc)
		pc.cell(gid, nc)
		cell = tmpcell
		tmpcell = nil
	return cell

objref spikevec, idvec
proc want_all_spikes() {local gid
	spikevec = new Vector(10000)
	idvec = new Vector(10000)
	// real cells
	for pcitr(&gid) {
		if (pc.gid_exists(gid) > 1) {
			pc.spike_record(gid, spikevec, idvec)
	// AlphaFF and AlphaFB --- see MuBurst_10.hoc
	for (gid =; gid < gidAlphaend; gid += pc.nhost) {
		pc.spike_record(gid, spikevec, idvec)
	// FF, FB, FF2
	if ( == 0) {
		pc.spike_record(FFgid, spikevec, idvec)
		pc.spike_record(FBgid, spikevec, idvec)
		pc.spike_record(FF2gid, spikevec, idvec)

proc spikeout() {local rank, i  localobj f
	f = new File("out.dat")
	if ( == 0) { f.wopen()  f.close() }
	for rank = 0, pc.nhost-1 {
		if (rank == {
			for i=0, idvec.size-1 {
				f.printf("%.3f %d\n", spikevec.x[i], idvec.x[i])

total_process_dipole = 0
objref dipole_record, t_record
dipole_record = new Vector()
t_record = new Vector()

proc calc_total_dipole() {
	// need to do a reduce sum on all the dipole_record
	pc.allreduce(dipole_record, 1)
proc print_total_dipole() {local i
	if ( == 0) {
		for i=0, t_record.size-1 {
			fprint("%f %f  \n", t_record.x[i], dipole_record.x[i])

proc psolve() {

Loading data, please wait...