initA.hoc

/* Extracellular stimulation of FH myelinated axon model.
Bipolar electrodes applied to surface of a semi-infinite conductive medium.
*/

load_file("nrngui.hoc")

///// parameters

RHOE = 300 // extracellular resistivity in ohm cm

///// anatomical and biophysical properties of the axon
///// assuming extracellular medium is perfect conductor

load_file("axon10.hoc") // external diameter is 10 um

v_init = -70 // from ModelDB entry 3507
// actual resting potential is closer to -69.77 mV

///// stimulation

// steps:
// 0.  insert extracellular mechanism and specify its parameters
// 1.  set up transfer resistances
// 2.  set up stimulus waveform
// 3.  couple stim waveform to xstim

// 0.  insert extracellular mechanism and specify its parameters
//     when using extracellular to implement extracelluar stimulation,
//     use extracellular's xg and xc to play the role of myelin--see axon.hoc

forall insert extracellular
forsec internodes {
  for (x,0) {
    cm(x) = CM // since extracellular's default xc is 0; CM is defined in axon.hoc
    for i=0,1 xg(x) = 1e-9 // "effectively a perfect insulator"
  }
}

// 1.  set up transfer resistances

forall {
  insert xstim
  for (x,0) setpointer ex_xstim(x), e_extracellular(x)
}
is_xstim = 0 // for development & testing
             // eventually is driven by a forcing function

/*
Electrodes:  point sources on the surface of a semi-infinite conductive medium
so a stimulus current with value istim is twice as effective
as if the entire extracellular volume were conductive.
Axonal geometry coordinates:  axon lies along the x axis, 
with the 0 end of the axon at (0,0,0) and the electrode at (xe,ye,ze).
Electrode geometry coordinates used in "test cases":
axon lies along the x axis,
but (0,0,0) corresponds to the middle of the axon.
The (xa,ya) and (xc,yc) used in the "test cases"
correspond to (xa + axlen/2, ya, 0) and (xc + axlen/2, yc, 0).
where axlen == total length of the axon
*/

load_file("interpxyz.hoc") // defines proc grindaway()

grindaway(all) // find xyz coords of all internal nodes

// $1   rho in ohm cm
// $2-4 xyz coords of point A in um
// $5-7                     B in um
// $8   lower limit to distance in um
//        i.e. A and B can never be closer than local neurite radius
// returns transfer resistance in megohms

// rxm returns transfer resistance between an electrode and a location along the axon
// includes factor of 2 because electrode is on surface of medium

