Deconstruction of cortical evoked potentials generated by subthalamic DBS (Kumaravelu et al 2018)

 Download zip file 
Help downloading and running models
"... High frequency deep brain stimulation (DBS) of the subthalamic nucleus (STN) suppresses parkinsonian motor symptoms and modulates cortical activity. ... Cortical evoked potentials (cEP) generated by STN DBS reflect the response of cortex to subcortical stimulation, and the goal was to determine the neural origin of cEP using a two-step approach. First, we recorded cEP over ipsilateral primary motor cortex during different frequencies of STN DBS in awake healthy and unilateral 6-OHDA lesioned parkinsonian rats. Second, we used a biophysically-based model of the thalamocortical network to deconstruct the neural origin of the cEP. The in vivo cEP included short (R1), intermediate (R2) and long-latency (R3) responses. Model-based cortical responses to simulated STN DBS matched remarkably well the in vivo responses. R1 was generated by antidromic activation of layer 5 pyramidal neurons, while recurrent activation of layer 5 pyramidal neurons via excitatory axon collaterals reproduced R2. R3 was generated by polysynaptic activation of layer 2/3 pyramidal neurons via the cortico-thalamic-cortical pathway. Antidromic activation of the hyperdirect pathway and subsequent intracortical and cortico-thalamo-cortical synaptic interactions were sufficient to generate cEP by STN DBS, and orthodromic activation through basal ganglia-thalamus-cortex pathways was not required. These results demonstrate the utility of cEP to determine the neural elements activated by STN DBS that might modulate cortical activity and contribute to the suppression of parkinsonian symptoms."
1 . Kumaravelu K, Oza CS, Behrend CE, Grill WM (2018) Model-based deconstruction of cortical evoked potentials generated by subthalamic nucleus deep brain stimulation. J Neurophysiol 120:662-680 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex; Thalamus;
Cell Type(s): Neocortex M1 L6 pyramidal corticothalamic GLU cell; Neocortex M1 L5B pyramidal pyramidal tract GLU cell; Neocortex M1 L4 stellate GLU cell; Hodgkin-Huxley neuron; Neocortex layer 4 neuron; Neocortex fast spiking (FS) interneuron; Neocortex primary motor area pyramidal layer 5 corticospinal cell;
Channel(s): I Na,p; I K; I Sodium; I_KD; I Calcium; I T low threshold; I L high threshold; I_AHP;
Gap Junctions: Gap junctions;
Receptor(s): AMPA; Gaba; NMDA;
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Deep brain stimulation; Evoked LFP;
Implementer(s): Kumaravelu, Karthik [kk192 at];
Search NeuronDB for information about:  Neocortex M1 L6 pyramidal corticothalamic GLU cell; Neocortex M1 L5B pyramidal pyramidal tract GLU cell; Neocortex M1 L4 stellate GLU cell; AMPA; NMDA; Gaba; I Na,p; I L high threshold; I T low threshold; I K; I Sodium; I Calcium; I_AHP; I_KD; Gaba; Glutamate;
balcomp.hoc *
defvar.hoc *
lbcreate.hoc *
mscreate.hoc *
parlib2.hoc *
traubcon.hoc *
traubcon_net.hoc *
if (name_declared("use_traubexact") == 0) {
	use_traubexact = 0
begintemplate TraubExactInfo
public tci, traub_parent
objref tci[14], traub_parent[14]
endtemplate TraubExactInfo
objref traubExactInfo
if (use_traubexact) {	// for network context traub exact connections
	traubExactInfo = new TraubExactInfo()

objref tci, traub_parent
strdef tstr
proc traub_connect() {local x, ip, ic  localobj child
	if (numarg() == 1) {
		tci = new Matrix($1, $1, 2)
		traub_parent = new Vector($1)
	if (use_traubexact) {
		if ([$o1.type] == nil) {[$o1.type] = tci
			traubExactInfo.traub_parent[$o1.type] = traub_parent
	if (numarg() == 5) {x = $5} else { x = 1 }
	ip = $2  ic = $3
	// which is the child and which is the parent?
	// assume a child is connected to the parent before it
	// makes a delta connection (to child connected to same parent)
	// if $3 is already connected assume the child is $2
	$o1.comp[$3] child = new SectionRef()
	if (child.has_parent) {
		$o1.comp[$2] child = new SectionRef()
		ip = $3  ic = $2
	// standard NEURON connection only if child not already connected
	if (!child.has_parent) {
		$o1.comp[ip] connect child.sec(0), x
		traub_parent.x[ic] = ip
	// save info for optional traub style connections
	tci.x[$2][$3] = $4
	tci.x[$3][$2] = $4
	// note that the above matrix has the topology of a tree unless
	// child to child connections form delta loops.

// exact() transforms the delta and pyramidal axial resistance branch points to Y's and X's.
// single_cell points at current cell object, .e.g TCR[0]
proc exact(){ 
	print "exact() has been called."
	execute("traubexact(single_cell, tci)")

proc traubexact() {
	//$o1 is the cell, $o2 is the connection matrix
	// assume the names of sections range from
	// $o1.comp[1] to $o1.comp[$o2.nrow-1]
	// the connection matrix

	// the traub connection rule is that a childs 0 end is connected
	// to the .5 location of the parent and the resistance is defined
	// by the connection matrix (uS) with one exception. Sometimes
	// pairs or triples of children are also connected to each other
	// and these children are moved to the 1 end of the parent and
	// the delta network is replaced by the equivalent wye network.

	traubexact_topology($o1, $o2, traub_parent)
	traubexact_coef($o1, $o2, traub_parent)

proc traubexact_topology() { local j, jc, ip, ic, gc, x
	// move all children to .5 location of parent except if
	// a group of children have a delta topology in which case move
	// them to the 1 location of the parent.
	for ic = 2, $o2.nrow-1 $o1.comp[ic] {
		// parent must be the least index
		ip = $o3.x[ic]
		x = .5
		for j = 0, $o2.sprowlen(ic)-1 {
			gc = $o2.spgetrowval(ic, j, &jc)
			// the question is whether jc is also connected
			// to ic's parent.
			if ($o2.getval(ip, jc)) { // yes
				x = 1
		if (x != parent_connection()) {
			$o1.comp[ip] connect $o1.comp[ic](0), x

proc traubexact_coef() { local i, j, id, ip, ic, gp, gc, jc, x localobj dtopol, dmat, dindx, wye
	// now that the topology is stable we can safely change the
	// connection coefficients so they reflect the desired axial
	// resistance
	dtopol = new Vector($o2.nrow)
	for ic = 2, $o2.nrow-1 $o1.comp[ic] {
		ip = $o3.x[ic]
		gp = $o2.getval(ip,ic)
		if (parent_connection() == 0.5) {
			rold = ri(0.5)
			scale_connection_coef(0.5, rold*gp)
		}else{ // this is delta-topology child.
			// so save info for further analysis
			dtopol.x[ic] = ip
	// examine the dtopol info. For each delta topology, construct a
	// conductance matrix and a section index vector (the parent is always
	// the first index) and then find the equivalent wye conductances and
	// use those to specify the connection coefficients
	for (id = dtopol.indwhere(">",0); id != -1; id = dtopol.indwhere(">",0)){
		// a delta exists and the parent is
		ip = dtopol.x[id]
		// note: error if ip != $o3.x[id]
		// so all the children are
		dindx = dtopol.c.indvwhere("==", ip)
		// prepare for next iteration over id
		for i=0, dindx.size-1 { dtopol.x[dindx.x[i]] = 0 }
		// and the complete delta index vector is
		dindx.insrt(0, ip)		
		// the delta conductance matrix
		dmat = new Matrix(dindx.size, dindx.size)
		if (dindx.size < 3) {
			execerror("not a delta", 0)
		for i = 0, dindx.size-1 {
			ic = dindx.x[i]
			for j = i, dindx.size-1 {
				jc = dindx.x[j]
				dmat.x[i][j] = $o2.getval(ic,jc)
				dmat.x[j][i] = dmat.x[i][j]
		// transform into an equivalent wye
		wye = delta2wye(dmat)
		// modify the connection coefficients
		for i = 0, dindx.size-1 $o1.comp[dindx.x[i]] {
			if (i == 0) { x = 1 } else { x = .5 }
			rold = ri(x)
			scale_connection_coef(x, rold/wye.x[i])

obfunc delta2wye() {localobj g, x, row0
	// $o1 is the conductance matrix. return vector is the equivalent
	// wye vector of resistances.
	// construct a full Gij*(Vi - Vj) = Ii matrix equation
	x = new Vector($o1.nrow)
	$o1.mulv(x, x)
	g = new Matrix($o1.nrow, $o1.ncol)
	g.setdiag(0, x)
	//presently indeterminate. save first row and replace by ground equation
	row0 = g.getrow(0)
	g.setrow(0, 0)
	g.x[0][0] = 1
	$o1.getdiag(0, x) // 1 to n-1 are the r0 + ri resistances and x.x[0]=1
	// need one more equation (not involving rp) and assume
	// numerical accuracy allows an arbitrary pair of nodes for
	// current injection and ground.
	// so put the row back, select node 1 as the ground point,
	// and inject into node 2.
	g.setrow(0, row0)
	g.setrow(1, 0)
	g.x[1][1] = 1
	//wasteful but...
	x.x[0] = $o1.x[2][2]
	// now the 0 element is r1+r2
	// set up the matrix equation to solve for ri
	g.setcol(0, 1)
	g.setdiag(0, 1)
	// but the first equation is r1 + r2 = x0
	g.x[0][0] = 0
	g.x[0][1] = 1
	g.x[0][2] = 1
	x = g.solv(x)

	return x

Loading data, please wait...