Peak duration

The basics of how to develop, test, and use models.
Post Reply
menica
Posts: 68
Joined: Wed Sep 02, 2015 11:02 am

Peak duration

Post by menica » Fri Feb 16, 2018 11:03 am

Hi,
I would like to modify this code in order to save the vmax_mhw(0.5), vhalf_mhw(0.5), hw_mhw(0.5) in vectors because I would like to consider the case of multiple spikes in one simulation.

Code: Select all

COMMENT
Measures peak depol
and calculates spike half width from the times at which v crosses a (depolarized) threshold.
Threshold may be specified by the user, or determined in the previous run.

USAGE EXAMPLES

//////////////////////////
// User-specified threshold
forall insert mhw
forall for (x,0) vhalf_mhw(x) = THRESH // must assign value everywhere
mode_mhw = 0 // determine half width from fixed threshold
run()
printf(" base \t peak \t vhalf \thalf width\n")
printf("%6.2f \t%6.2f \t%6.2f \t%6.2f\n", \
       vinit_mhw(0.5), vmax_mhw(0.5), vhalf_mhw(0.5), hw_mhw(0.5))
//////////////////////////

//////////////////////////
// Dynamically-determined threshold
// run two simulations, first time with parameter mode_meas = 1
//   and second time with mode_meas = 2
// at end of first run, tmax and vmax will equal time and value of peak depol
// at end of second run, vhalf will be threshold for measurement of spike half width,
//   t0 and t1 will be threshold crossing times, and hw will be spike half width
forall insert mhw
mode_mhw = 1
run() // find local vmax and tmax
mode_mhw = 2
run() // find spike hw
printf(" base \t peak \t vhalf \thalf width\n")
printf("%6.2f \t%6.2f \t%6.2f \t%6.2f\n", \
       vinit_mhw(0.5), vmax_mhw(0.5), vhalf_mhw(0.5), hw_mhw(0.5))
//////////////////////////

hw will be -1 if there is no max, or if simulation ends before v falls below vhalf

Be cautious when using with adaptive integration--if the integrator uses long dt, 
t0 or t1 may be missed by a wide margin.
ENDCOMMENT

NEURON {
  SUFFIX mhw
  : mode values--
  : fixed threshold--0 use user-specified vhalf
  : dynamic threshold--1 measure vmax, 2 calc vhalf and measure halfwidth
  GLOBAL mode
  RANGE vinit, vmax, tmax
  RANGE vhalf, t0, t1, hw
}

UNITS {
  (mA) = (milliamp)
  (mV) = (millivolt)
  (mM) = (milli/liter)
}

PARAMETER {
  : mode values--
  : fixed threshold--0 use user-specified vhalf
  : dynamic threshold--1 measure vmax, 2 calc vhalf and measure halfwidth
  mode = 0 (1) : default is fixed (user-specified) threshold
}

ASSIGNED {
  v (mV)     : local v
  vinit (mV) : initial local v
  vmax (mV)  : max local v during previous run
  tmax (ms)  : time at which vmax occurred
  vhalf (mV) : (vinit + vmax)/2
  t0 (ms)    : time in rising phase of spike when v first > vhalf
  t1 (ms)    : time in falling phase of spike when v first < vhalf
  hw (ms)    : t1-t0
  findwhich (1) : 0 to find t0, 1 to find t1, 2 to find neither
}

INITIAL {
  if (mode==1) { : measure peak v then calc vhalf
: printf("Finding vmax\n")
    vinit = v
    vmax = v
    tmax = -1 (ms) : nonsense values for tmax, t0, t1, hw
    vhalf = v
    t0 = -1 (ms)
    t1 = -1 (ms)
    hw = -1 (ms)
  } else if (mode==2) { : calc vhalf from vinit and vmax in order to determine halfwidth
: printf("Determining depolarization halfwidth\n")
    vhalf = (vinit + vmax)/2
    findwhich = 0 : 0 to find t0, 1 to find t1
  } else if (mode==0) {
    vinit = v
    vmax = v
    tmax = -1 (ms) : nonsense values for tmax, t0, t1, hw
    t0 = -1 (ms)
    t1 = -1 (ms)
    hw = -1 (ms)
    findwhich = 0 : 0 to find t0, 1 to find t1
  }
}

PROCEDURE findmax() {
  if (v>vmax) {
    vmax = v
    tmax = t
  }
}

