From fd6cb7fe759dd6a709294366e34fbf905eb83530 Mon Sep 17 00:00:00 2001 From: Subhasis Ray Date: Sat, 7 Dec 2024 10:24:28 +0530 Subject: [PATCH 1/6] Fixes in docs and HDF5 support - Fixed issue #495 - Enabled HDF5 supported build and fixed HDF5 linking issue - Updated installation instructions, README. - Updated packaging to upload wheels only on master branch --- .github/workflows/package.yml | 3 ++- WindowsBuild.md | 8 ++++---- meson.build | 3 ++- meson_options.txt | 2 +- 4 files changed, 9 insertions(+), 7 deletions(-) diff --git a/.github/workflows/package.yml b/.github/workflows/package.yml index a752800cee..611061e316 100644 --- a/.github/workflows/package.yml +++ b/.github/workflows/package.yml @@ -8,7 +8,7 @@ jobs: strategy: fail-fast: false matrix: - os: [ubuntu-latest, macos-14, windows-latest] + os: [ubuntu-latest, macos-14, macos-15, windows-latest] build_type: [Release] c_compiler: [clang] python-version: ['3.11'] @@ -105,6 +105,7 @@ jobs: python -m cibuildwheel --output-dir wheelhouse dir wheelhouse - name: Upload packages + if: github.ref == 'refs/heads/main' uses: xresloader/upload-to-github-release@v1 env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/WindowsBuild.md b/WindowsBuild.md index bda4c5dff5..377c14a06e 100644 --- a/WindowsBuild.md +++ b/WindowsBuild.md @@ -53,8 +53,8 @@ $env:PATH="$env:PATH;C:\Program Files\LLVM\bin" ``` cd moose-core -meson setup --wipe _build -Duse_mpi=false --buildtype=release -ninja -v -C _build +meson setup --wipe _build -Duse_mpi=false --buildtype=release --vsenv +meson compile -vC _build meson install -C _build ``` *Note: If you want to keep the installation in local folder replace the third command with: `meson setup --wipe _build --prefix={absolute-path} -Duse_mpi=false --buildtype=release` where `absolute-path` is the full path of the folder you want to install it in. You will have to add the `Lib\site-packages` directory within this to your PATH environment variable to make `moose` module visible to Python.* @@ -68,7 +68,7 @@ Meson provides many builtin options: https://mesonbuild.com/Builtin-options.html ``` - meson setup --wipe _build --prefix=%CD%\\_build_install -Duse_mpi=false -Dbuildtype=debug + meson setup --wipe _build --prefix=%CD%\\_build_install -Duse_mpi=false -Dbuildtype=debug --vsenv ``` You can either use `buildtype` option alone or use the two options `debug` and `optimization` for finer grained control over the build. According to `meson` documentation `-Dbuildtype=debug` will create a debug build with optimization level 0 (i.e., no optimization, passing `-O0 -g` to GCC), `-Dbuildtype=debugoptimized` will create a debug build with optimization level 2 (equivalent to `-Ddebug=true -Doptimization=2`), `-Dbuildtype=release` will create a release build with optimization level 3 (equivalent to `-Ddebug=false -Doptimization=3`), and `-Dbuildtype=minsize` will create a release build with space optimization (passing `-Os` to GCC). @@ -91,7 +91,7 @@ is to make a copy of `python3x.lib` named `python3x_d.lib` in the same directory (`libs`). After that, you can run meson setup as follows: ``` -meson setup --wipe _build --prefix=%CD%\\_build_install -Duse_mpi=false --buildtype=debug +meson setup --wipe _build --prefix=%CD%\\_build_install -Duse_mpi=false --buildtype=debug --vsenv ``` and then go through the rest of the steps. diff --git a/meson.build b/meson.build index 0f1075c51a..e0d617bed3 100644 --- a/meson.build +++ b/meson.build @@ -23,8 +23,9 @@ gsl_dep = dependency('gsl', version: '>=1.16') # HDF5 cannot be found via pkg-config. Skip for now. TODO: fix this in next release use_hdf5 = get_option('use_hdf5') if use_hdf5 - hdf5_dep = dependency('hdf5_cpp', version: '>=1.8') + hdf5_dep = dependency('hdf5', version: '>=1.8') add_global_arguments('-DUSE_HDF5', language: ['c', 'cpp']) + add_global_arguments('-DH5_BUILT_AS_DYNAMIC_LIB', language: ['c', 'cpp']) # to fix this issue: https://github.com/conda-forge/hdf5-feedstock/issues/58 else hdf5_dep = [] endif diff --git a/meson_options.txt b/meson_options.txt index c05bfed36e..e6cd2b2f30 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -1,4 +1,4 @@ option('use_mpi', type: 'boolean', value: false, description: 'If specified, build with MPI support') -option('use_hdf5', type: 'boolean', value: false, +option('use_hdf5', type: 'boolean', value: true, description: 'If specified, build with HDF5 support. Needed for NSDF.') From 19547cc66c4740571b4797baaa84aadaf64fd7fd Mon Sep 17 00:00:00 2001 From: Subhasis Ray Date: Tue, 10 Dec 2024 16:03:30 +0530 Subject: [PATCH 2/6] Fix for https://github.com/mesonbuild/meson/issues/13985 Fixes build issue on Ubuntu --- UbuntuBuild.md | 2 +- builtins/meson.build | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/UbuntuBuild.md b/UbuntuBuild.md index 6776de9225..5d3ccf9280 100644 --- a/UbuntuBuild.md +++ b/UbuntuBuild.md @@ -47,7 +47,7 @@ conda activate moose cd moose-core meson setup --wipe _build --prefix=`pwd`/_build_install -Duse_mpi=false -Dbuildtype=release -ninja -v -C _build +meson compile -v -C _build meson install -C _build ``` diff --git a/builtins/meson.build b/builtins/meson.build index 077c4fb6f6..31a4e18bba 100644 --- a/builtins/meson.build +++ b/builtins/meson.build @@ -31,7 +31,10 @@ if host_machine.system() != 'windows' builtins_src += files(['SocketStreamer.cpp']) endif include_directories = ['../external/fmt/include'] + if use_hdf5 - include_directories += hdf5_dep.get_variable(pkgconfig: 'includedir') + builtins_lib = static_library('builtins', builtins_src, include_directories: include_directories, dependencies: hdf5_dep) +else + builtins_lib = static_library('builtins', builtins_src, include_directories: include_directories) endif -builtins_lib = static_library('builtins', builtins_src, include_directories: include_directories) + From c123aaaa4cb8e6902d3b2b99cc5a617cd2aafd88 Mon Sep 17 00:00:00 2001 From: Subhasis Ray Date: Mon, 16 Dec 2024 10:50:36 +0530 Subject: [PATCH 3/6] Raise exception in `connect` and `element` when failed - The shell only generates an error message when an attempt to connect two objects via message fails. This caused hard to debug errors. Now the wrapper function raises RuntimeError in such case. - `element(path)` returned the roor element when `path` does not exist, again causing hard to find bugs. This now raises RuntimeError. --- README.md | 2 +- python/moose/__init__.py | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 63de8c02b0..a7ea97b82c 100644 --- a/README.md +++ b/README.md @@ -41,7 +41,7 @@ https://github.com/BhallaLab/moose-examples. # Build To build `pymoose`, follow instructions given in -[INSTALLATION.md](INSTALLATION.md) and for platform specific +[INSTALL.md](INSTALL.md) and for platform specific information see: - Linux: [UbuntuBuild.md](UbuntuBuild.md) - MacOSX: [AppleM1Build.md](AppleM1Build.md) diff --git a/python/moose/__init__.py b/python/moose/__init__.py index a0f726f55a..e90976c2e8 100644 --- a/python/moose/__init__.py +++ b/python/moose/__init__.py @@ -204,7 +204,10 @@ def connect(src, srcfield, dest, destfield, msgtype="Single"): """ src = _moose.element(src) dest = _moose.element(dest) - return src.connect(srcfield, dest, destfield, msgtype) + msg = src.connect(srcfield, dest, destfield, msgtype) + if msg.name == '/': + raise RuntimeError(f'Could not connect {src}.{srcfield} with {dest}.{destfield}') + return msg def delete(arg): @@ -239,6 +242,8 @@ def element(arg): MOOSE element (object) corresponding to the `arg` converted to write subclass. """ + if not _moose.exists(arg): + raise RuntimeError(f'{arg}: element at path does not exist') return _moose.element(arg) From affaca718dceab714220eab67616b020e58a69ad Mon Sep 17 00:00:00 2001 From: Subhasis Ray Date: Mon, 16 Dec 2024 11:07:21 +0530 Subject: [PATCH 4/6] Minor fix in installation instructions --- INSTALL.md | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/INSTALL.md b/INSTALL.md index e00e3a7508..4d6efb0fbd 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -43,9 +43,9 @@ To build `MOOSE` from source, you need `meson`, `ninja`, `pkg-config`, and `pyth recommend to use Python 3.9 or higher. To build and install with `pip`, you also need `meson-python`. For platform specific instructions, see: -- Linux: UbuntuBuild.md -- MacOSX: AppleM1Build.md -- Windows: WindowsBuild.md +- Linux: [UbuntuBuild.md](UbuntuBuild.md) +- MacOSX: [AppleM1Build.md](AppleM1Build.md) +- Windows: [WindowsBuild.md](WindowsBuild.md) Briefly, 1. Build environment: make sure that the following packages are installed on your system. @@ -68,6 +68,14 @@ Briefly, $ pip install git+https://github.com/BhallaLab/moose-core --user ``` +This will install the `master` branch of MOOSE. If you want a specific fork/branch, modify the address: + +``` +$ pip install git+https://github.com/subhacom/moose-core@fix495merge --user +``` + +will install the fix495merge branch from `subhacom`'s fork of `moose-core`. + ## Post installation You can check that moose is installed and initializes correctly by running: From 3c719ad5a29601b6bf01239e4c895d0c5ccb3780 Mon Sep 17 00:00:00 2001 From: Subhasis Ray Date: Mon, 6 Jan 2025 16:25:22 +0530 Subject: [PATCH 5/6] Updated moose.element() to raise error for nonexistent path --- python/moose/__init__.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/python/moose/__init__.py b/python/moose/__init__.py index e90976c2e8..f1ebf9ea57 100644 --- a/python/moose/__init__.py +++ b/python/moose/__init__.py @@ -12,7 +12,13 @@ # Notes # ----- -# Use these guidelines for docstring: https://numpydoc.readthedocs.io/en/latest/format.html +# +# 1. Use these guidelines for docstring: +# https://numpydoc.readthedocs.io/en/latest/format.html. +# +# 2. We redefine many functions defined in _moose just to add the +# docstring since Python C-API does not provide a way to add docstring +# to a function defined in the C/C++ extension import sys import pydoc @@ -242,7 +248,7 @@ def element(arg): MOOSE element (object) corresponding to the `arg` converted to write subclass. """ - if not _moose.exists(arg): + if isinstance(arg, str) and not _moose.exists(arg): raise RuntimeError(f'{arg}: element at path does not exist') return _moose.element(arg) From db0730facdec25df5ab44c51d341d953120209ac Mon Sep 17 00:00:00 2001 From: Subhasis Ray Date: Fri, 24 Jan 2025 14:05:03 +0530 Subject: [PATCH 6/6] Implemented HHChannelF and HHGateF for formula based evaluation - Refactored HHChannelBase and created a class hierarchy for HHGate - Added test scripts for HHChannelF (`test_hhchannelf.py`) and updated test for HHChannel (`test_hhchannel.py`) - Moved voltage clamp code for test fixture into `tests/core/ephys.py` - Removed old `Func` class from default scheduling --- biophysics/HHChannel.cpp | 333 ++++++------ biophysics/HHChannel.h | 124 ++--- biophysics/HHChannel2D.cpp | 991 +++++++++++++++-------------------- biophysics/HHChannel2D.h | 362 ++++++------- biophysics/HHChannelBase.cpp | 352 +++++++------ biophysics/HHChannelBase.h | 376 +++++++------ biophysics/HHChannelF.cpp | 295 +++++++++++ biophysics/HHChannelF.h | 124 +++++ biophysics/HHGate.cpp | 43 +- biophysics/HHGate.h | 34 +- biophysics/HHGate2D.cpp | 5 +- biophysics/HHGate2D.h | 4 +- biophysics/HHGateBase.cpp | 62 +++ biophysics/HHGateBase.h | 87 +++ biophysics/HHGateF.cpp | 227 ++++++++ biophysics/HHGateF.h | 119 +++++ biophysics/meson.build | 10 +- builtins/Function.cpp | 14 +- builtins/Function.h | 4 +- hsolve/ZombieHHChannel.cpp | 28 +- hsolve/ZombieHHChannel.h | 28 +- pybind11/Finfo.cpp | 2 +- scheduling/Clock.cpp | 6 +- tests/core/ephys.py | 113 ++++ tests/core/test_hhchan.py | 228 ++++++++ tests/core/test_hhchanf.py | 228 ++++++++ 26 files changed, 2742 insertions(+), 1457 deletions(-) create mode 100644 biophysics/HHChannelF.cpp create mode 100644 biophysics/HHChannelF.h create mode 100644 biophysics/HHGateBase.cpp create mode 100644 biophysics/HHGateBase.h create mode 100644 biophysics/HHGateF.cpp create mode 100644 biophysics/HHGateF.h create mode 100644 tests/core/ephys.py create mode 100644 tests/core/test_hhchan.py create mode 100644 tests/core/test_hhchanf.py diff --git a/biophysics/HHChannel.cpp b/biophysics/HHChannel.cpp index 942d974ac5..e87ff937ed 100644 --- a/biophysics/HHChannel.cpp +++ b/biophysics/HHChannel.cpp @@ -7,33 +7,59 @@ ** See the file COPYING.LIB for the full notice. **********************************************************************/ -#include "../basecode/header.h" -#include "../basecode/ElementValueFinfo.h" -#include "HHGate.h" -#include "ChanBase.h" +// #include "../basecode/ElementValueFinfo.h" +// #include "../basecode/header.h" +// #include "ChanBase.h" #include "ChanCommon.h" #include "HHChannelBase.h" #include "HHChannel.h" +#include "HHGate.h" +// const double HHChannel::EPSILON = 1.0e-10; +// const int HHChannel::INSTANT_X = 1; +// const int HHChannel::INSTANT_Y = 2; +// const int HHChannel::INSTANT_Z = 4; + +const Cinfo* HHChannel::initCinfo() { + static FieldElementFinfo gateX( + "gateX", "Sets up HHGate X for channel", HHGate::initCinfo(), + &HHChannel::getXgate, &HHChannel::setNumGates, &HHChannel::getNumXgates + // 1 + ); + static FieldElementFinfo gateY( + "gateY", "Sets up HHGate Y for channel", HHGate::initCinfo(), + &HHChannel::getYgate, &HHChannel::setNumGates, &HHChannel::getNumYgates + // 1 + ); + static FieldElementFinfo gateZ( + "gateZ", "Sets up HHGate Z for channel", HHGate::initCinfo(), + &HHChannel::getZgate, &HHChannel::setNumGates, &HHChannel::getNumZgates + // 1 + ); + /////////////////////////////////////////////////////// + static Finfo* HHChannelFinfos[] = { + &gateX, // FieldElement + &gateY, // FieldElement + &gateZ // FieldElement + }; -const double HHChannel::EPSILON = 1.0e-10; -const int HHChannel::INSTANT_X = 1; -const int HHChannel::INSTANT_Y = 2; -const int HHChannel::INSTANT_Z = 4; - -const Cinfo* HHChannel::initCinfo() -{ /////////////////////////////////////////////////////// static string doc[] = { - "Name", "HHChannel", "Author", "Upinder S. Bhalla, 2007, NCBS", + "Name", + "HHChannel", + "Author", + "Upinder S. Bhalla, 2007, NCBS", "Description", "HHChannel: Hodgkin-Huxley type voltage-gated Ion channel. Something " "like the old tabchannel from GENESIS, but also presents " - "a similar interface as hhchan from GENESIS. ", }; + "a similar interface as hhchan from GENESIS. ", + }; static Dinfo dinfo; - static Cinfo HHChannelCinfo("HHChannel", HHChannelBase::initCinfo(), 0, 0, + static Cinfo HHChannelCinfo("HHChannel", HHChannelBase::initCinfo(), + HHChannelFinfos, + sizeof(HHChannelFinfos) / sizeof(Finfo*), &dinfo, doc, sizeof(doc) / sizeof(string)); return &HHChannelCinfo; @@ -45,27 +71,10 @@ static const Cinfo* hhChannelCinfo = HHChannel::initCinfo(); /////////////////////////////////////////////////// // Constructor /////////////////////////////////////////////////// -HHChannel::HHChannel() - : conc_(0.0), - instant_(0), - X_(0.0), - Y_(0.0), - Z_(0.0), - xInited_(false), - yInited_(false), - zInited_(false), - g_(0.0), - xGate_(0), - yGate_(0), - zGate_(0), - myId_() -{ - ; -} +HHChannel::HHChannel() : conc_(0.0), xGate_(nullptr), yGate_(nullptr), zGate_(nullptr) { ; } -HHChannel::~HHChannel() -{ - ; +HHChannel::~HHChannel() { + ; // if ( xGate_ && reinterpret_cast< char* >( this ) == // ObjId( xGate_->originalChannelId(), 0 ).data() ) // delete xGate_; @@ -77,44 +86,46 @@ HHChannel::~HHChannel() // delete zGate_; } -bool HHChannel::setGatePower(const Eref& e, double power, double* assignee, - const string& gateType) -{ - if (doubleEq(power, *assignee)) return false; - - if (doubleEq(*assignee, 0.0) && power > 0) { - createGate(e, gateType); - } else if (doubleEq(power, 0.0)) { - // destroyGate( e, gateType ); - } - *assignee = power; - - return true; -} - -/** - * Assigns the Xpower for this gate. If the gate exists and has - * only this element for input, then change the gate value. - * If the gate exists and has multiple parents, then make a new gate. - * If the gate does not exist, make a new gate - */ -void HHChannel::vSetXpower(const Eref& e, double power) -{ - if (setGatePower(e, power, &Xpower_, "X")) takeXpower_ = selectPower(power); -} - -void HHChannel::vSetYpower(const Eref& e, double power) -{ - if (setGatePower(e, power, &Ypower_, "Y")) takeYpower_ = selectPower(power); -} - -void HHChannel::vSetZpower(const Eref& e, double power) -{ - if (setGatePower(e, power, &Zpower_, "Z")) { - takeZpower_ = selectPower(power); - useConcentration_ = 1; // Not sure about this. - } -} +// bool HHChannel::setGatePower(const Eref& e, double power, double* assignee, +// const string& gateType) +// { +// if (doubleEq(power, *assignee)) return false; + +// if (doubleEq(*assignee, 0.0) && power > 0) { +// createGate(e, gateType); +// } else if (doubleEq(power, 0.0)) { +// // destroyGate( e, gateType ); +// } +// *assignee = power; + +// return true; +// } + +// /** +// * Assigns the Xpower for this gate. If the gate exists and has +// * only this element for input, then change the gate value. +// * If the gate exists and has multiple parents, then make a new gate. +// * If the gate does not exist, make a new gate +// */ +// void HHChannel::vSetXpower(const Eref& e, double power) +// { +// if (setGatePower(e, power, &Xpower_, "X")) takeXpower_ = +// selectPower(power); +// } + +// void HHChannel::vSetYpower(const Eref& e, double power) +// { +// if (setGatePower(e, power, &Ypower_, "Y")) takeYpower_ = +// selectPower(power); +// } + +// void HHChannel::vSetZpower(const Eref& e, double power) +// { +// if (setGatePower(e, power, &Zpower_, "Z")) { +// takeZpower_ = selectPower(power); +// useConcentration_ = 1; // Not sure about this. +// } +// } /** * If the gate exists and has only this element for input, then change @@ -134,22 +145,21 @@ void HHChannel::vSetZpower(const Eref& e, double power) // Assuming that the elements are simple elements. Use Eref for // general case -bool HHChannel::checkOriginal(Id chanId) const -{ - bool isOriginal = 1; - if (xGate_) { +bool HHChannel::checkOriginal(Id chanId) const { + bool isOriginal = true; + // cerr << "# HHChannel::checkOriginal(Id chanId) chanId: " << chanId << ", xGate: " << xGate_ << endl; + if (xGate_ != nullptr) { isOriginal = xGate_->isOriginalChannel(chanId); - } else if (yGate_) { + } else if (yGate_ != nullptr) { isOriginal = yGate_->isOriginalChannel(chanId); - } else if (zGate_) { + } else if (zGate_ != nullptr) { isOriginal = zGate_->isOriginalChannel(chanId); } return isOriginal; } void HHChannel::innerCreateGate(const string& gateName, HHGate** gatePtr, - Id chanId, Id gateId) -{ + Id chanId, Id gateId) { // Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() ); if (*gatePtr) { cout << "Warning: HHChannel::createGate: '" << gateName @@ -159,11 +169,11 @@ void HHChannel::innerCreateGate(const string& gateName, HHGate** gatePtr, *gatePtr = new HHGate(chanId, gateId); } -void HHChannel::vCreateGate(const Eref& e, string gateType) -{ +void HHChannel::vCreateGate(const Eref& e, string gateType) { if (!checkOriginal(e.id())) { cout << "Warning: HHChannel::createGate: Not allowed from copied " - "channel:\n" << e.id().path() << "\n"; + "channel:\n" + << e.id().path() << "\n"; return; } @@ -179,22 +189,21 @@ void HHChannel::vCreateGate(const Eref& e, string gateType) } void HHChannel::innerDestroyGate(const string& gateName, HHGate** gatePtr, - Id chanId) -{ - if (*gatePtr == 0) { + Id chanId) { + if (*gatePtr == nullptr) { cout << "Warning: HHChannel::destroyGate: '" << gateName << "' on Element '" << chanId.path() << "' not present\n"; return; } delete (*gatePtr); - *gatePtr = 0; + *gatePtr = nullptr; } -void HHChannel::destroyGate(const Eref& e, string gateType) -{ +void HHChannel::destroyGate(const Eref& e, string gateType) { if (!checkOriginal(e.id())) { cout << "Warning: HHChannel::destroyGate: Not allowed from copied " - "channel:\n" << e.id().path() << "\n"; + "channel:\n" + << e.id().path() << "\n"; return; } @@ -208,75 +217,83 @@ void HHChannel::destroyGate(const Eref& e, string gateType) cout << "Warning: HHChannel::destroyGate: Unknown gate type '" << gateType << "'. Ignored\n"; } - /////////////////////////////////////////////////// -// Field function definitions +// HHGate functions +// +// These are breaking the design as the return type is HHGate for +// HHChannel but HHGate2D for HHChannel2D. Making a common HHGateBase +// also turns out to be problematic as the field element can no longer +// be accessed as an HHGate or HHGate2D. /////////////////////////////////////////////////// -void HHChannel::vSetInstant(const Eref& e, int instant) -{ - instant_ = instant; -} +HHGate* HHChannel::getXgate(unsigned int i) { return xGate_; } -int HHChannel::vGetInstant(const Eref& e) const -{ - return instant_; -} +HHGate* HHChannel::getYgate(unsigned int i) { return yGate_; } -void HHChannel::vSetX(const Eref& e, double X) -{ - X_ = X; - xInited_ = true; -} -double HHChannel::vGetX(const Eref& e) const -{ - return X_; -} +HHGate* HHChannel::getZgate(unsigned int i) { return zGate_; } -void HHChannel::vSetY(const Eref& e, double Y) -{ - Y_ = Y; - yInited_ = true; -} -double HHChannel::vGetY(const Eref& e) const -{ - return Y_; -} +void HHChannel::setNumGates(unsigned int num) { ; } -void HHChannel::vSetZ(const Eref& e, double Z) -{ - Z_ = Z; - zInited_ = true; -} -double HHChannel::vGetZ(const Eref& e) const -{ - return Z_; -} +unsigned int HHChannel::getNumXgates() const { return xGate_ != nullptr; } -void HHChannel::vSetUseConcentration(const Eref& e, int value) -{ - useConcentration_ = value; -} +unsigned int HHChannel::getNumYgates() const { return yGate_ != nullptr; } + +unsigned int HHChannel::getNumZgates() const { return zGate_ != nullptr; } /////////////////////////////////////////////////// -// Dest function definitions +// Field function definitions /////////////////////////////////////////////////// -/** - * Returns the state variable for the new timestep based on - * the internal variables A_ and B_ which were passed in from the gate. - */ -double HHChannel::integrate(double state, double dt, double A, double B) -{ - if (B > EPSILON) { - double x = exp(-B * dt); - return state * x + (A / B) * (1 - x); - } - return state + A * dt; -} +// void HHChannel::vSetInstant(const Eref& e, int instant) +// { +// instant_ = instant; +// } + +// int HHChannel::vGetInstant(const Eref& e) const +// { +// return instant_; +// } + +// void HHChannel::vSetX(const Eref& e, double X) +// { +// X_ = X; +// xInited_ = true; +// } +// double HHChannel::vGetX(const Eref& e) const +// { +// return X_; +// } + +// void HHChannel::vSetY(const Eref& e, double Y) +// { +// Y_ = Y; +// yInited_ = true; +// } +// double HHChannel::vGetY(const Eref& e) const +// { +// return Y_; +// } + +// void HHChannel::vSetZ(const Eref& e, double Z) +// { +// Z_ = Z; +// zInited_ = true; +// } +// double HHChannel::vGetZ(const Eref& e) const +// { +// return Z_; +// } + +// void HHChannel::vSetUseConcentration(const Eref& e, int value) +// { +// useConcentration_ = value; +// } + +/////////////////////////////////////////////////// +// Dest function definitions +/////////////////////////////////////////////////// -void HHChannel::vProcess(const Eref& e, ProcPtr info) -{ +void HHChannel::vProcess(const Eref& e, ProcPtr info) { g_ += ChanCommon::vGetGbar(e); double A = 0; double B = 0; @@ -327,8 +344,7 @@ void HHChannel::vProcess(const Eref& e, ProcPtr info) * Here we get the steady-state values for the gate (the 'instant' * calculation) as A_/B_. */ -void HHChannel::vReinit(const Eref& er, ProcPtr info) -{ +void HHChannel::vReinit(const Eref& er, ProcPtr info) { g_ = ChanCommon::vGetGbar(er); Element* e = er.element(); @@ -385,27 +401,14 @@ void HHChannel::vReinit(const Eref& er, ProcPtr info) g_ = 0.0; } -void HHChannel::vHandleConc(const Eref& e, double conc) -{ - conc_ = conc; -} - +void HHChannel::vHandleConc(const Eref& e, double conc) { conc_ = conc; } /////////////////////////////////////////////////// // HHGate functions /////////////////////////////////////////////////// -HHGate* HHChannel::vGetXgate(unsigned int i) const -{ - return xGate_; -} +// HHGate* HHChannel::vGetXgate(unsigned int i) const { return xGate_; } -HHGate* HHChannel::vGetYgate(unsigned int i) const -{ - return yGate_; -} +// HHGate* HHChannel::vGetYgate(unsigned int i) const { return yGate_; } -HHGate* HHChannel::vGetZgate(unsigned int i) const -{ - return zGate_; -} +// HHGate* HHChannel::vGetZgate(unsigned int i) const { return zGate_; } diff --git a/biophysics/HHChannel.h b/biophysics/HHChannel.h index 3ceb14985f..d33749e8f2 100644 --- a/biophysics/HHChannel.h +++ b/biophysics/HHChannel.h @@ -53,23 +53,22 @@ class HHGate; */ class HHChannel : public HHChannelBase { - #ifdef DO_UNIT_TESTS friend void testHHChannel(); friend void testHHGateCreation(); #endif // DO_UNIT_TESTS -public: + public: HHChannel(); ~HHChannel(); ////////////////////////////////////////////////////////////////////////////////////////////// - // Avoid warning C4250 from MSVC: + // Avoid warning C4250 from MSVC: // "'HHChannel': inherits 'ChanCommon::ChanCommon::vSetGbar' via dominance" - // Although vSetX where X is the field to be set is pure virtual in ChanBase, and is - // implemented in ChanCommon, MSVC seems to still resolve it by dominance. - // Probably the ambiguity between ChanBase::vSetEk and ChanCommon::vSetEk etc. cause Ek to be - // garbage when built with MSVC. - // Looks like this bug persists since 2005: + // Although vSetX where X is the field to be set is pure virtual in + // ChanBase, and is implemented in ChanCommon, MSVC seems to still resolve + // it by dominance. Probably the ambiguity between ChanBase::vSetEk and + // ChanCommon::vSetEk etc. cause Ek to be garbage when built with MSVC. + // Looks like this bug persists since 2005: // https://stackoverflow.com/questions/469508/visual-studio-compiler-warning-c4250-class1-inherits-class2member-via-d ////////////////////////////////////////////////////////////////////////////////////////////// // void vSetGbar(const Eref&e, double Gbar ); @@ -85,20 +84,21 @@ class HHChannel : public HHChannelBase { ///////////////////////////////////////////////////////////// // Value field access function definitions ///////////////////////////////////////////////////////////// - void vSetXpower(const Eref& e, double Xpower) override; - void vSetYpower(const Eref& e, double Ypower) override; - void vSetZpower(const Eref& e, double Zpower) override; - void vSetInstant(const Eref& e, int Instant) override; - int vGetInstant(const Eref& e) const override; - void vSetX(const Eref& e, double X) override; - double vGetX(const Eref& e) const override; - void vSetY(const Eref& e, double Y) override; - double vGetY(const Eref& e) const override; - void vSetZ(const Eref& e, double Z) override; - double vGetZ(const Eref& e) const override; - void vSetUseConcentration(const Eref& e, int value) override; - // void vSetModulation(const Eref& e, double modulation) override; // defined in ChanCommon - // double vGetModulation(const Eref& e) const override; // defined in ChanCommon + // void vSetXpower(const Eref& e, double Xpower) override; + // void vSetYpower(const Eref& e, double Ypower) override; + // void vSetZpower(const Eref& e, double Zpower) override; + // void vSetInstant(const Eref& e, int Instant) override; + // int vGetInstant(const Eref& e) const override; + // void vSetX(const Eref& e, double X) override; + // double vGetX(const Eref& e) const override; + // void vSetY(const Eref& e, double Y) override; + // double vGetY(const Eref& e) const override; + // void vSetZ(const Eref& e, double Z) override; + // double vGetZ(const Eref& e) const override; + // void vSetUseConcentration(const Eref& e, int value) override; + // void vSetModulation(const Eref& e, double modulation) override; // + // defined in ChanCommon double vGetModulation(const Eref& e) const + // override; // defined in ChanCommon void innerSetXpower(double Xpower); void innerSetYpower(double Ypower); @@ -120,7 +120,7 @@ class HHChannel : public HHChannelBase { * send back to the parent compartment through regular * messages. */ - void vProcess(const Eref& e, ProcPtr p); + void vProcess(const Eref& e, ProcPtr p) override; /** * Reinitializes the values for the channel. This involves @@ -129,7 +129,7 @@ class HHChannel : public HHChannelBase { * involves a similar cycle through the gates and then * updates to the parent compartment as for the processFunc. */ - void vReinit(const Eref& e, ProcPtr p); + void vReinit(const Eref& e, ProcPtr p) override; /** * Assign the local Vm_ to the incoming Vm from the compartment @@ -142,80 +142,48 @@ class HHChannel : public HHChannelBase { * the message source will be a CaConc object, but there * are other options for computing the conc. */ - void vHandleConc(const Eref& e, double conc); + void vHandleConc(const Eref& e, double conc) override; ///////////////////////////////////////////////////////////// // Gate handling functions ///////////////////////////////////////////////////////////// - /** - * Access function used for the X gate. The index is ignored. - */ - HHGate* vGetXgate(unsigned int i) const; - - /** - * Access function used for the Y gate. The index is ignored. - */ - HHGate* vGetYgate(unsigned int i) const; + HHGate* getXgate(unsigned int i); + HHGate* getYgate(unsigned int i); + HHGate* getZgate(unsigned int i); - /** - * Access function used for the Z gate. The index is ignored. - */ - HHGate* vGetZgate(unsigned int i) const; + void setNumGates(unsigned int num); + unsigned int getNumXgates() const; + unsigned int getNumYgates() const; + unsigned int getNumZgates() const; /// Inner utility function for creating the gate. void innerCreateGate(const string& gateName, HHGate** gatePtr, Id chanId, Id gateId); /// Returns true if channel is original, false if copy. - bool checkOriginal(Id chanId) const; + bool checkOriginal(Id chanId) const override; - void vCreateGate(const Eref& e, string gateType); - /** - * Utility function for destroying gate. Works only on original - * HHChannel. Somewhat dangerous, should never be used after a - * copy has been made as the pointer of the gate will be in use - * elsewhere. - */ - void destroyGate(const Eref& e, string gateType); + void vCreateGate(const Eref& e, string gateType) override; + void destroyGate(const Eref& e, string gateType) override; /** * Inner utility for destroying the gate */ void innerDestroyGate(const string& gateName, HHGate** gatePtr, Id chanId); - /** - * Utility for altering gate powers - */ - bool setGatePower(const Eref& e, double power, double* assignee, - const string& gateType); + // /** + // * Utility for altering gate powers + // */ + // bool setGatePower(const Eref& e, double power, double* assignee, + // const string& gateType); ///////////////////////////////////////////////////////////// static const Cinfo* initCinfo(); -private: + private: /// Conc_ is input variable for Ca-dependent channels. double conc_; - double (*takeXpower_)(double, double); - double (*takeYpower_)(double, double); - double (*takeZpower_)(double, double); - - /// bitmapped flag for X, Y, Z, to do equil calculation for gate - int instant_; - /// Channel actual conductance depending on opening of gates. - /// State variable for X gate - double X_; - /// State variable for Y gate - double Y_; - /// State variable for Z gate - double Z_; - - bool xInited_, yInited_, zInited_; // true when a state variable - // has been initialized - double g_; /// Internal variable used to calculate conductance - - double integrate(double state, double dt, double A, double B); - /** * HHGate data structure for the xGate. This is writable only * on the HHChannel that originally created the HHGate, for others @@ -229,12 +197,12 @@ class HHChannel : public HHChannelBase { /// HHGate data structure for the yGate. HHGate* zGate_; - Id myId_; + // Id myId_; - static const double EPSILON; - static const int INSTANT_X; - static const int INSTANT_Y; - static const int INSTANT_Z; + // static const double EPSILON; + // static const int INSTANT_X; + // static const int INSTANT_Y; + // static const int INSTANT_Z; }; #endif // _HHChannel_h diff --git a/biophysics/HHChannel2D.cpp b/biophysics/HHChannel2D.cpp index d0d2cf02d1..b3ce2d69f9 100644 --- a/biophysics/HHChannel2D.cpp +++ b/biophysics/HHChannel2D.cpp @@ -7,577 +7,429 @@ ** See the file COPYING.LIB for the full notice. **********************************************************************/ -#include "../basecode/header.h" +#include "HHChannel2D.h" + #include "../basecode/ElementValueFinfo.h" +#include "../basecode/header.h" #include "../builtins/Interpol2D.h" #include "ChanBase.h" #include "ChanCommon.h" +#include "HHChannelBase.h" #include "HHGate2D.h" -#include "HHChannel2D.h" - -const double HHChannel2D::EPSILON = 1.0e-10; -const int HHChannel2D::INSTANT_X = 1; -const int HHChannel2D::INSTANT_Y = 2; -const int HHChannel2D::INSTANT_Z = 4; -const Cinfo* HHChannel2D::initCinfo() -{ - ///////////////////////////////////////////////////////////////////// - // Shared messages - ///////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////// -// Field definitions -/////////////////////////////////////////////////////// - static ValueFinfo< HHChannel2D, string > Xindex( "Xindex", - "String for setting X index.", - &HHChannel2D::setXindex, - &HHChannel2D::getXindex - ); - static ValueFinfo< HHChannel2D, string > Yindex( "Yindex", - "String for setting Y index.", - &HHChannel2D::setYindex, - &HHChannel2D::getYindex - ); - static ValueFinfo< HHChannel2D, string > Zindex( "Zindex", - "String for setting Z index.", - &HHChannel2D::setZindex, - &HHChannel2D::getZindex - ); - static ElementValueFinfo< HHChannel2D, double > Xpower( "Xpower", - "Power for X gate", - &HHChannel2D::setXpower, - &HHChannel2D::getXpower - ); - static ElementValueFinfo< HHChannel2D, double > Ypower( "Ypower", - "Power for Y gate", - &HHChannel2D::setYpower, - &HHChannel2D::getYpower - ); - static ElementValueFinfo< HHChannel2D, double > Zpower( "Zpower", - "Power for Z gate", - &HHChannel2D::setZpower, - &HHChannel2D::getZpower - ); - static ValueFinfo< HHChannel2D, int > instant( "instant", - "Bitmapped flag: bit 0 = Xgate, bit 1 = Ygate, bit 2 = Zgate" - "When true, specifies that the lookup table value should be" - "used directly as the state of the channel, rather than used" - "as a rate term for numerical integration for the state", - &HHChannel2D::setInstant, - &HHChannel2D::getInstant - ); - static ValueFinfo< HHChannel2D, double > X( "X", - "State variable for X gate", - &HHChannel2D::setX, - &HHChannel2D::getX - ); - static ValueFinfo< HHChannel2D, double > Y( "Y", - "State variable for Y gate", - &HHChannel2D::setY, - &HHChannel2D::getY - ); - static ValueFinfo< HHChannel2D, double > Z( "Z", - "State variable for Y gate", - &HHChannel2D::setZ, - &HHChannel2D::getZ - ); - -/////////////////////////////////////////////////////// -// MsgSrc definitions -/////////////////////////////////////////////////////// - -/////////////////////////////////////////////////////// -// MsgDest definitions -/////////////////////////////////////////////////////// - static DestFinfo concen( "concen", - "Incoming message from Concen object to specific conc to use" - "as the first concen variable", - new OpFunc1< HHChannel2D, double >( &HHChannel2D::conc1 ) - ); - static DestFinfo concen2( "concen2", - "Incoming message from Concen object to specific conc to use" - "as the second concen variable", - new OpFunc1< HHChannel2D, double >( &HHChannel2D::conc2 ) - ); -/////////////////////////////////////////////////////// -// FieldElementFinfo definition for HHGates. Note that these are made -// with the deferCreate flag off, so that the HHGates are created -// right away even if they are empty. -// I assume that we only have a single HHGate entry for each one. -/////////////////////////////////////////////////////// - static FieldElementFinfo< HHChannel2D, HHGate2D > gateX( "gateX", - "Sets up HHGate X for channel", - HHGate2D::initCinfo(), - &HHChannel2D::getXgate, - &HHChannel2D::setNumGates, - &HHChannel2D::getNumXgates - ); - static FieldElementFinfo< HHChannel2D, HHGate2D > gateY( "gateY", - "Sets up HHGate Y for channel", - HHGate2D::initCinfo(), - &HHChannel2D::getYgate, - &HHChannel2D::setNumGates, - &HHChannel2D::getNumYgates - ); - static FieldElementFinfo< HHChannel2D, HHGate2D > gateZ( "gateZ", - "Sets up HHGate Z for channel", - HHGate2D::initCinfo(), - &HHChannel2D::getZgate, - &HHChannel2D::setNumGates, - &HHChannel2D::getNumZgates - ); - static Finfo* HHChannel2DFinfos[] = - { - &Xindex, // Value - &Yindex, // Value - &Zindex, // Value - &Xpower, // Value - &Ypower, // Value - &Zpower, // Value - &instant, // Value - &X, // Value - &Y, // Value - &Z, // Value - &concen, // Dest - &concen2, // Dest - &gateX, // FieldElement - &gateY, // FieldElement - &gateZ // FieldElement - }; - - static string doc[] = - { - "Name", "HHChannel2D", - "Author", "Niraj Dudani, 2009, NCBS, Updated Upi Bhalla, 2011", - "Description", "HHChannel2D: Hodgkin-Huxley type voltage-gated Ion channel. Something " - "like the old tabchannel from GENESIS, but also presents " - "a similar interface as hhchan from GENESIS. ", - }; - - static Dinfo< HHChannel2D > dinfo; - static Cinfo HHChannel2DCinfo( - "HHChannel2D", - ChanBase::initCinfo(), - HHChannel2DFinfos, - sizeof( HHChannel2DFinfos ) / sizeof(Finfo *), - &dinfo, - doc, - sizeof(doc) / sizeof(string) - ); - - return &HHChannel2DCinfo; -} - -static const Cinfo* HHChannel2DCinfo = HHChannel2D::initCinfo(); +// const double HHChannel2D::EPSILON = 1.0e-10; +// const int HHChannel2D::INSTANT_X = 1; +// const int HHChannel2D::INSTANT_Y = 2; +// const int HHChannel2D::INSTANT_Z = 4; + +const Cinfo *HHChannel2D::initCinfo() +{ + ///////////////////////////////////////////////////////////////////// + // Shared messages + ///////////////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////// + // Field definitions + /////////////////////////////////////////////////////// + static ValueFinfo Xindex( + "Xindex", "String for setting X index.", &HHChannel2D::setXindex, + &HHChannel2D::getXindex); + static ValueFinfo Yindex( + "Yindex", "String for setting Y index.", &HHChannel2D::setYindex, + &HHChannel2D::getYindex); + static ValueFinfo Zindex( + "Zindex", "String for setting Z index.", &HHChannel2D::setZindex, + &HHChannel2D::getZindex); + static ElementValueFinfo Xpower( + "Xpower", "Power for X gate", &HHChannel2D::setXpower, + &HHChannel2D::getXpower); + static ElementValueFinfo Ypower( + "Ypower", "Power for Y gate", &HHChannel2D::setYpower, + &HHChannel2D::getYpower); + static ElementValueFinfo Zpower( + "Zpower", "Power for Z gate", &HHChannel2D::setZpower, + &HHChannel2D::getZpower); + // static ValueFinfo< HHChannel2D, int > instant( "instant", + // "Bitmapped flag: bit 0 = Xgate, bit 1 = Ygate, bit 2 = Zgate" + // "When true, specifies that the lookup table value should be" + // "used directly as the state of the channel, rather than used" + // "as a rate term for numerical integration for the state", + // &HHChannel2D::setInstant, + // &HHChannel2D::getInstant + // ); + // static ValueFinfo< HHChannel2D, double > X( "X", + // "State variable for X gate", + // &HHChannel2D::setX, + // &HHChannel2D::getX + // ); + // static ValueFinfo< HHChannel2D, double > Y( "Y", + // "State variable for Y gate", + // &HHChannel2D::setY, + // &HHChannel2D::getY + // ); + // static ValueFinfo< HHChannel2D, double > Z( "Z", + // "State variable for Y gate", + // &HHChannel2D::setZ, + // &HHChannel2D::getZ + // ); + + /////////////////////////////////////////////////////// + // MsgSrc definitions + /////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////// + // MsgDest definitions + /////////////////////////////////////////////////////// + static DestFinfo concen( + "concen", + "Incoming message from Concen object to specific conc to use" + "as the first concen variable", + new OpFunc1(&HHChannel2D::conc1)); + static DestFinfo concen2( + "concen2", + "Incoming message from Concen object to specific conc to use" + "as the second concen variable", + new OpFunc1(&HHChannel2D::conc2)); + /////////////////////////////////////////////////////// + // FieldElementFinfo definition for HHGates. Note that these are made + // with the deferCreate flag off, so that the HHGates are created + // right away even if they are empty. + // I assume that we only have a single HHGate entry for each one. + /////////////////////////////////////////////////////// + static FieldElementFinfo gateX( + "gateX", "Sets up HHGate X for channel", HHGate2D::initCinfo(), + &HHChannel2D::getXgate, &HHChannel2D::setNumGates, + &HHChannel2D::getNumXgates); + static FieldElementFinfo gateY( + "gateY", "Sets up HHGate Y for channel", HHGate2D::initCinfo(), + &HHChannel2D::getYgate, &HHChannel2D::setNumGates, + &HHChannel2D::getNumYgates); + static FieldElementFinfo gateZ( + "gateZ", "Sets up HHGate Z for channel", HHGate2D::initCinfo(), + &HHChannel2D::getZgate, &HHChannel2D::setNumGates, + &HHChannel2D::getNumZgates); + static Finfo *HHChannel2DFinfos[] = { + &Xindex, // Value + &Yindex, // Value + &Zindex, // Value + &Xpower, // Value + &Ypower, // Value + &Zpower, // Value + // &instant, // Value + // &X, // Value + // &Y, // Value + // &Z, // Value + &concen, // Dest + &concen2, // Dest + &gateX, // FieldElement + &gateY, // FieldElement + &gateZ // FieldElement + }; + + static string doc[] = { + "Name", + "HHChannel2D", + "Author", + "Niraj Dudani, 2009, NCBS, Updated Upi Bhalla, 2011", + "Description", + "HHChannel2D: Hodgkin-Huxley type voltage-gated Ion channel. Something " + "like the old tabchannel from GENESIS, but also presents " + "a similar interface as hhchan from GENESIS. ", + }; + + static Dinfo dinfo; + static Cinfo HHChannel2DCinfo("HHChannel2D", ChanBase::initCinfo(), + HHChannel2DFinfos, + sizeof(HHChannel2DFinfos) / sizeof(Finfo *), + &dinfo, doc, sizeof(doc) / sizeof(string)); + + return &HHChannel2DCinfo; +} + +static const Cinfo *HHChannel2DCinfo = HHChannel2D::initCinfo(); HHChannel2D::HHChannel2D() - : ChanCommon(), - Xpower_( 0.0 ), - Ypower_( 0.0 ), - Zpower_( 0.0 ), - instant_(0), - X_(0.0), - Y_(0.0), - Z_(0.0), - xInited_(false), - yInited_(false), - zInited_(false), - conc1_(0.0), - conc2_(0.0), - Xdep0_( -1 ), - Xdep1_( -1 ), - Ydep0_( -1 ), - Ydep1_( -1 ), - Zdep0_( -1 ), - Zdep1_( -1 ), - xGate_( 0 ), - yGate_( 0 ), - zGate_( 0 ) -{;} - + : HHChannelBase(), + conc1_(0.0), + conc2_(0.0), + Xdep0_(-1), + Xdep1_(-1), + Ydep0_(-1), + Ydep1_(-1), + Zdep0_(-1), + Zdep1_(-1), + xGate_(0), + yGate_(0), + zGate_(0) +{ + ; +} /////////////////////////////////////////////////// // Field function definitions /////////////////////////////////////////////////// -/** - * Assigns the Xpower for this gate. If the gate exists and has - * only this element for input, then change the gate value. - * If the gate exists and has multiple parents, then make a new gate. - * If the gate does not exist, make a new gate - */ - -double HHChannel2D::getXpower( const Eref& e ) const -{ - return Xpower_; -} - -double HHChannel2D::getYpower( const Eref& e ) const -{ - return Ypower_; -} - -double HHChannel2D::getZpower( const Eref& e ) const -{ - return Zpower_; -} - -void HHChannel2D::setInstant( int instant ) -{ - instant_ = instant; -} -int HHChannel2D::getInstant() const -{ - return instant_; -} - -void HHChannel2D::setX( double X ) -{ - X_ = X; - xInited_ = true; -} -double HHChannel2D::getX() const -{ - return X_; -} - -void HHChannel2D::setY( double Y ) -{ - Y_ = Y; - yInited_ = true; -} -double HHChannel2D::getY() const -{ - return Y_; -} - -void HHChannel2D::setZ( double Z ) -{ - Z_ = Z; - zInited_ = true; -} -double HHChannel2D::getZ() const -{ - return Z_; -} string HHChannel2D::getXindex() const { - return Xindex_; + return Xindex_; } -void HHChannel2D::setXindex( string Xindex ) +void HHChannel2D::setXindex(string Xindex) { - if ( Xindex == Xindex_ ) - return; + if(Xindex == Xindex_) + return; - Xindex_ = Xindex; - Xdep0_ = dependency( Xindex, 0 ); - Xdep1_ = dependency( Xindex, 1 ); + Xindex_ = Xindex; + Xdep0_ = dependency(Xindex, 0); + Xdep1_ = dependency(Xindex, 1); - assert( Xdep0_ >= 0 ); + assert(Xdep0_ >= 0); } string HHChannel2D::getYindex() const { - return Yindex_; + return Yindex_; } -void HHChannel2D::setYindex( string Yindex ) +void HHChannel2D::setYindex(string Yindex) { - if ( Yindex == Yindex_ ) - return; + if(Yindex == Yindex_) + return; - Yindex_ = Yindex; - Ydep0_ = dependency( Yindex, 0 ); - Ydep1_ = dependency( Yindex, 1 ); + Yindex_ = Yindex; + Ydep0_ = dependency(Yindex, 0); + Ydep1_ = dependency(Yindex, 1); - assert( Ydep0_ >= 0 ); + assert(Ydep0_ >= 0); } string HHChannel2D::getZindex() const { - return Zindex_; + return Zindex_; } -void HHChannel2D::setZindex( string Zindex ) +void HHChannel2D::setZindex(string Zindex) { - if ( Zindex == Zindex_ ) - return; + if(Zindex == Zindex_) + return; - Zindex_ = Zindex; - Zdep0_ = dependency( Zindex, 0 ); - Zdep1_ = dependency( Zindex, 1 ); + Zindex_ = Zindex; + Zdep0_ = dependency(Zindex, 0); + Zdep1_ = dependency(Zindex, 1); - assert( Zdep0_ >= 0 ); + assert(Zdep0_ >= 0); } //////////////////////////////////////////////////////////////////// // HHGate2D access funcs //////////////////////////////////////////////////////////////////// -HHGate2D* HHChannel2D::getXgate( unsigned int i ) +HHGate2D *HHChannel2D::getXgate(unsigned int i) { - return xGate_; + return xGate_; } -HHGate2D* HHChannel2D::getYgate( unsigned int i ) +HHGate2D *HHChannel2D::getYgate(unsigned int i) { - return yGate_; + return yGate_; } -HHGate2D* HHChannel2D::getZgate( unsigned int i ) +HHGate2D *HHChannel2D::getZgate(unsigned int i) { - return zGate_; + return zGate_; } -void HHChannel2D::setNumGates( unsigned int num ) -{ ; } - -unsigned int HHChannel2D::getNumXgates() const +void HHChannel2D::setNumGates(unsigned int num) { - return ( xGate_ != 0 ); + ; } - -unsigned int HHChannel2D::getNumYgates() const +unsigned int HHChannel2D::getNumXgates() const { - return ( yGate_ != 0 ); + return xGate_ != nullptr; } - -unsigned int HHChannel2D::getNumZgates() const +unsigned int HHChannel2D::getNumYgates() const +{ + return yGate_ != nullptr; +} +unsigned int HHChannel2D::getNumZgates() const { - return ( zGate_ != 0 ); + return zGate_ != nullptr; } -double HHChannel2D::depValue( int dep ) +double HHChannel2D::depValue(int dep) { - switch( dep ) - { - case 0: return Vm_; - case 1: return conc1_; - case 2: return conc2_; - default: assert( 0 ); return 0.0; - } + switch(dep) { + case 0: + return Vm_; + case 1: + return conc1_; + case 2: + return conc2_; + default: + assert(0); + return 0.0; + } } -int HHChannel2D::dependency( string index, unsigned int dim ) +int HHChannel2D::dependency(string index, unsigned int dim) { - static vector< map< string, int > > dep; - if ( dep.empty() ) { - dep.resize( 2 ); - - dep[ 0 ][ "VOLT_INDEX" ] = 0; - dep[ 0 ][ "C1_INDEX" ] = 1; - dep[ 0 ][ "C2_INDEX" ] = 2; - - dep[ 0 ][ "VOLT_C1_INDEX" ] = 0; - dep[ 0 ][ "VOLT_C2_INDEX" ] = 0; - dep[ 0 ][ "C1_C2_INDEX" ] = 1; - - dep[ 1 ][ "VOLT_INDEX" ] = -1; - dep[ 1 ][ "C1_INDEX" ] = -1; - dep[ 1 ][ "C2_INDEX" ] = -1; - - dep[ 1 ][ "VOLT_C1_INDEX" ] = 1; - dep[ 1 ][ "VOLT_C2_INDEX" ] = 2; - dep[ 1 ][ "C1_C2_INDEX" ] = 2; - } - - if ( dep[ dim ].find( index ) == dep[ dim ].end() ) - return -1; - - if ( dep[ dim ][ index ] == 0 ) - return 0; - if ( dep[ dim ][ index ] == 1 ) - return 1; - if ( dep[ dim ][ index ] == 2 ) - return 2; - - return -1; + static vector> dep; + if(dep.empty()) { + dep.resize(2); + + dep[0]["VOLT_INDEX"] = 0; + dep[0]["C1_INDEX"] = 1; + dep[0]["C2_INDEX"] = 2; + + dep[0]["VOLT_C1_INDEX"] = 0; + dep[0]["VOLT_C2_INDEX"] = 0; + dep[0]["C1_C2_INDEX"] = 1; + + dep[1]["VOLT_INDEX"] = -1; + dep[1]["C1_INDEX"] = -1; + dep[1]["C2_INDEX"] = -1; + + dep[1]["VOLT_C1_INDEX"] = 1; + dep[1]["VOLT_C2_INDEX"] = 2; + dep[1]["C1_C2_INDEX"] = 2; + } + + if(dep[dim].find(index) == dep[dim].end()) + return -1; + + if(dep[dim][index] == 0) + return 0; + if(dep[dim][index] == 1) + return 1; + if(dep[dim][index] == 2) + return 2; + + return -1; } /////////////////////////////////////////////////// // Dest function definitions /////////////////////////////////////////////////// -void HHChannel2D::conc1( double conc ) +void HHChannel2D::conc1(double conc) { - conc1_ = conc; + conc1_ = conc; } -void HHChannel2D::conc2( double conc ) +void HHChannel2D::conc2(double conc) { - conc2_ = conc; + conc2_ = conc; } /////////////////////////////////////////////////// // utility function definitions /////////////////////////////////////////////////// -/** - * Returns the state variable for the new timestep based on - * the internal variables A_ and B_ which were passed in from the gate. - */ -double HHChannel2D::integrate( double state, double dt, double A, double B ) -{ - if ( B > EPSILON ) { - double x = exp( -B * dt ); - return state * x + ( A / B ) * ( 1 - x ); - } - return state + A * dt ; +void HHChannel2D::vProcess(const Eref &e, ProcPtr info) +{ + g_ += ChanBase::getGbar(e); + double A = 0; + double B = 0; + if(Xpower_ > 0) { + xGate_->lookupBoth(depValue(Xdep0_), depValue(Xdep1_), &A, &B); + if(instant_ & INSTANT_X) + X_ = A / B; + else + X_ = integrate(X_, info->dt, A, B); + g_ *= takeXpower_(X_, Xpower_); + } + + if(Ypower_ > 0) { + yGate_->lookupBoth(depValue(Ydep0_), depValue(Ydep1_), &A, &B); + if(instant_ & INSTANT_Y) + Y_ = A / B; + else + Y_ = integrate(Y_, info->dt, A, B); + + g_ *= takeYpower_(Y_, Ypower_); + } + + if(Zpower_ > 0) { + zGate_->lookupBoth(depValue(Zdep0_), depValue(Zdep1_), &A, &B); + if(instant_ & INSTANT_Z) + Z_ = A / B; + else + Z_ = integrate(Z_, info->dt, A, B); + + g_ *= takeZpower_(Z_, Zpower_); + } + + ChanBase::setGk(e, g_ * vGetModulation(e)); + updateIk(); + // Gk_ = g_; + // Ik_ = ( Ek_ - Vm_ ) * g_; + + // Send out the relevant channel messages. + sendProcessMsgs(e, info); + g_ = 0.0; } -void HHChannel2D::vProcess( const Eref& e, ProcPtr info ) -{ - g_ += ChanBase::getGbar( e ); - double A = 0; - double B = 0; - if ( Xpower_ > 0 ) { - xGate_->lookupBoth( depValue( Xdep0_ ), depValue( Xdep1_ ), &A, &B ); - if ( instant_ & INSTANT_X ) - X_ = A/B; - else - X_ = integrate( X_, info->dt, A, B ); - g_ *= takeXpower_( X_, Xpower_ ); - } - - if ( Ypower_ > 0 ) { - yGate_->lookupBoth( depValue( Ydep0_ ), depValue( Ydep1_ ), &A, &B ); - if ( instant_ & INSTANT_Y ) - Y_ = A/B; - else - Y_ = integrate( Y_, info->dt, A, B ); - - g_ *= takeYpower_( Y_, Ypower_ ); - } - - if ( Zpower_ > 0 ) { - zGate_->lookupBoth( depValue( Zdep0_ ), depValue( Zdep1_ ), &A, &B ); - if ( instant_ & INSTANT_Z ) - Z_ = A/B; - else - Z_ = integrate( Z_, info->dt, A, B ); - - g_ *= takeZpower_( Z_, Zpower_ ); - } - - ChanBase::setGk( e, g_ * vGetModulation( e ) ); - updateIk(); - // Gk_ = g_; - // Ik_ = ( Ek_ - Vm_ ) * g_; - - // Send out the relevant channel messages. - sendProcessMsgs( e, info ); +/** + * Here we get the steady-state values for the gate (the 'instant' + * calculation) as A_/B_. + */ +void HHChannel2D::vReinit(const Eref &er, ProcPtr info) +{ + g_ = ChanBase::getGbar(er); + Element *e = er.element(); + + double A = 0.0; + double B = 0.0; + if(Xpower_ > 0) { + xGate_->lookupBoth(depValue(Xdep0_), depValue(Xdep1_), &A, &B); + if(B < EPSILON) { + cout << "Warning: B_ value for " << e->getName() + << " is ~0. Check X table\n"; + return; + } + if(!xInited_) + X_ = A / B; + g_ *= takeXpower_(X_, Xpower_); + } + + if(Ypower_ > 0) { + yGate_->lookupBoth(depValue(Ydep0_), depValue(Ydep1_), &A, &B); + if(B < EPSILON) { + cout << "Warning: B value for " << e->getName() + << " is ~0. Check Y table\n"; + return; + } + if(!yInited_) + Y_ = A / B; + g_ *= takeYpower_(Y_, Ypower_); + } + + if(Zpower_ > 0) { + zGate_->lookupBoth(depValue(Zdep0_), depValue(Zdep1_), &A, &B); + if(B < EPSILON) { + cout << "Warning: B value for " << e->getName() + << " is ~0. Check Z table\n"; + return; + } + if(!zInited_) + Z_ = A / B; + g_ *= takeZpower_(Z_, Zpower_); + } + + ChanBase::setGk(er, g_ * vGetModulation(er)); + updateIk(); + // Gk_ = g_; + // Ik_ = ( Ek_ - Vm_ ) * g_; + + // Send out the relevant channel messages. + // Same for reinit as for process. + sendReinitMsgs(er, info); g_ = 0.0; - } - - /** - * Here we get the steady-state values for the gate (the 'instant' - * calculation) as A_/B_. - */ - void HHChannel2D::vReinit( const Eref& er, ProcPtr info ) - { - g_ = ChanBase::getGbar( er ); - Element* e = er.element(); - - double A = 0.0; - double B = 0.0; - if ( Xpower_ > 0 ) { - xGate_->lookupBoth( depValue( Xdep0_ ), depValue( Xdep1_ ), &A, &B ); - if ( B < EPSILON ) { - cout << "Warning: B_ value for " << e->getName() << - " is ~0. Check X table\n"; - return; - } - if (!xInited_) - X_ = A/B; - g_ *= takeXpower_( X_, Xpower_ ); - } - - if ( Ypower_ > 0 ) { - yGate_->lookupBoth( depValue( Ydep0_ ), depValue( Ydep1_ ), &A, &B ); - if ( B < EPSILON ) { - cout << "Warning: B value for " << e->getName() << - " is ~0. Check Y table\n"; - return; - } - if (!yInited_) - Y_ = A/B; - g_ *= takeYpower_( Y_, Ypower_ ); - } - - if ( Zpower_ > 0 ) { - zGate_->lookupBoth( depValue( Zdep0_ ), depValue( Zdep1_ ), &A, &B ); - if ( B < EPSILON ) { - cout << "Warning: B value for " << e->getName() << - " is ~0. Check Z table\n"; - return; - } - if (!zInited_) - Z_ = A/B; - g_ *= takeZpower_( Z_, Zpower_ ); - } - - ChanBase::setGk( er, g_ * vGetModulation( er ) ); - updateIk(); - // Gk_ = g_; - // Ik_ = ( Ek_ - Vm_ ) * g_; - - // Send out the relevant channel messages. - // Same for reinit as for process. - sendReinitMsgs( er, info ); - g_ = 0.0; } //////////////////////////////////////////////////////////////////////// // Gate management stuff. //////////////////////////////////////////////////////////////////////// -bool HHChannel2D::setGatePower( const Eref& e, double power, - double *assignee, const string& gateType ) -{ - if ( power < 0 ) { - cout << "Error: HHChannel2D::set" << gateType << - "power: Cannot use negative power: " << power << endl; - return 0; - } - - if ( doubleEq( power, *assignee ) ) - return 0; - - if ( doubleEq( *assignee, 0.0 ) && power > 0 ) { - createGate( e, gateType ); - } else if ( doubleEq( power, 0.0 ) ) { - destroyGate( e, gateType ); - } - - *assignee = power; - return 1; -} - -/** - * Assigns the Xpower for this gate. If the gate exists and has - * only this element for input, then change the gate value. - * If the gate exists and has multiple parents, then make a new gate. - * If the gate does not exist, make a new gate - */ -void HHChannel2D::setXpower( const Eref& e, double Xpower ) -{ - if ( setGatePower( e, Xpower, &Xpower_, "X" ) ) - takeXpower_ = selectPower( Xpower ); -} - -void HHChannel2D::setYpower( const Eref& e, double Ypower ) -{ - if ( setGatePower( e, Ypower, &Ypower_, "Y" ) ) - takeYpower_ = selectPower( Ypower ); -} - -void HHChannel2D::setZpower( const Eref& e, double Zpower ) -{ - if ( setGatePower( e, Zpower, &Zpower_, "Z" ) ) - takeZpower_ = selectPower( Zpower ); -} - /** * If the gate exists and has only this element for input, then change * the gate power. @@ -596,102 +448,83 @@ void HHChannel2D::setZpower( const Eref& e, double Zpower ) // Assuming that the elements are simple elements. Use Eref for // general case -bool HHChannel2D::checkOriginal( Id chanId ) const -{ - bool isOriginal = 1; - if ( xGate_ ) { - isOriginal = xGate_->isOriginalChannel( chanId ); - } else if ( yGate_ ) { - isOriginal = yGate_->isOriginalChannel( chanId ); - } else if ( zGate_ ) { - isOriginal = zGate_->isOriginalChannel( chanId ); - } - return isOriginal; -} - -void HHChannel2D::innerCreateGate( const string& gateName, - HHGate2D** gatePtr, Id chanId, Id gateId ) -{ - //Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() ); - if ( *gatePtr ) { - cout << "Warning: HHChannel2D::createGate: '" << gateName << - "' on Element '" << chanId.path() << "' already present\n"; - return; - } - *gatePtr = new HHGate2D( chanId, gateId ); -} - -void HHChannel2D::createGate( const Eref& e, - string gateType ) -{ - if ( !checkOriginal( e.id() ) ) { - cout << "Warning: HHChannel2D::createGate: Not allowed from copied channel:\n" << e.id().path() << "\n"; - return; - } - - if ( gateType == "X" ) - innerCreateGate( "xGate", &xGate_, e.id(), Id(e.id().value() + 1) ); - else if ( gateType == "Y" ) - innerCreateGate( "yGate", &yGate_, e.id(), Id(e.id().value() + 2) ); - else if ( gateType == "Z" ) - innerCreateGate( "zGate", &zGate_, e.id(), Id(e.id().value() + 3) ); - else - cout << "Warning: HHChannel2D::createGate: Unknown gate type '" << - gateType << "'. Ignored\n"; -} - -void HHChannel2D::innerDestroyGate( const string& gateName, - HHGate2D** gatePtr, Id chanId ) -{ - if ( *gatePtr == 0 ) { - cout << "Warning: HHChannel2D::destroyGate: '" << gateName << - "' on Element '" << chanId.path() << "' not present\n"; - return; - } - delete (*gatePtr); - *gatePtr = 0; -} - -void HHChannel2D::destroyGate( const Eref& e, - string gateType ) -{ - if ( !checkOriginal( e.id() ) ) { - cout << "Warning: HHChannel2D::destroyGate: Not allowed from copied channel:\n" << e.id().path() << "\n"; - return; - } - - if ( gateType == "X" ) - innerDestroyGate( "xGate", &xGate_, e.id() ); - else if ( gateType == "Y" ) - innerDestroyGate( "yGate", &yGate_, e.id() ); - else if ( gateType == "Z" ) - innerDestroyGate( "zGate", &zGate_, e.id() ); - else - cout << "Warning: HHChannel2D::destroyGate: Unknown gate type '" << - gateType << "'. Ignored\n"; -} - -double HHChannel2D::powerN( double x, double p ) -{ - if ( x > 0.0 ) - return exp( p * log( x ) ); - return 0.0; -} - -PFDD HHChannel2D::selectPower( double power ) -{ - if ( power == 0.0 ) - return powerN; - else if ( power == 1.0 ) - return power1; - else if ( power == 2.0 ) - return power2; - else if ( power == 3.0 ) - return power3; - else if ( power == 4.0 ) - return power4; - else - return powerN; +bool HHChannel2D::checkOriginal(Id chanId) const +{ + bool isOriginal = 1; + if(xGate_) { + isOriginal = xGate_->isOriginalChannel(chanId); + } + else if(yGate_) { + isOriginal = yGate_->isOriginalChannel(chanId); + } + else if(zGate_) { + isOriginal = zGate_->isOriginalChannel(chanId); + } + return isOriginal; +} + +void HHChannel2D::innerCreateGate(const string &gateName, HHGate2D **gatePtr, + Id chanId, Id gateId) +{ + // Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() ); + if(*gatePtr) { + cout << "Warning: HHChannel2D::createGate: '" << gateName + << "' on Element '" << chanId.path() << "' already present\n"; + return; + } + *gatePtr = new HHGate2D(chanId, gateId); +} + +void HHChannel2D::vCreateGate(const Eref &e, string gateType) +{ + if(!checkOriginal(e.id())) { + cout << "Warning: HHChannel2D::createGate: Not allowed from copied " + "channel:\n" + << e.id().path() << "\n"; + return; + } + + if(gateType == "X") + innerCreateGate("xGate", &xGate_, e.id(), Id(e.id().value() + 1)); + else if(gateType == "Y") + innerCreateGate("yGate", &yGate_, e.id(), Id(e.id().value() + 2)); + else if(gateType == "Z") + innerCreateGate("zGate", &zGate_, e.id(), Id(e.id().value() + 3)); + else + cout << "Warning: HHChannel2D::createGate: Unknown gate type '" + << gateType << "'. Ignored\n"; +} + +void HHChannel2D::innerDestroyGate(const string &gateName, HHGate2D **gatePtr, + Id chanId) +{ + if(*gatePtr == nullptr) { + cout << "Warning: HHChannel2D::destroyGate: '" << gateName + << "' on Element '" << chanId.path() << "' not present\n"; + return; + } + delete(*gatePtr); + *gatePtr = nullptr; +} + +void HHChannel2D::destroyGate(const Eref &e, string gateType) +{ + if(!checkOriginal(e.id())) { + cout << "Warning: HHChannel2D::destroyGate: Not allowed from copied " + "channel:\n" + << e.id().path() << "\n"; + return; + } + + if(gateType == "X") + innerDestroyGate("xGate", &xGate_, e.id()); + else if(gateType == "Y") + innerDestroyGate("yGate", &yGate_, e.id()); + else if(gateType == "Z") + innerDestroyGate("zGate", &zGate_, e.id()); + else + cout << "Warning: HHChannel2D::destroyGate: Unknown gate type '" + << gateType << "'. Ignored\n"; } /////////////////////////////////////////////////// @@ -702,6 +535,6 @@ PFDD HHChannel2D::selectPower( double power ) void testHHChannel2D2D() { - ; + ; } #endif diff --git a/biophysics/HHChannel2D.h b/biophysics/HHChannel2D.h index b29398274c..1f8b71d509 100644 --- a/biophysics/HHChannel2D.h +++ b/biophysics/HHChannel2D.h @@ -25,208 +25,172 @@ // Ported to asyn13 on 2014-05-30 by Subhasis Ray -typedef double ( *PFDD )( double, double ); - class HHGate2D; -class HHChannel2D: public ChanCommon +class HHChannel2D : public HHChannelBase { #ifdef DO_UNIT_TESTS - friend void testHHChannel2D(); -#endif // DO_UNIT_TESTS - - public: - HHChannel2D(); - - ///////////////////////////////////////////////////////////// - // Value field access function definitions - ///////////////////////////////////////////////////////////// - void setUseConcentration( int value ); - int getUseConcentration( ); - void setXindex( string index ); - string getXindex() const; - void setYindex( string index ); - string getYindex() const; - void setZindex( string index ); - string getZindex() const; - - void setXpower( const Eref& e, double Xpower ); - double getXpower( const Eref& e) const; - void setYpower( const Eref& e, double Ypower ); - double getYpower( const Eref& e ) const; - void setZpower( const Eref& e, double Zpower ); - double getZpower( const Eref& e ) const; - void setInstant( int Instant ); - int getInstant() const; - void setX( double X ); - double getX() const; - void setY( double Y ); - double getY() const; - void setZ( double Z ); - double getZ() const; - - /// Access function used for the X gate. The index is ignored. - HHGate2D* getXgate( unsigned int i ); - /// Access function used for the Y gate. The index is ignored. - HHGate2D* getYgate( unsigned int i ); - /// Access function used for the Z gate. The index is ignored. - HHGate2D* getZgate( unsigned int i ); - /// Dummy assignment function for the number of gates. - void setNumGates( unsigned int num ); - - /** - * Access function for the number of Xgates. Gives 1 if present, - * otherwise 0. - */ - unsigned int getNumXgates() const; - /// Returns 1 if Y gate present, otherwise 0 - unsigned int getNumYgates() const; - /// Returns 1 if Z gate present, otherwise 0 - unsigned int getNumZgates() const; - - ///////////////////////////////////////////////////////////// - // Dest function definitions - ///////////////////////////////////////////////////////////// - /** - * processFunc handles the update and calculations every - * clock tick. It first sends the request for evaluation of - * the gate variables to the respective gate objects and - * recieves their response immediately through a return - * message. This is done so that many channel instances can - * share the same gate lookup tables, but do so cleanly. - * Such messages should never go to a remote node. - * Then the function does its own little calculations to - * send back to the parent compartment through regular - * messages. - */ - void vProcess( const Eref& e, ProcPtr p ); - - /** - * Reinitializes the values for the channel. This involves - * computing the steady-state value for the channel gates - * using the provided Vm from the parent compartment. It - * involves a similar cycle through the gates and then - * updates to the parent compartment as for the processFunc. - */ - void vReinit( const Eref& e, ProcPtr p ); - /** - * Assign the local conc_ to the incoming conc from the - * concentration calculations for the compartment. Typically - * the message source will be a CaConc object, but there - * are other options for computing the conc. - */ - void conc1( double conc ); - void conc2( double conc ); - - double ( *takeXpower_ )( double, double ); - double ( *takeYpower_ )( double, double ); - double ( *takeZpower_ )( double, double ); - - /** - * Function for safely creating each gate, identified by strings - * as X, Y and Z. Will only work on a new channel, not on a - * copy. The idea is that the gates are always referred to the - * original 'library' channel, and their contents cannot be touched - * except by the original. - */ - void createGate( const Eref& e, string gateType ); - - /// Inner utility function for creating the gate. - void innerCreateGate( - const string& gateName, - HHGate2D** gatePtr, Id chanId, Id gateId ); - - /// Returns true if channel is original, false if copy. - bool checkOriginal( Id chanId ) const; - - /** - * Utility function for destroying gate. Works only on original - * HHChannel. Somewhat dangerous, should never be used after a - * copy has been made as the pointer of the gate will be in use - * elsewhere. - */ - void destroyGate( const Eref& e, string gateType ); - - /** - * Inner utility for destroying the gate - */ - void innerDestroyGate( const string& gateName, - HHGate2D** gatePtr, Id chanId ); - - /** - * Utility for altering gate powers - */ - bool setGatePower( const Eref& e, double power, - double* assignee, const string& gateType ); - - static PFDD selectPower( double power); - - static const Cinfo* initCinfo(); - - private: - int dependency( string index, unsigned int dim ); - double depValue( int dependency ); - double integrate( double state, double dt, double A, double B ); - - double Xpower_; /// Exponent for X gate - double Ypower_; /// Exponent for Y gate - double Zpower_; /// Exponent for Z gate - - /// bitmapped flag for X, Y, Z, to do equil calculation for gate - int instant_; - double X_; /// State variable for X gate - double Y_; /// State variable for Y gate - double Z_; /// State variable for Z gate - - /** - * true when the matching state variable has been initialized - */ - bool xInited_, yInited_, zInited_; - - double g_; /// Internal variable used to calculate conductance - - double conc1_; - double conc2_; - - string Xindex_; - string Yindex_; - string Zindex_; - - int Xdep0_; - int Xdep1_; - int Ydep0_; - int Ydep1_; - int Zdep0_; - int Zdep1_; - - /** - * HHGate2D data structure for the xGate. This is writable only - * on the HHChannel that originally created the HHGate, for others - * it must be treated as readonly. - */ - HHGate2D* xGate_; - HHGate2D* yGate_; /// HHGate2D for the yGate - HHGate2D* zGate_; /// HHGate2D for the zGate - - static const double EPSILON; - static const int INSTANT_X; - static const int INSTANT_Y; - static const int INSTANT_Z; - - static double power1( double x, double p ) { - return x; - } - static double power2( double x, double p ) { - return x * x; - } - static double power3( double x, double p ) { - return x * x * x; - } - static double power4( double x, double p ) { - return power2( x * x, p ); - } - static double powerN( double x, double p ); + friend void testHHChannel2D(); +#endif // DO_UNIT_TESTS + +public: + HHChannel2D(); + + ///////////////////////////////////////////////////////////// + // Value field access function definitions + ///////////////////////////////////////////////////////////// + // void setUseConcentration(int value); + // int getUseConcentration(); + void setXindex(string index); + string getXindex() const; + void setYindex(string index); + string getYindex() const; + void setZindex(string index); + string getZindex() const; + + // void setXpower(const Eref& e, double Xpower); + // double getXpower(const Eref& e) const; + // void setYpower(const Eref& e, double Ypower); + // double getYpower(const Eref& e) const; + // void setZpower(const Eref& e, double Zpower); + // double getZpower(const Eref& e) const; + // void setInstant(int Instant); + // int getInstant() const; + // void setX(double X); + // double getX() const; + // void setY(double Y); + // double getY() const; + // void setZ(double Z); + // double getZ() const; + + /// Access function used for the X gate. The index is ignored. + HHGate2D* getXgate(unsigned int i); + /// Access function used for the Y gate. The index is ignored. + HHGate2D* getYgate(unsigned int i); + /// Access function used for the Z gate. The index is ignored. + HHGate2D* getZgate(unsigned int i); + /// Dummy assignment function for the number of gates. + void setNumGates(unsigned int num); + + /** + * Access function for the number of Xgates. Gives 1 if present, + * otherwise 0. + */ + unsigned int getNumXgates() const; + /// Returns 1 if Y gate present, otherwise 0 + unsigned int getNumYgates() const; + /// Returns 1 if Z gate present, otherwise 0 + unsigned int getNumZgates() const; + + ///////////////////////////////////////////////////////////// + // Dest function definitions + ///////////////////////////////////////////////////////////// + /** + * processFunc handles the update and calculations every + * clock tick. It first sends the request for evaluation of + * the gate variables to the respective gate objects and + * recieves their response immediately through a return + * message. This is done so that many channel instances can + * share the same gate lookup tables, but do so cleanly. + * Such messages should never go to a remote node. + * Then the function does its own little calculations to + * send back to the parent compartment through regular + * messages. + */ + void vProcess(const Eref& e, ProcPtr p) override; + + /** + * Reinitializes the values for the channel. This involves + * computing the steady-state value for the channel gates + * using the provided Vm from the parent compartment. It + * involves a similar cycle through the gates and then + * updates to the parent compartment as for the processFunc. + */ + void vReinit(const Eref& e, ProcPtr p) override; + /** + * Assign the local conc_ to the incoming conc from the + * concentration calculations for the compartment. Typically + * the message source will be a CaConc object, but there + * are other options for computing the conc. + */ + void conc1(double conc); + void conc2(double conc); + + + /** + * Function for safely creating each gate, identified by strings + * as X, Y and Z. Will only work on a new channel, not on a + * copy. The idea is that the gates are always referred to the + * original 'library' channel, and their contents cannot be touched + * except by the original. + */ + void vCreateGate(const Eref& e, string gateType) override; + + /// Inner utility function for creating the gate. + void innerCreateGate(const string& gateName, HHGate2D** gatePtr, Id chanId, + Id gateId); + + /// Returns true if channel is original, false if copy. + bool checkOriginal(Id chanId) const override; + + /** + * Utility function for destroying gate. Works only on original + * HHChannel. Somewhat dangerous, should never be used after a + * copy has been made as the pointer of the gate will be in use + * elsewhere. + */ + void destroyGate(const Eref& e, string gateType) override; + + /** + * Inner utility for destroying the gate + */ + void innerDestroyGate(const string& gateName, HHGate2D** gatePtr, + Id chanId); + static const Cinfo* initCinfo(); + +private: + int dependency(string index, unsigned int dim); + double depValue(int dependency); + double conc1_; + double conc2_; + + string Xindex_; + string Yindex_; + string Zindex_; + + int Xdep0_; + int Xdep1_; + int Ydep0_; + int Ydep1_; + int Zdep0_; + int Zdep1_; + + /** + * HHGate2D data structure for the xGate. This is writable only + * on the HHChannel that originally created the HHGate, for others + * it must be treated as readonly. + */ + HHGate2D* xGate_; + HHGate2D* yGate_; /// HHGate2D for the yGate + HHGate2D* zGate_; /// HHGate2D for the zGate + + static double power1(double x, double p) + { + return x; + } + static double power2(double x, double p) + { + return x * x; + } + static double power3(double x, double p) + { + return x * x * x; + } + static double power4(double x, double p) + { + return power2(x * x, p); + } + static double powerN(double x, double p); }; - -#endif // _HHChannel2D_h +#endif // _HHChannel2D_h diff --git a/biophysics/HHChannelBase.cpp b/biophysics/HHChannelBase.cpp index e2926c5be7..3ca9543dd0 100644 --- a/biophysics/HHChannelBase.cpp +++ b/biophysics/HHChannelBase.cpp @@ -7,14 +7,19 @@ ** See the file COPYING.LIB for the full notice. **********************************************************************/ -#include "../basecode/header.h" +#include "HHChannelBase.h" + #include "../basecode/ElementValueFinfo.h" -#include "HHGate.h" +#include "../basecode/header.h" #include "ChanBase.h" -#include "HHChannelBase.h" +#include "HHGate.h" -const Cinfo* HHChannelBase::initCinfo() -{ +const double HHChannelBase::EPSILON = 1.0e-10; +const int HHChannelBase::INSTANT_X = 1; +const int HHChannelBase::INSTANT_Y = 2; +const int HHChannelBase::INSTANT_Z = 4; + +const Cinfo *HHChannelBase::initCinfo() { ///////////////////////////////////////////////////////////////////// // Shared messages ///////////////////////////////////////////////////////////////////// @@ -78,27 +83,9 @@ const Cinfo* HHChannelBase::initCinfo() // right away even if they are empty. // Assume only a single entry allocated in each gate. /////////////////////////////////////////////////////// - static FieldElementFinfo gateX( - "gateX", "Sets up HHGate X for channel", HHGate::initCinfo(), - &HHChannelBase::getXgate, &HHChannelBase::setNumGates, - &HHChannelBase::getNumXgates - // 1 - ); - static FieldElementFinfo gateY( - "gateY", "Sets up HHGate Y for channel", HHGate::initCinfo(), - &HHChannelBase::getYgate, &HHChannelBase::setNumGates, - &HHChannelBase::getNumYgates - // 1 - ); - static FieldElementFinfo gateZ( - "gateZ", "Sets up HHGate Z for channel", HHGate::initCinfo(), - &HHChannelBase::getZgate, &HHChannelBase::setNumGates, - &HHChannelBase::getNumZgates - // 1 - ); /////////////////////////////////////////////////////// - static Finfo* HHChannelBaseFinfos[] = { + static Finfo *HHChannelBaseFinfos[] = { &Xpower, // Value &Ypower, // Value &Zpower, // Value @@ -109,9 +96,9 @@ const Cinfo* HHChannelBase::initCinfo() &useConcentration, // Value &concen, // Dest &createGate, // Dest - &gateX, // FieldElement - &gateY, // FieldElement - &gateZ // FieldElement + // &gateX, // FieldElement + // &gateY, // FieldElement + // &gateZ // FieldElement }; static string doc[] = { @@ -130,39 +117,43 @@ const Cinfo* HHChannelBase::initCinfo() static Cinfo HHChannelBaseCinfo( "HHChannelBase", ChanBase::initCinfo(), HHChannelBaseFinfos, - sizeof(HHChannelBaseFinfos) / sizeof(Finfo*), &dinfo, doc, + sizeof(HHChannelBaseFinfos) / sizeof(Finfo *), &dinfo, doc, sizeof(doc) / sizeof(string)); return &HHChannelBaseCinfo; } -static const Cinfo* hhChannelCinfo = HHChannelBase::initCinfo(); +static const Cinfo *hhChannelCinfo = HHChannelBase::initCinfo(); ////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////// // Constructor /////////////////////////////////////////////////// HHChannelBase::HHChannelBase() - : Xpower_(0.0), + : instant_(0), + X_(0.0), + Y_(0.0), + Z_(0.0), + xInited_(false), + yInited_(false), + zInited_(false), + g_(0.0), + Xpower_(0.0), Ypower_(0.0), Zpower_(0.0), - useConcentration_(0) -{ - ; + useConcentration_(0), + myId_() { + ; } -HHChannelBase::~HHChannelBase() -{ - ; -} +HHChannelBase::~HHChannelBase() { ; } -bool checkPower(double power) -{ - if(power < 0.0) { +bool checkPower(double power) { + if (power < 0.0) { cout << "Warning: HHChannelBase::setPower: Cannot be negative\n"; return false; } - if(power > 5.0) { + if (power > 5.0) { cout << "Warning: HHChannelBase::setPower: unlikely to be > 5\n"; return false; } @@ -174,93 +165,91 @@ bool checkPower(double power) * If the gate exists and has multiple parents, then make a new gate. * If the gate does not exist, make a new gate */ -void HHChannelBase::setXpower(const Eref& e, double power) -{ - if(checkPower(power)) - vSetXpower(e, power); +void HHChannelBase::setXpower(const Eref &e, double power) { + if (checkPower(power)) vSetXpower(e, power); } -void HHChannelBase::setYpower(const Eref& e, double power) -{ - if(checkPower(power)) - vSetYpower(e, power); +void HHChannelBase::setYpower(const Eref &e, double power) { + if (checkPower(power)) vSetYpower(e, power); } -void HHChannelBase::setZpower(const Eref& e, double power) -{ - if(checkPower(power)) - vSetZpower(e, power); +void HHChannelBase::setZpower(const Eref &e, double power) { + if (checkPower(power)) vSetZpower(e, power); } -void HHChannelBase::createGate(const Eref& e, string gateType) -{ +void HHChannelBase::createGate(const Eref &e, string gateType) { vCreateGate(e, gateType); } +void HHChannelBase::destroyGate(const Eref &e, string gateType) { + cerr << "HHChannelBase::destroyGate: This should never have been reached." + << endl; +} + + +bool HHChannelBase::checkOriginal(Id id) const +{ + cout << "Warning: HHChannelBase::checkOriginal(Id id) should be overriden by all derived classes." << endl; + return true; +} /////////////////////////////////////////////////// // Field function definitions /////////////////////////////////////////////////// -double HHChannelBase::getXpower(const Eref& e) const -{ - return Xpower_; -} +double HHChannelBase::getXpower(const Eref &e) const { return Xpower_; } -double HHChannelBase::getYpower(const Eref& e) const -{ - return Ypower_; -} +double HHChannelBase::getYpower(const Eref &e) const { return Ypower_; } -double HHChannelBase::getZpower(const Eref& e) const -{ - return Zpower_; -} +double HHChannelBase::getZpower(const Eref &e) const { return Zpower_; } -void HHChannelBase::setInstant(const Eref& e, int instant) -{ +void HHChannelBase::setInstant(const Eref &e, int instant) { vSetInstant(e, instant); } -int HHChannelBase::getInstant(const Eref& e) const -{ - return vGetInstant(e); -} +int HHChannelBase::getInstant(const Eref &e) const { return vGetInstant(e); } -void HHChannelBase::setX(const Eref& e, double X) -{ - vSetX(e, X); -} -double HHChannelBase::getX(const Eref& e) const -{ - return vGetX(e); +void HHChannelBase::setX(const Eref &e, double X) { vSetX(e, X); } +double HHChannelBase::getX(const Eref &e) const { return vGetX(e); } + +void HHChannelBase::setY(const Eref &e, double Y) { vSetY(e, Y); } +double HHChannelBase::getY(const Eref &e) const { return vGetY(e); } + +void HHChannelBase::setZ(const Eref &e, double Z) { vSetZ(e, Z); } +double HHChannelBase::getZ(const Eref &e) const { return vGetZ(e); } + +void HHChannelBase::setUseConcentration(const Eref &e, int value) { + vSetUseConcentration(e, value); } -void HHChannelBase::setY(const Eref& e, double Y) -{ - vSetY(e, Y); +int HHChannelBase::getUseConcentration(const Eref &e) const { + return useConcentration_; } -double HHChannelBase::getY(const Eref& e) const -{ - return vGetY(e); + +void HHChannelBase::vSetInstant(const Eref &e, int instant) { + instant_ = instant; } -void HHChannelBase::setZ(const Eref& e, double Z) -{ - vSetZ(e, Z); +int HHChannelBase::vGetInstant(const Eref &e) const { return instant_; } + +void HHChannelBase::vSetX(const Eref &e, double X) { + X_ = X; + xInited_ = true; } -double HHChannelBase::getZ(const Eref& e) const -{ - return vGetZ(e); +double HHChannelBase::vGetX(const Eref &e) const { return X_; } + +void HHChannelBase::vSetY(const Eref &e, double Y) { + Y_ = Y; + yInited_ = true; } +double HHChannelBase::vGetY(const Eref &e) const { return Y_; } -void HHChannelBase::setUseConcentration(const Eref& e, int value) -{ - useConcentration_ = value; - vSetUseConcentration(e, value); +void HHChannelBase::vSetZ(const Eref &e, double Z) { + Z_ = Z; + zInited_ = true; } +double HHChannelBase::vGetZ(const Eref &e) const { return Z_; } -int HHChannelBase::getUseConcentration(const Eref& e) const -{ - return useConcentration_; +void HHChannelBase::vSetUseConcentration(const Eref &e, int value) { + useConcentration_ = value; } // double HHChannelBase::vGetModulation(const Eref& e) const @@ -272,99 +261,152 @@ int HHChannelBase::getUseConcentration(const Eref& e) const // Dest function definitions /////////////////////////////////////////////////// -void HHChannelBase::handleConc(const Eref& e, double conc) -{ +/** + * Returns the state variable for the new timestep based on + * the internal variables A_ and B_ which were passed in from the gate. + */ +double HHChannelBase::integrate(double state, double dt, double A, double B) { + if (B > EPSILON) { + double x = exp(-B * dt); + return state * x + (A / B) * (1 - x); + } + return state + A * dt; +} + +void HHChannelBase::vHandleConc(const Eref &e, double conc) { + cerr << "HHChannelBase::vHandleConc: This function should never be called." + << endl; +} + +void HHChannelBase::handleConc(const Eref &e, double conc) { vHandleConc(e, conc); } /////////////////////////////////////////////////// // HHGate functions +// +// These are breaking the design as the return type is HHGate for +// HHChannel but HHGate2D for HHChannel2D. Making a common HHGateBase +// also turns out to be problematic as the field element can no longer +// be accessed as an HHGate or HHGate2D. + + /////////////////////////////////////////////////// +// HHGate *HHChannelBase::vGetXgate(unsigned int i) const { +// cerr << "HHChannelBase::vGetXgate: This function should never be called." +// << endl; +// return nullptr; +// } -HHGate* HHChannelBase::getXgate(unsigned int i) -{ - return vGetXgate(i); -} +// HHGateBase *HHChannelBase::vGetYgate(unsigned int i) const { +// cerr << "HHChannelBase::vGetYgate: This function should never be called." +// << endl; +// return nullptr; +// } -HHGate* HHChannelBase::getYgate(unsigned int i) -{ - return vGetYgate(i); -} +// HHGateBase *HHChannelBase::vGetZgate(unsigned int i) const { +// cerr << "HHChannelBase::vGetZgate: This function should never be called." +// << endl; +// return nullptr; +// } -HHGate* HHChannelBase::getZgate(unsigned int i) -{ - return vGetZgate(i); -} +// HHGateBase* *HHChannelBase::getXgate(unsigned int i) { return vGetXgate(i); } -void HHChannelBase::setNumGates(unsigned int num) -{ - ; -} +// HHGateBase* *HHChannelBase::getYgate(unsigned int i) { return vGetYgate(i); } -unsigned int HHChannelBase::getNumXgates() const -{ - return (vGetXgate(0) != 0); -} +// HHGateBase* *HHChannelBase::getZgate(unsigned int i) { return vGetZgate(i); } -unsigned int HHChannelBase::getNumYgates() const -{ - return (vGetYgate(0) != 0); -} +// void HHChannelBase::setNumGates(unsigned int num) { ; } + +// unsigned int HHChannelBase::getNumXgates() const { return (vGetXgate(0) != nullptr); } + +// unsigned int HHChannelBase::getNumYgates() const { return (vGetYgate(0) != nullptr); } + +// unsigned int HHChannelBase::getNumZgates() const { return (vGetZgate(0) != nullptr); } -unsigned int HHChannelBase::getNumZgates() const -{ - return (vGetZgate(0) != 0); -} /////////////////////////////////////////////////// // Utility function /////////////////////////////////////////////////// -double HHChannelBase::powerN(double x, double p) -{ - if(x > 0.0) - return exp(p * log(x)); +double HHChannelBase::powerN(double x, double p) { + if (x > 0.0) return exp(p * log(x)); return 0.0; } -PFDD HHChannelBase::selectPower(double power) -{ - if(doubleEq(power, 0.0)) +PFDD HHChannelBase::selectPower(double power) { + if (doubleEq(power, 0.0)) return powerN; - else if(doubleEq(power, 1.0)) + else if (doubleEq(power, 1.0)) return power1; - else if(doubleEq(power, 2.0)) + else if (doubleEq(power, 2.0)) return power2; - else if(doubleEq(power, 3.0)) + else if (doubleEq(power, 3.0)) return power3; - else if(doubleEq(power, 4.0)) + else if (doubleEq(power, 4.0)) return power4; else return powerN; } + +bool HHChannelBase::setGatePower(const Eref &e, double power, double *assignee, + const string &gateType) { + if (power < 0) { + cout << "Error: HHChannel2D::set" << gateType + << "power: Cannot use negative power: " << power << endl; + return false; + } + + if (doubleEq(power, *assignee)) return false; + + if (doubleEq(*assignee, 0.0) && power > 0) { + createGate(e, gateType); + } else if (doubleEq(power, 0.0)) { + destroyGate(e, gateType); + } + + *assignee = power; + return true; +} + +/** + * Assigns the Xpower for this gate. If the gate exists and has + * only this element for input, then change the gate value. + * If the gate exists and has multiple parents, then make a new gate. + * If the gate does not exist, make a new gate + */ +void HHChannelBase::vSetXpower(const Eref &e, double power) { + if (setGatePower(e, power, &Xpower_, "X")) takeXpower_ = selectPower(power); +} + +void HHChannelBase::vSetYpower(const Eref &e, double power) { + if (setGatePower(e, power, &Ypower_, "Y")) takeYpower_ = selectPower(power); +} + +void HHChannelBase::vSetZpower(const Eref &e, double power) { + if (setGatePower(e, power, &Zpower_, "Z")) { + takeZpower_ = selectPower(power); + useConcentration_ = 1; // Not sure about this. + } +} + ///////////////////////////////////////////////////////////////////// // Dummy instantiation, the zombie derivatives make the real function -void HHChannelBase::vSetSolver(const Eref& e, Id hsolve) -{ - ; -} +void HHChannelBase::vSetSolver(const Eref &e, Id hsolve) { ; } -void HHChannelBase::zombify(Element* orig, const Cinfo* zClass, Id hsolve) -{ - if(orig->cinfo() == zClass) - return; +void HHChannelBase::zombify(Element *orig, const Cinfo *zClass, Id hsolve) { + if (orig->cinfo() == zClass) return; unsigned int start = orig->localDataStart(); unsigned int num = orig->numLocalData(); - if(num == 0) - return; + if (num == 0) return; // Parameters are Gbar, Ek, Xpower, Ypower, Zpower, useConcentration // We also want to haul the original gates over, this is done earlier // in the HSolve building process. So just six terms. vector chandata(num * 6, 0.0); vector::iterator j = chandata.begin(); - for(unsigned int i = 0; i < num; ++i) { + for (unsigned int i = 0; i < num; ++i) { Eref er(orig, i + start); - const HHChannelBase* hb = - reinterpret_cast(er.data()); + const HHChannelBase *hb = + reinterpret_cast(er.data()); *j = hb->vGetGbar(er); *(j + 1) = hb->vGetEk(er); *(j + 2) = hb->getXpower(er); @@ -375,9 +417,9 @@ void HHChannelBase::zombify(Element* orig, const Cinfo* zClass, Id hsolve) } orig->zombieSwap(zClass); j = chandata.begin(); - for(unsigned int i = 0; i < num; ++i) { + for (unsigned int i = 0; i < num; ++i) { Eref er(orig, i + start); - HHChannelBase* hb = reinterpret_cast(er.data()); + HHChannelBase *hb = reinterpret_cast(er.data()); hb->vSetSolver(er, hsolve); hb->vSetGbar(er, *j); hb->vSetEk(er, *(j + 1)); diff --git a/biophysics/HHChannelBase.h b/biophysics/HHChannelBase.h index d66896edd5..95a3e2a956 100644 --- a/biophysics/HHChannelBase.h +++ b/biophysics/HHChannelBase.h @@ -14,8 +14,7 @@ #include "ChanCommon.h" - -typedef double ( *PFDD )( double, double ); +typedef double (*PFDD)(double, double); class HHGate; @@ -27,172 +26,211 @@ class HHGate; * fields. */ -class HHChannelBase: public ChanCommon -{ - public: - HHChannelBase(); - virtual ~HHChannelBase() = 0; // this class is not to be instantiated - - ///////////////////////////////////////////////////////////// - // Value field access function definitions - ///////////////////////////////////////////////////////////// - - void setXpower( const Eref& e, double Xpower ); - double getXpower( const Eref& e) const; - void setYpower( const Eref& e, double Ypower ); - double getYpower( const Eref& e) const; - void setZpower( const Eref& e, double Zpower ); - double getZpower( const Eref& e) const; - void setInstant( const Eref& e, int Instant ); - int getInstant( const Eref& e ) const; - void setX( const Eref& e, double X ); - double getX( const Eref& e ) const; - void setY( const Eref& e, double Y ); - double getY( const Eref& e ) const; - void setZ( const Eref& e, double Z ); - double getZ( const Eref& e ) const; - void setUseConcentration( const Eref& e, int value ); - int getUseConcentration( const Eref& e ) const; - // double vGetModulation( const Eref& e ) const; // provided by ChanCommon - ///////////////////////////////////////////////////////////// - // Dest function definitions - ///////////////////////////////////////////////////////////// - /** - * Assign the local conc_ to the incoming conc from the - * concentration calculations for the compartment. Typically - * the message source will be a CaConc object, but there - * are other options for computing the conc. - */ - void handleConc( const Eref& e, double conc ); - - ///////////////////////////////////////////////////////////// - // Gate handling functions - ///////////////////////////////////////////////////////////// - /** - * Access function used for the X gate. The index is ignored. - */ - HHGate* getXgate( unsigned int i ); - - /** - * Access function used for the Y gate. The index is ignored. - */ - HHGate* getYgate( unsigned int i ); - - /** - * Access function used for the Z gate. The index is ignored. - */ - HHGate* getZgate( unsigned int i ); - - /** - * Dummy assignment function for the number of gates. - */ - void setNumGates( unsigned int num ); - - /** - * Access function for the number of Xgates. Gives 1 if present, - * otherwise 0. - */ - unsigned int getNumXgates() const; - /// Returns 1 if Y gate present, otherwise 0 - unsigned int getNumYgates() const; - /// Returns 1 if Z gate present, otherwise 0 - unsigned int getNumZgates() const; - - /** - * Function for safely creating each gate, identified by strings - * as X, Y and Z. Will only work on a new channel, not on a - * copy. The idea is that the gates are always referred to the - * original 'library' channel, and their contents cannot be touched - * except by the original. - */ - void createGate( const Eref& e, string gateType ); - ///////////////////////////////////////////////////////////// - // Virtual Value field access function definitions - ///////////////////////////////////////////////////////////// - virtual void vSetXpower( const Eref& e, double Xpower ) = 0; - virtual void vSetYpower( const Eref& e, double Ypower ) = 0; - virtual void vSetZpower( const Eref& e, double Zpower ) = 0; - // getXpower etc functions are implemented here in the baseclass. - virtual void vSetInstant( const Eref& e, int Instant ) = 0; - virtual int vGetInstant( const Eref& e ) const = 0; - virtual void vSetX( const Eref& e, double X ) = 0; - virtual double vGetX( const Eref& e ) const = 0; - virtual void vSetY( const Eref& e, double Y ) = 0; - virtual double vGetY( const Eref& e ) const = 0; - virtual void vSetZ( const Eref& e, double Z ) = 0; - virtual double vGetZ( const Eref& e ) const = 0; - virtual void vSetUseConcentration( const Eref& e, int value ) = 0; - - ///////////////////////////////////////////////////////////// - // Some more Virtual Value field functions from ChanBase, - // to be defined in derived classes. Listed here for clarity. - ///////////////////////////////////////////////////////////// - // void vSetGbar( double Gbar ); - // double vGetGbar() const; - // void vSetEk( double Ek ); - // double vGetEk() const; - // void vSetGk( double Gk ); - // double vGetGk() const; - // void vSetIk( double Ic ); - // double vGetIk() const; - // void vHandleVm( double Vm ); - - ///////////////////////////////////////////////////////////// - // Virtual Dest function definitions - ///////////////////////////////////////////////////////////// - // void vProcess( const Eref& e, ProcPtr p ); // Listed for clarity - // void vReinit( const Eref& e, ProcPtr p ); // Listed for clarity - - virtual void vHandleConc( const Eref& e, double conc ) = 0; - - ///////////////////////////////////////////////////////////// - // Virtual Gate handling functions - ///////////////////////////////////////////////////////////// - virtual HHGate* vGetXgate( unsigned int i ) const = 0; - virtual HHGate* vGetYgate( unsigned int i ) const = 0; - virtual HHGate* vGetZgate( unsigned int i ) const = 0; - virtual void vCreateGate( const Eref& e, string gateType ) = 0; - - ///////////////////////////////////////////////////////////// - // Utility functions for taking integer powers. - ///////////////////////////////////////////////////////////// - static double power1( double x, double p ) { - return x; - } - static double power2( double x, double p ) { - return x * x; - } - static double power3( double x, double p ) { - return x * x * x; - } - static double power4( double x, double p ) { - return power2( x * x, p ); - } - static double powerN( double x, double p ); - - static PFDD selectPower( double power); - - ///////////////////////////////////////////////////////////// - // Zombification functions. - ///////////////////////////////////////////////////////////// - virtual void vSetSolver( const Eref& e, Id hsolve ); - static void zombify( Element* orig, const Cinfo* zClass, Id hsolve); - - ///////////////////////////////////////////////////////////// - static const Cinfo* initCinfo(); - protected: - /// Exponent for X gate - double Xpower_; - /// Exponent for Y gate - double Ypower_; - /// Exponent for Z gate - double Zpower_; - /// Flag for use of conc for input to Z gate calculations. - bool useConcentration_; - - /// Value used to scale channel conductance up or down - // double modulation_; // this clashes with same field in ChanCommon +class HHChannelBase : public ChanCommon { +public: + HHChannelBase(); + virtual ~HHChannelBase() = 0; // this class is not to be instantiated + + ///////////////////////////////////////////////////////////// + // Value field access function definitions + ///////////////////////////////////////////////////////////// + + void setXpower(const Eref& e, double Xpower); + double getXpower(const Eref& e) const; + void setYpower(const Eref& e, double Ypower); + double getYpower(const Eref& e) const; + void setZpower(const Eref& e, double Zpower); + double getZpower(const Eref& e) const; + void setInstant(const Eref& e, int Instant); + int getInstant(const Eref& e) const; + void setX(const Eref& e, double X); + double getX(const Eref& e) const; + void setY(const Eref& e, double Y); + double getY(const Eref& e) const; + void setZ(const Eref& e, double Z); + double getZ(const Eref& e) const; + void setUseConcentration(const Eref& e, int value); + int getUseConcentration(const Eref& e) const; + // double vGetModulation( const Eref& e ) const; // provided by ChanCommon + ///////////////////////////////////////////////////////////// + // Dest function definitions + ///////////////////////////////////////////////////////////// + /** + * Assign the local conc_ to the incoming conc from the + * concentration calculations for the compartment. Typically + * the message source will be a CaConc object, but there + * are other options for computing the conc. + */ + void handleConc(const Eref& e, double conc); + + ///////////////////////////////////////////////////////////// + // Gate handling functions + ///////////////////////////////////////////////////////////// + /** + * Access function used for the X gate. The index is ignored. + */ + HHGate* getXgate(unsigned int i); + + /** + * Access function used for the Y gate. The index is ignored. + */ + HHGate* getYgate(unsigned int i); + + /** + * Access function used for the Z gate. The index is ignored. + */ + HHGate* getZgate(unsigned int i); + + /** + * Dummy assignment function for the number of gates. + */ + void setNumGates(unsigned int num); + + /** + * Access function for the number of Xgates. Gives 1 if present, + * otherwise 0. + */ + // unsigned int getNumXgates() const; + // /// Returns 1 if Y gate present, otherwise 0 + // unsigned int getNumYgates() const; + // /// Returns 1 if Z gate present, otherwise 0 + // unsigned int getNumZgates() const; + + /** + * Function for safely creating each gate, identified by strings + * as X, Y and Z. Will only work on a new channel, not on a + * copy. The idea is that the gates are always referred to the + * original 'library' channel, and their contents cannot be touched + * except by the original. + */ + void createGate(const Eref& e, string gateType); + /** + * Utility function for destroying gate. Works only on original + * HHChannel. Somewhat dangerous, should never be used after a + * copy has been made as the pointer of the gate will be in use + * elsewhere. + * This needs to be virtual as it needs to call subclass specific delete + * methods for the gate pointer. + */ + virtual void destroyGate(const Eref& e, string gateType); + /** + * Utility for altering gate powers + */ + bool setGatePower(const Eref& e, double power, double* assignee, + const string& gateType); + + ///////////////////////////////////////////////////////////// + // Virtual Value field access function definitions + ///////////////////////////////////////////////////////////// + virtual void vSetXpower(const Eref& e, double Xpower); + virtual void vSetYpower(const Eref& e, double Ypower); + virtual void vSetZpower(const Eref& e, double Zpower); + // getXpower etc functions are implemented here in the baseclass. + virtual void vSetInstant(const Eref& e, int Instant); + virtual int vGetInstant(const Eref& e) const; + virtual void vSetX(const Eref& e, double X); + virtual double vGetX(const Eref& e) const; + virtual void vSetY(const Eref& e, double Y); + virtual double vGetY(const Eref& e) const; + virtual void vSetZ(const Eref& e, double Z); + virtual double vGetZ(const Eref& e) const; + virtual void vSetUseConcentration(const Eref& e, int value); + virtual bool checkOriginal(Id chanId) const; + ///////////////////////////////////////////////////////////// + // Some more Virtual Value field functions from ChanBase, + // to be defined in derived classes. Listed here for clarity. + ///////////////////////////////////////////////////////////// + // void vSetGbar( double Gbar ); + // double vGetGbar() const; + // void vSetEk( double Ek ); + // double vGetEk() const; + // void vSetGk( double Gk ); + // double vGetGk() const; + // void vSetIk( double Ic ); + // double vGetIk() const; + // void vHandleVm( double Vm ); + + ///////////////////////////////////////////////////////////// + // Virtual Dest function definitions + ///////////////////////////////////////////////////////////// + // void vProcess( const Eref& e, ProcPtr p ); // Listed for clarity + // void vReinit( const Eref& e, ProcPtr p ); // Listed for clarity + + virtual void vHandleConc(const Eref& e, double conc); + + ///////////////////////////////////////////////////////////// + // Virtual Gate handling functions + ///////////////////////////////////////////////////////////// + // virtual HHGate* vGetXgate(unsigned int i) const; + // virtual HHGate* vGetYgate(unsigned int i) const; + // virtual HHGate* vGetZgate(unsigned int i) const; + virtual void vCreateGate(const Eref& e, string gateType) = 0; + ///////////////////////////////////////////////////////////// + // Utility functions for taking integer powers. + ///////////////////////////////////////////////////////////// + static double power1(double x, double p) + { + return x; + } + static double power2(double x, double p) + { + return x * x; + } + static double power3(double x, double p) + { + return x * x * x; + } + static double power4(double x, double p) + { + return power2(x * x, p); + } + static double powerN(double x, double p); + + static PFDD selectPower(double power); + + ///////////////////////////////////////////////////////////// + // Zombification functions. + ///////////////////////////////////////////////////////////// + virtual void vSetSolver(const Eref& e, Id hsolve); + static void zombify(Element* orig, const Cinfo* zClass, Id hsolve); + + ///////////////////////////////////////////////////////////// + static const Cinfo* initCinfo(); + +protected: + double integrate(double state, double dt, double A, double B); + /// Function pointers to specific power calculation function. This + /// is an optimization to avoid the full pow function for simple + /// commonly occurring integer exponents + double (*takeXpower_)(double, double); + double (*takeYpower_)(double, double); + double (*takeZpower_)(double, double); + /// Exponent for X gate + double Xpower_; + /// Exponent for Y gate + double Ypower_; + /// Exponent for Z gate + double Zpower_; + /// Flag for use of conc for input to Z gate calculations. + bool useConcentration_; + /// Flag to indicate if gates have been created + bool xInited_, yInited_, zInited_; + /// Bitmasked flag indicating which gates are to be computed as + /// instantaneous (i.e., no integration) + int instant_; + + double X_, Y_, Z_; + double g_; /// conductance without modulation component is different from + /// Gk_ in ChanCommon + + Id myId_; + + static const double EPSILON; + static const int INSTANT_X; + static const int INSTANT_Y; + static const int INSTANT_Z; + /// Value used to scale channel conductance up or down + // double modulation_; // this clashes with same field in ChanCommon }; - -#endif // _HHChannelBase_h +#endif // _HHChannelBase_h diff --git a/biophysics/HHChannelF.cpp b/biophysics/HHChannelF.cpp new file mode 100644 index 0000000000..508da161ed --- /dev/null +++ b/biophysics/HHChannelF.cpp @@ -0,0 +1,295 @@ +// Filename: HHChannelF.cpp +// Description: Formula evaluation based HHChannel +// Author: Subhasis Ray +// Created: Fri Jan 24 13:19:07 2025 (+0530) +// + +#include "ChanBase.h" +#include "ChanCommon.h" +#include "HHChannelBase.h" +#include "HHGateF.h" +#include "HHChannelF.h" + +const Cinfo* HHChannelF::initCinfo() { + static FieldElementFinfo gateX( + "gateX", "Sets up HHGate X for channel", HHGateF::initCinfo(), + &HHChannelF::getXgate, &HHChannelF::setNumGates, &HHChannelF::getNumXgates + // 1 + ); + static FieldElementFinfo gateY( + "gateY", "Sets up HHGateF Y for channel", HHGateF::initCinfo(), + &HHChannelF::getYgate, &HHChannelF::setNumGates, &HHChannelF::getNumYgates + // 1 + ); + static FieldElementFinfo gateZ( + "gateZ", "Sets up HHGateF Z for channel", HHGateF::initCinfo(), + &HHChannelF::getZgate, &HHChannelF::setNumGates, &HHChannelF::getNumZgates + // 1 + ); + /////////////////////////////////////////////////////// + static Finfo* HHChannelFFinfos[] = { + &gateX, // FieldElement + &gateY, // FieldElement + &gateZ // FieldElement + }; + + /////////////////////////////////////////////////////// + static string doc[] = { + "Name", + "HHChannelF", + "Author", + "Subhasis Ray, 2025, CHINTA", + "Description", + "HHChannelF: Hodgkin-Huxley type voltage-gated Ion channel. Unlike " + "HHChannel, which uses table lookup for speed, this version evaluates " + "an expression to compute the gate variables for better accuracy.", + }; + + static Dinfo dinfo; + + static Cinfo HHChannelFCinfo("HHChannelF", HHChannelBase::initCinfo(), + HHChannelFFinfos, + sizeof(HHChannelFFinfos) / sizeof(Finfo*), + &dinfo, doc, sizeof(doc) / sizeof(string)); + + return &HHChannelFCinfo; +} + +static const Cinfo* hhChannelCinfo = HHChannelF::initCinfo(); +////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////// +// Constructor +/////////////////////////////////////////////////// +HHChannelF::HHChannelF() : conc_(0.0), xGate_(0), yGate_(0), zGate_(0) { ; } + +HHChannelF::~HHChannelF() { + ; +} +/** + * If the gate exists and has only this element for input, then change + * the gate power. + * If the gate exists and has multiple parents, then make a new gate, + * set its power. + * If the gate does not exist, make a new gate, set its power. + * + * The function is designed with the idea that if copies of this + * channel are made, then they all point back to the original HHGateF. + * (Unless they are cross-node copies). + * It is only if we subsequently alter the HHGateF of this channel that + * we need to make our own variant of the HHGateF, or disconnect from + * an existing one. + * \todo: May need to convert to handling arrays and Erefs. + */ +// Assuming that the elements are simple elements. Use Eref for +// general case + +bool HHChannelF::checkOriginal(Id chanId) const { + bool isOriginal = true; + if (xGate_) { + isOriginal = xGate_->isOriginalChannel(chanId); + } else if (yGate_) { + isOriginal = yGate_->isOriginalChannel(chanId); + } else if (zGate_) { + isOriginal = zGate_->isOriginalChannel(chanId); + } + return isOriginal; +} + +void HHChannelF::innerCreateGate(const string& gateName, HHGateF** gatePtr, + Id chanId, Id gateId) { + // Shell* shell = reinterpret_cast< Shell* >( ObjId( Id(), 0 ).data() ); + if (*gatePtr) { + cout << "Warning: HHChannelF::createGate: '" << gateName + << "' on Element '" << chanId.path() << "' already present\n"; + return; + } + *gatePtr = new HHGateF(chanId, gateId); +} + +void HHChannelF::vCreateGate(const Eref& e, string gateType) { + if (!checkOriginal(e.id())) { + cout << "Warning: HHChannelF::createGate: Not allowed from copied " + "channel:\n" + << e.id().path() << "\n"; + return; + } + + if (gateType == "X") + innerCreateGate("xGate", &xGate_, e.id(), Id(e.id().value() + 1)); + else if (gateType == "Y") + innerCreateGate("yGate", &yGate_, e.id(), Id(e.id().value() + 2)); + else if (gateType == "Z") + innerCreateGate("zGate", &zGate_, e.id(), Id(e.id().value() + 3)); + else + cout << "Warning: HHChannelF::createGate: Unknown gate type '" + << gateType << "'. Ignored\n"; +} + +void HHChannelF::innerDestroyGate(const string& gateName, HHGateF** gatePtr, + Id chanId) { + if (*gatePtr == nullptr) { + cout << "Warning: HHChannelF::destroyGate: '" << gateName + << "' on Element '" << chanId.path() << "' not present\n"; + return; + } + delete (*gatePtr); + *gatePtr = nullptr; +} + +void HHChannelF::destroyGate(const Eref& e, string gateType) { + if (!checkOriginal(e.id())) { + cout << "Warning: HHChannelF::destroyGate: Not allowed from copied " + "channel:\n" + << e.id().path() << "\n"; + return; + } + + if (gateType == "X") + innerDestroyGate("xGate", &xGate_, e.id()); + else if (gateType == "Y") + innerDestroyGate("yGate", &yGate_, e.id()); + else if (gateType == "Z") + innerDestroyGate("zGate", &zGate_, e.id()); + else + cout << "Warning: HHChannelF::destroyGate: Unknown gate type '" + << gateType << "'. Ignored\n"; +} +/////////////////////////////////////////////////// +// HHGateF functions +// +// These are breaking the design as the return type is HHGateF for +// HHChannelF but HHGateF2D for HHChannelF2D. Making a common HHGateBase +// also turns out to be problematic as the field element can no longer +// be accessed as an HHGateF or HHGateF2D. +/////////////////////////////////////////////////// + +HHGateF* HHChannelF::getXgate(unsigned int i) { return xGate_; } + +HHGateF* HHChannelF::getYgate(unsigned int i) { return yGate_; } + +HHGateF* HHChannelF::getZgate(unsigned int i) { return zGate_; } + +void HHChannelF::setNumGates(unsigned int num) { ; } + +unsigned int HHChannelF::getNumXgates() const { return xGate_ != nullptr; } + +unsigned int HHChannelF::getNumYgates() const { return yGate_ != nullptr; } + +unsigned int HHChannelF::getNumZgates() const { return zGate_ != nullptr; } + + +/////////////////////////////////////////////////// +// Dest function definitions +/////////////////////////////////////////////////// + +void HHChannelF::vProcess(const Eref& e, ProcPtr info) { + g_ += ChanCommon::vGetGbar(e); + double A = 0; + double B = 0; + if (Xpower_ > 0) { + xGate_->lookupBoth(Vm_, &A, &B); + if (instant_ & INSTANT_X) + X_ = A / B; + else + X_ = integrate(X_, info->dt, A, B); + g_ *= takeXpower_(X_, Xpower_); + } + + if (Ypower_ > 0) { + yGate_->lookupBoth(Vm_, &A, &B); + if (instant_ & INSTANT_Y) + Y_ = A / B; + else + Y_ = integrate(Y_, info->dt, A, B); + + g_ *= takeYpower_(Y_, Ypower_); + } + + if (Zpower_ > 0) { + if (useConcentration_) + zGate_->lookupBoth(conc_, &A, &B); + else + zGate_->lookupBoth(Vm_, &A, &B); + if (instant_ & INSTANT_Z) + Z_ = A / B; + else + Z_ = integrate(Z_, info->dt, A, B); + + g_ *= takeZpower_(Z_, Zpower_); + } + ChanCommon::vSetGk(e, g_ * ChanCommon::vGetModulation(e)); + ChanCommon::updateIk(); + // Gk_ = g_; + // Ik_ = ( Ek_ - Vm_ ) * g_; + + // Send out the relevant channel messages. + ChanCommon::sendProcessMsgs(e, info); + + g_ = 0.0; +} + +/** + * Here we get the steady-state values for the gate (the 'instant' + * calculation) as A_/B_. + */ +void HHChannelF::vReinit(const Eref& er, ProcPtr info) { + g_ = ChanCommon::vGetGbar(er); + Element* e = er.element(); + + double A = 0.0; + double B = 0.0; + if (Xpower_ > 0) { + assert(xGate_); + xGate_->lookupBoth(Vm_, &A, &B); + if (B < EPSILON) { + cout << "Warning: B_ value for " << e->getName() + << " is ~0. Check X table\n"; + return; + } + if (!xInited_) X_ = A / B; + g_ *= takeXpower_(X_, Xpower_); + } + + if (Ypower_ > 0) { + assert(yGate_); + yGate_->lookupBoth(Vm_, &A, &B); + if (B < EPSILON) { + cout << "Warning: B value for " << e->getName() + << " is ~0. Check Y table\n"; + return; + } + if (!yInited_) Y_ = A / B; + g_ *= takeYpower_(Y_, Ypower_); + } + + if (Zpower_ > 0) { + assert(zGate_); + if (useConcentration_) + zGate_->lookupBoth(conc_, &A, &B); + else + zGate_->lookupBoth(Vm_, &A, &B); + if (B < EPSILON) { + cout << "Warning: B value for " << e->getName() + << " is ~0. Check Z table\n"; + return; + } + if (!zInited_) Z_ = A / B; + g_ *= takeZpower_(Z_, Zpower_); + } + ChanCommon::vSetGk(er, g_ * ChanCommon::vGetModulation(er)); + ChanCommon::updateIk(); + // Gk_ = g_; + // Ik_ = ( Ek_ - Vm_ ) * g_; + + // Send out the relevant channel messages. + // Same for reinit as for process. + ChanCommon::sendReinitMsgs(er, info); + + g_ = 0.0; +} + +void HHChannelF::vHandleConc(const Eref& e, double conc) { conc_ = conc; } + +// +// HHChannelF.cpp ends here diff --git a/biophysics/HHChannelF.h b/biophysics/HHChannelF.h new file mode 100644 index 0000000000..74bdc89f6e --- /dev/null +++ b/biophysics/HHChannelF.h @@ -0,0 +1,124 @@ +/* HHChannelF.h --- + * + * Filename: HHChannelF.h + * Description: + * Author: Subhasis Ray + * Created: Fri Jan 24 13:09:34 2025 (+0530) + */ + +#ifndef _HHChannelF_h +#define _HHChannelF_h + +#include "HHChannelBase.h" + +class HHGateF; + +class HHChannelF : public HHChannelBase { +#ifdef DO_UNIT_TESTS + /* friend void testHHChannelF(); */ + /* friend void testHHGateFCreation(); */ +#endif +public: + HHChannelF(); + ~HHChannelF(); + void innerSetXpower(double Xpower); + void innerSetYpower(double Ypower); + void innerSetZpower(double Zpower); + ///////////////////////////////////////////////////////////// + // Dest function definitions + ///////////////////////////////////////////////////////////// + + /** + * processFunc handles the update and calculations every + * clock tick. It first sends the request for evaluation of + * the gate variables to the respective gate objects and + * recieves their response immediately through a return + * message. This is done so that many channel instances can + * share the same gate lookup tables, but do so cleanly. + * Such messages should never go to a remote node. + * Then the function does its own little calculations to + * send back to the parent compartment through regular + * messages. + */ + void vProcess(const Eref& e, ProcPtr p) override; + + /** + * Reinitializes the values for the channel. This involves + * computing the steady-state value for the channel gates + * using the provided Vm from the parent compartment. It + * involves a similar cycle through the gates and then + * updates to the parent compartment as for the processFunc. + */ + void vReinit(const Eref& e, ProcPtr p) override; + + /** + * Assign the local Vm_ to the incoming Vm from the compartment + void handleVm( double Vm ); + */ + + /** + * Assign the local conc_ to the incoming conc from the + * concentration calculations for the compartment. Typically + * the message source will be a CaConc object, but there + * are other options for computing the conc. + */ + void vHandleConc(const Eref& e, double conc) override; + + ///////////////////////////////////////////////////////////// + // Gate handling functions + ///////////////////////////////////////////////////////////// + HHGateF* getXgate(unsigned int i); + HHGateF* getYgate(unsigned int i); + HHGateF* getZgate(unsigned int i); + + void setNumGates(unsigned int num); + unsigned int getNumXgates() const; + unsigned int getNumYgates() const; + unsigned int getNumZgates() const; + + /// Inner utility function for creating the gate. + void innerCreateGate(const string& gateName, HHGateF** gatePtr, Id chanId, + Id gateId); + + /// Returns true if channel is original, false if copy. + bool checkOriginal(Id chanId) const override; + + void vCreateGate(const Eref& e, string gateType) override; + void destroyGate(const Eref& e, string gateType) override; + + /** + * Inner utility for destroying the gate + */ + void innerDestroyGate(const string& gateName, HHGateF** gatePtr, Id chanId); + + // /** + // * Utility for altering gate powers + // */ + // bool setGatePower(const Eref& e, double power, double* assignee, + // const string& gateType); + + ///////////////////////////////////////////////////////////// + static const Cinfo* initCinfo(); + + private: + /// Conc_ is input variable for Ca-dependent channels. + double conc_; + + /** + * HHGate data structure for the xGate. This is writable only + * on the HHChannel that originally created the HHGate, for others + * it must be treated as readonly. + */ + HHGateF* xGate_; + + /// HHGate data structure for the yGate. + HHGateF* yGate_; + + /// HHGate data structure for the yGate. + HHGateF* zGate_; + +}; + +#endif + +/* HHChannelF.h ends here */ diff --git a/biophysics/HHGate.cpp b/biophysics/HHGate.cpp index 930a66e592..7bd4d4e439 100644 --- a/biophysics/HHGate.cpp +++ b/biophysics/HHGate.cpp @@ -9,6 +9,8 @@ #include "../basecode/header.h" #include "../basecode/ElementValueFinfo.h" +#include "../builtins/MooseParser.h" +#include "HHGateBase.h" #include "HHGate.h" static const double SINGULARITY = 1.0e-6; @@ -199,25 +201,23 @@ static const Cinfo* hhGateCinfo = HHGate::initCinfo(); // Core class functions /////////////////////////////////////////////////// HHGate::HHGate() - : xmin_(0), + : HHGateBase(0, 0), + xmin_(0), xmax_(1), invDx_(1), - originalChanId_(0), - originalGateId_(0), lookupByInterpolation_(0), isDirectTable_(0) { - ; + cerr << "# HHGate::HHGate() should never be called" << endl; } HHGate::HHGate(Id originalChanId, Id originalGateId) - : A_(1, 0.0), + : HHGateBase(originalChanId, originalGateId), + A_(1, 0.0), B_(1, 0.0), xmin_(0), xmax_(1), invDx_(1), - originalChanId_(originalChanId), - originalGateId_(originalGateId), lookupByInterpolation_(0), isDirectTable_(0) { @@ -772,35 +772,6 @@ void HHGate::tabFill(vector& table, unsigned int newXdivs, lookupByInterpolation_ = origLookupMode; } -bool HHGate::checkOriginal(Id id, const string& field) const -{ - if(id == originalGateId_) - return 1; - - cout << "Warning: HHGate: attempt to set field '" << field << "' on " - << id.path() << ", which is not the original Gate element. Ignored.\n"; - return 0; -} - -bool HHGate::isOriginalChannel(Id id) const -{ - return (id == originalChanId_); -} - -bool HHGate::isOriginalGate(Id id) const -{ - return (id == originalGateId_); -} - -Id HHGate::originalChannelId() const -{ - return originalChanId_; -} - -Id HHGate::originalGateId() const -{ - return originalGateId_; -} void HHGate::updateAlphaBeta() { diff --git a/biophysics/HHGate.h b/biophysics/HHGate.h index f9f4ce7ea2..f274128d8d 100644 --- a/biophysics/HHGate.h +++ b/biophysics/HHGate.h @@ -9,6 +9,9 @@ #ifndef _HHGate_h #define _HHGate_h + +#include "HHGateBase.h" + /** * This class handles a single gate on an HHChannel. It is equivalent to the * m and h terms on the Hodgkin-Huxley Na channel, or the n term on the @@ -29,7 +32,7 @@ * original HHChannel, but all the others do have read permission. */ -class HHGate { +class HHGate: public HHGateBase { friend void testHHGateLookup(); friend void testHHGateSetup(); @@ -120,35 +123,6 @@ class HHGate { */ double lookupTable(const vector& tab, double v) const; - /** - * Checks if the provided Id is the one that the HHGate was created - * on. If true, fine, otherwise complains about trying to set the - * field. - */ - bool checkOriginal(Id id, const string& field) const; - - /** - * isOriginalChannel returns true if the provided Id is the Id of - * the channel on which the HHGate was created. - */ - bool isOriginalChannel(Id id) const; - - /** - * isOriginalChannel returns true if the provided Id is the Id of - * the Gate created at the same time as the original channel. - */ - bool isOriginalGate(Id id) const; - - /** - * Returns the Id of the original Channel. - */ - Id originalChannelId() const; - - /** - * Returns the Id of the original Gate. - */ - Id originalGateId() const; - /** * tabFill does interpolation and range resizing for * a table representing a lookup function. diff --git a/biophysics/HHGate2D.cpp b/biophysics/HHGate2D.cpp index 29fff06462..5185660162 100644 --- a/biophysics/HHGate2D.cpp +++ b/biophysics/HHGate2D.cpp @@ -10,6 +10,7 @@ #include "../basecode/header.h" #include "../basecode/ElementValueFinfo.h" #include "../builtins/Interpol2D.h" +#include "HHGateBase.h" #include "HHGate2D.h" static const double SINGULARITY = 1.0e-6; @@ -129,8 +130,8 @@ HHGate2D::HHGate2D() HHGate2D::HHGate2D( Id originalChanId, Id originalGateId ) : - originalChanId_( originalChanId ), - originalGateId_( originalGateId ) + HHGateBase(originalChanId, + originalGateId ) {;} /////////////////////////////////////////////////// diff --git a/biophysics/HHGate2D.h b/biophysics/HHGate2D.h index fe920c9c04..f33732ccb2 100644 --- a/biophysics/HHGate2D.h +++ b/biophysics/HHGate2D.h @@ -10,7 +10,9 @@ #ifndef _HHGate2D_h #define _HHGate2D_h -class HHGate2D +#include "HHGateBase.h" + +class HHGate2D: public HHGateBase { public: HHGate2D(); diff --git a/biophysics/HHGateBase.cpp b/biophysics/HHGateBase.cpp new file mode 100644 index 0000000000..ec9d3f5ed9 --- /dev/null +++ b/biophysics/HHGateBase.cpp @@ -0,0 +1,62 @@ +/********************************************************************** + ** This program is part of 'MOOSE', the + ** Messaging Object Oriented Simulation Environment. + ** Copyright (C) 2003-2007 Upinder S. Bhalla. and NCBS + ** It is made available under the terms of the + ** GNU Lesser General Public License version 2.1 + ** See the file COPYING.LIB for the full notice. + **********************************************************************/ + +#include "../basecode/header.h" +#include "../basecode/ElementValueFinfo.h" +#include "HHGateBase.h" + +/////////////////////////////////////////////////// +// Core class functions +/////////////////////////////////////////////////// +HHGateBase::HHGateBase() : originalChanId_(0), originalGateId_(0) +{ + cerr << "# HHGateBase::HHGateBase() should never be called" << endl; +} + +HHGateBase::HHGateBase(Id originalChanId, Id originalGateId) + : originalChanId_(originalChanId), originalGateId_(originalGateId) +{ + // cerr << "# HHGateBase::HHGateBase(): originalChanId:" << originalChanId << ", originalGateId: " << originalGateId << endl; + ; +} + +/////////////////////////////////////////////////////////////////////// +// Utility funcs +/////////////////////////////////////////////////////////////////////// + +bool HHGateBase::checkOriginal(Id id, const string& field) const +{ + if(id == originalGateId_) + return true; + + cout << "Warning: HHGateBase: attempt to set field '" << field << "' on " + << id.path() << ", which is not the original Gate element. Ignored.\n"; + return false; +} + +bool HHGateBase::isOriginalChannel(Id id) const +{ + // cerr << "# Received: " << id << ", original: " << originalChanId_ << endl; + return (id == originalChanId_); +} + +bool HHGateBase::isOriginalGate(Id id) const +{ + return (id == originalGateId_); +} + +Id HHGateBase::originalChannelId() const +{ + return originalChanId_; +} + +Id HHGateBase::originalGateId() const +{ + return originalGateId_; +} diff --git a/biophysics/HHGateBase.h b/biophysics/HHGateBase.h new file mode 100644 index 0000000000..afb24786d6 --- /dev/null +++ b/biophysics/HHGateBase.h @@ -0,0 +1,87 @@ +/********************************************************************** +** This program is part of 'MOOSE', the +** Messaging Object Oriented Simulation Environment. +** copyright (C) 2003-2011 Upinder S. Bhalla. and NCBS +** It is made available under the terms of the +** GNU Lesser General Public License version 2.1 +** See the file COPYING.LIB for the full notice. +**********************************************************************/ +#ifndef _HHGateBase_h +#define _HHGateBase_h + +/** + * This class is the base for gates in Hodgkin-Huxley type channels + * (HHGates). HHGates are typically shared. This means that when you + * make a copy or a vector of an HHChannel, there is only a single + * HHGate created, and its pointer is used by all the copies. The + * lookup functions are thread-safe. Field assignment to the HHGate + * should be possible only from the original HHChannel, but all the + * others do have read permission. + * This base class only factors out the tracking of the ownership. + */ + +class HHGateBase { +public: + /** + * Dummy constructor, to keep Dinfo happy. Should never be used + */ + HHGateBase(); + + /** + * This constructor is the one meant to be used. It takes the + * originalId of the parent HHChannel as a required argument, + * so that any subsequent 'write' functions can be checked to + * see if they are legal. Also tracks its own Id. + */ + HHGateBase(Id originalChanId, Id originalGateId); + + ///////////////////////////////////////////////////////////////// + // Utility funcs + ///////////////////////////////////////////////////////////////// + /** + * Checks if the provided Id is the one that the HHGate was created + * on. If true, fine, otherwise complains about trying to set the + * field. + */ + bool checkOriginal(Id id, const string& field) const; + + /** + * isOriginalChannel returns true if the provided Id is the Id of + * the channel on which the HHGate was created. + */ + bool isOriginalChannel(Id id) const; + + /** + * isOriginalChannel returns true if the provided Id is the Id of + * the Gate created at the same time as the original channel. + */ + bool isOriginalGate(Id id) const; + + /** + * Returns the Id of the original Channel. + */ + Id originalChannelId() const; + + /** + * Returns the Id of the original Gate. + */ + Id originalGateId() const; + + static const Cinfo* initCinfo(); + +protected: + /** + * Id of original channel, the one which has actually allocated it, + * All other Elements have to treat the values as readonly. + */ + Id originalChanId_; + + /** + * Id of original Gate, the one which was created with the original + * channel. + * All other Elements have to treat the values as readonly. + */ + Id originalGateId_; +}; + +#endif // _HHGateBase_h diff --git a/biophysics/HHGateF.cpp b/biophysics/HHGateF.cpp new file mode 100644 index 0000000000..392e468fcc --- /dev/null +++ b/biophysics/HHGateF.cpp @@ -0,0 +1,227 @@ +/********************************************************************** + ** This program is part of 'MOOSE', the + ** Messaging Object Oriented Simulation Environment. + ** Copyright (C) 2003-2007 Upinder S. Bhalla. and NCBS + ** It is made available under the terms of the + ** GNU Lesser General Public License version 2.1 + ** See the file COPYING.LIB for the full notice. + **********************************************************************/ + +#include "../basecode/header.h" +#include "../basecode/ElementValueFinfo.h" +#include "HHGateF.h" + +const Cinfo* HHGateF::initCinfo() +{ + /////////////////////////////////////////////////////// + // Field definitions. + /////////////////////////////////////////////////////// + static ReadOnlyLookupValueFinfo A( + "A", + "lookupA: Compute the A gate value from a double. " + "This is done by evaluating the expressions for alpha/beta" + " or tau/inf.", + &HHGateF::lookupA); + static ReadOnlyLookupValueFinfo B( + "B", + "lookupB: Look up the B gate value from a double." + "This is done by evaluating the expressions for alpha/beta" + " or tau/inf.", + &HHGateF::lookupB); + + static ElementValueFinfo alpha( + "alpha", + "Expression for voltage-dependent rates, alpha. " + "This requires the expression for beta to be defined as well.", + &HHGateF::setAlpha, &HHGateF::getAlpha); + + static ElementValueFinfo beta( + "beta", + "Expression for voltage-dependent rates, beta. " + "This requires the expression for alpha to be defined as well.", + &HHGateF::setBeta, &HHGateF::getBeta); + + static ElementValueFinfo tau( + "tau", + "Expression for voltage-dependent rates, tau. " + "This requires the expression for mInfinity to be defined as well.", + &HHGateF::setTau, &HHGateF::getTau); + + static ElementValueFinfo mInfinity( + "mInfinity", + "Expression for voltage-dependent rates, mInfinity. " + "This requires the expression for tau to be defined as well.", + &HHGateF::setMinfinity, &HHGateF::getMinfinity); + + /////////////////////////////////////////////////////// + // DestFinfos + /////////////////////////////////////////////////////// + static Finfo* HHGateFFinfos[] = { + &A, // ReadOnlyLookupValue + &B, // ReadOnlyLookupValue + &alpha, // Value + &beta, // Value + &tau, // Value + &mInfinity, // Value + }; + + static string doc[] = { + "Name", + "HHGateF", + "Author", + "Subhasis Ray, 2025, CHINTA", + "Description", + "HHGateF: Gate for Hodkgin-Huxley type channels, equivalent to the " + "m and h terms on the Na squid channel and the n term on K. " + "This takes the voltage and state variable from the channel, " + "computes the new value of the state variable and a scaling, " + "depending on gate power, for the conductance. As opposed to HHGate, " + "which uses lookup tables for speed, this evaluates explicit " + "expressions for accuracy. This is a single variable gate, either " + "voltage or concentration. So the expression also allows only one " + "indpendent variable, which is assumed `v`. See the documentation of " + "``Function`` class for details on the praser.", + }; + + static Dinfo dinfo; + static Cinfo HHGateFCinfo("HHGateF", Neutral::initCinfo(), HHGateFFinfos, + sizeof(HHGateFFinfos) / sizeof(Finfo*), &dinfo, + doc, sizeof(doc) / sizeof(string)); + + return &HHGateFCinfo; +} + +static const Cinfo* hhGateCinfo = HHGateF::initCinfo(); +/////////////////////////////////////////////////// +// Core class functions +/////////////////////////////////////////////////// +HHGateF::HHGateF() : HHGateBase(0, 0) +{ + cerr << "Warning: HHGateF::HHGateF(): this should never be called" << endl; +} + +HHGateF::HHGateF(Id originalChanId, Id originalGateId) + : HHGateBase(originalChanId, originalGateId) + +{ + symTab_.add_variable("v", v_); + symTab_.add_constants(); + alpha_.register_symbol_table(symTab_); + beta_.register_symbol_table(symTab_); +} + +HHGateF& HHGateF::operator=(const HHGateF& rhs) +{ + // protect from self-assignment. + if( this == &rhs) + return *this; + + v_ = rhs.v_; + symTab_.add_variable("v_", v_); + symTab_.add_constants(); + alpha_.register_symbol_table(symTab_); + beta_.register_symbol_table(symTab_); + alphaExpr_ = rhs.alphaExpr_; + betaExpr_ = rhs.betaExpr_; + parser_.compile(alphaExpr_, alpha_); + parser_.compile(betaExpr_, beta_); + tauInf_ = rhs.tauInf_; + return *this; +} + +/////////////////////////////////////////////////// +// Field function definitions +/////////////////////////////////////////////////// + +double HHGateF::lookupA(double v) const +{ + // TODO: check for divide by zero? + v_ = v; + return tauInf_ ? beta_.value() / alpha_.value() : alpha_.value(); +} + +double HHGateF::lookupB(double v) const +{ + // TODO: check for divide by zero? + v_ = v; + return tauInf_ ? 1.0 / alpha_.value() : alpha_.value() + beta_.value(); +} + +void HHGateF::lookupBoth(double v, double* A, double* B) const +{ + *A = lookupA(v); + *B = lookupB(v); + cerr << "# HHGateF::lookupBoth: v=" << v << ", A=" << *A << ", B="<< *B << endl; +} + +void HHGateF::setAlpha(const Eref& e, const string expr) +{ + if(checkOriginal(e.id(), "alpha")) { + if(!parser_.compile(expr, alpha_)) { + cerr << "Error: HHGateF::setAlpha: cannot compile expression!\n" + << parser_.error() << endl; + return; + } + tauInf_ = false; + alphaExpr_ = expr; + } +} + +string HHGateF::getAlpha(const Eref& e) const +{ + return tauInf_ ? "" : alphaExpr_; +} + +void HHGateF::setBeta(const Eref& e, const string expr) +{ + if(checkOriginal(e.id(), "beta")) { + if(!parser_.compile(expr, beta_)) { + cerr << "Error: HHGateF::setBeta: cannot compile expression!\n" + << parser_.error() << endl; + return; + } + tauInf_ = false; + betaExpr_ = expr; + } +} + +string HHGateF::getBeta(const Eref& e) const +{ + return tauInf_ ? "" : betaExpr_; +} + +void HHGateF::setTau(const Eref& e, const string expr) +{ + if(checkOriginal(e.id(), "alpha")) { + if(!parser_.compile(expr, alpha_)) { + cerr << "Error: HHGateF::setTau: cannot compile expression!\n" + << parser_.error() << endl; + return; + } + tauInf_ = true; + alphaExpr_ = expr; + } +} + +string HHGateF::getTau(const Eref& e) const +{ + return tauInf_ ? alphaExpr_ : ""; +} + +void HHGateF::setMinfinity(const Eref& e, const string expr) +{ + if(checkOriginal(e.id(), "beta")) { + if(!parser_.compile(expr, beta_)) { + cerr << "Error: HHGateF::setMinfinity: cannot compile expression!\n" + << parser_.error() << endl; + return; + } + tauInf_ = true; + betaExpr_ = expr; + } +} + +string HHGateF::getMinfinity(const Eref& e) const +{ + return tauInf_ ? betaExpr_ : ""; +} diff --git a/biophysics/HHGateF.h b/biophysics/HHGateF.h new file mode 100644 index 0000000000..ff2763d5f6 --- /dev/null +++ b/biophysics/HHGateF.h @@ -0,0 +1,119 @@ +/********************************************************************** +** This program is part of 'MOOSE', the +** Messaging Object Oriented Simulation Environment. +** copyright (C) 2003-2011 Upinder S. Bhalla. and NCBS +** It is made available under the terms of the +** GNU Lesser General Public License version 2.1 +** See the file COPYING.LIB for the full notice. +**********************************************************************/ +#ifndef _HHGateF_h +#define _HHGateF_h + +#include "exprtk.hpp" +#include "../basecode/global.h" + +#include "../basecode/header.h" +#include "../basecode/ElementValueFinfo.h" +#include "HHGateBase.h" + +/** + * This class handles a single gate on an HHChannel. It is equivalent to the + * m and h terms on the Hodgkin-Huxley Na channel, or the n term on the + * K channel. It stores the + * voltage-dependence (sometimes concentration-dependence) of the + * gating variables for opening the channel. It does so in a tabular form + * which can be directly filled using experimental data points. + * It also provides a set of + * utility functions for defining the gate in functional forms, and + * accessing those original functional forms. + * The HHGateF is + * accessed as a FieldElement, which means that it is available as a + * pointer on the HHChannel. HHGateFs are typically shared. This means that + * when you make a copy or a vector of an HHChannel, there is only a single + * HHGateF created, and its pointer is used by all the copies. + * The lookup functions are thread-safe. + * Field assignment to the HHGateF should be possible only from the + * original HHChannel, but all the others do have read permission. + * Whereas HHGate uses interpolation tables, HHGateF uses direct + * formula evaluation, hence slower but possibly more accurate. + */ + +class HHGateF : public HHGateBase { +public: + /** + * Dummy constructor, to keep Dinfo happy. Should never be used + */ + HHGateF(); + + /** + * This constructor is the one meant to be used. It takes the + * originalId of the parent HHChannel as a required argument, + * so that any subsequent 'write' functions can be checked to + * see if they are legal. Also tracks its own Id. + */ + HHGateF(Id originalChanId, Id originalGateId); + /// HHGates remain shared between copies of a channel, so it + /// should never be copied. Yet we need to define this because + /// eprtk parser deletes its copy assignment, which deletes + /// compiler generated copy assignment operator for HHGateF, which + /// raises error with Dinfo, which tries to reference the copy + /// operator. + HHGateF& operator=(const HHGateF&); + ////////////////////////////////////////////////////////// + // LookupValueFinfos + ////////////////////////////////////////////////////////// + /** + * lookupA: Look up the A vector from a double. Typically does + * so by direct scaling and offset to an integer lookup, using + * a fine enough table granularity that there is little error. + * Alternatively uses linear interpolation. + * The range of the double is predefined based on knowledge of + * voltage or conc ranges, and the granularity is specified by + * the xmin, xmax, and invDx fields. + */ + double lookupA(double v) const; + + /** + * lookupB: Look up the B vector from a double, similar to lookupA. + */ + double lookupB(double v) const; + ////////////////////////////////////////////////////////// + // DestFinfos + ////////////////////////////////////////////////////////// + /** + * Single call to get both A and B values by lookup + */ + void lookupBoth(double v, double* A, double* B) const; + /// Set the expression for evaluating alpha + void setAlpha(const Eref& e, const string expr); + string getAlpha(const Eref& e) const; + /// Set the expression for evaluating beta + void setBeta(const Eref& e, const string expr); + string getBeta(const Eref& e) const; + /// Set the expression for evaluating tau + void setTau(const Eref& e, const string expr); + string getTau(const Eref& e) const; + /// Set the expression for evaluating mInfinity + void setMinfinity(const Eref& e, const string expr); + string getMinfinity(const Eref& e) const; + + ///////////////////////////////////////////////////////////////// + // Utility funcs + ///////////////////////////////////////////////////////////////// + static const Cinfo* initCinfo(); + +private: + /// Whether the gate is expressed in tau-inf form. If false, it is + /// alpha-beta form + bool tauInf_; + exprtk::symbol_table symTab_; + exprtk::expression alpha_; + exprtk::expression beta_; + exprtk::parser parser_; + mutable double v_; + /// Store the user-specified expression strings + string alphaExpr_; + string betaExpr_; +}; + +#endif // _HHGateF_h diff --git a/biophysics/meson.build b/biophysics/meson.build index 268cf2097e..7dfd5436f9 100644 --- a/biophysics/meson.build +++ b/biophysics/meson.build @@ -11,12 +11,14 @@ biophysics_src = ['IntFire.cpp', 'GapJunction.cpp', 'ChanBase.cpp', 'ChanCommon.cpp', - 'HHChannel.cpp', 'HHChannelBase.cpp', + 'HHChannel.cpp', 'HHChannel2D.cpp', + 'HHChannelF.cpp', + 'HHGateBase.cpp', 'HHGate.cpp', 'HHGate2D.cpp', - 'HHChannel2D.cpp', + 'HHGateF.cpp', 'CaConcBase.cpp', 'CaConc.cpp', 'MgBlock.cpp', @@ -45,4 +47,6 @@ biophysics_src = ['IntFire.cpp', 'MarkovOdeSolver.cpp', 'testBiophysics.cpp'] -biophysics_lib = static_library('biophysics', biophysics_src, dependencies: gsl_dep, include_directories: gsl_dep.get_variable(pkgconfig:'includedir')) +include_directories = ['../external/exprtk'] +include_directories += gsl_dep.get_variable(pkgconfig:'includedir') +biophysics_lib = static_library('biophysics', biophysics_src, dependencies: gsl_dep, include_directories: include_directories) diff --git a/builtins/Function.cpp b/builtins/Function.cpp index 4f9b5a492c..a9df69f32b 100644 --- a/builtins/Function.cpp +++ b/builtins/Function.cpp @@ -295,12 +295,11 @@ const Cinfo * Function::initCinfo() derivativeOut(), }; - static string doc[] = - { - "Name", "Function", - "Author", "Subhasis Ray/Dilawar Singh", + static string doc[] = { + "Name", "Function", "Author", "Subhasis Ray/Dilawar Singh", "Description", - R""""(General purpose function calculator using real numbers. + R"( +General purpose function calculator using real numbers. It can parse mathematical expression defining a function and evaluate it and/or its derivative for specified variable values. You can assign expressions of @@ -322,7 +321,8 @@ specified in the expression as y{i} and connecting the messages should happen in the same order as the y indices. This class handles only real numbers (C-double). Predefined constants - are: pi=3.141592..., e=2.718281...)"""" + are: pi=3.141592..., e=2.718281... +)" }; static Dinfo< Function > dinfo; @@ -575,7 +575,7 @@ void Function::setExpr(const Eref& eref, const string expression) * * @Param eref * @Param expr Expression to set. - * @Param dynamicLookup Weather to allow unknown variables in the expression. + * @Param dynamicLookup Whether to allow unknown variables in the expression. * (default to true in moose>=4.0.0) * * @Returns True if compilation was successful. diff --git a/builtins/Function.h b/builtins/Function.h index 618cd59281..57f43ba0b7 100644 --- a/builtins/Function.h +++ b/builtins/Function.h @@ -41,9 +41,9 @@ class Function // get a list of variable identifiers. vector getVars() const; - void setVarValues(vector vars, vector vals); + void setVarValues(vector vars, vector vals); // subha: where is this defined? - // get/set the value of variable `name` + /// set the value of `index`-th variable void setVar(unsigned int index, double value); Variable* getX(unsigned int ii); diff --git a/hsolve/ZombieHHChannel.cpp b/hsolve/ZombieHHChannel.cpp index d58fc3641e..2488cb4cfe 100644 --- a/hsolve/ZombieHHChannel.cpp +++ b/hsolve/ZombieHHChannel.cpp @@ -191,20 +191,20 @@ void ZombieHHChannel::vCreateGate(const Eref& e, string name) // HHGate functions /////////////////////////////////////////////////// -HHGate* ZombieHHChannel::vGetXgate( unsigned int i ) const -{ - return 0; -} - -HHGate* ZombieHHChannel::vGetYgate( unsigned int i ) const -{ - return 0; -} - -HHGate* ZombieHHChannel::vGetZgate( unsigned int i ) const -{ - return 0; -} +// HHGate* ZombieHHChannel::vGetXgate( unsigned int i ) const +// { +// return 0; +// } + +// HHGate* ZombieHHChannel::vGetYgate( unsigned int i ) const +// { +// return 0; +// } + +// HHGate* ZombieHHChannel::vGetZgate( unsigned int i ) const +// { +// return 0; +// } /////////////////////////////////////////////////// // Assign solver diff --git a/hsolve/ZombieHHChannel.h b/hsolve/ZombieHHChannel.h index fc70fe8dbb..fcbd65f725 100644 --- a/hsolve/ZombieHHChannel.h +++ b/hsolve/ZombieHHChannel.h @@ -99,20 +99,20 @@ class ZombieHHChannel: public HHChannelBase ///////////////////////////////////////////////////////////// // Gate handling functions ///////////////////////////////////////////////////////////// - /** - * Access function used for the X gate. The index is ignored. - */ - HHGate* vGetXgate( unsigned int i ) const override; - - /** - * Access function used for the Y gate. The index is ignored. - */ - HHGate* vGetYgate( unsigned int i ) const override; - - /** - * Access function used for the Z gate. The index is ignored. - */ - HHGate* vGetZgate( unsigned int i ) const override; + // /** + // * Access function used for the X gate. The index is ignored. + // */ + // HHGate* vGetXgate( unsigned int i ) const override; + + // /** + // * Access function used for the Y gate. The index is ignored. + // */ + // HHGate* vGetYgate( unsigned int i ) const override; + + // /** + // * Access function used for the Z gate. The index is ignored. + // */ + // HHGate* vGetZgate( unsigned int i ) const override; ///////////////////////////////////////////////////////////// void vSetSolver( const Eref& e , Id hsolve ) override; diff --git a/pybind11/Finfo.cpp b/pybind11/Finfo.cpp index 5b40c38690..9bad327ac2 100644 --- a/pybind11/Finfo.cpp +++ b/pybind11/Finfo.cpp @@ -29,7 +29,7 @@ std::function getSetGetFunc1(const ObjId &oid, const string &fname) std::function func = [oid, fname](const T &val) { return SetGet1::set(oid, fname, val); }; - std::cout << "getSetGet1Func" << std::endl; + // std::cout << "getSetGet1Func" << std::endl; return func; } diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp index db1553864b..89fa8560c4 100644 --- a/scheduling/Clock.cpp +++ b/scheduling/Clock.cpp @@ -396,6 +396,7 @@ const Cinfo* Clock::initCinfo() " GapJunction 4 50e-6\n" " HHChannel 4 50e-6\n" " HHChannel2D 4 50e-6\n" + " HHChannelF 4 50e-6\n" " Leakage 4 50e-6\n" " MarkovChannel 4 50e-6\n" " SpikeGen 5 50e-6\n" @@ -406,7 +407,7 @@ const Cinfo* Clock::initCinfo() " Dsolve 10 0.01\n" " Adaptor 11 0.1\n" - " Func 12 0.1\n" + // " Func 12 0.1\n" " Function 12 0.1\n" " Arith 12 0.1\n" " Gsolve (init) 15 0.1\n" @@ -913,6 +914,7 @@ void Clock::buildDefaultTick() defaultTick_["GapJunction"] = 4; defaultTick_["HHChannel"] = 4; defaultTick_["HHChannel2D"] = 4; + defaultTick_["HHChannelF"] = 4; defaultTick_["Leakage"] = 4; defaultTick_["MarkovChannel"] = 4; defaultTick_["SpikeGen"] = 5; @@ -922,7 +924,7 @@ void Clock::buildDefaultTick() defaultTick_["TimeTable"] = 8; defaultTick_["Dsolve"] = 10; defaultTick_["Adaptor"] = 11; - defaultTick_["Func"] = 12; + // defaultTick_["Func"] = 12; // as of 2025 this class has been removed defaultTick_["Function"] = 12; defaultTick_["Arith"] = 12; defaultTick_["Gsolve"] = 16; // Note this uses an 'init' at t-1 diff --git a/tests/core/ephys.py b/tests/core/ephys.py new file mode 100644 index 0000000000..b0a3cbe30d --- /dev/null +++ b/tests/core/ephys.py @@ -0,0 +1,113 @@ +# ephys.py --- +# +# Filename: ephys.py +# Description: +# Author: Subhasis Ray +# Created: Wed Feb 5 16:00:48 2025 (+0530) +# + +# Code: +"""Utility functions for ephys tests""" +import moose + + +def create_voltage_clamp( + compartment, + modelpath='./elec', + datapath=None, + vclampname='vclamp', + cmdname='command', +): + """Creates a voltage clamp object under `modelpath` and + a table under `datapath` to record the command voltage. + + Parameters + ---------- + compartment: moose.Compartment + Compartment to be voltage clamped + modelpath: str (default: './elec') + Path to container for the electrical circuit + datapath: moose.melement (default: None) + Container for data recorders, if `None` the command voltage is not + recorded, and the returned `command_tab` is also `None`. + vclampname: str (default: 'vclamp' + Name of the voltage clamp object. + cmdname: str (default: 'command') + Name of the command pulse generator object. + + Returns + ------- + (vclamp, command, commandtab): `vclamp` is the VoltageClamp object, + `command` the `moose.PulseGen` that sets the command value to the + voltage clamp, `command_tab` a `moose.Table` that records the + command value during simulation. If `datapath` is `None`, then this + will be `None`. + + """ + # create the elec container if it does not exist + _ = moose.Neutral(modelpath) + vclamp_path = f'{modelpath}/{vclampname}' + existent = False + if moose.exists(vclamp_path): + # avoid duplicate connect + existent = True + print( + f'{modelpath}: Object already exists. ' + 'Returning references to existing components.' + ) + _ = moose.Neutral(modelpath) + vclamp = moose.VClamp(vclamp_path) + command = moose.PulseGen(f'{modelpath}/{cmdname}') + # Also setup a table to record the command voltage of the VClamp directly + if datapath is not None: + commandtab = moose.Table(f'{datapath}/command') + else: + commandtab = None + if not existent: + # The voltage clamp's output is `currentOut` which will be + # injected into the compartment + moose.connect(vclamp, 'currentOut', compartment, 'injectMsg') + # The voltage clamp object senses the voltage of the compartment + moose.connect(compartment, 'VmOut', vclamp, 'sensedIn') + # Connect the output of the command pulse to the command input of + # the voltage clamp circuit + moose.connect(command, 'output', vclamp, 'commandIn') + if commandtab is not None: + moose.connect(commandtab, 'requestOut', command, 'getOutputValue') + # ==================================================== + # set the parameters for voltage clamp circuit + # ---------------------------------------------------- + # compartment.dt is the integration time step for the + # compartment. `tau` is the time constant of the RC filter in + # the circuit. 5 times the integration timestep value is a + # good starting point for tau + vclamp.tau = 5 * compartment.dt + # `gain` is the proportional gain of the PID + # controller. `Cm/dt` is a good value + vclamp.gain = compartment.Cm / compartment.dt + # `ti` is the integral time of the PID controller, `dt` is a good value + vclamp.ti = compartment.dt + # `td` is the derivative time of the PID controller. We can + # keep it 0, the default value + return vclamp, command, commandtab + + +def setup_step_command(command, base, delay, level): + """Set up an existing pulse generator `command` to output `base` + as initial value and `level` after `delay` time + + This provides a pulse that is a single step function. If you want + repeated pulses directly modify the vectors of delay, width, and + level. + + """ + command.baseLevel = base + command.firstDelay = delay + command.secondDelay = 1e9 # Never stop + command.firstWidth = 1e9 # Never stop + command.firstLevel = level + + + +# +# ephys.py ends here diff --git a/tests/core/test_hhchan.py b/tests/core/test_hhchan.py new file mode 100644 index 0000000000..1b788cca96 --- /dev/null +++ b/tests/core/test_hhchan.py @@ -0,0 +1,228 @@ +# Filename: test_hhchan.py +# Description: +# Author: Subhasis Ray +# Created: Thu Jan 23 13:33:22 2025 (+0530) +# + +# Code: +"""Tests for HHChannel class""" + +import moose +import math + +from ephys import create_voltage_clamp, setup_step_command + + +def test_hhgate_creation(container='test'): + _ = moose.Neutral(container) + moose.ce(container) + ch = moose.HHChannel('ch') + ch.Xpower = 1 + ch.Ypower = 2 + assert math.isclose(ch.Xpower, 1) + assert moose.exists(f'ch/gateX') + assert math.isclose(ch.Ypower, 2), 'Ypower not set' + assert moose.exists(f'ch/gateY'), 'gateY object does not exist' + moose.ce('..') + moose.delete(container) + +def test_hh_k_vclamp(container='test', steptime=5.0): + """Simulate a voltage clamp with Hodhkin and Huxley's K+ channel. + + Parameters + ---------- + steptime: float + Time of the voltage step + """ + _ = moose.Neutral(container) + moose.ce(container) + comp = moose.Compartment('comp0') + chan = moose.HHChannel(f'{comp.path}/K') + moose.connect(chan, 'channel', comp, 'channel') + chan.Gbar = 36.0 + chan.Ek = -12.0 + chan.Xpower = 4 + n_gate = moose.element(f'{chan.path}/gateX') + vdivs = 150 + vmin = -30.0 + vmax = 120.0 + x_alpha_params = [0.1, -0.01, -1.0, -10.0, -10.0] + x_beta_params = [0.125, 0, 0, 0, 80.0] + # Note that `+` operator with lists as operands concatenates them + x_params = x_alpha_params + x_beta_params + [vdivs, vmin, vmax] + n_gate.setupAlpha(x_params) + n_gate.useInterpolation = True + comp.Em = 0 # Hodgkin and Huxley used resting voltage as 0 + comp.Vm = 0 + comp.initVm = 0 + comp.Cm = 1 + comp.Rm = 1 / 0.3 # G_leak is 0.3 mS/cm^2 + dt = 0.01 + for tick in range(8): + moose.setClock(tick, dt) + # This must be after setting clock dt because parameters depend on dt + vclamp, command, commandtab = create_voltage_clamp(comp) + # Hodgkin and Huxley 1952, "Currents carried by sodium and + # potassium ions ...", find maximum conductance + simtime = 100.0 + steptime + # Precalculated steady state K conductance + vm_gk = { + 0: 0.3666444556069115, + # 10: 1.8401, # this is a discontinuity + 20: 5.287062775435651, + 30: 10.176967265218654, + 40: 15.220230442170081, + 50: 19.596740114488156, + 60: 23.1009379487905, + 70: 25.821065780738117, + 80: 27.91721528657564, + 90: 29.537324370047777, + 100: 30.79812392487768, + } + + for vstep, gk in vm_gk.items(): + setup_step_command(command, 0.0, delay=steptime, level=vstep) + moose.reinit() + moose.start(simtime) + assert math.isclose( + gk, chan.Gk, abs_tol=1e-6 + ), f'Vm={vstep} Gk={chan.Gk}, expected={gk}' + moose.ce('..') + moose.delete(container) + + +def test_hh_na_vclamp(container='test', steptime=5.0): + """Simulate a voltage clamp with Hodhkin and Huxley's K+ channel. + + Parameters + ---------- + steptime: float + Time of the voltage step + """ + _ = moose.Neutral(container) + moose.ce(container) + comp = moose.Compartment('comp0') + chan = moose.HHChannel(f'{comp.path}/Na') + moose.connect(chan, 'channel', comp, 'channel') + chan.Gbar = 120.0 + chan.Ek = 115.0 + chan.Xpower = 3 + chan.Ypower = 1 + vdivs = 150 + vmin = -30.0 + vmax = 120.0 + m_gate = moose.element(f'{chan.path}/gateX') + x_alpha_params = [2.5, -0.1, -1.0, -25.0, -10.0] + x_beta_params = [4, 0, 0, 0, 18.0] + # Note that `+` operator with lists as operands concatenates them + x_params = x_alpha_params + x_beta_params + [vdivs, vmin, vmax] + m_gate.setupAlpha(x_params) + m_gate.useInterpolation = True + h_gate = moose.element(f'{chan.path}/gateY') + y_alpha_params = [0.07, 0, 0, 0, 20.0] + y_beta_params = [1, 0, 1, -30, -10.0] + # Note that `+` operator with lists as operands concatenates them + y_params = y_alpha_params + y_beta_params + [vdivs, vmin, vmax] + h_gate.setupAlpha(y_params) + h_gate.useInterpolation = True + comp.Em = 0 # Hodgkin and Huxley used resting voltage as 0 + comp.Vm = 0 + comp.initVm = 0 + comp.Cm = 1 + comp.Rm = 1 / 0.3 # G_leak is 0.3 mS/cm^2 + dt = 0.01 + for tick in range(8): + moose.setClock(tick, dt) + # This must be after setting clock dt because parameters depend on dt + vclamp, command, commandtab = create_voltage_clamp(comp) + # Hodgkin and Huxley 1952, "Currents carried by sodium and + # potassium ions ...", find maximum conductance + simtime = 100.0 + steptime + # Precalculated steady state Na conductance + vm_gna = { + 0: 0.010609192838829854, + 10: 0.12443211458767735, + 20: 0.527787746095522, + 30: 0.8966173682113173, + 40: 0.8361209504788413, + 50: 0.5983994566094422, + 60: 0.3893933499995894, + 70: 0.24432304663465443, + 80: 0.15080726000137287, + 90: 0.09232255449417581, + 100: 0.0562750224683144, + } + + for vstep, gna in vm_gna.items(): + setup_step_command(command, 0.0, delay=steptime, level=vstep) + moose.reinit() + moose.start(simtime) + assert math.isclose( + gna, chan.Gk, abs_tol=1e-6 + ), f'Vm={vstep} Gk={chan.Gk}, expected={gna}' + moose.ce('..') + moose.delete(container) + + +def test_hhchan_eval(container='test'): + """Test the evaluation of hhchannel conductance using Hodgkin and + Huxleys Na channel model + + """ + _ = moose.Neutral(container) + moose.ce(container) + comp = moose.Compartment('comp0') + chan = moose.HHChannel(f'{comp.path}/Na') + moose.connect(chan, 'channel', comp, 'channel') + chan.Gbar = 120.0 + chan.Ek = 115.0 + chan.Xpower = 3 + chan.Ypower = 1 + m_gate = moose.element(f'{chan.path}/gateX') + h_gate = moose.element(f'{chan.path}/gateY') + vdivs = 150 + vmin = -30.0 + vmax = 120.0 + m_gate = moose.element(f'{chan.path}/gateX') + x_alpha_params = [2.5, -0.1, -1.0, -25.0, -10.0] + x_beta_params = [4, 0, 0, 0, 18.0] + # Note that `+` operator with lists as operands concatenates them + x_params = x_alpha_params + x_beta_params + [vdivs, vmin, vmax] + m_gate.setupAlpha(x_params) + m_gate.useInterpolation = True + h_gate = moose.element(f'{chan.path}/gateY') + y_alpha_params = [0.07, 0, 0, 0, 20.0] + y_beta_params = [1, 0, 1, -30, -10.0] + # Note that `+` operator with lists as operands concatenates them + y_params = y_alpha_params + y_beta_params + [vdivs, vmin, vmax] + h_gate.setupAlpha(y_params) + h_gate.useInterpolation = True + comp.Em = 0 # Hodgkin and Huxley used resting voltage as 0 + comp.Vm = 0 + comp.initVm = 0 + comp.Cm = 1 + comp.Rm = 1 / 0.3 # G_leak is 0.3 mS/cm^2 + moose.reinit() + t = 0 + tdelta = comp.dt * 1000 + for ii in range(10): + moose.start(tdelta) + t += tdelta + print(f'time={t} Gk={chan.Gk} Vm={comp.Vm}') + comp.inject = 1.0 + print('Starting current injection...') + for ii in range(10): + moose.start(tdelta) + t += tdelta + print(f'time={t} Gk={chan.Gk} Vm={comp.Vm}') + moose.ce('..') + moose.delete(container) + + +if __name__ == '__main__': + test_hhgate_creation() + test_hh_k_vclamp() + test_hh_na_vclamp() + test_hhchan_eval() +# +# test_hhchan.py ends here diff --git a/tests/core/test_hhchanf.py b/tests/core/test_hhchanf.py new file mode 100644 index 0000000000..86a908557d --- /dev/null +++ b/tests/core/test_hhchanf.py @@ -0,0 +1,228 @@ +# Filename: test_hhchanf.py +# Description: +# Author: Subhasis Ray +# Created: Wed Jan 29 13:15:18 2025 (+0530) +# + +"""Tests HHChannelF class""" +import moose +import math +import pytest +from ephys import create_voltage_clamp, setup_step_command + + +def test_hhgatef_creation(container='test'): + _ = moose.Neutral(container) + moose.ce(container) + ch = moose.HHChannelF('ch0') + ch.Xpower = 1 + ch.Ypower = 2 + assert math.isclose(ch.Xpower, 1), 'Xpower not set' + assert moose.exists('ch0/gateX'), 'gateX object does not exist' + assert math.isclose(ch.Ypower, 2), 'Ypower not set' + assert moose.exists('ch0/gateY'), 'gateY object does not exist' + moose.ce('..') + moose.delete(container) + + +def test_hhgatef_alpha_beta(container='test'): + """Test set/get alpha and beta expressions""" + _ = moose.Neutral(container) + moose.ce(container) + ch = moose.HHChannelF('ch1') + ch.Xpower = 1 + xgate = moose.element('ch1/gateX') + alpha = '0.01 * (10 - v) / exp((10 - v)/10 - 1)' + beta = '0.125 * exp( - v/80)' + xgate.alpha = alpha + xgate.beta = beta + assert xgate.alpha == alpha, 'alpha not set' + assert xgate.beta == beta, 'beta not set' + assert xgate.tau == '', 'tau not reset' + assert xgate.mInfinity == '', 'mInfinity not reset' + moose.ce('..') + moose.delete(container) + + +def test_hhgatef_tau_inf(container='test'): + """Test set/get tau and mInfinity expressions""" + _ = moose.Neutral(container) + moose.ce(container) + ch = moose.HHChannelF('ch2') + ch.Xpower = 1 + xgate = moose.element('ch2/gateX') + alpha = '0.01 * (10 - v) / exp((10 - v)/10 - 1)' + beta = '0.125 * exp( - v/80)' + xgate.tau = alpha + xgate.mInfinity = beta + assert xgate.tau == alpha, 'tau not set' + assert xgate.mInfinity == beta, 'mInfinity not set' + assert xgate.alpha == '', 'alpha not reset' + assert xgate.beta == '', 'beta not reset' + moose.ce('..') + moose.delete(container) + +def test_hh_k_vclamp(container='test', steptime=5.0): + """Simulate a voltage clamp with Hodhkin and Huxley's K+ channel. + + Parameters + ---------- + steptime: float + Time of the voltage step + """ + _ = moose.Neutral(container) + moose.ce(container) + comp = moose.Compartment('comp0') + chan = moose.HHChannelF(f'{comp.path}/K') + moose.connect(chan, 'channel', comp, 'channel') + chan.Gbar = 36.0 + chan.Ek = -12.0 + chan.Xpower = 4 + n_gate = moose.element(f'{chan.path}/gateX') + n_gate.alpha = '0.01 * (10 - v) / (exp((10-v)/10) - 1)' + n_gate.beta = '0.125 * exp(-v/80)' + comp.Em = 0 # Hodgkin and Huxley used resting voltage as 0 + comp.Vm = 0 + comp.initVm = 0 + comp.Cm = 1 + comp.Rm = 1 / 0.3 # G_leak is 0.3 mS/cm^2 + dt = 0.01 + for tick in range(8): + moose.setClock(tick, dt) + # vclamp needs to be setup after clock dt because it uses dt + vclamp, command, commandtab = create_voltage_clamp(comp) + # Hodgkin and Huxley 1952, "Currents carried by sodium and + # potassium ions ...", find maximum conductance try 0-100 mV at 10 mV steps + simtime = 100.0 + steptime + # Precalculated steady state K conductance + vm_gk = { + 0: 0.3666444556069115, + # 10: 1.8401, + 20: 5.287062775435651, + 30: 10.176967265218654, + 40: 15.220230442170081, + 50: 19.596740114488156, + 60: 23.1009379487905, + 70: 25.821065780738117, + 80: 27.91721528657564, + 90: 29.537324370047777, + 100: 30.79812392487768, + } + + for vstep, gk in vm_gk.items(): + setup_step_command(command, 0.0, delay=steptime, level=vstep) + moose.reinit() + moose.start(simtime) + assert math.isclose( + gk, chan.Gk, abs_tol=1e-6 + ), f'Vm={vstep} Gk={chan.Gk}, expected={gk}' + moose.ce('..') + moose.delete(container) + + +def test_hh_na_vclamp(container='test', steptime=5.0): + """Test the evaluation of hhchannel conductance using Hodgkin and + Huxleys Na channel model""" + _ = moose.Neutral(container) + moose.ce(container) + comp = moose.Compartment('comp0') + chan = moose.HHChannelF(f'{comp.path}/Na') + moose.connect(chan, 'channel', comp, 'channel') + chan.Gbar = 120.0 + chan.Ek = 115.0 + chan.Xpower = 3 + chan.Ypower = 1 + m_gate = moose.element(f'{chan.path}/gateX') + h_gate = moose.element(f'{chan.path}/gateY') + m_gate.alpha = '0.1 * (25 - v) / (exp((25 - v) / 10) - 1)' + m_gate.beta = '4 * exp(- v / 18)' + h_gate.alpha = '0.07 * exp(- v / 20)' + h_gate.beta = '1 / (exp((30 - v) / 10) + 1)' + comp.Em = 0 # Hodgkin and Huxley used resting voltage as 0 + comp.Vm = 0 + comp.initVm = 0 + comp.Cm = 1 + comp.Rm = 1 / 0.3 # G_leak is 0.3 mS/cm^2 + dt = 0.01 + for tick in range(8): + moose.setClock(tick, dt) + # This must be after setting clock dt because parameters depend on dt + vclamp, command, commandtab = create_voltage_clamp(comp) + # Hodgkin and Huxley 1952, "Currents carried by sodium and + # potassium ions ...", find maximum conductance + simtime = 100.0 + steptime + # Precalculated steady state Na conductance + vm_gna = { + 0: 0.010609192838829854, + 10: 0.12443211458767735, + 20: 0.527787746095522, + 30: 0.8966173682113173, + 40: 0.8361209504788413, + 50: 0.5983994566094422, + 60: 0.3893933499995894, + 70: 0.24432304663465443, + 80: 0.15080726000137287, + 90: 0.09232255449417581, + 100: 0.0562750224683144, + } + + for vstep, gna in vm_gna.items(): + setup_step_command(command, 0.0, delay=steptime, level=vstep) + moose.reinit() + moose.start(simtime) + assert math.isclose( + gna, chan.Gk, abs_tol=1e-6 + ), f'Vm={vstep} Gk={chan.Gk}, expected={gna}' + moose.ce('..') + moose.delete(container) + + +def test_hhchanf_eval(container='test'): + """Test the evaluation of hhchannel conductance using Hodgkin and + Huxleys Na channel model""" + _ = moose.Neutral(container) + moose.ce(container) + comp = moose.Compartment('comp0') + chan = moose.HHChannelF(f'{comp.path}/Na') + moose.connect(chan, 'channel', comp, 'channel') + chan.Gbar = 120.0 + chan.Ek = 115.0 + chan.Xpower = 3 + chan.Ypower = 1 + m_gate = moose.element(f'{chan.path}/gateX') + h_gate = moose.element(f'{chan.path}/gateY') + m_gate.alpha = '0.1 * (25 - v) / (exp((25 - v) / 10) - 1)' + m_gate.beta = '4 * exp(- v / 18)' + h_gate.alpha = '0.07 * exp(- v / 20)' + h_gate.beta = '1 / (exp((30 - v) / 10) + 1)' + comp.Em = 0 # Hodgkin and Huxley used resting voltage as 0 + comp.Vm = 0 + comp.initVm = 0 + comp.Cm = 1 + comp.Rm = 1 / 0.3 # G_leak is 0.3 mS/cm^2 + moose.reinit() + t = 0 + tdelta = comp.dt * 1000 + for ii in range(10): + moose.start(tdelta) + t += tdelta + print(f'time={t} Gk={chan.Gk} Vm={comp.Vm}') + comp.inject = 1.0 + print('Starting current injection...') + for ii in range(10): + moose.start(tdelta) + t += tdelta + print(f'time={t} Gk={chan.Gk} Vm={comp.Vm}') + moose.ce('..') + moose.delete(container) + + +if __name__ == '__main__': + test_hhgatef_creation() + test_hhgatef_alpha_beta() + test_hhgatef_tau_inf() + test_hh_k_vclamp() + test_hh_na_vclamp() + test_hhchanf_eval() +# +# test_hhchanf.py ends here