From 483f9d93e5771aaf7c6ab849ca8c5718d9e87ea2 Mon Sep 17 00:00:00 2001 From: Pramod Kumbhar Date: Sat, 10 Oct 2020 03:04:04 +0530 Subject: [PATCH] Add SEClamp support in CoreNEURON (#381) * Add SEClamp support in CoreNEURON - add svclmp.mod mod file from NEURON - update registeration callbacks - expose at_time in the header file * Fixed compilation issue for GPU with at_time and updated NMODL commit * Added definition of at_time in nrnoc_ml.ispc * Small fix for at_time declaration outside ispc Co-authored-by: Ioannis Magkanaris --- coreneuron/mechanism/mech/cfile/cabvars.h | 4 +- coreneuron/mechanism/mech/modfile/svclmp.mod | 107 +++++++++++++++++++ coreneuron/mechanism/membfunc.hpp | 8 ++ coreneuron/mechanism/nrnoc_ml.ispc | 5 +- coreneuron/network/cvodestb.cpp | 5 +- 5 files changed, 124 insertions(+), 5 deletions(-) create mode 100755 coreneuron/mechanism/mech/modfile/svclmp.mod diff --git a/coreneuron/mechanism/mech/cfile/cabvars.h b/coreneuron/mechanism/mech/cfile/cabvars.h index ce6a511ca..fa1ee2c05 100644 --- a/coreneuron/mechanism/mech/cfile/cabvars.h +++ b/coreneuron/mechanism/mech/cfile/cabvars.h @@ -32,7 +32,7 @@ extern void capacitance_reg(void), _passive_reg(void), #if EXTRACELLULAR extracell_reg_(void), #endif - _stim_reg(void), _hh_reg(void), _netstim_reg(void), _expsyn_reg(void), _exp2syn_reg(void); + _stim_reg(void), _hh_reg(void), _netstim_reg(void), _expsyn_reg(void), _exp2syn_reg(void), _svclmp_reg(void); static void (*mechanism[])(void) = {/* type will start at 3 */ capacitance_reg, _passive_reg, @@ -40,6 +40,6 @@ static void (*mechanism[])(void) = {/* type will start at 3 */ /* extracellular requires special handling and must be type 5 */ extracell_reg_, #endif - _stim_reg, _hh_reg, _expsyn_reg, _netstim_reg, _exp2syn_reg, 0}; + _stim_reg, _hh_reg, _expsyn_reg, _netstim_reg, _exp2syn_reg, _svclmp_reg, 0}; } // namespace coreneuron diff --git a/coreneuron/mechanism/mech/modfile/svclmp.mod b/coreneuron/mechanism/mech/modfile/svclmp.mod new file mode 100755 index 000000000..daa5de2d4 --- /dev/null +++ b/coreneuron/mechanism/mech/modfile/svclmp.mod @@ -0,0 +1,107 @@ +TITLE svclmp.mod +COMMENT +Single electrode Voltage clamp with three levels. +Clamp is on at time 0, and off at time +dur1+dur2+dur3. When clamp is off the injected current is 0. +The clamp levels are amp1, amp2, amp3. +i is the injected current, vc measures the control voltage) +Do not insert several instances of this model at the same location in order to +make level changes. That is equivalent to independent clamps and they will +have incompatible internal state values. +The electrical circuit for the clamp is exceedingly simple: +vc ---'\/\/`--- cell + rs + +Note that since this is an electrode current model v refers to the +internal potential which is equivalent to the membrane potential v when +there is no extracellular membrane mechanism present but is v+vext when +one is present. +Also since i is an electrode current, +positive values of i depolarize the cell. (Normally, positive membrane currents +are outward and thus hyperpolarize the cell) +ENDCOMMENT + +INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} + +DEFINE NSTEP 3 + +NEURON { + POINT_PROCESS SEClamp + ELECTRODE_CURRENT i + RANGE dur1, amp1, dur2, amp2, dur3, amp3, rs, vc, i +} + +UNITS { + (nA) = (nanoamp) + (mV) = (millivolt) + (uS) = (microsiemens) +} + + +PARAMETER { + rs = 1 (megohm) <1e-9, 1e9> + dur1 (ms) amp1 (mV) + dur2 (ms) <0,1e9> amp2 (mV) + dur3 (ms) <0,1e9> amp3 (mV) +} + +ASSIGNED { + v (mV) : automatically v + vext when extracellular is present + i (nA) + vc (mV) + tc2 (ms) + tc3 (ms) + on +} + +INITIAL { + tc2 = dur1 + dur2 + tc3 = tc2 + dur3 + on = 0 +} + +BREAKPOINT { + SOLVE icur METHOD after_cvode + vstim() +} + +PROCEDURE icur() { + if (on) { + i = (vc - v)/rs + }else{ + i = 0 + } +} + +COMMENT +The SOLVE of icur() in the BREAKPOINT block is necessary to compute +i=(vc - v(t))/rs instead of i=(vc - v(t-dt))/rs +This is important for time varying vc because the actual i used in +the implicit method is equivalent to (vc - v(t)/rs due to the +calculation of di/dv from the BREAKPOINT block. +The reason this works is because the SOLVE statement in the BREAKPOINT block +is executed after the membrane potential is advanced. + +It is a shame that vstim has to be called twice but putting the call +in a SOLVE block would cause playing a Vector into vc to be off by one +time step. +ENDCOMMENT + +PROCEDURE vstim() { + on = 1 + if (dur1) {at_time(dur1)} + if (dur2) {at_time(tc2)} + if (dur3) {at_time(tc3)} + if (t < dur1) { + vc = amp1 + }else if (t < tc2) { + vc = amp2 + }else if (t < tc3) { + vc = amp3 + }else { + vc = 0 + on = 0 + } + icur() +} + diff --git a/coreneuron/mechanism/membfunc.hpp b/coreneuron/mechanism/membfunc.hpp index 36204bd15..ea073e59b 100644 --- a/coreneuron/mechanism/membfunc.hpp +++ b/coreneuron/mechanism/membfunc.hpp @@ -196,3 +196,11 @@ extern void hoc_malchk(void); /* just a stub */ extern void* hoc_Emalloc(size_t); } // namespace coreneuron + +/* + * at_time function declared outside of coreneuron namespace + * because it can be also called from ISPC code. That's also + * the reason for extern "C". + */ +#pragma acc routine seq +extern "C" int at_time(coreneuron::NrnThread*, double); diff --git a/coreneuron/mechanism/nrnoc_ml.ispc b/coreneuron/mechanism/nrnoc_ml.ispc index a2956cb80..1036dc935 100644 --- a/coreneuron/mechanism/nrnoc_ml.ispc +++ b/coreneuron/mechanism/nrnoc_ml.ispc @@ -166,4 +166,7 @@ struct NrnThread { void* mapping; }; -extern uniform double ispc_celsius; \ No newline at end of file +extern uniform double ispc_celsius; + +extern "C" +int at_time(NrnThread* nt, varying double te); diff --git a/coreneuron/network/cvodestb.cpp b/coreneuron/network/cvodestb.cpp index d84515e58..372bbea89 100644 --- a/coreneuron/network/cvodestb.cpp +++ b/coreneuron/network/cvodestb.cpp @@ -105,11 +105,12 @@ void fixed_play_continuous(NrnThread* nt) { } } -int at_time(NrnThread* nt, double te) { +} // namespace coreneuron + +extern "C" int at_time(coreneuron::NrnThread* nt, double te) { double x = te - 1e-11; if (x <= nt->_t && x > (nt->_t - nt->_dt)) { return 1; } return 0; } -} // namespace coreneuron