: find threshold crossings
PROCEDURE findx() {
  if (findwhich==0) {
    if (v > vhalf) {
      t0 = t
      findwhich = 1
    }
  } else if (findwhich==1) {
    if (v < vhalf) { 
      t1 = t
      hw = t1-t0
      findwhich = 2 : stop looking already
    }
  }
}

: BREAKPOINT {
AFTER SOLVE { : works as well, executed half as many times
  if (mode==1) { : measure peak v
    findmax()
  } else if (mode==2) {
    findx()
  } else if (mode==0) {
    findmax() : might as well, even if we don't use it to find threshold
    findx()
  }
}


now I am saving the vmax_mhw(0.5), vhalf_mhw(0.5), hw_mhw(0.5) in vectors but in this way I have a value for each time point.

Best
Menica

ted
Site Admin
Posts: 5334
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Peak duration

Post by ted » Tue Feb 20, 2018 5:23 pm

Don't even bother. NMODL is not well suited to that task. Develop your own code in Python or hoc that identifies each spike and determines its half width. That will involve many design decisions (especially if the spikes vary in amplitude, or if there is much subthreshold fluctuation of membrane potential) , and many cycles of revision and refinement.

menica
Posts: 68
Joined: Wed Sep 02, 2015 11:02 am

Re: Peak duration

Post by menica » Wed Feb 21, 2018 6:39 am

dear Ted,
thanks for your reply, I tried to implement this in hoc

//Determine the baseline
t2= tstop
j2 = tvec.indwhere(">=", t2-1)
V0baseline = vvec.mean(j2-stim[0].del, j2)

//find the max and the corresponding index
vmax = vvecE.max(j0start, j2)
vmax_index = vvecE.max_ind(j0start, j2)

//Calculate half width
vhalf = (vmax + V0baseline)/2

//find index at which v first cross vhalf in rise phase
jup = vvec.indwhere(">=", vhalf)

//find index at which v first cross vhalf in fall phase
jbelow = vvec.indwhere("<=", vhalf)

//Half width
deltaT = jbelow - jup


but something is wrong. how to properly calculate the time at which the voltage cross the threshold in the falling phase?
and in the case of multiple spikes, how to save the values of the half width of each maxima?

I thought I can use NetCon in this way:

// to record spike times
objref ncc,spvec, nilc
somaI ncc = new NetCon(&v(0.5), nilc)
ncc.threshold = -20
spvecc = new Vector()
ncc.record(spvecc)

and calculate the maximum in each interval between two time points recorded my netcon
for i=0,spveccI.size()-1{
j1 = tvec.indwhere("==", spvecc.x)
j2 = tvec.indwhere("==", spvecc.x[i+1])
vmax = vvecI.max(j1, j2)
print vmax
}

but indwhere is returning j1 and j2 =-1

Best
Menica

menica
Posts: 68
Joined: Wed Sep 02, 2015 11:02 am

Re: Peak duration

Post by menica » Thu Feb 22, 2018 10:33 am

Hi,
I implemented the code in hoc:

// to record spike times
objref vvec, tvec, inac
tvec = new Vector()
tvec.record(&t) // record time
vvec = new Vector()
vvec.record(&soma.v(0.5))

objref ncc,spvecc,nilc
soma ncc = new NetCon(&v(0.5), nilc)
ncc.threshold = -20
spvecc = new Vector()
ncc.record(spvecc)

objref vbx, hwV, maxV
hwV = new Vector()
maxV = new Vector()

proc calcola() {
nspikes = spvecc.size()
tend = tstop
jend = tvec.indwhere(">=", tend-1)
V0baseline = vvec.mean(jend-stim[0].del, jend)
if (nspikes >= 1) {
for i=0,spvecc.size()-1{
j1 = tvec.indwhere(">=", spvecc.x)
if ((i+1) == spvecc.size() ) {
j2 = jend
}else{
j2 = tvec.indwhere(">=", spvecc.x[i+1])
}
maxv= vvec.max(j1, j2) // amplitude maxima
maxV.insrt(i, maxv)

//vhalf = vvec.x[j1]
vhalf = (maxv+V0baseline)/2
//print vhalf

vbx = vvec.c(j1, j2)
b = vbx.indwhere("<", vhalf)
b1= j1+b
td = tvec.x[b1]
hw = td - spvecc.x //half width duration
hW.insrt(i, hwI)
}
}
}

It seems to work properly.
Now, my question is how to calculate in a correct way the vhalf?
vhalf = vvec.x[j1]
vhalf = (maxv+V0baseline)/2

Best
Menica

Post Reply