diff --git a/spike_travel/coreneuron/expsyn2.cpp b/spike_travel/coreneuron/expsyn2.cpp new file mode 100644 index 0000000..fd70b44 --- /dev/null +++ b/spike_travel/coreneuron/expsyn2.cpp @@ -0,0 +1,424 @@ +/********************************************************* +Model Name : ExpSyn2 +Filename : expsyn2.mod +NMODL Version : 7.7.0 +Vectorized : true +Threadsafe : true +Created : DATE +Simulator : CoreNEURON +Backend : C++ (api-compatibility) +NMODL Compiler : VERSION +*********************************************************/ + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace coreneuron { + #ifndef NRN_PRCELLSTATE + #define NRN_PRCELLSTATE 0 + #endif + + + /** channel information */ + static const char *mechanism_info[] = { + "7.7.0", + "ExpSyn2", + "tau", + "e", + 0, + "i", + 0, + "g", + 0, + 0 + }; + + + /** all global variables */ + struct ExpSyn2_Store { + int point_type{}; + double g0{}; + int reset{}; + int mech_type{}; + int slist1[1]{3}; + int dlist1[1]{4}; + }; + static_assert(std::is_trivially_copy_constructible_v); + static_assert(std::is_trivially_move_constructible_v); + static_assert(std::is_trivially_copy_assignable_v); + static_assert(std::is_trivially_move_assignable_v); + static_assert(std::is_trivially_destructible_v); + ExpSyn2_Store ExpSyn2_global; + + + /** all mechanism instance variables and global variables */ + struct ExpSyn2_Instance { + const double* tau{}; + const double* e{}; + double* i{}; + double* g{}; + double* Dg{}; + double* v_unused{}; + double* g_unused{}; + double* tsave{}; + const double* node_area{}; + const int* point_process{}; + ExpSyn2_Store* global{&ExpSyn2_global}; + }; + + + /** connect global (scalar) variables to hoc -- */ + static DoubScal hoc_scalar_double[] = { + {nullptr, nullptr} + }; + + + /** connect global (array) variables to hoc -- */ + static DoubVec hoc_vector_double[] = { + {nullptr, nullptr, 0} + }; + + + static inline int first_pointer_var_index() { + return -1; + } + + + static inline int first_random_var_index() { + return -1; + } + + + static inline int num_net_receive_args() { + return 1; + } + + + static inline int float_variables_size() { + return 8; + } + + + static inline int int_variables_size() { + return 2; + } + + + static inline int get_mech_type() { + return ExpSyn2_global.mech_type; + } + + + static inline Memb_list* get_memb_list(NrnThread* nt) { + if (!nt->_ml_list) { + return nullptr; + } + return nt->_ml_list[get_mech_type()]; + } + + + static inline void* mem_alloc(size_t num, size_t size, size_t alignment = 16) { + void* ptr; + posix_memalign(&ptr, alignment, num*size); + memset(ptr, 0, size); + return ptr; + } + + + static inline void mem_free(void* ptr) { + free(ptr); + } + + + static inline void coreneuron_abort() { + abort(); + } + + // Allocate instance structure + static void nrn_private_constructor_ExpSyn2(NrnThread* nt, Memb_list* ml, int type) { + assert(!ml->instance); + assert(!ml->global_variables); + assert(ml->global_variables_size == 0); + auto* const inst = new ExpSyn2_Instance{}; + assert(inst->global == &ExpSyn2_global); + ml->instance = inst; + ml->global_variables = inst->global; + ml->global_variables_size = sizeof(ExpSyn2_Store); + } + + // Deallocate the instance structure + static void nrn_private_destructor_ExpSyn2(NrnThread* nt, Memb_list* ml, int type) { + auto* const inst = static_cast(ml->instance); + assert(inst); + assert(inst->global); + assert(inst->global == &ExpSyn2_global); + assert(inst->global == ml->global_variables); + assert(ml->global_variables_size == sizeof(ExpSyn2_Store)); + delete inst; + ml->instance = nullptr; + ml->global_variables = nullptr; + ml->global_variables_size = 0; + } + + /** initialize mechanism instance variables */ + static inline void setup_instance(NrnThread* nt, Memb_list* ml) { + auto* const inst = static_cast(ml->instance); + assert(inst); + assert(inst->global); + assert(inst->global == &ExpSyn2_global); + assert(inst->global == ml->global_variables); + assert(ml->global_variables_size == sizeof(ExpSyn2_Store)); + int pnodecount = ml->_nodecount_padded; + Datum* indexes = ml->pdata; + inst->tau = ml->data+0*pnodecount; + inst->e = ml->data+1*pnodecount; + inst->i = ml->data+2*pnodecount; + inst->g = ml->data+3*pnodecount; + inst->Dg = ml->data+4*pnodecount; + inst->v_unused = ml->data+5*pnodecount; + inst->g_unused = ml->data+6*pnodecount; + inst->tsave = ml->data+7*pnodecount; + inst->node_area = nt->_data; + inst->point_process = ml->pdata; + } + + + + static void nrn_alloc_ExpSyn2(double* data, Datum* indexes, int type) { + // do nothing + } + + + void nrn_constructor_ExpSyn2(NrnThread* nt, Memb_list* ml, int type) { + #ifndef CORENEURON_BUILD + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #endif + } + + + void nrn_destructor_ExpSyn2(NrnThread* nt, Memb_list* ml, int type) { + #ifndef CORENEURON_BUILD + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #endif + } + + + static inline void net_receive_kernel_ExpSyn2(double t, Point_process* pnt, ExpSyn2_Instance* inst, NrnThread* nt, Memb_list* ml, int weight_index, double flag) { + int tid = pnt->_tid; + int id = pnt->_i_instance; + double v = 0; + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + double* data = ml->data; + double* weights = nt->weights; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + + double* weight = weights + weight_index + 0; + inst->tsave[id] = t; + { + inst->g[id] = inst->g[id] + (*weight); + } + } + + + static void net_receive_ExpSyn2(Point_process* pnt, int weight_index, double flag) { + NrnThread* nt = nrn_threads + pnt->_tid; + Memb_list* ml = get_memb_list(nt); + NetReceiveBuffer_t* nrb = ml->_net_receive_buffer; + if (nrb->_cnt >= nrb->_size) { + realloc_net_receive_buffer(nt, ml); + } + int id = nrb->_cnt; + nrb->_pnt_index[id] = pnt-nt->pntprocs; + nrb->_weight_index[id] = weight_index; + nrb->_nrb_t[id] = nt->_t; + nrb->_nrb_flag[id] = flag; + nrb->_cnt++; + } + + + void net_buf_receive_ExpSyn2(NrnThread* nt) { + Memb_list* ml = get_memb_list(nt); + if (!ml) { + return; + } + + NetReceiveBuffer_t* nrb = ml->_net_receive_buffer; + auto* const inst = static_cast(ml->instance); + int count = nrb->_displ_cnt; + #pragma omp simd + #pragma ivdep + for (int i = 0; i < count; i++) { + int start = nrb->_displ[i]; + int end = nrb->_displ[i+1]; + for (int j = start; j < end; j++) { + int index = nrb->_nrb_index[j]; + int offset = nrb->_pnt_index[index]; + double t = nrb->_nrb_t[index]; + int weight_index = nrb->_weight_index[index]; + double flag = nrb->_nrb_flag[index]; + Point_process* point_process = nt->pntprocs + offset; + net_receive_kernel_ExpSyn2(t, point_process, inst, nt, ml, weight_index, flag); + } + } + nrb->_displ_cnt = 0; + nrb->_cnt = 0; + } + + + /** initialize channel */ + void nrn_init_ExpSyn2(NrnThread* nt, Memb_list* ml, int type) { + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + + setup_instance(nt, ml); + auto* const inst = static_cast(ml->instance); + + if (_nrn_skip_initmodel == 0) { + #pragma omp simd + #pragma ivdep + for (int id = 0; id < nodecount; id++) { + inst->tsave[id] = -1e20; + int node_id = node_index[id]; + double v = voltage[node_id]; + #if NRN_PRCELLSTATE + inst->v_unused[id] = v; + #endif + inst->g[id] = inst->global->g0; + inst->g[id] = 0.0; + } + } + } + + + inline double nrn_current_ExpSyn2(int id, int pnodecount, ExpSyn2_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v) { + double current = 0.0; + inst->i[id] = inst->g[id] * (v - inst->e[id]); + current += inst->i[id]; + return current; + } + + + /** update current */ + void nrn_cur_ExpSyn2(NrnThread* nt, Memb_list* ml, int type) { + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + double* vec_rhs = nt->_actual_rhs; + double* vec_d = nt->_actual_d; + double* shadow_rhs = nt->_shadow_rhs; + double* shadow_d = nt->_shadow_d; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #pragma omp simd + #pragma ivdep + for (int id = 0; id < nodecount; id++) { + int node_id = node_index[id]; + double v = voltage[node_id]; + #if NRN_PRCELLSTATE + inst->v_unused[id] = v; + #endif + double g = nrn_current_ExpSyn2(id, pnodecount, inst, data, indexes, thread, nt, v+0.001); + double rhs = nrn_current_ExpSyn2(id, pnodecount, inst, data, indexes, thread, nt, v); + g = (g-rhs)/0.001; + double mfactor = 1.e2/inst->node_area[indexes[0*pnodecount + id]]; + g = g*mfactor; + rhs = rhs*mfactor; + #if NRN_PRCELLSTATE + inst->g_unused[id] = g; + #endif + shadow_rhs[id] = rhs; + shadow_d[id] = g; + } + for (int id = 0; id < nodecount; id++) { + int node_id = node_index[id]; + vec_rhs[node_id] -= shadow_rhs[id]; + vec_d[node_id] += shadow_d[id]; + } + } + + + /** update state */ + void nrn_state_ExpSyn2(NrnThread* nt, Memb_list* ml, int type) { + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #pragma omp simd + #pragma ivdep + for (int id = 0; id < nodecount; id++) { + int node_id = node_index[id]; + double v = voltage[node_id]; + #if NRN_PRCELLSTATE + inst->v_unused[id] = v; + #endif + inst->g[id] = inst->g[id] + (1.0 - exp(nt->_dt * (( -1.0) / inst->tau[id]))) * ( -(0.0) / (( -1.0) / inst->tau[id]) - inst->g[id]); + } + } + + + /** register channel with the simulator */ + void _expsyn2_reg() { + + int mech_type = nrn_get_mechtype("ExpSyn2"); + ExpSyn2_global.mech_type = mech_type; + if (mech_type == -1) { + return; + } + + _nrn_layout_reg(mech_type, 0); + point_register_mech(mechanism_info, nrn_alloc_ExpSyn2, nrn_cur_ExpSyn2, nullptr, nrn_state_ExpSyn2, nrn_init_ExpSyn2, nrn_private_constructor_ExpSyn2, nrn_private_destructor_ExpSyn2, first_pointer_var_index(), nullptr, nullptr, 1); + + hoc_register_prop_size(mech_type, float_variables_size(), int_variables_size()); + hoc_register_dparam_semantics(mech_type, 0, "area"); + hoc_register_dparam_semantics(mech_type, 1, "pntproc"); + hoc_register_net_receive_buffering(net_buf_receive_ExpSyn2, mech_type); + set_pnt_receive(mech_type, net_receive_ExpSyn2, nullptr, num_net_receive_args()); + hoc_register_var(hoc_scalar_double, hoc_vector_double, NULL); + } +} diff --git a/spike_travel/coreneuron/hodhux.cpp b/spike_travel/coreneuron/hodhux.cpp new file mode 100644 index 0000000..1ac5ac6 --- /dev/null +++ b/spike_travel/coreneuron/hodhux.cpp @@ -0,0 +1,472 @@ +/********************************************************* +Model Name : hodhux +Filename : hodhux.mod +NMODL Version : 7.7.0 +Vectorized : true +Threadsafe : true +Created : DATE +Simulator : CoreNEURON +Backend : C++ (api-compatibility) +NMODL Compiler : VERSION +*********************************************************/ + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace coreneuron { + #ifndef NRN_PRCELLSTATE + #define NRN_PRCELLSTATE 0 + #endif + + + /** channel information */ + static const char *mechanism_info[] = { + "7.7.0", + "hodhux", + "gnabar_hodhux", + "gkbar_hodhux", + "gl_hodhux", + "el_hodhux", + 0, + "il_hodhux", + "minf_hodhux", + "hinf_hodhux", + "ninf_hodhux", + "mexp_hodhux", + "hexp_hodhux", + "nexp_hodhux", + 0, + "m_hodhux", + "h_hodhux", + "n_hodhux", + 0, + 0 + }; + + + /** all global variables */ + struct hodhux_Store { + int na_type{}; + int k_type{}; + double m0{}; + double h0{}; + double n0{}; + int reset{}; + int mech_type{}; + }; + static_assert(std::is_trivially_copy_constructible_v); + static_assert(std::is_trivially_move_constructible_v); + static_assert(std::is_trivially_copy_assignable_v); + static_assert(std::is_trivially_move_assignable_v); + static_assert(std::is_trivially_destructible_v); + hodhux_Store hodhux_global; + + + /** all mechanism instance variables and global variables */ + struct hodhux_Instance { + double* celsius{&coreneuron::celsius}; + const double* gnabar{}; + const double* gkbar{}; + const double* gl{}; + const double* el{}; + double* il{}; + double* minf{}; + double* hinf{}; + double* ninf{}; + double* mexp{}; + double* hexp{}; + double* nexp{}; + double* m{}; + double* h{}; + double* n{}; + double* ena{}; + double* ek{}; + double* Dm{}; + double* Dh{}; + double* Dn{}; + double* ina{}; + double* ik{}; + double* v_unused{}; + double* g_unused{}; + const double* ion_ena{}; + double* ion_ina{}; + double* ion_dinadv{}; + const double* ion_ek{}; + double* ion_ik{}; + double* ion_dikdv{}; + hodhux_Store* global{&hodhux_global}; + }; + + + /** connect global (scalar) variables to hoc -- */ + static DoubScal hoc_scalar_double[] = { + {nullptr, nullptr} + }; + + + /** connect global (array) variables to hoc -- */ + static DoubVec hoc_vector_double[] = { + {nullptr, nullptr, 0} + }; + + + static inline int first_pointer_var_index() { + return -1; + } + + + static inline int first_random_var_index() { + return -1; + } + + + static inline int float_variables_size() { + return 23; + } + + + static inline int int_variables_size() { + return 6; + } + + + static inline int get_mech_type() { + return hodhux_global.mech_type; + } + + + static inline Memb_list* get_memb_list(NrnThread* nt) { + if (!nt->_ml_list) { + return nullptr; + } + return nt->_ml_list[get_mech_type()]; + } + + + static inline void* mem_alloc(size_t num, size_t size, size_t alignment = 16) { + void* ptr; + posix_memalign(&ptr, alignment, num*size); + memset(ptr, 0, size); + return ptr; + } + + + static inline void mem_free(void* ptr) { + free(ptr); + } + + + static inline void coreneuron_abort() { + abort(); + } + + // Allocate instance structure + static void nrn_private_constructor_hodhux(NrnThread* nt, Memb_list* ml, int type) { + assert(!ml->instance); + assert(!ml->global_variables); + assert(ml->global_variables_size == 0); + auto* const inst = new hodhux_Instance{}; + assert(inst->global == &hodhux_global); + ml->instance = inst; + ml->global_variables = inst->global; + ml->global_variables_size = sizeof(hodhux_Store); + } + + // Deallocate the instance structure + static void nrn_private_destructor_hodhux(NrnThread* nt, Memb_list* ml, int type) { + auto* const inst = static_cast(ml->instance); + assert(inst); + assert(inst->global); + assert(inst->global == &hodhux_global); + assert(inst->global == ml->global_variables); + assert(ml->global_variables_size == sizeof(hodhux_Store)); + delete inst; + ml->instance = nullptr; + ml->global_variables = nullptr; + ml->global_variables_size = 0; + } + + /** initialize mechanism instance variables */ + static inline void setup_instance(NrnThread* nt, Memb_list* ml) { + auto* const inst = static_cast(ml->instance); + assert(inst); + assert(inst->global); + assert(inst->global == &hodhux_global); + assert(inst->global == ml->global_variables); + assert(ml->global_variables_size == sizeof(hodhux_Store)); + int pnodecount = ml->_nodecount_padded; + Datum* indexes = ml->pdata; + inst->gnabar = ml->data+0*pnodecount; + inst->gkbar = ml->data+1*pnodecount; + inst->gl = ml->data+2*pnodecount; + inst->el = ml->data+3*pnodecount; + inst->il = ml->data+4*pnodecount; + inst->minf = ml->data+5*pnodecount; + inst->hinf = ml->data+6*pnodecount; + inst->ninf = ml->data+7*pnodecount; + inst->mexp = ml->data+8*pnodecount; + inst->hexp = ml->data+9*pnodecount; + inst->nexp = ml->data+10*pnodecount; + inst->m = ml->data+11*pnodecount; + inst->h = ml->data+12*pnodecount; + inst->n = ml->data+13*pnodecount; + inst->ena = ml->data+14*pnodecount; + inst->ek = ml->data+15*pnodecount; + inst->Dm = ml->data+16*pnodecount; + inst->Dh = ml->data+17*pnodecount; + inst->Dn = ml->data+18*pnodecount; + inst->ina = ml->data+19*pnodecount; + inst->ik = ml->data+20*pnodecount; + inst->v_unused = ml->data+21*pnodecount; + inst->g_unused = ml->data+22*pnodecount; + inst->ion_ena = nt->_data; + inst->ion_ina = nt->_data; + inst->ion_dinadv = nt->_data; + inst->ion_ek = nt->_data; + inst->ion_ik = nt->_data; + inst->ion_dikdv = nt->_data; + } + + + + static void nrn_alloc_hodhux(double* data, Datum* indexes, int type) { + // do nothing + } + + + void nrn_constructor_hodhux(NrnThread* nt, Memb_list* ml, int type) { + #ifndef CORENEURON_BUILD + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #endif + } + + + void nrn_destructor_hodhux(NrnThread* nt, Memb_list* ml, int type) { + #ifndef CORENEURON_BUILD + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #endif + } + + + inline double vtrap_hodhux(int id, int pnodecount, hodhux_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v, double x, double y); + inline int states_hodhux(int id, int pnodecount, hodhux_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v); + inline int rates_hodhux(int id, int pnodecount, hodhux_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v, double arg_v); + + + inline int states_hodhux(int id, int pnodecount, hodhux_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v) { + int ret_states = 0; + rates_hodhux(id, pnodecount, inst, data, indexes, thread, nt, v, v); + inst->m[id] = inst->m[id] + inst->mexp[id] * (inst->minf[id] - inst->m[id]); + inst->h[id] = inst->h[id] + inst->hexp[id] * (inst->hinf[id] - inst->h[id]); + inst->n[id] = inst->n[id] + inst->nexp[id] * (inst->ninf[id] - inst->n[id]); + return ret_states; + } + + + inline int rates_hodhux(int id, int pnodecount, hodhux_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v, double arg_v) { + int ret_rates = 0; + double q10, tinc, alpha, beta, sum; + q10 = pow(3.0, ((*(inst->celsius) - 6.3) / 10.0)); + tinc = -nt->_dt * q10; + alpha = .1 * vtrap_hodhux(id, pnodecount, inst, data, indexes, thread, nt, v, -(arg_v + 40.0), 10.0); + beta = 4.0 * exp( -(arg_v + 65.0) / 18.0); + sum = alpha + beta; + inst->minf[id] = alpha / sum; + inst->mexp[id] = 1.0 - exp(tinc * sum); + alpha = .07 * exp( -(arg_v + 65.0) / 20.0); + beta = 1.0 / (exp( -(arg_v + 35.0) / 10.0) + 1.0); + sum = alpha + beta; + inst->hinf[id] = alpha / sum; + inst->hexp[id] = 1.0 - exp(tinc * sum); + alpha = .01 * vtrap_hodhux(id, pnodecount, inst, data, indexes, thread, nt, v, -(arg_v + 55.0), 10.0); + beta = .125 * exp( -(arg_v + 65.0) / 80.0); + sum = alpha + beta; + inst->ninf[id] = alpha / sum; + inst->nexp[id] = 1.0 - exp(tinc * sum); + return ret_rates; + } + + + inline double vtrap_hodhux(int id, int pnodecount, hodhux_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v, double x, double y) { + double ret_vtrap = 0.0; + if (fabs(x / y) < 1e-6) { + ret_vtrap = y * (1.0 - x / y / 2.0); + } else { + ret_vtrap = x / (exp(x / y) - 1.0); + } + return ret_vtrap; + } + + + /** initialize channel */ + void nrn_init_hodhux(NrnThread* nt, Memb_list* ml, int type) { + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + + setup_instance(nt, ml); + auto* const inst = static_cast(ml->instance); + + if (_nrn_skip_initmodel == 0) { + #pragma omp simd + #pragma ivdep + for (int id = 0; id < nodecount; id++) { + int node_id = node_index[id]; + double v = voltage[node_id]; + #if NRN_PRCELLSTATE + inst->v_unused[id] = v; + #endif + inst->ena[id] = inst->ion_ena[indexes[0*pnodecount + id]]; + inst->ek[id] = inst->ion_ek[indexes[3*pnodecount + id]]; + inst->m[id] = inst->global->m0; + inst->h[id] = inst->global->h0; + inst->n[id] = inst->global->n0; + rates_hodhux(id, pnodecount, inst, data, indexes, thread, nt, v, v); + inst->m[id] = inst->minf[id]; + inst->h[id] = inst->hinf[id]; + inst->n[id] = inst->ninf[id]; + } + } + } + + + inline double nrn_current_hodhux(int id, int pnodecount, hodhux_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v) { + double current = 0.0; + inst->ina[id] = inst->gnabar[id] * inst->m[id] * inst->m[id] * inst->m[id] * inst->h[id] * (v - inst->ena[id]); + inst->ik[id] = inst->gkbar[id] * inst->n[id] * inst->n[id] * inst->n[id] * inst->n[id] * (v - inst->ek[id]); + inst->il[id] = inst->gl[id] * (v - inst->el[id]); + current += inst->il[id]; + current += inst->ina[id]; + current += inst->ik[id]; + return current; + } + + + /** update current */ + void nrn_cur_hodhux(NrnThread* nt, Memb_list* ml, int type) { + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + double* vec_rhs = nt->_actual_rhs; + double* vec_d = nt->_actual_d; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #pragma omp simd + #pragma ivdep + for (int id = 0; id < nodecount; id++) { + int node_id = node_index[id]; + double v = voltage[node_id]; + #if NRN_PRCELLSTATE + inst->v_unused[id] = v; + #endif + inst->ena[id] = inst->ion_ena[indexes[0*pnodecount + id]]; + inst->ek[id] = inst->ion_ek[indexes[3*pnodecount + id]]; + double g = nrn_current_hodhux(id, pnodecount, inst, data, indexes, thread, nt, v+0.001); + double dina = inst->ina[id]; + double dik = inst->ik[id]; + double rhs = nrn_current_hodhux(id, pnodecount, inst, data, indexes, thread, nt, v); + g = (g-rhs)/0.001; + inst->ion_dinadv[indexes[2*pnodecount + id]] += (dina-inst->ina[id])/0.001; + inst->ion_dikdv[indexes[5*pnodecount + id]] += (dik-inst->ik[id])/0.001; + inst->ion_ina[indexes[1*pnodecount + id]] += inst->ina[id]; + inst->ion_ik[indexes[4*pnodecount + id]] += inst->ik[id]; + #if NRN_PRCELLSTATE + inst->g_unused[id] = g; + #endif + vec_rhs[node_id] -= rhs; + vec_d[node_id] += g; + } + } + + + /** update state */ + void nrn_state_hodhux(NrnThread* nt, Memb_list* ml, int type) { + int nodecount = ml->nodecount; + int pnodecount = ml->_nodecount_padded; + const int* node_index = ml->nodeindices; + double* data = ml->data; + const double* voltage = nt->_actual_v; + Datum* indexes = ml->pdata; + ThreadDatum* thread = ml->_thread; + auto* const inst = static_cast(ml->instance); + + #pragma omp simd + #pragma ivdep + for (int id = 0; id < nodecount; id++) { + int node_id = node_index[id]; + double v = voltage[node_id]; + #if NRN_PRCELLSTATE + inst->v_unused[id] = v; + #endif + inst->ena[id] = inst->ion_ena[indexes[0*pnodecount + id]]; + inst->ek[id] = inst->ion_ek[indexes[3*pnodecount + id]]; + rates_hodhux(id, pnodecount, inst, data, indexes, thread, nt, v, v); + inst->m[id] = inst->m[id] + inst->mexp[id] * (inst->minf[id] - inst->m[id]); + inst->h[id] = inst->h[id] + inst->hexp[id] * (inst->hinf[id] - inst->h[id]); + inst->n[id] = inst->n[id] + inst->nexp[id] * (inst->ninf[id] - inst->n[id]); + } + } + + + /** register channel with the simulator */ + void _hodhux_reg() { + + int mech_type = nrn_get_mechtype("hodhux"); + hodhux_global.mech_type = mech_type; + if (mech_type == -1) { + return; + } + + _nrn_layout_reg(mech_type, 0); + register_mech(mechanism_info, nrn_alloc_hodhux, nrn_cur_hodhux, nullptr, nrn_state_hodhux, nrn_init_hodhux, nrn_private_constructor_hodhux, nrn_private_destructor_hodhux, first_pointer_var_index(), 1); + hodhux_global.na_type = nrn_get_mechtype("na_ion"); + hodhux_global.k_type = nrn_get_mechtype("k_ion"); + + hoc_register_prop_size(mech_type, float_variables_size(), int_variables_size()); + hoc_register_dparam_semantics(mech_type, 0, "na_ion"); + hoc_register_dparam_semantics(mech_type, 1, "na_ion"); + hoc_register_dparam_semantics(mech_type, 2, "na_ion"); + hoc_register_dparam_semantics(mech_type, 3, "k_ion"); + hoc_register_dparam_semantics(mech_type, 4, "k_ion"); + hoc_register_dparam_semantics(mech_type, 5, "k_ion"); + hoc_register_var(hoc_scalar_double, hoc_vector_double, NULL); + } +} diff --git a/spike_travel/neuron/expsyn2.cpp b/spike_travel/neuron/expsyn2.cpp new file mode 100644 index 0000000..74b73c8 --- /dev/null +++ b/spike_travel/neuron/expsyn2.cpp @@ -0,0 +1,343 @@ +/********************************************************* +Model Name : ExpSyn2 +Filename : expsyn2.mod +NMODL Version : 7.7.0 +Vectorized : true +Threadsafe : true +Created : DATE +Simulator : NEURON +Backend : C++ (api-compatibility) +NMODL Compiler : VERSION +*********************************************************/ + +#include +#include +#include + +#include "mech_api.h" +#include "neuron/cache/mechanism_range.hpp" +#include "nrniv_mf.h" +#include "section_fwd.hpp" + +/* NEURON global macro definitions */ +/* VECTORIZED */ +#define NRN_VECTORIZED 1 + +static constexpr auto number_of_datum_variables = 2; +static constexpr auto number_of_floating_point_variables = 8; + +namespace { +template +using _nrn_mechanism_std_vector = std::vector; +using _nrn_model_sorted_token = neuron::model_sorted_token; +using _nrn_mechanism_cache_range = neuron::cache::MechanismRange; +using _nrn_mechanism_cache_instance = neuron::cache::MechanismInstance; +using _nrn_non_owning_id_without_container = neuron::container::non_owning_identifier_without_container; +template +using _nrn_mechanism_field = neuron::mechanism::field; +template +void _nrn_mechanism_register_data_fields(Args&&... args) { + neuron::mechanism::register_data_fields(std::forward(args)...); +} +} // namespace + +extern Prop* nrn_point_prop_; + + +namespace neuron { + #ifndef NRN_PRCELLSTATE + #define NRN_PRCELLSTATE 0 + #endif + + + /** channel information */ + static const char *mechanism_info[] = { + "7.7.0", + "ExpSyn2", + "tau", + "e", + 0, + "i", + 0, + "g", + 0, + 0 + }; + + + /* NEURON global variables */ + static neuron::container::field_index _slist1[1], _dlist1[1]; + static int mech_type; + static int _pointtype; + static int hoc_nrnpointerindex = -1; + static _nrn_mechanism_std_vector _extcall_thread; + + + /** all global variables */ + struct ExpSyn2_Store { + double g0{}; + }; + static_assert(std::is_trivially_copy_constructible_v); + static_assert(std::is_trivially_move_constructible_v); + static_assert(std::is_trivially_copy_assignable_v); + static_assert(std::is_trivially_move_assignable_v); + static_assert(std::is_trivially_destructible_v); + ExpSyn2_Store ExpSyn2_global; + + + /** all mechanism instance variables and global variables */ + struct ExpSyn2_Instance { + double* tau{}; + double* e{}; + double* i{}; + double* g{}; + double* Dg{}; + double* v_unused{}; + double* g_unused{}; + double* tsave{}; + const double* const* node_area{}; + ExpSyn2_Store* global{&ExpSyn2_global}; + }; + + + struct ExpSyn2_NodeData { + int const * nodeindices; + double const * node_voltages; + double * node_diagonal; + double * node_rhs; + int nodecount; + }; + + + static ExpSyn2_Instance make_instance_ExpSyn2(_nrn_mechanism_cache_range& _ml) { + return ExpSyn2_Instance { + _ml.template fpfield_ptr<0>(), + _ml.template fpfield_ptr<1>(), + _ml.template fpfield_ptr<2>(), + _ml.template fpfield_ptr<3>(), + _ml.template fpfield_ptr<4>(), + _ml.template fpfield_ptr<5>(), + _ml.template fpfield_ptr<6>(), + _ml.template fpfield_ptr<7>(), + _ml.template dptr_field_ptr<0>() + }; + } + + + static ExpSyn2_NodeData make_node_data_ExpSyn2(NrnThread& _nt, Memb_list& _ml_arg) { + return ExpSyn2_NodeData { + _ml_arg.nodeindices, + _nt.node_voltage_storage(), + _nt.node_d_storage(), + _nt.node_rhs_storage(), + _ml_arg.nodecount + }; + } + + + static void nrn_alloc_ExpSyn2(Prop* _prop) { + Prop *prop_ion{}; + Datum *_ppvar{}; + if (nrn_point_prop_) { + _nrn_mechanism_access_alloc_seq(_prop) = _nrn_mechanism_access_alloc_seq(nrn_point_prop_); + _ppvar = _nrn_mechanism_access_dparam(nrn_point_prop_); + } else { + _ppvar = nrn_prop_datum_alloc(mech_type, 2, _prop); + _nrn_mechanism_access_dparam(_prop) = _ppvar; + _nrn_mechanism_cache_instance _ml_real{_prop}; + auto* const _ml = &_ml_real; + size_t const _iml{}; + assert(_nrn_mechanism_get_num_vars(_prop) == 8); + /*initialize range parameters*/ + _ml->template fpfield<0>(_iml) = 0.1; /* tau */ + _ml->template fpfield<1>(_iml) = 0; /* e */ + } + _nrn_mechanism_access_dparam(_prop) = _ppvar; + } + + + /* Point Process specific functions */ + static void* _hoc_create_pnt(Object* _ho) { + return create_point_process(_pointtype, _ho); + } + static void _hoc_destroy_pnt(void* _vptr) { + destroy_point_process(_vptr); + } + static double _hoc_loc_pnt(void* _vptr) { + return loc_point_process(_pointtype, _vptr); + } + static double _hoc_has_loc(void* _vptr) { + return has_loc_point(_vptr); + } + static double _hoc_get_loc_pnt(void* _vptr) { + return (get_loc_point_process(_vptr)); + } + /* Neuron setdata functions */ + extern void _nrn_setdata_reg(int, void(*)(Prop*)); + static void _setdata(Prop* _prop) { + } + static void _hoc_setdata(void* _vptr) { + Prop* _prop; + _prop = ((Point_process*)_vptr)->prop; + _setdata(_prop); + } + /* Mechanism procedures and functions */ + + + /** connect global (scalar) variables to hoc -- */ + static DoubScal hoc_scalar_double[] = { + {nullptr, nullptr} + }; + + + /** connect global (array) variables to hoc -- */ + static DoubVec hoc_vector_double[] = { + {nullptr, nullptr, 0} + }; + + + /* declaration of user functions */ + + + /* connect user functions to hoc names */ + static VoidFunc hoc_intfunc[] = { + {0, 0} + }; + static Member_func _member_func[] = { + {"loc", _hoc_loc_pnt}, + {"has_loc", _hoc_has_loc}, + {"get_loc", _hoc_get_loc_pnt}, + {0, 0} + }; + + + void nrn_init_ExpSyn2(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_ExpSyn2(_lmr); + auto node_data = make_node_data_ExpSyn2(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + auto* const _ml = &_lmr; + auto* _thread = _ml_arg->_thread; + for (int id = 0; id < nodecount; id++) { + + int node_id = node_data.nodeindices[id]; + auto* _ppvar = _ml_arg->pdata[id]; + auto v = node_data.node_voltages[node_id]; + inst.v_unused[id] = v; + inst.g[id] = 0.0; + } + } + + + inline double nrn_current_ExpSyn2(size_t id, ExpSyn2_Instance& inst, ExpSyn2_NodeData& node_data, double v) { + double current = 0.0; + inst.i[id] = inst.g[id] * (v - inst.e[id]); + current += inst.i[id]; + return current; + } + + + /** update current */ + void nrn_cur_ExpSyn2(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_ExpSyn2(_lmr); + auto node_data = make_node_data_ExpSyn2(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + auto* const _ml = &_lmr; + auto* _thread = _ml_arg->_thread; + for (int id = 0; id < nodecount; id++) { + int node_id = node_data.nodeindices[id]; + double v = node_data.node_voltages[node_id]; + double I1 = nrn_current_ExpSyn2(id, inst, node_data, v+0.001); + double I0 = nrn_current_ExpSyn2(id, inst, node_data, v); + double rhs = I0; + double g = (I1-I0)/0.001; + double mfactor = 1.e2/(*inst.node_area[id]); + g = g*mfactor; + rhs = rhs*mfactor; + node_data.node_rhs[node_id] -= rhs; + // remember the conductances so we can set them later + inst.g_unused[id] = g; + } + } + + + void nrn_state_ExpSyn2(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_ExpSyn2(_lmr); + auto node_data = make_node_data_ExpSyn2(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + auto* const _ml = &_lmr; + auto* _thread = _ml_arg->_thread; + for (int id = 0; id < nodecount; id++) { + + int node_id = node_data.nodeindices[id]; + auto* _ppvar = _ml_arg->pdata[id]; + auto v = node_data.node_voltages[node_id]; + inst.g[id] = inst.g[id] + (1.0 - exp(_nt->_dt * (( -1.0) / inst.tau[id]))) * ( -(0.0) / (( -1.0) / inst.tau[id]) - inst.g[id]); + } + } + + + /** nrn_jacob function */ + static void nrn_jacob_ExpSyn2(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_ExpSyn2(_lmr); + auto node_data = make_node_data_ExpSyn2(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + for (int id = 0; id < nodecount; id++) { + // set conductances properly + int node_id = node_data.nodeindices[id]; + node_data.node_diagonal[node_id] += inst.g_unused[id]; + } + } + static void nrn_net_receive_ExpSyn2(Point_process* _pnt, double* _args, double _lflag) { + _nrn_mechanism_cache_instance _ml_obj{_pnt->prop}; + auto * _nt = static_cast(_pnt->_vnt); + auto * _ml = &_ml_obj; + auto inst = make_instance_ExpSyn2(_ml_obj); + size_t id = 0; + double t = _nt->_t; + inst.g[id] = inst.g[id] + _args[0]; + + } + + + static void _initlists() { + /* g */ + _slist1[0] = {3, 0}; + /* Dg */ + _dlist1[0] = {4, 0}; + } + + + /** register channel with the simulator */ + extern "C" void _expsyn2_reg() { + _initlists(); + + + + _pointtype = point_register_mech(mechanism_info, nrn_alloc_ExpSyn2, nrn_cur_ExpSyn2, nrn_jacob_ExpSyn2, nrn_state_ExpSyn2, nrn_init_ExpSyn2, hoc_nrnpointerindex, 1, _hoc_create_pnt, _hoc_destroy_pnt, _member_func); + + mech_type = nrn_get_mechtype(mechanism_info[1]); + _nrn_mechanism_register_data_fields(mech_type, + _nrn_mechanism_field{"tau"} /* 0 */, + _nrn_mechanism_field{"e"} /* 1 */, + _nrn_mechanism_field{"i"} /* 2 */, + _nrn_mechanism_field{"g"} /* 3 */, + _nrn_mechanism_field{"Dg"} /* 4 */, + _nrn_mechanism_field{"v_unused"} /* 5 */, + _nrn_mechanism_field{"g_unused"} /* 6 */, + _nrn_mechanism_field{"tsave"} /* 7 */, + _nrn_mechanism_field{"node_area", "area"} /* 0 */, + _nrn_mechanism_field{"point_process", "pntproc"} /* 1 */ + ); + + hoc_register_prop_size(mech_type, 8, 2); + hoc_register_dparam_semantics(mech_type, 0, "area"); + hoc_register_dparam_semantics(mech_type, 1, "pntproc"); + hoc_register_var(hoc_scalar_double, hoc_vector_double, hoc_intfunc); + pnt_receive[mech_type] = nrn_net_receive_ExpSyn2; + pnt_receive_size[mech_type] = 1; + } +} diff --git a/spike_travel/neuron/hodhux.cpp b/spike_travel/neuron/hodhux.cpp new file mode 100644 index 0000000..05b2831 --- /dev/null +++ b/spike_travel/neuron/hodhux.cpp @@ -0,0 +1,587 @@ +/********************************************************* +Model Name : hodhux +Filename : hodhux.mod +NMODL Version : 7.7.0 +Vectorized : true +Threadsafe : true +Created : DATE +Simulator : NEURON +Backend : C++ (api-compatibility) +NMODL Compiler : VERSION +*********************************************************/ + +#include +#include +#include + +#include "mech_api.h" +#include "neuron/cache/mechanism_range.hpp" +#include "nrniv_mf.h" +#include "section_fwd.hpp" + +/* NEURON global macro definitions */ +/* VECTORIZED */ +#define NRN_VECTORIZED 1 + +static constexpr auto number_of_datum_variables = 6; +static constexpr auto number_of_floating_point_variables = 23; + +namespace { +template +using _nrn_mechanism_std_vector = std::vector; +using _nrn_model_sorted_token = neuron::model_sorted_token; +using _nrn_mechanism_cache_range = neuron::cache::MechanismRange; +using _nrn_mechanism_cache_instance = neuron::cache::MechanismInstance; +using _nrn_non_owning_id_without_container = neuron::container::non_owning_identifier_without_container; +template +using _nrn_mechanism_field = neuron::mechanism::field; +template +void _nrn_mechanism_register_data_fields(Args&&... args) { + neuron::mechanism::register_data_fields(std::forward(args)...); +} +} // namespace + +Prop* hoc_getdata_range(int type); +extern double celsius; + + +namespace neuron { + #ifndef NRN_PRCELLSTATE + #define NRN_PRCELLSTATE 0 + #endif + + + /** channel information */ + static const char *mechanism_info[] = { + "7.7.0", + "hodhux", + "gnabar_hodhux", + "gkbar_hodhux", + "gl_hodhux", + "el_hodhux", + 0, + "il_hodhux", + "minf_hodhux", + "hinf_hodhux", + "ninf_hodhux", + "mexp_hodhux", + "hexp_hodhux", + "nexp_hodhux", + 0, + "m_hodhux", + "h_hodhux", + "n_hodhux", + 0, + 0 + }; + + + /* NEURON global variables */ + static Symbol* _na_sym; + static Symbol* _k_sym; + static int mech_type; + static Prop* _extcall_prop; + /* _prop_id kind of shadows _extcall_prop to allow validity checking. */ + static _nrn_non_owning_id_without_container _prop_id{}; + static int hoc_nrnpointerindex = -1; + static _nrn_mechanism_std_vector _extcall_thread; + + + /** all global variables */ + struct hodhux_Store { + double m0{}; + double h0{}; + double n0{}; + }; + static_assert(std::is_trivially_copy_constructible_v); + static_assert(std::is_trivially_move_constructible_v); + static_assert(std::is_trivially_copy_assignable_v); + static_assert(std::is_trivially_move_assignable_v); + static_assert(std::is_trivially_destructible_v); + hodhux_Store hodhux_global; + + + /** all mechanism instance variables and global variables */ + struct hodhux_Instance { + double* celsius{&::celsius}; + double* gnabar{}; + double* gkbar{}; + double* gl{}; + double* el{}; + double* il{}; + double* minf{}; + double* hinf{}; + double* ninf{}; + double* mexp{}; + double* hexp{}; + double* nexp{}; + double* m{}; + double* h{}; + double* n{}; + double* ena{}; + double* ek{}; + double* Dm{}; + double* Dh{}; + double* Dn{}; + double* ina{}; + double* ik{}; + double* v_unused{}; + double* g_unused{}; + const double* const* ion_ena{}; + double* const* ion_ina{}; + double* const* ion_dinadv{}; + const double* const* ion_ek{}; + double* const* ion_ik{}; + double* const* ion_dikdv{}; + hodhux_Store* global{&hodhux_global}; + }; + + + struct hodhux_NodeData { + int const * nodeindices; + double const * node_voltages; + double * node_diagonal; + double * node_rhs; + int nodecount; + }; + + + static hodhux_Instance make_instance_hodhux(_nrn_mechanism_cache_range& _ml) { + return hodhux_Instance { + &::celsius, + _ml.template fpfield_ptr<0>(), + _ml.template fpfield_ptr<1>(), + _ml.template fpfield_ptr<2>(), + _ml.template fpfield_ptr<3>(), + _ml.template fpfield_ptr<4>(), + _ml.template fpfield_ptr<5>(), + _ml.template fpfield_ptr<6>(), + _ml.template fpfield_ptr<7>(), + _ml.template fpfield_ptr<8>(), + _ml.template fpfield_ptr<9>(), + _ml.template fpfield_ptr<10>(), + _ml.template fpfield_ptr<11>(), + _ml.template fpfield_ptr<12>(), + _ml.template fpfield_ptr<13>(), + _ml.template fpfield_ptr<14>(), + _ml.template fpfield_ptr<15>(), + _ml.template fpfield_ptr<16>(), + _ml.template fpfield_ptr<17>(), + _ml.template fpfield_ptr<18>(), + _ml.template fpfield_ptr<19>(), + _ml.template fpfield_ptr<20>(), + _ml.template fpfield_ptr<21>(), + _ml.template fpfield_ptr<22>(), + _ml.template dptr_field_ptr<0>(), + _ml.template dptr_field_ptr<1>(), + _ml.template dptr_field_ptr<2>(), + _ml.template dptr_field_ptr<3>(), + _ml.template dptr_field_ptr<4>(), + _ml.template dptr_field_ptr<5>() + }; + } + + + static hodhux_NodeData make_node_data_hodhux(NrnThread& _nt, Memb_list& _ml_arg) { + return hodhux_NodeData { + _ml_arg.nodeindices, + _nt.node_voltage_storage(), + _nt.node_d_storage(), + _nt.node_rhs_storage(), + _ml_arg.nodecount + }; + } + + + static void nrn_alloc_hodhux(Prop* _prop) { + Prop *prop_ion{}; + Datum *_ppvar{}; + _ppvar = nrn_prop_datum_alloc(mech_type, 6, _prop); + _nrn_mechanism_access_dparam(_prop) = _ppvar; + _nrn_mechanism_cache_instance _ml_real{_prop}; + auto* const _ml = &_ml_real; + size_t const _iml{}; + assert(_nrn_mechanism_get_num_vars(_prop) == 23); + /*initialize range parameters*/ + _ml->template fpfield<0>(_iml) = 0.12; /* gnabar */ + _ml->template fpfield<1>(_iml) = 0.036; /* gkbar */ + _ml->template fpfield<2>(_iml) = 0.0003; /* gl */ + _ml->template fpfield<3>(_iml) = -54.3; /* el */ + _nrn_mechanism_access_dparam(_prop) = _ppvar; + Symbol * na_sym = hoc_lookup("na_ion"); + Prop * na_prop = need_memb(na_sym); + _ppvar[0] = _nrn_mechanism_get_param_handle(na_prop, 0); + _ppvar[1] = _nrn_mechanism_get_param_handle(na_prop, 3); + _ppvar[2] = _nrn_mechanism_get_param_handle(na_prop, 4); + Symbol * k_sym = hoc_lookup("k_ion"); + Prop * k_prop = need_memb(k_sym); + _ppvar[3] = _nrn_mechanism_get_param_handle(k_prop, 0); + _ppvar[4] = _nrn_mechanism_get_param_handle(k_prop, 3); + _ppvar[5] = _nrn_mechanism_get_param_handle(k_prop, 4); + } + + + /* Neuron setdata functions */ + extern void _nrn_setdata_reg(int, void(*)(Prop*)); + static void _setdata(Prop* _prop) { + _extcall_prop = _prop; + _prop_id = _nrn_get_prop_id(_prop); + } + static void _hoc_setdata() { + Prop *_prop = hoc_getdata_range(mech_type); + _setdata(_prop); + hoc_retpushx(1.); + } + /* Mechanism procedures and functions */ + inline double vtrap_hodhux(_nrn_mechanism_cache_range* _ml, hodhux_Instance& inst, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* _nt, double x, double y); + inline int states_hodhux(_nrn_mechanism_cache_range* _ml, hodhux_Instance& inst, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* _nt); + inline int rates_hodhux(_nrn_mechanism_cache_range* _ml, hodhux_Instance& inst, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* _nt, double v); + + + /** connect global (scalar) variables to hoc -- */ + static DoubScal hoc_scalar_double[] = { + {nullptr, nullptr} + }; + + + /** connect global (array) variables to hoc -- */ + static DoubVec hoc_vector_double[] = { + {nullptr, nullptr, 0} + }; + + + /* declaration of user functions */ + static void _hoc_states(void); + static void _hoc_rates(void); + static void _hoc_vtrap(void); + static double _npy_states(Prop*); + static double _npy_rates(Prop*); + static double _npy_vtrap(Prop*); + + + /* connect user functions to hoc names */ + static VoidFunc hoc_intfunc[] = { + {"setdata_hodhux", _hoc_setdata}, + {"states_hodhux", _hoc_states}, + {"rates_hodhux", _hoc_rates}, + {"vtrap_hodhux", _hoc_vtrap}, + {0, 0} + }; + static NPyDirectMechFunc npy_direct_func_proc[] = { + {"states", _npy_states}, + {"rates", _npy_rates}, + {"vtrap", _npy_vtrap}, + }; + static void _hoc_states(void) { + double _r{}; + Datum* _ppvar; + Datum* _thread; + NrnThread* _nt; + if (!_prop_id) { + hoc_execerror("No data for states_hodhux. Requires prior call to setdata_hodhux and that the specified mechanism instance still be in existence.", NULL); + } + Prop* _local_prop = _extcall_prop; + _nrn_mechanism_cache_instance _ml_real{_local_prop}; + auto* const _ml = &_ml_real; + size_t const id{}; + _ppvar = _local_prop ? _nrn_mechanism_access_dparam(_local_prop) : nullptr; + _thread = _extcall_thread.data(); + _nt = nrn_threads; + auto inst = make_instance_hodhux(_ml_real); + _r = 1.; + states_hodhux(_ml, inst, id, _ppvar, _thread, _nt); + hoc_retpushx(_r); + } + static double _npy_states(Prop* _prop) { + double _r{}; + Datum* _ppvar; + Datum* _thread; + NrnThread* _nt; + _nrn_mechanism_cache_instance _ml_real{_prop}; + auto* const _ml = &_ml_real; + size_t const id{}; + _ppvar = _nrn_mechanism_access_dparam(_prop); + _thread = _extcall_thread.data(); + _nt = nrn_threads; + auto inst = make_instance_hodhux(_ml_real); + _r = 1.; + states_hodhux(_ml, inst, id, _ppvar, _thread, _nt); + return(_r); + } + static void _hoc_rates(void) { + double _r{}; + Datum* _ppvar; + Datum* _thread; + NrnThread* _nt; + if (!_prop_id) { + hoc_execerror("No data for rates_hodhux. Requires prior call to setdata_hodhux and that the specified mechanism instance still be in existence.", NULL); + } + Prop* _local_prop = _extcall_prop; + _nrn_mechanism_cache_instance _ml_real{_local_prop}; + auto* const _ml = &_ml_real; + size_t const id{}; + _ppvar = _local_prop ? _nrn_mechanism_access_dparam(_local_prop) : nullptr; + _thread = _extcall_thread.data(); + _nt = nrn_threads; + auto inst = make_instance_hodhux(_ml_real); + _r = 1.; + rates_hodhux(_ml, inst, id, _ppvar, _thread, _nt, *getarg(1)); + hoc_retpushx(_r); + } + static double _npy_rates(Prop* _prop) { + double _r{}; + Datum* _ppvar; + Datum* _thread; + NrnThread* _nt; + _nrn_mechanism_cache_instance _ml_real{_prop}; + auto* const _ml = &_ml_real; + size_t const id{}; + _ppvar = _nrn_mechanism_access_dparam(_prop); + _thread = _extcall_thread.data(); + _nt = nrn_threads; + auto inst = make_instance_hodhux(_ml_real); + _r = 1.; + rates_hodhux(_ml, inst, id, _ppvar, _thread, _nt, *getarg(1)); + return(_r); + } + static void _hoc_vtrap(void) { + double _r{}; + Datum* _ppvar; + Datum* _thread; + NrnThread* _nt; + Prop* _local_prop = _prop_id ? _extcall_prop : nullptr; + _nrn_mechanism_cache_instance _ml_real{_local_prop}; + auto* const _ml = &_ml_real; + size_t const id{}; + _ppvar = _local_prop ? _nrn_mechanism_access_dparam(_local_prop) : nullptr; + _thread = _extcall_thread.data(); + _nt = nrn_threads; + auto inst = make_instance_hodhux(_ml_real); + _r = vtrap_hodhux(_ml, inst, id, _ppvar, _thread, _nt, *getarg(1), *getarg(2)); + hoc_retpushx(_r); + } + static double _npy_vtrap(Prop* _prop) { + double _r{}; + Datum* _ppvar; + Datum* _thread; + NrnThread* _nt; + _nrn_mechanism_cache_instance _ml_real{_prop}; + auto* const _ml = &_ml_real; + size_t const id{}; + _ppvar = _nrn_mechanism_access_dparam(_prop); + _thread = _extcall_thread.data(); + _nt = nrn_threads; + auto inst = make_instance_hodhux(_ml_real); + _r = vtrap_hodhux(_ml, inst, id, _ppvar, _thread, _nt, *getarg(1), *getarg(2)); + return(_r); + } + + + inline int states_hodhux(_nrn_mechanism_cache_range* _ml, hodhux_Instance& inst, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* _nt) { + int ret_states = 0; + auto v = inst.v_unused[id]; + rates_hodhux(_ml, inst, id, _ppvar, _thread, _nt, v); + inst.m[id] = inst.m[id] + inst.mexp[id] * (inst.minf[id] - inst.m[id]); + inst.h[id] = inst.h[id] + inst.hexp[id] * (inst.hinf[id] - inst.h[id]); + inst.n[id] = inst.n[id] + inst.nexp[id] * (inst.ninf[id] - inst.n[id]); + return ret_states; + } + + + inline int rates_hodhux(_nrn_mechanism_cache_range* _ml, hodhux_Instance& inst, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* _nt, double v) { + int ret_rates = 0; + double q10, tinc, alpha, beta, sum; + q10 = pow(3.0, ((*(inst.celsius) - 6.3) / 10.0)); + tinc = -_nt->_dt * q10; + alpha = .1 * vtrap_hodhux(_ml, inst, id, _ppvar, _thread, _nt, -(v + 40.0), 10.0); + beta = 4.0 * exp( -(v + 65.0) / 18.0); + sum = alpha + beta; + inst.minf[id] = alpha / sum; + inst.mexp[id] = 1.0 - exp(tinc * sum); + alpha = .07 * exp( -(v + 65.0) / 20.0); + beta = 1.0 / (exp( -(v + 35.0) / 10.0) + 1.0); + sum = alpha + beta; + inst.hinf[id] = alpha / sum; + inst.hexp[id] = 1.0 - exp(tinc * sum); + alpha = .01 * vtrap_hodhux(_ml, inst, id, _ppvar, _thread, _nt, -(v + 55.0), 10.0); + beta = .125 * exp( -(v + 65.0) / 80.0); + sum = alpha + beta; + inst.ninf[id] = alpha / sum; + inst.nexp[id] = 1.0 - exp(tinc * sum); + return ret_rates; + } + + + inline double vtrap_hodhux(_nrn_mechanism_cache_range* _ml, hodhux_Instance& inst, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* _nt, double x, double y) { + double ret_vtrap = 0.0; + auto v = inst.v_unused[id]; + if (fabs(x / y) < 1e-6) { + ret_vtrap = y * (1.0 - x / y / 2.0); + } else { + ret_vtrap = x / (exp(x / y) - 1.0); + } + return ret_vtrap; + } + + + void nrn_init_hodhux(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_hodhux(_lmr); + auto node_data = make_node_data_hodhux(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + auto* const _ml = &_lmr; + auto* _thread = _ml_arg->_thread; + for (int id = 0; id < nodecount; id++) { + + int node_id = node_data.nodeindices[id]; + auto* _ppvar = _ml_arg->pdata[id]; + auto v = node_data.node_voltages[node_id]; + inst.v_unused[id] = v; + inst.ena[id] = (*inst.ion_ena[id]); + inst.ek[id] = (*inst.ion_ek[id]); + rates_hodhux(_ml, inst, id, _ppvar, _thread, _nt, v); + inst.m[id] = inst.minf[id]; + inst.h[id] = inst.hinf[id]; + inst.n[id] = inst.ninf[id]; + } + } + + + inline double nrn_current_hodhux(size_t id, hodhux_Instance& inst, hodhux_NodeData& node_data, double v) { + double current = 0.0; + inst.ina[id] = inst.gnabar[id] * inst.m[id] * inst.m[id] * inst.m[id] * inst.h[id] * (v - inst.ena[id]); + inst.ik[id] = inst.gkbar[id] * inst.n[id] * inst.n[id] * inst.n[id] * inst.n[id] * (v - inst.ek[id]); + inst.il[id] = inst.gl[id] * (v - inst.el[id]); + current += inst.il[id]; + current += inst.ina[id]; + current += inst.ik[id]; + return current; + } + + + /** update current */ + void nrn_cur_hodhux(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_hodhux(_lmr); + auto node_data = make_node_data_hodhux(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + auto* const _ml = &_lmr; + auto* _thread = _ml_arg->_thread; + for (int id = 0; id < nodecount; id++) { + int node_id = node_data.nodeindices[id]; + double v = node_data.node_voltages[node_id]; + inst.ena[id] = (*inst.ion_ena[id]); + inst.ek[id] = (*inst.ion_ek[id]); + double I1 = nrn_current_hodhux(id, inst, node_data, v+0.001); + double dina = inst.ina[id]; + double dik = inst.ik[id]; + double I0 = nrn_current_hodhux(id, inst, node_data, v); + double rhs = I0; + double g = (I1-I0)/0.001; + (*inst.ion_dinadv[id]) += (dina-inst.ina[id])/0.001; + (*inst.ion_dikdv[id]) += (dik-inst.ik[id])/0.001; + (*inst.ion_ina[id]) += inst.ina[id]; + (*inst.ion_ik[id]) += inst.ik[id]; + node_data.node_rhs[node_id] -= rhs; + // remember the conductances so we can set them later + inst.g_unused[id] = g; + } + } + + + void nrn_state_hodhux(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_hodhux(_lmr); + auto node_data = make_node_data_hodhux(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + auto* const _ml = &_lmr; + auto* _thread = _ml_arg->_thread; + for (int id = 0; id < nodecount; id++) { + + int node_id = node_data.nodeindices[id]; + auto* _ppvar = _ml_arg->pdata[id]; + auto v = node_data.node_voltages[node_id]; + inst.ena[id] = (*inst.ion_ena[id]); + inst.ek[id] = (*inst.ion_ek[id]); + rates_hodhux(_ml, inst, id, _ppvar, _thread, _nt, v); + inst.m[id] = inst.m[id] + inst.mexp[id] * (inst.minf[id] - inst.m[id]); + inst.h[id] = inst.h[id] + inst.hexp[id] * (inst.hinf[id] - inst.h[id]); + inst.n[id] = inst.n[id] + inst.nexp[id] * (inst.ninf[id] - inst.n[id]); + } + } + + + /** nrn_jacob function */ + static void nrn_jacob_hodhux(_nrn_model_sorted_token const& _sorted_token, NrnThread* _nt, Memb_list* _ml_arg, int _type) { + _nrn_mechanism_cache_range _lmr{_sorted_token, *_nt, *_ml_arg, _type}; + auto inst = make_instance_hodhux(_lmr); + auto node_data = make_node_data_hodhux(*_nt, *_ml_arg); + auto nodecount = _ml_arg->nodecount; + for (int id = 0; id < nodecount; id++) { + // set conductances properly + int node_id = node_data.nodeindices[id]; + node_data.node_diagonal[node_id] += inst.g_unused[id]; + } + } + + + static void _initlists() { + } + + + /** register channel with the simulator */ + extern "C" void _hodhux_reg() { + _initlists(); + + ion_reg("na", -10000.); + ion_reg("k", -10000.); + + _na_sym = hoc_lookup("na_ion"); + _k_sym = hoc_lookup("k_ion"); + + register_mech(mechanism_info, nrn_alloc_hodhux, nrn_cur_hodhux, nrn_jacob_hodhux, nrn_state_hodhux, nrn_init_hodhux, hoc_nrnpointerindex, 1); + + mech_type = nrn_get_mechtype(mechanism_info[1]); + _nrn_mechanism_register_data_fields(mech_type, + _nrn_mechanism_field{"gnabar"} /* 0 */, + _nrn_mechanism_field{"gkbar"} /* 1 */, + _nrn_mechanism_field{"gl"} /* 2 */, + _nrn_mechanism_field{"el"} /* 3 */, + _nrn_mechanism_field{"il"} /* 4 */, + _nrn_mechanism_field{"minf"} /* 5 */, + _nrn_mechanism_field{"hinf"} /* 6 */, + _nrn_mechanism_field{"ninf"} /* 7 */, + _nrn_mechanism_field{"mexp"} /* 8 */, + _nrn_mechanism_field{"hexp"} /* 9 */, + _nrn_mechanism_field{"nexp"} /* 10 */, + _nrn_mechanism_field{"m"} /* 11 */, + _nrn_mechanism_field{"h"} /* 12 */, + _nrn_mechanism_field{"n"} /* 13 */, + _nrn_mechanism_field{"ena"} /* 14 */, + _nrn_mechanism_field{"ek"} /* 15 */, + _nrn_mechanism_field{"Dm"} /* 16 */, + _nrn_mechanism_field{"Dh"} /* 17 */, + _nrn_mechanism_field{"Dn"} /* 18 */, + _nrn_mechanism_field{"ina"} /* 19 */, + _nrn_mechanism_field{"ik"} /* 20 */, + _nrn_mechanism_field{"v_unused"} /* 21 */, + _nrn_mechanism_field{"g_unused"} /* 22 */, + _nrn_mechanism_field{"ion_ena", "na_ion"} /* 0 */, + _nrn_mechanism_field{"ion_ina", "na_ion"} /* 1 */, + _nrn_mechanism_field{"ion_dinadv", "na_ion"} /* 2 */, + _nrn_mechanism_field{"ion_ek", "k_ion"} /* 3 */, + _nrn_mechanism_field{"ion_ik", "k_ion"} /* 4 */, + _nrn_mechanism_field{"ion_dikdv", "k_ion"} /* 5 */ + ); + + hoc_register_prop_size(mech_type, 23, 6); + hoc_register_dparam_semantics(mech_type, 0, "na_ion"); + hoc_register_dparam_semantics(mech_type, 1, "na_ion"); + hoc_register_dparam_semantics(mech_type, 2, "na_ion"); + hoc_register_dparam_semantics(mech_type, 3, "k_ion"); + hoc_register_dparam_semantics(mech_type, 4, "k_ion"); + hoc_register_dparam_semantics(mech_type, 5, "k_ion"); + hoc_register_var(hoc_scalar_double, hoc_vector_double, hoc_intfunc); + hoc_register_npy_direct(mech_type, npy_direct_func_proc); + } +}