func rxm() { local rho,x1,y1,z1,x2,y2,z2,dmin,r
  rho = $1
  x1 = $2
  y1 = $3
  z1 = $4
  x2 = $5
  y2 = $6
  z2 = $7
  dmin = $8
  r = sqrt((x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2)   
  if (r<dmin) r=dmin
  // calculate the transfer resistance between the node and the grid point
  return 2*0.01*(rho / 4 / PI)*(1/r)
}

axlen = 0
forall axlen+=L

// bipolar stimulating electrodes
// the conductive medium is linear so the net effect of bipolar stimulation
// is the sum of the anodal and cathodal stimuli

XA = 50   // cm, must convert to um
YA = 0.25
XC = 0
YC = 0.25

proc calcrx() {
  forall for (x,0) rx_xstim(x) = rxm(RHOE, x_xstim(x), y_xstim(x), z_xstim(x), \
                                         XA*1e4 + axlen/2, YA*1e4, 0, diam(x)/2) \
                               - rxm(RHOE, x_xstim(x), y_xstim(x), z_xstim(x), \
                                         XC*1e4 + axlen/2, YC*1e4, 0, diam(x)/2)
}

calcrx()

// 2.  set up stimulus waveform
// and
// 3.  couple stim waveform to xstim

objref fsq
fsq = new Fsquare(0.5) // square wave generator
setpointer fsq.x, is_xstim

dummy = 0
objref fzap
fzap = new Fzap(0.5) // swept sine wave generator
                     // used here to produce a fixed frequency sine wave
setpointer fzap.x, dummy

///// graphical user interface

load_file("basicrig.ses") // RunControl
  // IClamp for direct intracellular stim at node[0](0.5) (for testing)
  // v vs. t
IClamp[0].amp = 0 // no intracellularly injected current

load_file("varstep.ses") // variable dt tool

load_file("vvsx.ses") // Movie Run, v vs. x

// additional graphs
// these have little effect on run time

load_file("vext_eext.ses") // vext and e_extracellular vs. distance along axon

// plot of rx vs distance along axon
// not updated during simulation!

objref xval, rxval, grx
xval = new Vector()
rxval = new Vector()
grx = new Graph(0)

proc plotrx() { localobj rvp
  grx = new Graph(0)
  grx.size(0,100000,0,0.0002)
  grx.view(0, 0, 100000, 0.0002, 327, 534, 300.48, 200.32)
  rvp = new RangeVarPlot("rx_xstim")
  node[0] rvp.begin(0)
  node[100] rvp.end(1)
  rvp.origin(0)
  grx.addobject(rvp)
  grx.exec_menu("View = plot")
}

plotrx()


///// automatic detection of threshold stimulus intensity

// spike detection
// for stimulation that triggers spike onset in middle of axon
// attach an APCount to the node at the proximal end of the axon
// this assumes spike initiation occurs far from node[0]
// so that stim artifact doesn't trigger spike detector

objref apc
node[0] apc = new APCount(0.5) // -20 mV is default thresh for spike detection

stimamp = 0

load_file("thresh4.hoc") // determines spike threshold to 4 significant figures

// This function requires an existing APCount[0] at the user desired location
// returns 1 if the voltage passed the APCount[0].thresh
func thresh_excited() {
  if (waveform==PULSE) {
    fsq.amp1 = stimamp
    fsq.amp2 = 0
  }
  if (waveform==SQUARE) {
    fsq.amp1 = stimamp
    fsq.amp2 = -stimamp
  }
  if (waveform==SINE) {
    fzap.amp = stimamp
  }
	run()
	return APCount[0].n > 0
}

// print "spike threshold is ", threshold(&stimamp), "mA"

///// define stimulus protocol and find corresponding threshold

load_file("protocolsA.hoc")

proc doprotocols() { local i, thresh, t0
  t0 = startsw()
  printf("test \twform \ttp \tnc \txa,ya \t \txc,yc \t\t   thresh\n")  
  for i=$1,$2 if (protocol(i)) {
    if (waveform==PULSE) {
      setpointer fsq.x, is_xstim
      fsq.del = 1
      fsq.num = NC
      fsq.dp = TP // update pulse or square wave stimulus waveform

      setpointer fzap.x, dummy
      fzap.del = 0
      fzap.dur = 0
    }
    if (waveform==SQUARE) {
      setpointer fsq.x, is_xstim
      fsq.del = 1
      fsq.dp = TP // update pulse or square wave stimulus waveform
      fsq.num = NC

      setpointer fzap.x, dummy
      fzap.del = 0
      fzap.dur = 0
    }
    if (waveform==SINE) {
      setpointer fsq.x, dummy
      fsq.del = 0
      fsq.dp = 0
      fsq.num = 0
      fsq.amp1 = 0
      fsq.amp2 = 0

      setpointer fzap.x, is_xstim
      fzap.del = 1
      fzap.dur = TP*2*NC
      fzap.f0 = 1000/2/TP
      fzap.f1 = fzap.f0
    }

    calcrx() // update rx, in case the new protocol changed electrode locations
    plotrx()
    stimamp = 0
    thresh = threshold(&stimamp)
    printf("%d \t%s \t%5.3f \t%d \t%5.2f,%5.2f \t%5.2f,%5.2f \t%11.5f\n", \
            i, wstr, TP, NC, XA, YA, XC, YC, thresh)
  }
  print " "
  print "run time: ", startsw()-t0
}

// doprotocols(1,2) // minimal test
// doprotocols(1,12) // the full battery

print "To run through cases A1-12, enter"
print "doprotocols(1,12)"
print "at the oc> prompt, then press return."

/*
Results obtained with adaptive integration, 
finding threshold to 4 place accuracy.
doprotocols(12)
test 	wform 	tp 	nc 	xa,ya 	 	xc,yc 		   thresh
1 	pls 	0.005 	1 	50.00, 0.25 	 0.00, 0.25 	   11.08643
2 	pls 	2.000 	1 	50.00, 0.25 	 0.00, 0.25 	    0.47032
3 	pls 	0.005 	1 	50.00, 0.25 	 0.00, 1.00 	  409.95312
4 	pls 	2.000 	1 	50.00, 0.25 	 0.00, 1.00 	   12.83545
5 	pls 	2.000 	1 	 0.00, 0.25 	50.00, 0.25 	    2.10657
6 	pls 	2.000 	1 	 0.00, 1.00 	 1.00, 1.00 	   11.00342
7 	sqr 	0.005 	1 	50.00, 0.25 	 0.00, 0.25 	   32.57227
8 	sqr 	2.000 	1 	50.00, 0.25 	 0.00, 0.25 	    0.47038
9 	sin 	0.005 	1 	50.00, 0.25 	 0.00, 0.25 	   48.41992
10 	sin 	0.100 	1 	50.00, 0.25 	 0.00, 0.25 	    1.44220
11 	sin 	0.005 	20000 	50.00, 0.25 	 0.00, 0.25 	   14.86279
12 	sin 	0.100 	10 	50.00, 0.25 	 0.00, 0.25 	    1.30255
*/