From a5122cfb6dd82f1a831f64d1fb1a29a7d1392334 Mon Sep 17 00:00:00 2001
From: Dilawar Singh <dilawars@ncbs.res.in>
Date: Mon, 23 Jan 2017 13:35:27 +0530
Subject: [PATCH] Squashed 'moose-core/' changes from 1b5b9a9..14e2d78
14e2d78 Fixes broken Makefile. Object MMPump.o not found error is resolved by commenting out the object. Either source is missing for it was introduced by recent PR.
aeddce7 Fixes to python3 test command in travis script. Previous fixes passed with python2.
00ad9e2 Fixes to recent Travis build failure. (#170)
3ebc0a5 Merge commit '5bceacfe2bebf128999e0565521b7e49ff2f6f25'
df5ff48 Merge branch 'difshell' of https://github.com/asiaszmek/moose into asiaszmek-difshell
e4a4bcd Change order of operations to slightly improve accuracy
aacebcc Small changes to make the code more alike genesis and other moose objects
4246932 Difshells and difbuffers initialize close to equilibrium
40a200f Add forgotten setting of helper variables to zero after calculations
dad35e6 Add state in t(i-1)
45c2286 Fix wrong definitions
567444a Fix errors in volume calculation
4e0a857 Indentation
cb40719 Add MMPump to scheduling and makefiles
bfc456c Add MMPump to DifShell
be1ceac Implementation of MMPump
c721713 Add ticks for DifBufferBase and DifShellBase
a3bee38 Fix bugs: DifShell now compiles and seems to give correct results
b131c44 Fix a typo in zombify
f1cbfb0 Rewrite DifShell. Inherit from DifShellBase
1486cdd Write base class for DifShell
336cbd1 Indentation
d8ee906 Implement exponential Euler for DifShell
f72d6a1 Move source messages to class. Correctly initialize Faraday's constant
b309dc4 Move source messages to class
5d3b7d2 Make DifBuffer look more like other moose classes and fix ticks
9a1e166 Add DifBuffer to clocks
f949b27 Add DifBuffer to compilation
5418e0b Fix no ticks when initializing difshells warning
5fba26e Remove useless variable
7e29b64 Reintroduce an old version of DifBuffer from moose
8f544e6 Moved to cmake folder.
ae1a5ba Merge commit '0d12ebb03bf7c3efce4b08bcbe40ed82137fa202'
3ac6c58 master branch.
65c371c Fixes to BhallaLab/moose#204.
git-subtree-dir: moose-core
git-subtree-split: 14e2d78a7b6e41e67bdce6e3a07fa0d65565744b
---
.travis/travis_build_linux.sh | 4 +-
CMakeLists.txt | 1 +
VERSION | 1 +
biophysics/CMakeLists.txt | 8 +-
biophysics/CaConcBase.cpp | 486 +++++------
biophysics/DifBuffer.cpp | 467 ++++++++++
biophysics/DifBuffer.h | 96 +++
biophysics/DifBufferBase.cpp | 461 ++++++++++
biophysics/DifBufferBase.h | 134 +++
biophysics/DifShell.cpp | 907 +++++++-------------
biophysics/DifShell.h | 253 +++---
biophysics/DifShellBase.cpp | 513 +++++++++++
biophysics/DifShellBase.h | 144 ++++
biophysics/HHGate.cpp | 1222 +++++++++++++--------------
biophysics/MMPump.cpp | 141 ++++
biophysics/MMPump.h | 45 +
biophysics/Makefile | 10 +-
python/moose/chemUtil/graphUtils.py | 16 +-
python/moose/moose.py | 8 +-
scheduling/Clock.cpp | 8 +
20 files changed, 3321 insertions(+), 1604 deletions(-)
create mode 100644 VERSION
create mode 100644 biophysics/DifBuffer.cpp
create mode 100644 biophysics/DifBuffer.h
create mode 100644 biophysics/DifBufferBase.cpp
create mode 100644 biophysics/DifBufferBase.h
create mode 100644 biophysics/DifShellBase.cpp
create mode 100644 biophysics/DifShellBase.h
create mode 100644 biophysics/MMPump.cpp
create mode 100644 biophysics/MMPump.h
diff --git a/.travis/travis_build_linux.sh b/.travis/travis_build_linux.sh
index 302dd129..cf339969 100755
--- a/.travis/travis_build_linux.sh
+++ b/.travis/travis_build_linux.sh
@@ -31,7 +31,7 @@ PYTHON3="/usr/bin/python3"
(
# Old makefile based flow.
python2 -m compileall -q .
- if type python3 -c 'import sys' > /dev/null; then python3 -m compileall -q . ; fi
+ if type $PYTHON3 > /dev/null; then python3 -m compileall -q . ; fi
# Traditional make.
make
## CMAKE based flow
@@ -46,7 +46,7 @@ PYTHON3="/usr/bin/python3"
# This is only applicable on linux build.
echo "Python3 support. Removed python2-networkx and install python3"
- if type $PYTHON3 -c 'import os' > /dev/null; then
+ if type $PYTHON3 > /dev/null; then
sudo apt-get remove -qq python-networkx
sudo apt-get install -qq python3-networkx
mkdir -p _GSL_BUILD2 && cd _GSL_BUILD2 && \
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 75f76d73..aa1c243d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -219,6 +219,7 @@ if(HDF5_FOUND)
)
set(HDF5_LIBRARIES ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES})
endif()
+
# Make sure, HDF5_HL_LIBRARIES are set. The COMPONENTS in find_package may
# or may not work. See BhallaLab/moose-core#163.
diff --git a/VERSION b/VERSION
new file mode 100644
index 00000000..a7ea0627
--- /dev/null
+++ b/VERSION
@@ -0,0 +1 @@
+3.1.1-73-g3ebc0a5
\ No newline at end of file
diff --git a/biophysics/CMakeLists.txt b/biophysics/CMakeLists.txt
index 27c34e12..c36def8e 100644
--- a/biophysics/CMakeLists.txt
+++ b/biophysics/CMakeLists.txt
@@ -33,8 +33,12 @@ set(BIOPHYSICS_SRCS
SynChan.cpp
NMDAChan.cpp
testBiophysics.cpp
- IzhikevichNrn.cpp
- DifShell.cpp
+ IzhikevichNrn.cpp
+ DifShellBase.cpp
+ DifShell.cpp
+ DifBufferBase.cpp
+ DifBuffer.cpp
+ MMPump.cpp
Leakage.cpp
VectorTable.cpp
MarkovRateTable.cpp
diff --git a/biophysics/CaConcBase.cpp b/biophysics/CaConcBase.cpp
index 2d66b709..c484b248 100644
--- a/biophysics/CaConcBase.cpp
+++ b/biophysics/CaConcBase.cpp
@@ -1,11 +1,11 @@
/**********************************************************************
-** This program is part of 'MOOSE', the
-** Messaging Object Oriented Simulation Environment.
-** Copyright (C) 2003-2014 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.
-**********************************************************************/
+ ** This program is part of 'MOOSE', the
+ ** Messaging Object Oriented Simulation Environment.
+ ** Copyright (C) 2003-2014 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 <cfloat>
#include "header.h"
@@ -23,178 +23,178 @@
*/
// Static function.
SrcFinfo1< double >* CaConcBase::concOut() {
- static SrcFinfo1< double > concOut( "concOut",
- "Concentration of Ca in pool" );
- return &concOut;
+ static SrcFinfo1< double > concOut( "concOut",
+ "Concentration of Ca in pool" );
+ return &concOut;
}
const Cinfo* CaConcBase::initCinfo()
{
- ///////////////////////////////////////////////////////
- // Shared message definitions
- ///////////////////////////////////////////////////////
- static DestFinfo process( "process",
- "Handles process call",
- new ProcOpFunc< CaConcBase >( &CaConcBase::process ) );
- static DestFinfo reinit( "reinit",
- "Handles reinit call",
- new ProcOpFunc< CaConcBase >( &CaConcBase::reinit ) );
-
- static Finfo* processShared[] =
- {
- &process, &reinit
- };
-
- static SharedFinfo proc( "proc",
- "Shared message to receive Process message from scheduler",
- processShared, sizeof( processShared ) / sizeof( Finfo* ) );
+ ///////////////////////////////////////////////////////
+ // Shared message definitions
+ ///////////////////////////////////////////////////////
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< CaConcBase >( &CaConcBase::process ) );
+ static DestFinfo reinit( "reinit",
+ "Handles reinit call",
+ new ProcOpFunc< CaConcBase >( &CaConcBase::reinit ) );
+
+ static Finfo* processShared[] =
+ {
+ &process, &reinit
+ };
+
+ static SharedFinfo proc( "proc",
+ "Shared message to receive Process message from scheduler",
+ processShared, sizeof( processShared ) / sizeof( Finfo* ) );
-///////////////////////////////////////////////////////
-// Field definitions
-///////////////////////////////////////////////////////
- static ElementValueFinfo< CaConcBase, double > Ca( "Ca",
- "Calcium concentration.",
- &CaConcBase::setCa,
- &CaConcBase::getCa
- );
- static ElementValueFinfo< CaConcBase, double > CaBasal( "CaBasal",
- "Basal Calcium concentration.",
- &CaConcBase::setCaBasal,
- &CaConcBase::getCaBasal
- );
- static ElementValueFinfo< CaConcBase, double > Ca_base( "Ca_base",
- "Basal Calcium concentration, synonym for CaBasal",
- &CaConcBase::setCaBasal,
- &CaConcBase::getCaBasal
- );
- static ElementValueFinfo< CaConcBase, double > tau( "tau",
- "Settling time for Ca concentration",
- &CaConcBase::setTau,
- &CaConcBase::getTau
- );
- static ElementValueFinfo< CaConcBase, double > B( "B",
- "Volume scaling factor. "
- "Deprecated. This is a legacy field from GENESIS and exposes "
- "internal calculations. Please do not use. \n"
- "B = 1/(vol * F) so that it obeys:\n"
- "dC/dt = B*I_Ca - C/tau",
- &CaConcBase::setB,
- &CaConcBase::getB
- );
- static ElementValueFinfo< CaConcBase, double > thick( "thick",
- "Thickness of Ca shell, assumed cylindrical. Legal range is "
- "between zero and the radius. If outside this range it is "
- "taken as the radius. Default zero, ie, the shell is the entire "
- "thickness of the cylinder",
- &CaConcBase::setThickness,
- &CaConcBase::getThickness
- );
- static ElementValueFinfo< CaConcBase, double > length( "length",
- "Length of Ca shell, assumed cylindrical",
- &CaConcBase::setLength,
- &CaConcBase::getLength
- );
- static ElementValueFinfo< CaConcBase, double > diameter( "diameter",
- "Diameter of Ca shell, assumed cylindrical",
- &CaConcBase::setDiameter,
- &CaConcBase::getDiameter
- );
- static ElementValueFinfo< CaConcBase, double > ceiling( "ceiling",
- "Ceiling value for Ca concentration. If Ca > ceiling, Ca = ceiling. If ceiling <= 0.0, there is no upper limit on Ca concentration value.",
- &CaConcBase::setCeiling,
- &CaConcBase::getCeiling
- );
- static ElementValueFinfo< CaConcBase, double > floor( "floor",
- "Floor value for Ca concentration. If Ca < floor, Ca = floor",
- &CaConcBase::setFloor,
- &CaConcBase::getFloor
- );
-
-///////////////////////////////////////////////////////
-// MsgDest definitions
-///////////////////////////////////////////////////////
-
- static DestFinfo current( "current",
- "Calcium Ion current, due to be converted to conc.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::current )
- );
-
- static DestFinfo currentFraction( "currentFraction",
- "Fraction of total Ion current, that is carried by Ca2+.",
- new EpFunc2< CaConcBase, double, double >( &CaConcBase::currentFraction )
- );
-
- static DestFinfo increase( "increase",
- "Any input current that increases the concentration.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::increase )
- );
-
- static DestFinfo decrease( "decrease",
- "Any input current that decreases the concentration.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::decrease )
- );
-
- static DestFinfo basal( "basal",
- "Synonym for assignment of basal conc.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::setCaBasal )
- );
-
- static Finfo* CaConcBaseFinfos[] =
- {
- &proc, // Shared
- concOut(), // Src
- &Ca, // Value
- &CaBasal, // Value
- &Ca_base, // Value
- &tau, // Value
- &B, // Value
- &thick, // Value
- &diameter, // Value
- &length, // Value
- &ceiling, // Value
- &floor, // Value
- ¤t, // Dest
- ¤tFraction, // Dest
- &increase, // Dest
- &decrease, // Dest
- &basal, // Dest
- };
-
- // We want the Ca updates before channel updates, so along with compts.
- // static SchedInfo schedInfo[] = { { process, 0, 0 } };
-
- static string doc[] =
- {
- "Name", "CaConcBase",
- "Author", "Upinder S. Bhalla, 2014, NCBS",
- "Description",
- "CaConcBase: Base class for Calcium concentration "
- "pool. Takes current from a channel and keeps track of "
- "calcium buildup and depletion by a "
- "single exponential process. ",
- };
+ ///////////////////////////////////////////////////////
+ // Field definitions
+ ///////////////////////////////////////////////////////
+ static ElementValueFinfo< CaConcBase, double > Ca( "Ca",
+ "Calcium concentration.",
+ &CaConcBase::setCa,
+ &CaConcBase::getCa
+ );
+ static ElementValueFinfo< CaConcBase, double > CaBasal( "CaBasal",
+ "Basal Calcium concentration.",
+ &CaConcBase::setCaBasal,
+ &CaConcBase::getCaBasal
+ );
+ static ElementValueFinfo< CaConcBase, double > Ca_base( "Ca_base",
+ "Basal Calcium concentration, synonym for CaBasal",
+ &CaConcBase::setCaBasal,
+ &CaConcBase::getCaBasal
+ );
+ static ElementValueFinfo< CaConcBase, double > tau( "tau",
+ "Settling time for Ca concentration",
+ &CaConcBase::setTau,
+ &CaConcBase::getTau
+ );
+ static ElementValueFinfo< CaConcBase, double > B( "B",
+ "Volume scaling factor. "
+ "Deprecated. This is a legacy field from GENESIS and exposes "
+ "internal calculations. Please do not use. \n"
+ "B = 1/(vol * F) so that it obeys:\n"
+ "dC/dt = B*I_Ca - C/tau",
+ &CaConcBase::setB,
+ &CaConcBase::getB
+ );
+ static ElementValueFinfo< CaConcBase, double > thick( "thick",
+ "Thickness of Ca shell, assumed cylindrical. Legal range is "
+ "between zero and the radius. If outside this range it is "
+ "taken as the radius. Default zero, ie, the shell is the entire "
+ "thickness of the cylinder",
+ &CaConcBase::setThickness,
+ &CaConcBase::getThickness
+ );
+ static ElementValueFinfo< CaConcBase, double > length( "length",
+ "Length of Ca shell, assumed cylindrical",
+ &CaConcBase::setLength,
+ &CaConcBase::getLength
+ );
+ static ElementValueFinfo< CaConcBase, double > diameter( "diameter",
+ "Diameter of Ca shell, assumed cylindrical",
+ &CaConcBase::setDiameter,
+ &CaConcBase::getDiameter
+ );
+ static ElementValueFinfo< CaConcBase, double > ceiling( "ceiling",
+ "Ceiling value for Ca concentration. If Ca > ceiling, Ca = ceiling. If ceiling <= 0.0, there is no upper limit on Ca concentration value.",
+ &CaConcBase::setCeiling,
+ &CaConcBase::getCeiling
+ );
+ static ElementValueFinfo< CaConcBase, double > floor( "floor",
+ "Floor value for Ca concentration. If Ca < floor, Ca = floor",
+ &CaConcBase::setFloor,
+ &CaConcBase::getFloor
+ );
+
+ ///////////////////////////////////////////////////////
+ // MsgDest definitions
+ ///////////////////////////////////////////////////////
+
+ static DestFinfo current( "current",
+ "Calcium Ion current, due to be converted to conc.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::current )
+ );
+
+ static DestFinfo currentFraction( "currentFraction",
+ "Fraction of total Ion current, that is carried by Ca2+.",
+ new EpFunc2< CaConcBase, double, double >( &CaConcBase::currentFraction )
+ );
+
+ static DestFinfo increase( "increase",
+ "Any input current that increases the concentration.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::increase )
+ );
+
+ static DestFinfo decrease( "decrease",
+ "Any input current that decreases the concentration.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::decrease )
+ );
+
+ static DestFinfo basal( "basal",
+ "Synonym for assignment of basal conc.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::setCaBasal )
+ );
+
+ static Finfo* CaConcBaseFinfos[] =
+ {
+ &proc, // Shared
+ concOut(), // Src
+ &Ca, // Value
+ &CaBasal, // Value
+ &Ca_base, // Value
+ &tau, // Value
+ &B, // Value
+ &thick, // Value
+ &diameter, // Value
+ &length, // Value
+ &ceiling, // Value
+ &floor, // Value
+ ¤t, // Dest
+ ¤tFraction, // Dest
+ &increase, // Dest
+ &decrease, // Dest
+ &basal, // Dest
+ };
+
+ // We want the Ca updates before channel updates, so along with compts.
+ // static SchedInfo schedInfo[] = { { process, 0, 0 } };
+
+ static string doc[] =
+ {
+ "Name", "CaConcBase",
+ "Author", "Upinder S. Bhalla, 2014, NCBS",
+ "Description",
+ "CaConcBase: Base class for Calcium concentration "
+ "pool. Takes current from a channel and keeps track of "
+ "calcium buildup and depletion by a "
+ "single exponential process. ",
+ };
- static ZeroSizeDinfo< int > dinfo;
-
- static Cinfo CaConcBaseCinfo(
- "CaConcBase",
- Neutral::initCinfo(),
- CaConcBaseFinfos,
- sizeof( CaConcBaseFinfos )/sizeof(Finfo *),
- &dinfo,
- doc,
- sizeof(doc)/sizeof(string)
- );
-
- return &CaConcBaseCinfo;
+ static ZeroSizeDinfo< int > dinfo;
+
+ static Cinfo CaConcBaseCinfo(
+ "CaConcBase",
+ Neutral::initCinfo(),
+ CaConcBaseFinfos,
+ sizeof( CaConcBaseFinfos )/sizeof(Finfo *),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string)
+ );
+
+ return &CaConcBaseCinfo;
}
///////////////////////////////////////////////////
static const Cinfo* caConcCinfo = CaConcBase::initCinfo();
CaConcBase::CaConcBase()
- :
- thickness_( 0.0 )
+ :
+ thickness_( 0.0 )
{;}
///////////////////////////////////////////////////
@@ -203,98 +203,98 @@ CaConcBase::CaConcBase()
void CaConcBase::setCa( const Eref& e, double Ca )
{
- vSetCa( e, Ca );
+ vSetCa( e, Ca );
}
double CaConcBase::getCa( const Eref& e ) const
{
- return vGetCa( e );
+ return vGetCa( e );
}
void CaConcBase::setCaBasal( const Eref& e, double CaBasal )
{
- vSetCaBasal( e, CaBasal );
+ vSetCaBasal( e, CaBasal );
}
double CaConcBase::getCaBasal( const Eref& e ) const
{
- return vGetCaBasal( e );
+ return vGetCaBasal( e );
}
void CaConcBase::setTau( const Eref& e, double tau )
{
- vSetTau( e, tau );
+ vSetTau( e, tau );
}
double CaConcBase::getTau( const Eref& e ) const
{
- return vGetTau( e );
+ return vGetTau( e );
}
void CaConcBase::setB( const Eref& e, double B )
{
- vSetB( e, B );
+ vSetB( e, B );
}
double CaConcBase::getB( const Eref& e ) const
{
- return vGetB( e );
+ return vGetB( e );
}
void CaConcBase::setCeiling( const Eref& e, double ceiling )
{
- vSetCeiling( e, ceiling );
+ vSetCeiling( e, ceiling );
}
double CaConcBase::getCeiling( const Eref& e ) const
{
- return vGetCeiling( e );
+ return vGetCeiling( e );
}
void CaConcBase::setFloor( const Eref& e, double floor )
{
- vSetFloor( e, floor );
+ vSetFloor( e, floor );
}
double CaConcBase::getFloor( const Eref& e ) const
{
- return vGetFloor( e );
+ return vGetFloor( e );
}
void CaConcBase::updateDimensions( const Eref& e )
{
- double vol = PI * diameter_ * diameter_ * length_ * 0.25;
- if ( thickness_ > 0 && thickness_ < diameter_/2.0 ) {
- double coreRadius = diameter_ / 2.0 - thickness_;
- vol -= PI * coreRadius * coreRadius * length_;
- }
- double B = 1.0 / ( FaradayConst * vol );
- vSetB( e, B );
+ double vol = PI * diameter_ * diameter_ * length_ * 0.25;
+ if ( thickness_ > 0 && thickness_ < diameter_/2.0 ) {
+ double coreRadius = diameter_ / 2.0 - thickness_;
+ vol -= PI * coreRadius * coreRadius * length_;
+ }
+ double B = 1.0 / ( FaradayConst * vol );
+ vSetB( e, B );
}
void CaConcBase::setThickness( const Eref& e, double thickness )
{
- thickness_ = thickness;
- updateDimensions( e );
+ thickness_ = thickness;
+ updateDimensions( e );
}
double CaConcBase::getThickness( const Eref& e ) const
{
- return thickness_;
+ return thickness_;
}
void CaConcBase::setDiameter( const Eref& e, double diameter )
{
- diameter_ = diameter;
- updateDimensions( e );
+ diameter_ = diameter;
+ updateDimensions( e );
}
double CaConcBase::getDiameter( const Eref& e ) const
{
- return diameter_;
+ return diameter_;
}
void CaConcBase::setLength( const Eref& e, double length )
{
- length_ = length;
- updateDimensions( e );
+ length_ = length;
+ updateDimensions( e );
}
double CaConcBase::getLength( const Eref& e ) const
{
- return length_;
+ return length_;
}
///////////////////////////////////////////////////
@@ -303,32 +303,32 @@ double CaConcBase::getLength( const Eref& e ) const
void CaConcBase::reinit( const Eref& e, ProcPtr p )
{
- vReinit( e, p );
+ vReinit( e, p );
}
void CaConcBase::process( const Eref& e, ProcPtr p )
{
- vProcess( e, p );
+ vProcess( e, p );
}
void CaConcBase::current( const Eref& e, double I )
{
- vCurrent( e, I );
+ vCurrent( e, I );
}
void CaConcBase::currentFraction( const Eref& e, double I, double fraction )
{
- vCurrentFraction( e, I, fraction );
+ vCurrentFraction( e, I, fraction );
}
void CaConcBase::increase( const Eref& e, double I )
{
- vIncrease( e, I );
+ vIncrease( e, I );
}
void CaConcBase::decrease( const Eref& e, double I )
{
- vDecrease( e, I );
+ vDecrease( e, I );
}
/////////////////////////////////////////////////////
@@ -342,48 +342,48 @@ void CaConcBase::vSetSolver( const Eref& e, Id hsolve )
// static func
void CaConcBase::zombify( Element* orig, const Cinfo* zClass,
- Id hsolve )
+ Id hsolve )
{
- if ( orig->cinfo() == zClass )
- return;
- unsigned int start = orig->localDataStart();
- unsigned int num = orig->numLocalData();
- if ( num == 0 )
- return;
- vector< double > data( num * 9 );
-
- unsigned int j = 0;
- for ( unsigned int i = 0; i < num; ++i ) {
- Eref er( orig, i + start );
- const CaConcBase* cb =
- reinterpret_cast< const CaConcBase* >( er.data() );
- data[j + 0] = cb->getCa( er );
- data[j + 1] = cb->getCaBasal( er );
- data[j + 2] = cb->getTau( er );
- data[j + 3] = cb->getB( er );
- data[j + 4] = cb->getCeiling( er );
- data[j + 5] = cb->getFloor( er );
- data[j + 6] = cb->getThickness( er );
- data[j + 7] = cb->getLength( er );
- data[j + 8] = cb->getDiameter( er );
- j += 9;
- }
- orig->zombieSwap( zClass );
- j = 0;
- for ( unsigned int i = 0; i < num; ++i ) {
- Eref er( orig, i + start );
- CaConcBase* cb = reinterpret_cast< CaConcBase* >( er.data() );
- cb->vSetSolver( er, hsolve );
- cb->setCa( er, data[j + 0] );
- cb->setCaBasal( er, data[j + 1] );
- cb->setTau( er, data[j + 2] );
- cb->setB( er, data[j + 3] );
- cb->setCeiling( er, data[j + 4] );
- cb->setFloor( er, data[j + 5] );
- cb->setThickness( er, data[j + 6] );
- cb->setLength( er, data[j + 7] );
- cb->setDiameter( er, data[j + 8] );
- j += 7;
- }
+ if ( orig->cinfo() == zClass )
+ return;
+ unsigned int start = orig->localDataStart();
+ unsigned int num = orig->numLocalData();
+ if ( num == 0 )
+ return;
+ vector< double > data( num * 9 );
+
+ unsigned int j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ const CaConcBase* cb =
+ reinterpret_cast< const CaConcBase* >( er.data() );
+ data[j + 0] = cb->getCa( er );
+ data[j + 1] = cb->getCaBasal( er );
+ data[j + 2] = cb->getTau( er );
+ data[j + 3] = cb->getB( er );
+ data[j + 4] = cb->getCeiling( er );
+ data[j + 5] = cb->getFloor( er );
+ data[j + 6] = cb->getThickness( er );
+ data[j + 7] = cb->getLength( er );
+ data[j + 8] = cb->getDiameter( er );
+ j += 9;
+ }
+ orig->zombieSwap( zClass );
+ j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ CaConcBase* cb = reinterpret_cast< CaConcBase* >( er.data() );
+ cb->vSetSolver( er, hsolve );
+ cb->setCa( er, data[j + 0] );
+ cb->setCaBasal( er, data[j + 1] );
+ cb->setTau( er, data[j + 2] );
+ cb->setB( er, data[j + 3] );
+ cb->setCeiling( er, data[j + 4] );
+ cb->setFloor( er, data[j + 5] );
+ cb->setThickness( er, data[j + 6] );
+ cb->setLength( er, data[j + 7] );
+ cb->setDiameter( er, data[j + 8] );
+ j += 9; //was 7?
+ }
}
diff --git a/biophysics/DifBuffer.cpp b/biophysics/DifBuffer.cpp
new file mode 100644
index 00000000..60763e82
--- /dev/null
+++ b/biophysics/DifBuffer.cpp
@@ -0,0 +1,467 @@
+// DifBuffer.cpp ---
+//
+// Filename: DifBuffer.cpp
+// Description:
+// Author: Subhasis Ray
+// Maintainer:
+// Created: Mon Feb 16 12:02:11 2015 (-0500)
+// Version:
+// Package-Requires: ()
+// Last-Updated: Mon Feb 23 13:07:56 2015 (-0500)
+// By: Subhasis Ray
+// Update #: 130
+// URL:
+// Doc URL:
+// Keywords:
+// Compatibility:
+//
+//
+
+// Commentary:
+//
+//
+//
+//
+
+// Change Log:
+// 5/25/16 Completing DifBuffer -- Asia J-Szmek (GMU)
+// 9/21/16 rewrote DifBuffer to account for DifBufferBase (AJS)
+//
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or (at
+// your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
+//
+//
+
+// Code:
+
+#include "header.h"
+#include "DifBufferBase.h"
+#include "DifBuffer.h"
+#include "ElementValueFinfo.h"
+#include "../utility/numutil.h"
+#include <cmath>
+const double DifBuffer::EPSILON = 1.0e-10;
+
+
+const Cinfo * DifBuffer::initCinfo()
+{
+
+ static string doc[] = {
+ "Name", "DifBuffer",
+ "Author", "Subhasis Ray (ported from GENESIS2)",
+ "Description", "Models diffusible buffer where total concentration is constant. It is"
+ " coupled with a DifShell.",
+ };
+ static Dinfo<DifBuffer> dinfo;
+ static Cinfo difBufferCinfo(
+ "DifBuffer",
+ DifBufferBase::initCinfo(),
+ 0,
+ 0,
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string));
+
+ return &difBufferCinfo;
+}
+
+static const Cinfo * difBufferCinfo = DifBuffer::initCinfo();
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Class functions
+////////////////////////////////////////////////////////////////////////////////
+
+DifBuffer::DifBuffer() :
+ activation_(0),
+ Af_(0),
+ Bf_(0),
+ bFree_(0),
+ bBound_(0),
+ prevFree_(0),
+ prevBound_(0),
+ bTot_(0),
+ kf_(0),
+ kb_(0),
+ D_(0),
+ shapeMode_(0),
+ length_(0),
+ diameter_(0),
+ thickness_(0),
+ volume_(0),
+ outerArea_(0),
+ innerArea_(0)
+
+{}
+
+////////////////////////////////////////////////////////////////////////////////
+// Field access functions
+////////////////////////////////////////////////////////////////////////////////
+double DifBuffer::vGetActivation(const Eref& e) const
+{
+ return activation_;
+}
+
+void DifBuffer::vSetActivation(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Activation cannot be negative!\n";
+ return;
+ }
+ activation_ = value;
+}
+
+
+double DifBuffer::vGetBFree(const Eref& e) const
+{
+ return bFree_;
+}
+void DifBuffer::vSetBFree(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Free Buffer cannot be negative!\n";
+ return;
+ }
+ if (value > bTot_){
+ cerr << "Error: DifBuffer: Free Buffer cannot exceed total buffer!\n";
+ return;
+ }
+ bFree_ = value;
+ bBound_ = bTot_ - bFree_;
+ prevFree_= bFree_;
+ prevBound_ = bBound_;
+
+}
+
+double DifBuffer::vGetBBound(const Eref& e) const
+{
+ return bBound_;
+}
+
+void DifBuffer::vSetBBound(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Bound Buffer cannot be negative!\n";
+ return;
+ }
+ if (value > bTot_){
+ cerr << "Error: DifBuffer: Bound Buffer cannot exceed total buffer!\n";
+ return;
+ }
+ bBound_ = value;
+ bFree_ = bTot_ - bBound_;
+ prevFree_= bFree_;
+ prevBound_ = bBound_;
+}
+
+
+double DifBuffer::vGetBTot(const Eref& e) const
+{
+ return bTot_;
+}
+
+void DifBuffer::vSetBTot(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Total buffer concentration cannot be negative!\n";
+ return;
+ }
+ bTot_ = value;
+ bFree_ = bTot_;
+ bBound_ = 0;
+}
+
+
+double DifBuffer::vGetKf(const Eref& e) const
+{
+ return kf_;
+}
+
+void DifBuffer::vSetKf(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Kf cannot be negative!\n";
+ return;
+ }
+ kf_ = value;
+}
+
+
+double DifBuffer::vGetKb(const Eref& e) const
+{
+ return kb_;
+}
+
+void DifBuffer::vSetKb(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Kb cannot be negative!\n";
+ return;
+ }
+ kb_ = value;
+}
+
+double DifBuffer::vGetD(const Eref& e) const
+{
+ return D_;
+}
+
+void DifBuffer::vSetD(const Eref& e,double value)
+{
+
+ if ( value < 0.0 ) {
+ cerr << " Error: DifBuffer: Diffusion constant, D, cannot be negative!\n";
+ return;
+ }
+ D_ = value;
+}
+
+
+void DifBuffer::vSetShapeMode(const Eref& e, unsigned int shapeMode )
+{
+ if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) {
+ cerr << "Error: DifBuffer: I only understand shapeModes 0, 1 and 3.\n";
+ return;
+ }
+ shapeMode_ = shapeMode;
+}
+
+unsigned int DifBuffer::vGetShapeMode(const Eref& e) const
+{
+ return shapeMode_;
+}
+
+void DifBuffer::vSetLength(const Eref& e, double length )
+{
+ if ( length < 0.0 ) {
+ cerr << "Error: DifBuffer: length cannot be negative!\n";
+ return;
+ }
+
+ length_ = length;
+}
+
+double DifBuffer::vGetLength(const Eref& e ) const
+{
+ return length_;
+}
+
+void DifBuffer::vSetDiameter(const Eref& e, double diameter )
+{
+ if ( diameter < 0.0 ) {
+ cerr << "Error: DifBuffer: diameter cannot be negative!\n";
+ return;
+ }
+
+ diameter_ = diameter;
+}
+
+double DifBuffer::vGetDiameter(const Eref& e ) const
+{
+ return diameter_;
+}
+
+void DifBuffer::vSetThickness( const Eref& e, double thickness )
+{
+ if ( thickness < 0.0 ) {
+ cerr << "Error: DifBuffer: thickness cannot be negative!\n";
+ return;
+ }
+
+ thickness_ = thickness;
+}
+
+double DifBuffer::vGetThickness(const Eref& e) const
+{
+ return thickness_;
+}
+
+void DifBuffer::vSetVolume(const Eref& e, double volume )
+{
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifBuffer: Trying to set volume, when shapeMode is not USER-DEFINED\n";
+
+ if ( volume < 0.0 ) {
+ cerr << "Error: DifBuffer: volume cannot be negative!\n";
+ return;
+ }
+
+ volume_ = volume;
+}
+
+double DifBuffer::vGetVolume(const Eref& e ) const
+{
+ return volume_;
+}
+
+void DifBuffer::vSetOuterArea(const Eref& e, double outerArea )
+{
+ if (shapeMode_ != 3 )
+ cerr << "Warning: DifBuffer: Trying to set outerArea, when shapeMode is not USER-DEFINED\n";
+
+ if ( outerArea < 0.0 ) {
+ cerr << "Error: DifBuffer: outerArea cannot be negative!\n";
+ return;
+ }
+
+ outerArea_ = outerArea;
+}
+
+double DifBuffer::vGetOuterArea(const Eref& e ) const
+{
+ return outerArea_;
+}
+
+void DifBuffer::vSetInnerArea(const Eref& e, double innerArea )
+{
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifBuffer: Trying to set innerArea, when shapeMode is not USER-DEFINED\n";
+
+ if ( innerArea < 0.0 ) {
+ cerr << "Error: DifBuffer: innerArea cannot be negative!\n";
+ return;
+ }
+
+ innerArea_ = innerArea;
+}
+
+double DifBuffer::vGetInnerArea(const Eref& e) const
+{
+ return innerArea_;
+}
+
+
+
+
+void DifBuffer::vBuffer(const Eref& e,
+ double C )
+{
+ activation_ = C;
+ //cout<<"Buffer"<< C<<" ";
+}
+
+
+double DifBuffer::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 DifBuffer::vProcess( const Eref & e, ProcPtr p )
+{
+ /**
+ * Send ion concentration and thickness to adjacent DifShells. They will
+ * then compute their incoming fluxes.
+ */
+
+
+
+ Af_ += kb_ * bBound_;
+ Bf_ += kf_ * activation_;
+
+ bFree_ = integrate(bFree_,p->dt,Af_,Bf_);
+ bBound_ = bTot_ - bFree_;
+ prevFree_ = bFree_;
+ prevBound_ = bBound_;
+
+ /**
+ * Send ion concentration to ion buffers. They will send back information on
+ * the reaction (forward / backward rates ; free / bound buffer concentration)
+ * immediately, which this DifShell will use to find amount of ion captured
+ * or released in the current time-step.
+ */
+ innerDifSourceOut()->send( e, prevFree_, thickness_ );
+ outerDifSourceOut()->send( e, prevFree_, thickness_ );
+ reactionOut()->send(e,kf_,kb_,bFree_,bBound_);
+ Af_ = 0;
+ Bf_= 0;
+
+}
+
+void DifBuffer::vReinit( const Eref& e, ProcPtr p )
+{
+
+ Af_ = 0;
+ Bf_= 0;
+ double rOut = diameter_/2.;
+
+ double rIn = rOut - thickness_;
+
+ if (rIn<0)
+ rIn = 0.;
+ switch ( shapeMode_ )
+ {
+ /*
+ * Onion Shell
+ */
+ case 0:
+ if ( length_ == 0.0 ) { // Spherical shell
+ volume_ = 4./3.* M_PI * ( rOut * rOut * rOut - rIn * rIn * rIn );
+ outerArea_ = 4*M_PI * rOut * rOut;
+ innerArea_ = 4*M_PI * rIn * rIn;
+ } else { // Cylindrical shell
+ volume_ = ( M_PI * length_ ) * ( rOut * rOut - rIn * rIn );
+ outerArea_ = 2*M_PI * rOut * length_;
+ innerArea_ = 2*M_PI * rIn * length_;
+ }
+
+ break;
+
+ /*
+ * Cylindrical Slice
+ */
+ case 1:
+ volume_ = M_PI * diameter_ * diameter_ * thickness_ / 4.0;
+ outerArea_ = M_PI * diameter_ * diameter_ / 4.0;
+ innerArea_ = outerArea_;
+ break;
+
+ /*
+ * User defined
+ */
+ case 3:
+ // Nothing to be done here. Volume and inner-, outer areas specified by
+ // user.
+ break;
+
+ default:
+ assert( 0 );
+ }
+
+ bFree_ = bTot_/(1+activation_*kf_/kb_);
+ prevFree_ = bFree_;
+ bBound_ = bTot_ - bFree_;
+ prevBound_ = bBound_;
+ innerDifSourceOut()->send( e, prevFree_, thickness_ );
+ outerDifSourceOut()->send( e, prevFree_, thickness_ );
+
+}
+
+void DifBuffer::vFluxFromIn(const Eref& e,double innerC, double innerThickness)
+{
+ double dif = 2 * D_ / volume_* innerArea_ / (thickness_ + innerThickness);
+ Af_ += dif * innerC;
+ Bf_ += dif;
+}
+
+void DifBuffer::vFluxFromOut(const Eref& e,double outerC, double outerThickness)
+{
+ double dif = 2 * D_ / volume_* outerArea_ / (thickness_ + outerThickness);
+ Af_ += dif * outerC;
+ Bf_ += dif;
+}
diff --git a/biophysics/DifBuffer.h b/biophysics/DifBuffer.h
new file mode 100644
index 00000000..60f57445
--- /dev/null
+++ b/biophysics/DifBuffer.h
@@ -0,0 +1,96 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani 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 _DifBuffer_h
+#define _DifBuffer_h
+
+class DifBuffer: public DifBufferBase{
+ public:
+ DifBuffer();
+
+ void vBuffer(const Eref& e,double C);
+ void vReinit( const Eref & e, ProcPtr p );
+ void vProcess(const Eref & e, ProcPtr p );
+ void vFluxFromOut(const Eref& e,double outerC, double outerThickness );
+ void vFluxFromIn( const Eref& e,double innerC, double innerThickness );
+ //Field access functions
+
+ double vGetActivation(const Eref& e) const;
+ void vSetActivation(const Eref& e,double value);
+
+ double vGetBFree(const Eref& e) const;
+ void vSetBFree(const Eref& e,double value);
+
+ double vGetBBound(const Eref& e) const;
+ void vSetBBound(const Eref& e,double value);
+
+ double vGetBTot(const Eref& e) const; // total buffer concentration in mM (free + bound)
+ void vSetBTot(const Eref& e,double value);
+
+ double vGetKf(const Eref& e) const; // forward rate constant in 1/(mM*sec)
+ void vSetKf(const Eref& e,double value);
+
+ double vGetKb(const Eref& e) const; // backward rate constant in 1/sec
+ void vSetKb(const Eref& e,double value);
+
+ double vGetD(const Eref& e) const; // diffusion constant of buffer molecules, m^2/sec
+ void vSetD(const Eref& e,double value);
+ void vSetShapeMode(const Eref& e, unsigned int shapeMode );
+ unsigned int vGetShapeMode(const Eref& e) const;
+
+ void vSetLength(const Eref& e, double length );
+ double vGetLength(const Eref& e) const;
+
+ void vSetDiameter(const Eref& e, double diameter );
+ double vGetDiameter(const Eref& e) const;
+
+ void vSetThickness(const Eref& e, double thickness );
+ double vGetThickness(const Eref& e) const;
+
+ void vSetVolume(const Eref& e, double volume );
+ double vGetVolume(const Eref& e) const;
+
+ void vSetOuterArea(const Eref& e, double outerArea );
+ double vGetOuterArea(const Eref& e) const;
+
+ void vSetInnerArea(const Eref& e, double innerArea );
+ double vGetInnerArea(const Eref& e) const;
+
+ static const Cinfo * initCinfo();
+
+ private:
+
+
+ double integrate(double state, double dt, double A, double B );
+
+ double activation_; //ion concentration from incoming CONCEN message
+ double Af_;
+ double Bf_;
+ double bFree_; //free buffer concentration
+ double bBound_; //bound buffer concentration
+ double prevFree_;
+ double prevBound_;
+ //double prevFree_; //bFree at previous time
+ //double prevBound_; //bBound at previous time
+ double bTot_; //bFree+bBound
+ double kf_; //forward rate constant
+ double kb_; //backward rate constant
+ double D_; //diffusion constant
+ unsigned int shapeMode_;
+ double length_;
+ double diameter_;
+ double thickness_;
+ double volume_;
+ double outerArea_;
+ double innerArea_;
+ static const double EPSILON;
+};
+
+#endif // _DifBuffer_h
+/* DifBuffer.h ends here */
diff --git a/biophysics/DifBufferBase.cpp b/biophysics/DifBufferBase.cpp
new file mode 100644
index 00000000..40628a21
--- /dev/null
+++ b/biophysics/DifBufferBase.cpp
@@ -0,0 +1,461 @@
+#include "header.h"
+#include "DifBufferBase.h"
+#include "ElementValueFinfo.h"
+#include "../utility/numutil.h"
+#include <cmath>
+
+SrcFinfo4< double, double, double, double >* DifBufferBase::reactionOut()
+{
+ static SrcFinfo4< double, double, double, double > reactionOut(
+ "reactionOut",
+ "Sends out reaction rates (forward and backward), and concentrations"
+ " (free-buffer and bound-buffer molecules).");
+ return &reactionOut;
+}
+
+
+SrcFinfo2< double, double >* DifBufferBase::innerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("innerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+SrcFinfo2< double, double >* DifBufferBase::outerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("outerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+const Cinfo * DifBufferBase::initCinfo()
+{
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< DifBufferBase >( &DifBufferBase::process) );
+ static DestFinfo reinit( "reinit",
+ "Reinit happens only in stage 0",
+ new ProcOpFunc< DifBufferBase >( &DifBufferBase::reinit));
+
+ static Finfo* processShared[] = {
+ &process,
+ &reinit
+ };
+
+ static SharedFinfo proc(
+ "proc",
+ "Here we create 2 shared finfos to attach with the Ticks. This is because we want to perform DifBufferBase "
+ "computations in 2 stages, much as in the Compartment object. "
+ "In the first stage we send out the concentration value to other DifBufferBases and Buffer elements. We also",
+ processShared,
+ sizeof( processShared ) / sizeof( Finfo* ));
+
+ static DestFinfo concentration("concentration",
+ "Receives concentration (from DifShell).",
+ new EpFunc1<DifBufferBase, double>(&DifBufferBase::buffer));
+ static Finfo* bufferShared[] = {
+ &concentration, DifBufferBase::reactionOut()
+ };
+ static SharedFinfo buffer( "buffer",
+ "This is a shared message with DifShell. "
+ "During stage 0:\n "
+ " - DifBufferBase sends ion concentration\n"
+ " - Buffer updates buffer concentration and sends it back immediately using a call-back.\n"
+ " - DifShell updates the time-derivative ( dC / dt ) \n"
+ "\n"
+ "During stage 1: \n"
+ " - DifShell advances concentration C \n\n"
+ "This scheme means that the Buffer does not need to be scheduled, and it does its computations when "
+ "it receives a cue from the DifShell. May not be the best idea, but it saves us from doing the above "
+ "computations in 3 stages instead of 2." ,
+ bufferShared,
+ sizeof( bufferShared ) / sizeof( Finfo* ));
+
+ static DestFinfo fluxFromOut( "fluxFromOut",
+ "Destination message",
+ new EpFunc2< DifBufferBase, double, double > ( &DifBufferBase::fluxFromOut ));
+
+ static Finfo* innerDifShared[] = {
+ &fluxFromOut,
+ DifBufferBase::innerDifSourceOut()
+
+ };
+
+ static SharedFinfo innerDif( "innerDif",
+ "This shared message (and the next) is between DifBufferBases: adjoining shells exchange information to "
+ "find out the flux between them. "
+ "Using this message, an inner shell sends to, and receives from its outer shell." ,
+ innerDifShared,
+ sizeof( innerDifShared ) / sizeof( Finfo* ));
+
+ static DestFinfo fluxFromIn( "fluxFromIn", "",
+ new EpFunc2< DifBufferBase, double, double> ( &DifBufferBase::fluxFromIn) );
+
+ static Finfo* outerDifShared[] = {
+ &fluxFromIn,
+ DifBufferBase::outerDifSourceOut(),
+
+ };
+
+ static SharedFinfo outerDif( "outerDif",
+ "Using this message, an outer shell sends to, and receives from its inner shell." ,
+ outerDifShared,
+ sizeof( outerDifShared ) / sizeof( Finfo* ));
+
+ ////////////////////////////
+ // Field defs
+ ////////////////////////////
+ static ElementValueFinfo<DifBufferBase, double> activation("activation",
+ "Ion concentration from incoming conc message.",
+ &DifBufferBase::setActivation,
+ &DifBufferBase::getActivation);
+ static ElementValueFinfo<DifBufferBase, double> kf("kf",
+ "Forward rate constant of buffer molecules 1/mM/s (?)",
+ &DifBufferBase::setKf,
+ &DifBufferBase::getKf);
+ static ElementValueFinfo<DifBufferBase, double> kb("kb",
+ "Backward rate constant of buffer molecules. 1/s",
+ &DifBufferBase::setKb,
+ &DifBufferBase::getKb);
+ static ElementValueFinfo<DifBufferBase, double> D("D",
+ "Diffusion constant of buffer molecules. m^2/s",
+ &DifBufferBase::setD,
+ &DifBufferBase::getD);
+ static ElementValueFinfo<DifBufferBase, double> bFree("bFree",
+ "Free buffer concentration",
+ &DifBufferBase::setBFree,
+ &DifBufferBase::getBFree);
+ static ElementValueFinfo<DifBufferBase, double> bBound("bBound",
+ "Bound buffer concentration",
+ &DifBufferBase::setBBound,
+ &DifBufferBase::getBBound);
+ static ElementValueFinfo<DifBufferBase, double> bTot("bTot",
+ "Total buffer concentration.",
+ &DifBufferBase::setBTot,
+ &DifBufferBase::getBTot);
+ static ElementValueFinfo<DifBufferBase, double> length("length",
+ "Length of shell",
+ &DifBufferBase::setLength,
+ &DifBufferBase::getLength);
+ static ElementValueFinfo<DifBufferBase, double> diameter("diameter",
+ "Diameter of shell",
+ &DifBufferBase::setDiameter,
+ &DifBufferBase::getDiameter);
+ static ElementValueFinfo<DifBufferBase, unsigned int> shapeMode("shapeMode",
+ "shape of the shell: SHELL=0, SLICE=SLAB=1, USERDEF=3",
+ &DifBufferBase::setShapeMode,
+ &DifBufferBase::getShapeMode);
+
+ static ElementValueFinfo<DifBufferBase, double> thickness("thickness",
+ "Thickness of shell",
+ &DifBufferBase::setThickness,
+ &DifBufferBase::getThickness);
+
+ static ElementValueFinfo<DifBufferBase, double> innerArea("innerArea",
+ "Inner area of shell",
+ &DifBufferBase::setInnerArea,
+ &DifBufferBase::getInnerArea);
+ static ElementValueFinfo<DifBufferBase, double> outerArea("outerArea",
+ "Outer area of shell",
+ &DifBufferBase::setOuterArea,
+ &DifBufferBase::getOuterArea);
+ static ElementValueFinfo< DifBufferBase, double> volume( "volume", "",
+ &DifBufferBase::setVolume,
+ &DifBufferBase::getVolume );
+
+ ////
+ // DestFinfo
+ ////
+ static Finfo * difBufferFinfos[] = {
+ //////////////////////////////////////////////////////////////////
+ // Field definitions
+ //////////////////////////////////////////////////////////////////
+
+ &activation,
+ &D,
+ &bFree,
+ &bBound,
+ &bTot,
+ &kf,
+ &kb,
+ //&prevFree,
+ //&prevBound,
+ &length,
+ &diameter,
+ &shapeMode,
+ &thickness,
+ &innerArea,
+ &outerArea,
+ &volume,
+ //////////////////////////////////////////////////////////////////
+ // SharedFinfo definitions
+ /////////////////////////////////////////////////////////////////
+ &proc,
+ &buffer,
+ &innerDif,
+ &outerDif,
+ //
+ reactionOut(),
+ innerDifSourceOut(),
+ outerDifSourceOut(),
+ //////////////////////////////////////////////////////////////////
+ // DestFinfo definitions
+ //////////////////////////////////////////////////////////////////
+ &concentration,
+ };
+
+ static string doc[] = {
+ "Name", "DifBufferBase",
+ "Author", "Subhasis Ray (ported from GENESIS2)",
+ "Description", "Models diffusible buffer where total concentration is constant. It is"
+ " coupled with a DifShell.",
+ };
+ static ZeroSizeDinfo<int> dinfo;
+ static Cinfo difBufferCinfo(
+ "DifBufferBase",
+ Neutral::initCinfo(),
+ difBufferFinfos,
+ sizeof(difBufferFinfos)/sizeof(Finfo*),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string));
+
+ return &difBufferCinfo;
+}
+
+static const Cinfo * difBufferCinfo = DifBufferBase::initCinfo();
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Class functions
+////////////////////////////////////////////////////////////////////////////////
+
+DifBufferBase::DifBufferBase()
+{ ; }
+
+
+double DifBufferBase::getActivation(const Eref& e) const
+{
+ return vGetActivation(e);
+}
+
+void DifBufferBase::setActivation(const Eref& e,double value)
+{
+ vSetActivation(e,value);
+}
+
+
+double DifBufferBase::getBFree(const Eref& e) const
+{
+ return vGetBFree(e);
+}
+
+void DifBufferBase::setBFree(const Eref& e,double value)
+{
+ vSetBFree(e,value);
+}
+
+double DifBufferBase::getBBound(const Eref& e) const
+{
+ return vGetBBound(e);
+}
+void DifBufferBase::setBBound(const Eref& e,double value)
+{
+ vSetBBound(e,value);
+}
+
+double DifBufferBase::getBTot(const Eref& e) const
+{
+ return vGetBTot(e);
+}
+
+void DifBufferBase::setBTot(const Eref& e,double value)
+{
+ vSetBTot(e,value);
+}
+
+
+double DifBufferBase::getKf(const Eref& e) const
+{
+ return vGetKf(e);
+}
+
+void DifBufferBase::setKf(const Eref& e,double value)
+{
+ vSetKf(e,value);
+}
+
+double DifBufferBase::getKb(const Eref& e) const
+{
+ return vGetKb(e);
+}
+
+void DifBufferBase::setKb(const Eref& e,double value)
+{
+ vSetKb(e,value);
+}
+
+double DifBufferBase::getD(const Eref& e) const
+{
+ return vGetD(e);
+}
+
+void DifBufferBase::setD(const Eref& e,double value)
+{
+ vSetD(e,value);
+}
+
+void DifBufferBase::setShapeMode(const Eref& e, unsigned int shapeMode )
+{
+ vSetShapeMode(e,shapeMode);
+}
+
+unsigned int DifBufferBase::getShapeMode(const Eref& e) const
+{
+ return vGetShapeMode(e);
+}
+
+void DifBufferBase::setLength(const Eref& e, double length )
+{
+ vSetLength(e,length);
+}
+
+double DifBufferBase::getLength(const Eref& e ) const
+{
+ return vGetLength(e);
+}
+
+void DifBufferBase::setDiameter(const Eref& e, double diameter )
+{
+ vSetDiameter(e,diameter);
+}
+
+double DifBufferBase::getDiameter(const Eref& e ) const
+{
+ return vGetDiameter(e);
+}
+
+void DifBufferBase::setThickness( const Eref& e, double thickness )
+{
+ vSetThickness(e,thickness);
+}
+
+double DifBufferBase::getThickness(const Eref& e) const
+{
+ return vGetThickness(e);
+}
+
+void DifBufferBase::setVolume(const Eref& e, double volume )
+{
+ vSetVolume(e,volume);
+}
+
+double DifBufferBase::getVolume(const Eref& e ) const
+{
+ return vGetVolume(e);
+}
+
+void DifBufferBase::setOuterArea(const Eref& e, double outerArea )
+{
+ vSetOuterArea(e,outerArea);
+}
+
+double DifBufferBase::getOuterArea(const Eref& e ) const
+{
+ return vGetOuterArea(e);
+}
+
+void DifBufferBase::setInnerArea(const Eref& e, double innerArea )
+{
+ vSetInnerArea(e,innerArea);
+}
+
+double DifBufferBase::getInnerArea(const Eref& e) const
+{
+ return vGetInnerArea(e);
+}
+
+
+
+void DifBufferBase::buffer(const Eref& e,double C)
+{
+ vBuffer(e,C);
+}
+
+void DifBufferBase::reinit( const Eref& e, ProcPtr p )
+{
+ vReinit( e, p );
+}
+
+void DifBufferBase::process( const Eref& e, ProcPtr p )
+{
+ vProcess( e, p );
+}
+void DifBufferBase:: fluxFromOut(const Eref& e,double outerC, double outerThickness )
+{
+ vFluxFromOut(e,outerC,outerThickness);
+}
+void DifBufferBase:: fluxFromIn(const Eref& e,double innerC, double innerThickness )
+{
+ vFluxFromIn(e,innerC,innerThickness);
+}
+void DifBufferBase::vSetSolver( const Eref& e, Id hsolve )
+{;}
+
+void DifBufferBase::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;
+ unsigned int len = 14;
+ vector< double > data( num * len );
+
+ unsigned int j = 0;
+
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ const DifBufferBase* ds =
+ reinterpret_cast< const DifBufferBase* >( er.data() );
+ data[j + 0] = ds->getActivation( er );
+ data[j + 1] = ds->getBFree( er );
+ data[j + 2] = ds->getBBound( er );
+ data[j + 3] = ds->getBTot( er );
+ data[j + 4] = ds->getKf( er );
+ data[j + 5] = ds->getKb( er );
+ data[j + 6] = ds->getD( er );
+ data[j + 7] = ds->getShapeMode( er );
+ data[j + 8] = ds->getLength( er );
+ data[j + 9] = ds->getDiameter( er );
+ data[j + 10] = ds->getThickness( er );
+ data[j + 11] = ds->getVolume( er );
+ data[j + 12] = ds->getOuterArea( er );
+ data[j + 13] = ds->getInnerArea( er );
+ j += len;
+ }
+ orig->zombieSwap( zClass );
+ j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ DifBufferBase* ds =
+ reinterpret_cast< DifBufferBase* >( er.data() );
+ ds->vSetSolver(er,hsolve);
+ ds->setActivation(er, data[j+0]);
+ ds->setBFree(er, data[j + 1]);
+ ds->setBBound(er, data[j + 2]);
+ ds->setBTot(er, data[j + 3]);
+ ds->setKf(er, data[j + 4]);
+ ds->setKb(er, data[j + 5]);
+ ds->setD(er, data[j + 6]);
+ ds->setShapeMode(er, data[j + 7]);
+ ds->setLength(er, data[j + 8]);
+ ds->setDiameter(er, data[j + 9]);
+ ds->setThickness(er, data[j + 10]);
+ ds->setVolume(er, data[j + 11]);
+ ds->setOuterArea(er, data[j + 12]);
+ ds->setInnerArea(er, data[j + 13]);
+ j += len; //??
+ }
+
+}
diff --git a/biophysics/DifBufferBase.h b/biophysics/DifBufferBase.h
new file mode 100644
index 00000000..e1bc23b6
--- /dev/null
+++ b/biophysics/DifBufferBase.h
@@ -0,0 +1,134 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani 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 _DIFBUFFER_BASE_H
+#define _DIFBUFFER_BASE_H
+/*This is base class for DifBuffer*/
+class DifBufferBase
+{
+public:
+ DifBufferBase();
+
+ void buffer(const Eref& e,double C);
+ void reinit( const Eref & e, ProcPtr p );
+ void process(const Eref & e, ProcPtr p );
+ void fluxFromOut(const Eref& e,double outerC, double outerThickness );
+ void fluxFromIn( const Eref& e,double innerC, double innerThickness );
+ virtual void vBuffer(const Eref& e,double C) = 0;
+ virtual void vReinit( const Eref & e, ProcPtr p ) = 0;
+ virtual void vProcess(const Eref & e, ProcPtr p ) = 0;
+ virtual void vFluxFromOut(const Eref& e,double outerC, double outerThickness ) = 0;
+ virtual void vFluxFromIn( const Eref& e,double innerC, double innerThickness ) = 0;
+
+
+ double getActivation(const Eref& e) const;
+ void setActivation(const Eref& e,double value);
+
+ double getBFree(const Eref& e) const;
+ void setBFree(const Eref& e,double value);
+
+ double getBBound(const Eref& e) const;
+ void setBBound(const Eref& e,double value);
+
+ double getBTot(const Eref& e) const; // total buffer concentration in mM (free + bound)
+ void setBTot(const Eref& e,double value);
+
+ double getKf(const Eref& e) const; // forward rate constant in 1/(mM*sec)
+ void setKf(const Eref& e,double value);
+
+ double getKb(const Eref& e) const; // backward rate constant in 1/sec
+ void setKb(const Eref& e,double value);
+
+ double getD(const Eref& e) const; // diffusion constant of buffer molecules, m^2/sec
+ void setD(const Eref& e,double value);
+
+
+ unsigned int getShapeMode(const Eref& e) const;
+ void setShapeMode(const Eref& e,unsigned int value); // variables SHELL=0, SLICE=SLAB=1, USERDEF=3.
+
+ double getLength(const Eref& e) const; // shell length
+ void setLength(const Eref& e,double value);
+
+ double getDiameter(const Eref& e) const; // shell diameter
+ void setDiameter(const Eref& e,double value);
+
+ double getThickness(const Eref& e) const; // shell thickness
+ void setThickness(const Eref& e,double value);
+
+ void setOuterArea( const Eref& e,double outerArea );
+ double getOuterArea(const Eref& e) const; // area of upper (outer) shell surface
+
+ void setInnerArea( const Eref& e,double innerArea );
+ double getInnerArea(const Eref& e) const; // area of lower (inner) shell surface
+
+ double getVolume(const Eref& e) const; // shell volume
+ void setVolume(const Eref& e,double volume); //
+
+
+ virtual double vGetActivation(const Eref& e) const = 0;
+ virtual void vSetActivation(const Eref& e,double value) = 0;
+
+ virtual double vGetBFree(const Eref& e) const = 0;
+ virtual void vSetBFree(const Eref& e,double value) = 0;
+
+ virtual double vGetBBound(const Eref& e) const = 0;
+ virtual void vSetBBound(const Eref& e,double value) = 0;
+
+
+ virtual double vGetBTot(const Eref& e) const = 0; // total buffer concentration in mM (free + bound)
+ virtual void vSetBTot(const Eref& e,double value) = 0;
+
+ virtual double vGetKf(const Eref& e) const = 0; // forward rate constant in 1/(mM*sec)
+ virtual void vSetKf(const Eref& e,double value) = 0;
+
+ virtual double vGetKb(const Eref& e) const = 0; // backward rate constant in 1/sec
+ virtual void vSetKb(const Eref& e,double value) = 0;
+
+ virtual double vGetD(const Eref& e) const = 0; // diffusion constant of buffer molecules, m^2/sec
+ virtual void vSetD(const Eref& e,double value) = 0;
+
+ virtual void vSetShapeMode(const Eref& e, unsigned int shapeMode ) = 0;
+ virtual unsigned int vGetShapeMode(const Eref& e) const = 0;
+
+ virtual void vSetLength(const Eref& e, double length ) = 0;
+ virtual double vGetLength(const Eref& e) const = 0;
+
+ virtual void vSetDiameter(const Eref& e, double diameter ) = 0;
+ virtual double vGetDiameter(const Eref& e) const = 0;
+
+ virtual void vSetThickness(const Eref& e, double thickness ) = 0;
+ virtual double vGetThickness(const Eref& e) const = 0;
+
+ virtual void vSetVolume(const Eref& e, double volume ) = 0;
+ virtual double vGetVolume(const Eref& e) const = 0;
+
+ virtual void vSetOuterArea(const Eref& e, double outerArea ) = 0;
+ virtual double vGetOuterArea(const Eref& e) const = 0;
+
+ virtual void vSetInnerArea(const Eref& e, double innerArea ) = 0;
+ virtual double vGetInnerArea(const Eref& e) const = 0;
+
+ static SrcFinfo4< double, double, double, double >* reactionOut();
+ static SrcFinfo2< double, double >* innerDifSourceOut();
+ static SrcFinfo2< double, double >* outerDifSourceOut();
+ static const Cinfo * initCinfo();
+ virtual void vSetSolver( const Eref& e, Id hsolve );
+ static void zombify( Element* orig, const Cinfo* zClass,
+ Id hsolve );
+private:
+
+
+
+};
+
+
+
+
+#endif //_DIFBUFFER_BASE_H
diff --git a/biophysics/DifShell.cpp b/biophysics/DifShell.cpp
index 98fe9a2b..35368325 100644
--- a/biophysics/DifShell.cpp
+++ b/biophysics/DifShell.cpp
@@ -1,253 +1,49 @@
/**********************************************************************
-** This program is part of 'MOOSE', the
-** Multiscale Object Oriented Simulation Environment.
-** copyright (C) 2003-2008
-** Upinder S. Bhalla, Niraj Dudani 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.
-**********************************************************************/
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani 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 "header.h"
+#include "DifShellBase.h"
#include "DifShell.h"
-#include "../utility/numutil.h"
-#include <cmath>
-static SrcFinfo1< double >* concentrationOut()
-{
- static SrcFinfo1< double > concentrationOut("concentrationOut",
- "Sends out concentration");
- return &concentrationOut;
-}
-static SrcFinfo2< double, double >* innerDifSourceOut(){
- static SrcFinfo2< double, double > sourceOut("innerDifSourceOut",
- "Sends out source information.");
- return & sourceOut;
-}
-
-static SrcFinfo2< double, double >* outerDifSourceOut(){
- static SrcFinfo2< double, double > sourceOut("outerDifSourceOut",
- "Sends out source information.");
- return & sourceOut;
-}
+const double DifShell::EPSILON = 1.0e-10;
+const double DifShell::F = 96485.3415; /* C / mol like in genesis */
const Cinfo* DifShell::initCinfo()
{
- static DestFinfo process( "process",
- "Handles process call",
- new ProcOpFunc< DifShell >( &DifShell::process_0 ) );
- static DestFinfo reinit( "reinit",
- "Reinit happens only in stage 0",
- new ProcOpFunc< DifShell >( &DifShell::reinit_0 ));
-
- static Finfo* processShared_0[] = {
- &process,
- &reinit
- };
-
- static SharedFinfo process_0(
- "process_0",
- "Here we create 2 shared finfos to attach with the Ticks. This is because we want to perform DifShell "
- "computations in 2 stages, much as in the Compartment object. "
- "In the first stage we send out the concentration value to other DifShells and Buffer elements. We also "
- "receive fluxes and currents and sum them up to compute ( dC / dt ). "
- "In the second stage we find the new C value using an explicit integration method. "
- "This 2-stage procedure eliminates the need to store and send prev_C values, as was common in GENESIS.",
- processShared_0,
- sizeof( processShared_0 ) / sizeof( Finfo* ));
-
- static DestFinfo process1( "process",
- "Handle process call",
- new ProcOpFunc< DifShell >( &DifShell::process_1 ) );
- static DestFinfo reinit1( "reinit",
- "Reinit happens only in stage 0",
- new ProcOpFunc< DifShell >( &DifShell::reinit_1)
- );
- static Finfo* processShared_1[] = {
- &process1, &reinit1
- };
-
- static SharedFinfo process_1( "process_1",
- "Second process call",
- processShared_1,
- sizeof( processShared_1 ) / sizeof( Finfo* ) );
-
-
- static DestFinfo reaction( "reaction",
- "Here the DifShell receives reaction rates (forward and backward), and concentrations for the "
- "free-buffer and bound-buffer molecules.",
- new OpFunc4< DifShell, double, double, double, double >( &DifShell::buffer ));
- static Finfo* bufferShared[] = {
- concentrationOut(), &reaction
+
+ static string doc[] =
+ {
+ "Name", "DifShell",
+ "Author", "Niraj Dudani. Ported to async13 by Subhasis Ray. Rewritten by Asia Jedrzejewska-Szmek",
+ "Description", "DifShell object: Models diffusion of an ion (typically calcium) within an "
+ "electric compartment. A DifShell is an iso-concentration region with respect to "
+ "the ion. Adjoining DifShells exchange flux of this ion, and also keep track of "
+ "changes in concentration due to pumping, buffering and channel currents, by "
+ "talking to the appropriate objects.",
};
- static SharedFinfo buffer( "buffer",
- "This is a shared message from a DifShell to a Buffer (FixBuffer or DifBuffer). "
- "During stage 0:\n "
- " - DifShell sends ion concentration\n"
- " - Buffer updates buffer concentration and sends it back immediately using a call-back.\n"
- " - DifShell updates the time-derivative ( dC / dt ) \n"
- "\n"
- "During stage 1: \n"
- " - DifShell advances concentration C \n\n"
- "This scheme means that the Buffer does not need to be scheduled, and it does its computations when "
- "it receives a cue from the DifShell. May not be the best idea, but it saves us from doing the above "
- "computations in 3 stages instead of 2." ,
- bufferShared,
- sizeof( bufferShared ) / sizeof( Finfo* ));
-
- static DestFinfo fluxFromOut( "fluxFromOut",
- "Destination message",
- new OpFunc2< DifShell, double, double > ( &DifShell::fluxFromOut ));
-
- static Finfo* innerDifShared[] = {
- innerDifSourceOut(),
- &fluxFromOut
- };
- static SharedFinfo innerDif( "innerDif",
- "This shared message (and the next) is between DifShells: adjoining shells exchange information to "
- "find out the flux between them. "
- "Using this message, an inner shell sends to, and receives from its outer shell." ,
- innerDifShared,
- sizeof( innerDifShared ) / sizeof( Finfo* ));
-
- static DestFinfo fluxFromIn( "fluxFromIn", "",
- new OpFunc2< DifShell, double, double> ( &DifShell::fluxFromIn ) );
- static Finfo* outerDifShared[] = {
- &fluxFromIn,
- outerDifSourceOut(),
- };
-
- static SharedFinfo outerDif( "outerDif",
- "Using this message, an outer shell sends to, and receives from its inner shell." ,
- outerDifShared,
- sizeof( outerDifShared ) / sizeof( Finfo* ));
-
- static ReadOnlyValueFinfo< DifShell, double> C( "C",
- "Concentration C is computed by the DifShell and is read-only",
- &DifShell::getC);
- static ValueFinfo< DifShell, double> Ceq( "Ceq", "",
- &DifShell::setCeq,
- &DifShell::getCeq);
- static ValueFinfo< DifShell, double> D( "D", "",
- &DifShell::setD,
- &DifShell::getD);
- static ValueFinfo< DifShell, double> valence( "valence", "",
- &DifShell::setValence,
- &DifShell::getValence);
- static ValueFinfo< DifShell, double> leak( "leak", "",
- &DifShell::setLeak,
- &DifShell::getLeak);
- static ValueFinfo< DifShell, unsigned int> shapeMode( "shapeMode", "",
- &DifShell::setShapeMode,
- &DifShell::getShapeMode);
- static ValueFinfo< DifShell, double> length( "length", "",
- &DifShell::setLength,
- &DifShell::getLength);
- static ValueFinfo< DifShell, double> diameter( "diameter", "",
- &DifShell::setDiameter,
- &DifShell::getDiameter );
- static ValueFinfo< DifShell, double> thickness( "thickness", "",
- &DifShell::setThickness,
- &DifShell::getThickness );
- static ValueFinfo< DifShell, double> volume( "volume", "",
- &DifShell::setVolume,
- &DifShell::getVolume );
- static ValueFinfo< DifShell, double> outerArea( "outerArea", "",
- &DifShell::setOuterArea,
- &DifShell::getOuterArea);
- static ValueFinfo< DifShell, double> innerArea( "innerArea", "",
- &DifShell::setInnerArea,
- &DifShell::getInnerArea );
+ static Dinfo< DifShell > dinfo;
+ static Cinfo difShellCinfo(
+ "DifShell",
+ DifShellBase::initCinfo(),
+ 0,
+ 0,
+ &dinfo,
+ doc,
+ sizeof( doc ) / sizeof( string ));
-
- static DestFinfo influx( "influx", "",
- new OpFunc1< DifShell, double > (&DifShell::influx ));
- static DestFinfo outflux( "outflux", "",
- new OpFunc1< DifShell, double >( &DifShell::influx ));
- static DestFinfo fInflux( "fInflux", "",
- new OpFunc2< DifShell, double, double >( &DifShell::fInflux ));
- static DestFinfo fOutflux( "fOutflux", "",
- new OpFunc2< DifShell, double, double> (&DifShell::fOutflux ));
- static DestFinfo storeInflux( "storeInflux", "",
- new OpFunc1< DifShell, double >( &DifShell::storeInflux ) );
- static DestFinfo storeOutflux( "storeOutflux", "",
- new OpFunc1< DifShell, double > ( &DifShell::storeOutflux ) );
- static DestFinfo tauPump( "tauPump","",
- new OpFunc2< DifShell, double, double >( &DifShell::tauPump ) );
- static DestFinfo eqTauPump( "eqTauPump", "",
- new OpFunc1< DifShell, double >( &DifShell::eqTauPump ) );
- static DestFinfo mmPump( "mmPump", "",
- new OpFunc2< DifShell, double, double >( &DifShell::mmPump ) );
- static DestFinfo hillPump( "hillPump", "",
- new OpFunc3< DifShell, double, double, unsigned int >( &DifShell::hillPump ) );
- static Finfo* difShellFinfos[] = {
- //////////////////////////////////////////////////////////////////
- // Field definitions
- //////////////////////////////////////////////////////////////////
- &C,
- &Ceq,
- &D,
- &valence,
- &leak,
- &shapeMode,
- &length,
- &diameter,
- &thickness,
- &volume,
- &outerArea,
- &innerArea,
- //////////////////////////////////////////////////////////////////
- // MsgSrc definitions
- //////////////////////////////////////////////////////////////////
-
- //////////////////////////////////////////////////////////////////
- // SharedFinfo definitions
- //////////////////////////////////////////////////////////////////
- &process_0,
- &process_1,
- &buffer,
- &innerDif,
- &outerDif,
- //////////////////////////////////////////////////////////////////
- // DestFinfo definitions
- //////////////////////////////////////////////////////////////////
- &influx,
- &outflux,
- &fInflux,
- &fOutflux,
- &storeInflux,
- &storeOutflux,
- &tauPump,
- &eqTauPump,
- &mmPump,
- &hillPump,
- };
-
- static string doc[] =
- {
- "Name", "DifShell",
- "Author", "Niraj Dudani. Ported to async13 by Subhasis Ray.",
- "Description", "DifShell object: Models diffusion of an ion (typically calcium) within an "
- "electric compartment. A DifShell is an iso-concentration region with respect to "
- "the ion. Adjoining DifShells exchange flux of this ion, and also keep track of "
- "changes in concentration due to pumping, buffering and channel currents, by "
- "talking to the appropriate objects.",
- };
- static Dinfo< DifShell > dinfo;
- static Cinfo difShellCinfo(
- "DifShell",
- Neutral::initCinfo(),
- difShellFinfos,
- sizeof( difShellFinfos ) / sizeof( Finfo* ),
- &dinfo,
- doc,
- sizeof( doc ) / sizeof( string ));
-
- return &difShellCinfo;
+ return &difShellCinfo;
}
+//Cinfo *object* corresponding to the class.
static const Cinfo* difShellCinfo = DifShell::initCinfo();
////////////////////////////////////////////////////////////////////////////////
@@ -255,499 +51,432 @@ static const Cinfo* difShellCinfo = DifShell::initCinfo();
////////////////////////////////////////////////////////////////////////////////
/// Faraday's constant (Coulomb / Mole)
-const double DifShell::F_ = 0.0;
+
DifShell::DifShell() :
- dCbyDt_( 0.0 ),
- C_( 0.0 ),
- Ceq_( 0.0 ),
- D_( 0.0 ),
- valence_( 0.0 ),
- leak_( 0.0 ),
- shapeMode_( 0 ),
- length_( 0.0 ),
- diameter_( 0.0 ),
- thickness_( 0.0 ),
- volume_( 0.0 ),
- outerArea_( 0.0 ),
- innerArea_( 0.0 )
+ dCbyDt_( 0.0 ),
+ Cmultiplier_(0.0),
+ C_( 0.0 ),
+ Ceq_( 0.0 ),
+ prevC_(0.0),
+ D_( 0.0 ),
+ valence_( 0.0 ),
+ leak_( 0.0 ),
+ shapeMode_(0),
+ length_(0),
+ diameter_(0),
+ thickness_(0),
+ volume_(0),
+ outerArea_(0),
+ innerArea_(0)
{ ; }
////////////////////////////////////////////////////////////////////////////////
// Field access functions
////////////////////////////////////////////////////////////////////////////////
/// C is a read-only field
-double DifShell::getC() const
+
+
+void DifShell::vSetC( const Eref& e, double C)
+{
+ if ( C < 0.0 ) {
+ cerr << "Error: DifShell: C cannot be negative!\n";
+ return;
+ }
+
+ C_ = C;
+ prevC_ = C_;
+}
+double DifShell::vGetC(const Eref& e) const
{
- return C_;
+ return C_;
}
-void DifShell::setCeq( double Ceq )
+void DifShell::vSetCeq( const Eref& e, double Ceq )
{
- if ( Ceq < 0.0 ) {
- cerr << "Error: DifShell: Ceq cannot be negative!\n";
- return;
- }
+ if ( Ceq < 0.0 ) {
+ cerr << "Error: DifShell: Ceq cannot be negative!\n";
+ return;
+ }
- Ceq_ = Ceq;
+ Ceq_ = Ceq;
}
-double DifShell::getCeq( ) const
+double DifShell::vGetCeq(const Eref& e ) const
{
- return Ceq_;
+ return Ceq_;
}
-void DifShell::setD( double D )
+void DifShell::vSetD(const Eref& e, double D )
{
- if ( D < 0.0 ) {
- cerr << "Error: DifShell: D cannot be negative!\n";
- return;
- }
+ if ( D < 0.0 ) {
+ cerr << "Error: DifShell: D cannot be negative!\n";
+ return;
+ }
- D_ = D;
+ D_ = D;
}
-double DifShell::getD( ) const
+double DifShell::vGetD(const Eref& e ) const
{
- return D_;
+ return D_;
}
-void DifShell::setValence( double valence )
+void DifShell::vSetValence(const Eref& e, double valence )
{
- if ( valence < 0.0 ) {
- cerr << "Error: DifShell: valence cannot be negative!\n";
- return;
- }
+ if ( valence < 0.0 ) {
+ cerr << "Error: DifShell: valence cannot be negative!\n";
+ return;
+ }
- valence_ = valence;
+ valence_ = valence;
}
-double DifShell::getValence( ) const
+double DifShell::vGetValence(const Eref& e ) const
{
- return valence_;
+ return valence_;
}
-void DifShell::setLeak( double leak )
+void DifShell::vSetLeak(const Eref& e, double leak )
{
- leak_ = leak;
+ leak_ = leak;
}
-double DifShell::getLeak( ) const
+double DifShell::vGetLeak(const Eref& e ) const
{
- return leak_;
+ return leak_;
}
-void DifShell::setShapeMode( unsigned int shapeMode )
+
+void DifShell::vSetShapeMode(const Eref& e, unsigned int shapeMode )
{
- if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) {
- cerr << "Error: DifShell: I only understand shapeModes 0, 1 and 3.\n";
- return;
- }
- shapeMode_ = shapeMode;
+ if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) {
+ cerr << "Error: DifShell: I only understand shapeModes 0, 1 and 3.\n";
+ return;
+ }
+ shapeMode_ = shapeMode;
}
-unsigned int DifShell::getShapeMode() const
+unsigned int DifShell::vGetShapeMode(const Eref& e) const
{
- return shapeMode_;
+ return shapeMode_;
}
-void DifShell::setLength( double length )
+void DifShell::vSetLength(const Eref& e, double length )
{
- if ( length < 0.0 ) {
- cerr << "Error: DifShell: length cannot be negative!\n";
- return;
- }
+ if ( length < 0.0 ) {
+ cerr << "Error: DifShell: length cannot be negative!\n";
+ return;
+ }
- length_ = length;
+ length_ = length;
}
-double DifShell::getLength( ) const
+double DifShell::vGetLength(const Eref& e ) const
{
- return length_;
+ return length_;
}
-void DifShell::setDiameter( double diameter )
+void DifShell::vSetDiameter(const Eref& e, double diameter )
{
- if ( diameter < 0.0 ) {
- cerr << "Error: DifShell: diameter cannot be negative!\n";
- return;
- }
+ if ( diameter < 0.0 ) {
+ cerr << "Error: DifShell: diameter cannot be negative!\n";
+ return;
+ }
- diameter_ = diameter;
+ diameter_ = diameter;
}
-double DifShell::getDiameter( ) const
+double DifShell::vGetDiameter(const Eref& e ) const
{
- return diameter_;
+ return diameter_;
}
-void DifShell::setThickness( double thickness )
+void DifShell::vSetThickness( const Eref& e, double thickness )
{
- if ( thickness < 0.0 ) {
- cerr << "Error: DifShell: thickness cannot be negative!\n";
- return;
- }
+ if ( thickness < 0.0 ) {
+ cerr << "Error: DifShell: thickness cannot be negative!\n";
+ return;
+ }
- thickness_ = thickness;
+ thickness_ = thickness;
}
-double DifShell::getThickness() const
+double DifShell::vGetThickness(const Eref& e) const
{
- return thickness_;
+ return thickness_;
}
-void DifShell::setVolume( double volume )
+void DifShell::vSetVolume(const Eref& e, double volume )
{
- if ( shapeMode_ != 3 )
- cerr << "Warning: DifShell: Trying to set volume, when shapeMode is not USER-DEFINED\n";
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifShell: Trying to set volume, when shapeMode is not USER-DEFINED\n";
- if ( volume < 0.0 ) {
- cerr << "Error: DifShell: volume cannot be negative!\n";
- return;
- }
+ if ( volume < 0.0 ) {
+ cerr << "Error: DifShell: volume cannot be negative!\n";
+ return;
+ }
- volume_ = volume;
+ volume_ = volume;
}
-double DifShell::getVolume( ) const
+double DifShell::vGetVolume(const Eref& e ) const
{
- return volume_;
+ return volume_;
}
-void DifShell::setOuterArea( double outerArea )
+void DifShell::vSetOuterArea(const Eref& e, double outerArea )
{
- if (shapeMode_ != 3 )
- cerr << "Warning: DifShell: Trying to set outerArea, when shapeMode is not USER-DEFINED\n";
+ if (shapeMode_ != 3 )
+ cerr << "Warning: DifShell: Trying to set outerArea, when shapeMode is not USER-DEFINED\n";
- if ( outerArea < 0.0 ) {
- cerr << "Error: DifShell: outerArea cannot be negative!\n";
- return;
- }
+ if ( outerArea < 0.0 ) {
+ cerr << "Error: DifShell: outerArea cannot be negative!\n";
+ return;
+ }
- outerArea_ = outerArea;
+ outerArea_ = outerArea;
}
-double DifShell::getOuterArea( ) const
+double DifShell::vGetOuterArea(const Eref& e ) const
{
- return outerArea_;
+ return outerArea_;
}
-void DifShell::setInnerArea( double innerArea )
+void DifShell::vSetInnerArea(const Eref& e, double innerArea )
{
- if ( shapeMode_ != 3 )
- cerr << "Warning: DifShell: Trying to set innerArea, when shapeMode is not USER-DEFINED\n";
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifShell: Trying to set innerArea, when shapeMode is not USER-DEFINED\n";
- if ( innerArea < 0.0 ) {
- cerr << "Error: DifShell: innerArea cannot be negative!\n";
- return;
- }
+ if ( innerArea < 0.0 ) {
+ cerr << "Error: DifShell: innerArea cannot be negative!\n";
+ return;
+ }
- innerArea_ = innerArea;
-}
-
-double DifShell::getInnerArea() const
-{
- return innerArea_;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-// Dest functions
-////////////////////////////////////////////////////////////////////////////////
-
-void DifShell::reinit_0( const Eref& e, ProcPtr p )
-{
- localReinit_0( e, p );
-}
-
-void DifShell::process_0( const Eref& e, ProcPtr p )
-{
- localProcess_0( e, p );
-}
-
-void DifShell::process_1(const Eref& e, ProcPtr p )
-{
- localProcess_1( e, p );
-}
-
-void DifShell::reinit_1(const Eref& e, ProcPtr p )
-{
- ;
-}
-
-void DifShell::buffer(
- double kf,
- double kb,
- double bFree,
- double bBound )
-{
- localBuffer( kf, kb, bFree, bBound );
-}
-
-void DifShell::fluxFromOut(
- double outerC,
- double outerThickness )
-{
- localFluxFromOut( outerC, outerThickness );
-}
-
-void DifShell::fluxFromIn(
- double innerC,
- double innerThickness )
-{
- localFluxFromIn( innerC, innerThickness );
-}
-
-void DifShell::influx(
- double I )
-{
- localInflux( I );
+ innerArea_ = innerArea;
}
-void DifShell::outflux(
- double I )
+double DifShell::vGetInnerArea(const Eref& e) const
{
- localOutflux( I );
+ return innerArea_;
}
-void DifShell::fInflux(
- double I,
- double fraction )
-{
-localFInflux( I, fraction );
-}
-void DifShell::fOutflux(
- double I,
- double fraction )
-{
-localFOutflux( I, fraction );
-}
-
-void DifShell::storeInflux(
- double flux )
-{
-localStoreInflux( flux );
-}
-
-void DifShell::storeOutflux(
- double flux )
-{
-localStoreOutflux( flux );
-}
-
-void DifShell::tauPump(
- double kP,
- double Ceq )
-{
-localTauPump( kP, Ceq );
-}
-
-void DifShell::eqTauPump(
- double kP )
-{
-localEqTauPump( kP );
-}
-
-void DifShell::mmPump(
- double vMax,
- double Kd )
-{
-localMMPump( vMax, Kd );
-}
-
-void DifShell::hillPump(
- double vMax,
- double Kd,
- unsigned int hill )
-{
-localHillPump( vMax, Kd, hill );
-}
////////////////////////////////////////////////////////////////////////////////
// Local dest functions
////////////////////////////////////////////////////////////////////////////////
-
-void DifShell::localReinit_0( const Eref& e, ProcPtr p )
+double DifShell::integrate( double state, double dt, double A, double B )
{
- dCbyDt_ = leak_;
-
- double Pi = M_PI;
- double dOut = diameter_;
- double dIn = diameter_ - thickness_;
-
- switch ( shapeMode_ )
- {
- /*
- * Onion Shell
- */
- case 0:
- if ( length_ == 0.0 ) { // Spherical shell
- volume_ = ( Pi / 6.0 ) * ( dOut * dOut * dOut - dIn * dIn * dIn );
- outerArea_ = Pi * dOut * dOut;
- innerArea_ = Pi * dIn * dIn;
- } else { // Cylindrical shell
- volume_ = ( Pi * length_ / 4.0 ) * ( dOut * dOut - dIn * dIn );
- outerArea_ = Pi * dOut * length_;
- innerArea_ = Pi * dIn * length_;
- }
+ if ( B > EPSILON ) {
+ double x = exp( -B * dt );
+ return state * x + ( A / B ) * ( 1 - x );
+ }
+
+ return state + A * dt ;
+}
+
+void DifShell::vReinit( const Eref& e, ProcPtr p )
+{
+ dCbyDt_ = leak_;
+ Cmultiplier_ = 0;
+
+ double rOut = diameter_/2.;
+
+ double rIn = rOut - thickness_;
+
+ if (rIn <0)
+ rIn = 0.;
+
+ switch ( shapeMode_ )
+ {
+ /*
+ * Onion Shell
+ */
+ case 0:
+ if ( length_ == 0.0 ) { // Spherical shell
+ volume_ = 4./3.* M_PI * ( rOut * rOut * rOut - rIn * rIn * rIn );
+ outerArea_ = 4*M_PI * rOut * rOut;
+ innerArea_ = 4*M_PI * rIn * rIn;
+ } else { // Cylindrical shell
+ volume_ = ( M_PI * length_ ) * ( rOut * rOut - rIn * rIn );
+ outerArea_ = 2*M_PI * rOut * length_;
+ innerArea_ = 2*M_PI * rIn * length_;
+ }
- break;
+ break;
- /*
- * Cylindrical Slice
- */
- case 1:
- volume_ = Pi * diameter_ * diameter_ * thickness_ / 4.0;
- outerArea_ = Pi * diameter_ * diameter_ / 4.0;
- innerArea_ = outerArea_;
- break;
+ /*
+ * Cylindrical Slice
+ */
+ case 1:
+ volume_ = M_PI * diameter_ * diameter_ * thickness_ / 4.0;
+ outerArea_ = M_PI * diameter_ * diameter_ / 4.0;
+ innerArea_ = outerArea_;
+ break;
- /*
- * User defined
- */
- case 3:
- // Nothing to be done here. Volume and inner-, outer areas specified by
- // user.
- break;
+ /*
+ * User defined
+ */
+ case 3:
+ // Nothing to be done here. Volume and inner-, outer areas specified by
+ // user.
+ break;
- default:
- assert( 0 );
- }
+ default:
+ assert( 0 );
+ }
+ C_= Ceq_;
+ prevC_ = Ceq_;
+ concentrationOut()->send( e, C_ );
+ innerDifSourceOut()->send( e, prevC_, thickness_ );
+ outerDifSourceOut()->send( e, prevC_, thickness_ );
}
-void DifShell::localProcess_0( const Eref & e, ProcPtr p )
+void DifShell::vProcess( const Eref & e, ProcPtr p )
{
- /**
- * Send ion concentration and thickness to adjacent DifShells. They will
- * then compute their incoming fluxes.
- */
- innerDifSourceOut()->send( e, C_, thickness_ );
- outerDifSourceOut()->send( e, C_, thickness_ );
-
- /**
- * Send ion concentration to ion buffers. They will send back information on
- * the reaction (forward / backward rates ; free / bound buffer concentration)
- * immediately, which this DifShell will use to find amount of ion captured
- * or released in the current time-step.
- */
- concentrationOut()->send( e, C_ );
-}
+ /**
+ * Send ion concentration and thickness to adjacent DifShells. They will
+ * then compute their incoming fluxes.
+ */
+
+
+ C_ = integrate(C_,p->dt,dCbyDt_,Cmultiplier_);
+
+ /**
+ * Send ion concentration to ion buffers. They will send back information on
+ * the reaction (forward / backward rates ; free / bound buffer concentration)
+ * immediately, which this DifShell will use to find amount of ion captured
+ * or released in the current time-step.
+ */
+ prevC_ = C_;
-void DifShell::localProcess_1( const Eref& e, ProcPtr p )
-{
- C_ += dCbyDt_ * p->dt;
- dCbyDt_ = leak_;
-}
+ //cout<<"Shell "<< C_<<" ";
+ dCbyDt_ = leak_;
+ Cmultiplier_ = 0;
+ innerDifSourceOut()->send( e, prevC_, thickness_ );
+ outerDifSourceOut()->send( e, prevC_, thickness_ );
+ concentrationOut()->send( e, C_ );
-void DifShell::localBuffer(
- double kf,
- double kb,
- double bFree,
- double bBound )
+}
+void DifShell::vBuffer(const Eref& e,
+ double kf,
+ double kb,
+ double bFree,
+ double bBound )
{
- dCbyDt_ += -kf * bFree * C_ + kb * bBound;
+ dCbyDt_ += kb * bBound;
+ Cmultiplier_ += kf * bFree;
}
-void DifShell::localFluxFromOut( double outerC, double outerThickness )
+void DifShell::vFluxFromOut(const Eref& e, double outerC, double outerThickness )
{
- double dx = ( outerThickness + thickness_ ) / 2.0;
-
- /**
- * We could pre-compute ( D / Volume ), but let us leave the optimizations
- * for the solver.
- */
- dCbyDt_ += ( D_ / volume_ ) * ( outerArea_ / dx ) * ( outerC - C_ );
+ double diff =2.* D_ /volume_ * outerArea_ / (outerThickness + thickness_) ;
+ //influx from outer shell
+ /**
+ * We could pre-compute ( D / Volume ), but let us leave the optimizations
+ * for the solver.
+ */
+
+ dCbyDt_ += diff * outerC;
+ Cmultiplier_ += diff ;
}
-void DifShell::localFluxFromIn( double innerC, double innerThickness )
+void DifShell::vFluxFromIn(const Eref& e, double innerC, double innerThickness )
{
- double dx = ( innerThickness + thickness_ ) / 2.0;
-
- dCbyDt_ += ( D_ / volume_ ) * ( innerArea_ / dx ) * ( innerC - C_ );
+ //influx from inner shell
+ //double dx = ( innerThickness + thickness_ ) / 2.0;
+ double diff = 2.* D_/volume_ * innerArea_ / (innerThickness + thickness_);
+ dCbyDt_ += diff * innerC ;
+ Cmultiplier_ += diff ;
}
-void DifShell::localInflux( double I )
+void DifShell::vInflux(const Eref& e, double I )
{
- /**
- * I: Amperes
- * F_: Faraday's constant: Coulomb / mole
- * valence_: charge on ion: dimensionless
- */
- dCbyDt_ += I / ( F_ * valence_ * volume_ );
+ /**
+ * I: Amperes
+ * F_: Faraday's constant: Coulomb / mole
+ * valence_: charge on ion: dimensionless
+ */
+ dCbyDt_ += I / ( F * valence_ * volume_ );
}
/**
* Same as influx, except subtracting.
*/
-void DifShell::localOutflux( double I )
+void DifShell::vOutflux(const Eref& e, double I )
{
- dCbyDt_ -= I / ( F_ * valence_ * volume_ );
+ dCbyDt_ -= I / ( F * valence_ * volume_ );
}
-void DifShell::localFInflux( double I, double fraction )
+void DifShell::vFInflux(const Eref& e, double I, double fraction )
{
- dCbyDt_ += fraction * I / ( F_ * valence_ * volume_ );
+ dCbyDt_ += fraction * I / ( F * valence_ * volume_ );
}
-void DifShell::localFOutflux( double I, double fraction )
+void DifShell::vFOutflux(const Eref& e, double I, double fraction )
{
- dCbyDt_ -= fraction * I / ( F_ * valence_ * volume_ );
+ dCbyDt_ -= fraction * I / ( F * valence_ * volume_ );
}
-void DifShell::localStoreInflux( double flux )
+void DifShell::vStoreInflux(const Eref& e, double flux )
{
- dCbyDt_ += flux / volume_;
+ dCbyDt_ += flux / volume_;
}
-void DifShell::localStoreOutflux( double flux )
+void DifShell::vStoreOutflux(const Eref& e, double flux )
{
- dCbyDt_ -= flux / volume_;
+ dCbyDt_ -= flux / volume_;
}
-void DifShell::localTauPump( double kP, double Ceq )
+void DifShell::vTauPump(const Eref& e, double kP, double Ceq )
{
- dCbyDt_ += -kP * ( C_ - Ceq );
+ //dCbyDt_ += -kP * ( C_ - Ceq );
+ dCbyDt_ += kP*Ceq;
+ Cmultiplier_ -= kP;
}
/**
- * Same as tauPump, except that we use the local value of Ceq.
+ * Same as tauPump, except that we use the v value of Ceq.
*/
-void DifShell::localEqTauPump( double kP )
-{
- dCbyDt_ += -kP * ( C_ - Ceq_ );
-}
-
-void DifShell::localMMPump( double vMax, double Kd )
-{
- dCbyDt_ += -( vMax / volume_ ) * ( C_ / ( C_ + Kd ) );
-}
-
-void DifShell::localHillPump( double vMax, double Kd, unsigned int hill )
-{
- double ch;
- switch( hill )
- {
- case 0:
- ch = 1.0;
- break;
- case 1:
- ch = C_;
- break;
- case 2:
- ch = C_ * C_;
- break;
- case 3:
- ch = C_ * C_ * C_;
- break;
- case 4:
- ch = C_ * C_;
- ch = ch * ch;
- break;
- default:
- ch = pow( C_, static_cast< double >( hill ) );
- };
+void DifShell::vEqTauPump(const Eref& e, double kP )
+{
+ dCbyDt_ += kP*Ceq_;
+ Cmultiplier_ -= kP;
+}
+
+void DifShell::vMMPump(const Eref& e, double vMax, double Kd )
+{
+ Cmultiplier_ += ( vMax / volume_ ) / ( C_ + Kd ) ;
+}
+
+void DifShell::vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill )
+{
+ //This needs fixing
+ double ch;
+ switch( hill )
+ {
+ case 0:
+ ch = 1.0;
+ break;
+ case 1:
+ ch = C_;
+ break;
+ case 2:
+ ch = C_ * C_;
+ break;
+ case 3:
+ ch = C_ * C_ * C_;
+ break;
+ case 4:
+ ch = C_ * C_;
+ ch = ch * ch;
+ break;
+ default:
+ ch = pow( C_, static_cast< double >( hill ) );
+ };
- dCbyDt_ += -( vMax / volume_ ) * ( ch / ( ch + Kd ) );
+ dCbyDt_ += -( vMax / volume_ ) * ( ch / ( ch + Kd ) );
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/biophysics/DifShell.h b/biophysics/DifShell.h
index 42f70600..42cf643d 100644
--- a/biophysics/DifShell.h
+++ b/biophysics/DifShell.h
@@ -1,156 +1,105 @@
/**********************************************************************
-** This program is part of 'MOOSE', the
-** Multiscale Object Oriented Simulation Environment.
-** copyright (C) 2003-2008
-** Upinder S. Bhalla, Niraj Dudani 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 _DifShell_h
-#define _DifShell_h
-
-class DifShell
-{
- public:
- DifShell();
-
- /////////////////////////////////////////////////////////////
- // Field access functions
- /////////////////////////////////////////////////////////////
- double getC( ) const;
-
- void setCeq(double Ceq );
- double getCeq() const;
-
- void setD( double D );
- double getD() const;
-
- void setValence( double valence );
- double getValence() const;
-
- void setLeak( double leak );
- double getLeak() const;
-
- void setShapeMode( unsigned int shapeMode );
- unsigned int getShapeMode() const;
-
- void setLength( double length );
- double getLength() const;
-
- void setDiameter( double diameter );
- double getDiameter() const;
-
- void setThickness( double thickness );
- double getThickness() const;
-
- void setVolume( double volume );
- double getVolume() const;
-
- void setOuterArea( double outerArea );
- double getOuterArea() const;
-
- void setInnerArea( double innerArea );
- double getInnerArea() const;
-
- /////////////////////////////////////////////////////////////
- // Dest functions
- /////////////////////////////////////////////////////////////
- void reinit_0( const Eref & e, ProcPtr p );
-
- void process_0(const Eref & e, ProcPtr p );
-
- void process_1(const Eref & e, ProcPtr p );
-
- void reinit_1(const Eref & e, ProcPtr p ); // dummyFunc
-
- void buffer(
- double kf,
- double kb,
- double bFree,
- double bBound );
-
- void fluxFromOut(
- double outerC,
- double outerThickness );
-
- void fluxFromIn(
- double innerC,
- double innerThickness );
-
- void influx(
- double I );
-
- void outflux(
- double I );
-
- void fInflux(
- double I,
- double fraction );
-
- void fOutflux(
- double I,
- double fraction );
-
- void storeInflux(
- double flux );
-
- void storeOutflux(
- double flux );
-
- void tauPump(
- double kP,
- double Ceq );
-
- void eqTauPump(
- double kP );
-
- void mmPump(
- double vMax,
- double Kd );
-
- void hillPump(
- double vMax,
- double Kd,
- unsigned int hill );
-
- static const Cinfo * initCinfo();
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani 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 _DIFSHELL_H
+#define _DIFSHELL_H
+
+class DifShell: public DifShellBase{
+ public:
+ DifShell();
+
+ /////////////////////////////////////////////////////////////
+ // Dest functions
+ /////////////////////////////////////////////////////////////
+ void vReinit( const Eref & e, ProcPtr p );
+ void vProcess(const Eref & e, ProcPtr p );
+ void vBuffer(const Eref& e, double kf, double kb, double bFree, double bBound );
+ void vFluxFromOut(const Eref& e, double outerC, double outerThickness );
+ void vFluxFromIn(const Eref& e, double innerC, double innerThickness );
+ void vInflux(const Eref& e, double I );
+ void vOutflux(const Eref& e, double I );
+ void vFInflux(const Eref& e, double I, double fraction );
+ void vFOutflux(const Eref& e, double I, double fraction );
+ void vStoreInflux(const Eref& e, double flux );
+ void vStoreOutflux(const Eref& e, double flux );
+ void vTauPump(const Eref& e, double kP, double Ceq );
+ void vEqTauPump(const Eref& e, double kP );
+ void vMMPump(const Eref& e, double vMax, double Kd );
+ void vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill );
+
+ /////////////////////////////////////////////////////////////
+ // Field access functions
+ /////////////////////////////////////////////////////////////
+
+ void vSetC(const Eref& e,double C);
+ double vGetC( const Eref& e) const;
+
+ void vSetCeq(const Eref& e,double Ceq );
+ double vGetCeq(const Eref& e) const;
+
+ void vSetD(const Eref& e, double D );
+ double vGetD(const Eref& e) const;
+
+ void vSetValence(const Eref& e, double valence );
+ double vGetValence(const Eref& e) const;
+
+ void vSetLeak(const Eref& e, double leak );
+ double vGetLeak(const Eref& e) const;
+
+ void vSetShapeMode(const Eref& e, unsigned int shapeMode );
+ unsigned int vGetShapeMode(const Eref& e) const;
+
+ void vSetLength(const Eref& e, double length );
+ double vGetLength(const Eref& e) const;
+
+ void vSetDiameter(const Eref& e, double diameter );
+ double vGetDiameter(const Eref& e) const;
+
+ void vSetThickness(const Eref& e, double thickness );
+ double vGetThickness(const Eref& e) const;
+
+ void vSetVolume(const Eref& e, double volume );
+ double vGetVolume(const Eref& e) const;
+
+ void vSetOuterArea(const Eref& e, double outerArea );
+ double vGetOuterArea(const Eref& e) const;
+
+ void vSetInnerArea(const Eref& e, double innerArea );
+ double vGetInnerArea(const Eref& e) const;
+ static const Cinfo * initCinfo();
+
- private:
- void localReinit_0( const Eref & e, ProcPtr p );
- void localProcess_0( const Eref & e, ProcPtr p );
- void localProcess_1( const Eref & e, ProcPtr p );
- void localBuffer( double kf, double kb, double bFree, double bBound );
- void localFluxFromOut( double outerC, double outerThickness );
- void localFluxFromIn( double innerC, double innerThickness );
- void localInflux( double I );
- void localOutflux( double I );
- void localFInflux( double I, double fraction );
- void localFOutflux( double I, double fraction );
- void localStoreInflux( double flux );
- void localStoreOutflux( double flux );
- void localTauPump( double kP, double Ceq );
- void localEqTauPump( double kP );
- void localMMPump( double vMax, double Kd );
- void localHillPump( double vMax, double Kd, unsigned int hill );
-
- double dCbyDt_;
- double C_;
- double Ceq_;
- double D_;
- double valence_;
- double leak_;
- unsigned int shapeMode_;
- double length_;
- double diameter_;
- double thickness_;
- double volume_;
- double outerArea_;
- double innerArea_;
-
- /// Faraday's constant (Coulomb / Mole)
- static const double F_;
+ private:
+
+ double integrate( double state, double dt, double A, double B );
+
+ double dCbyDt_;
+ double Cmultiplier_;
+ double C_;
+ double Ceq_;
+ double prevC_;
+ double D_;
+ double valence_;
+ double leak_;
+ unsigned int shapeMode_;
+ double length_;
+ double diameter_;
+ double thickness_;
+ double volume_;
+ double outerArea_;
+ double innerArea_;
+
+ static const double EPSILON;
+ /// Faraday's constant (Coulomb / Mole)
+ static const double F;
+
};
-#endif // _DifShell_h
+#endif // _DIFSHELL_H
diff --git a/biophysics/DifShellBase.cpp b/biophysics/DifShellBase.cpp
new file mode 100644
index 00000000..e6d1e533
--- /dev/null
+++ b/biophysics/DifShellBase.cpp
@@ -0,0 +1,513 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani 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 "header.h"
+#include "DifShellBase.h"
+#include "ElementValueFinfo.h"
+
+SrcFinfo1< double >* DifShellBase::concentrationOut()
+{
+ static SrcFinfo1< double > concentrationOut("concentrationOut",
+ "Sends out concentration");
+ return &concentrationOut;
+}
+
+SrcFinfo2< double, double >* DifShellBase::innerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("innerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+SrcFinfo2< double, double >* DifShellBase::outerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("outerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+const Cinfo* DifShellBase::initCinfo()
+{
+
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< DifShellBase>(&DifShellBase::process ) );
+ static DestFinfo reinit( "reinit",
+ "Reinit happens only in stage 0",
+ new ProcOpFunc< DifShellBase>( &DifShellBase::reinit ));
+
+ static Finfo* processShared[] = {
+ &process, &reinit
+ };
+
+ static SharedFinfo proc(
+ "proc",
+ "Shared message to receive Process message from scheduler",
+ processShared, sizeof( processShared ) / sizeof( Finfo* ));
+
+
+
+
+ static DestFinfo reaction( "reaction",
+ "Here the DifShell receives reaction rates (forward and backward), and concentrations for the "
+ "free-buffer and bound-buffer molecules.",
+ new EpFunc4< DifShellBase, double, double, double, double >( &DifShellBase::buffer ));
+
+ static Finfo* bufferShared[] = {
+ DifShellBase::concentrationOut(), &reaction
+ };
+
+ static SharedFinfo buffer( "buffer",
+ "This is a shared message from a DifShell to a Buffer (FixBuffer or DifBuffer). " ,
+ bufferShared,
+ sizeof( bufferShared ) / sizeof( Finfo* ));
+ /////
+
+
+
+
+ static DestFinfo fluxFromOut( "fluxFromOut",
+ "Destination message",
+ new EpFunc2< DifShellBase, double, double > ( &DifShellBase::fluxFromOut ));
+
+ static Finfo* innerDifShared[] = {
+ &fluxFromOut, DifShellBase::innerDifSourceOut(),
+ };
+ static SharedFinfo innerDif( "innerDif",
+ "This shared message (and the next) is between DifShells: adjoining shells exchange information to "
+ "find out the flux between them. "
+ "Using this message, an inner shell sends to, and receives from its outer shell." ,
+ innerDifShared,
+ sizeof( innerDifShared ) / sizeof( Finfo* ));
+
+ static DestFinfo fluxFromIn( "fluxFromIn", "",
+ new EpFunc2< DifShellBase, double, double> ( &DifShellBase::fluxFromIn ) );
+
+ static Finfo* outerDifShared[] = {
+ &fluxFromIn, DifShellBase::outerDifSourceOut(),
+ };
+
+ static SharedFinfo outerDif( "outerDif",
+ "Using this message, an outer shell sends to, and receives from its inner shell." ,
+ outerDifShared,
+ sizeof( outerDifShared ) / sizeof( Finfo* ));
+
+ static ElementValueFinfo< DifShellBase, double> C( "C",
+ "Concentration C",// is computed by the DifShell",
+ &DifShellBase::setC,
+ &DifShellBase::getC);
+ static ElementValueFinfo< DifShellBase, double> Ceq( "Ceq", "",
+ &DifShellBase::setCeq,
+ &DifShellBase::getCeq);
+ static ElementValueFinfo< DifShellBase, double> D( "D", "",
+ &DifShellBase::setD,
+ &DifShellBase::getD);
+ static ElementValueFinfo< DifShellBase, double> valence( "valence", "",
+ &DifShellBase::setValence,
+ &DifShellBase::getValence);
+ static ElementValueFinfo< DifShellBase, double> leak( "leak", "",
+ &DifShellBase::setLeak,
+ &DifShellBase::getLeak);
+ static ElementValueFinfo< DifShellBase, unsigned int> shapeMode( "shapeMode", "",
+ &DifShellBase::setShapeMode,
+ &DifShellBase::getShapeMode);
+ static ElementValueFinfo< DifShellBase, double> length( "length", "",
+ &DifShellBase::setLength,
+ &DifShellBase::getLength);
+ static ElementValueFinfo< DifShellBase, double> diameter( "diameter", "",
+ &DifShellBase::setDiameter,
+ &DifShellBase::getDiameter );
+ static ElementValueFinfo< DifShellBase, double> thickness( "thickness", "",
+ &DifShellBase::setThickness,
+ &DifShellBase::getThickness );
+ static ElementValueFinfo< DifShellBase, double> volume( "volume", "",
+ &DifShellBase::setVolume,
+ &DifShellBase::getVolume );
+ static ElementValueFinfo< DifShellBase, double> outerArea( "outerArea", "",
+ &DifShellBase::setOuterArea,
+ &DifShellBase::getOuterArea);
+ static ElementValueFinfo< DifShellBase, double> innerArea( "innerArea", "",
+ &DifShellBase::setInnerArea,
+ &DifShellBase::getInnerArea );
+
+ static DestFinfo mmPump( "mmPump", "Here DifShell receives pump outflux",
+ new EpFunc2< DifShellBase, double, double >( &DifShellBase::mmPump ) );
+ static DestFinfo influx( "influx", "",
+ new EpFunc1< DifShellBase, double > (&DifShellBase::influx ));
+ static DestFinfo outflux( "outflux", "",
+ new EpFunc1< DifShellBase, double >( &DifShellBase::influx ));
+ static DestFinfo fInflux( "fInflux", "",
+ new EpFunc2< DifShellBase, double, double >( &DifShellBase::fInflux ));
+ static DestFinfo fOutflux( "fOutflux", "",
+ new EpFunc2< DifShellBase, double, double> (&DifShellBase::fOutflux ));
+ static DestFinfo storeInflux( "storeInflux", "",
+ new EpFunc1< DifShellBase, double >( &DifShellBase::storeInflux ) );
+ static DestFinfo storeOutflux( "storeOutflux", "",
+ new EpFunc1< DifShellBase, double > ( &DifShellBase::storeOutflux ) );
+ static DestFinfo tauPump( "tauPump","",
+ new EpFunc2< DifShellBase, double, double >( &DifShellBase::tauPump ) );
+ static DestFinfo eqTauPump( "eqTauPump", "",
+ new EpFunc1< DifShellBase, double >( &DifShellBase::eqTauPump ) );
+ static DestFinfo hillPump( "hillPump", "",
+ new EpFunc3< DifShellBase, double, double, unsigned int >( &DifShellBase::hillPump ) );
+ static Finfo* difShellBaseFinfos[] = {
+ //////////////////////////////////////////////////////////////////
+ // Field definitions
+ //////////////////////////////////////////////////////////////////
+ &C,
+ &Ceq,
+ &D,
+ &valence,
+ &leak,
+ &shapeMode,
+ &length,
+ &diameter,
+ &thickness,
+ &volume,
+ &outerArea,
+ &innerArea,
+ //////////////////////////////////////////////////////////////////
+ // MsgSrc definitions
+ //////////////////////////////////////////////////////////////////
+
+ //////////////////////////////////////////////////////////////////
+ // SharedFinfo definitions
+ //////////////////////////////////////////////////////////////////
+ &proc,
+ &buffer,
+ concentrationOut(),
+ innerDifSourceOut(),
+ outerDifSourceOut(),
+ &innerDif,
+ &outerDif,
+ //////////////////////////////////////////////////////////////////
+ // DestFinfo definitions
+ //////////////////////////////////////////////////////////////////
+ &influx,
+ &outflux,
+ &fInflux,
+ &fOutflux,
+ &storeInflux,
+ &storeOutflux,
+ &tauPump,
+ &eqTauPump,
+ &mmPump,
+ &hillPump,
+ };
+
+ static string doc[] =
+ {
+ "Name", "DifShellBase",
+ "Author", "Niraj Dudani. Ported to async13 by Subhasis Ray/Asia Jedrzejewska-Szmek",
+ "Description", "DifShell object: Models diffusion of an ion (typically calcium) within an "
+ "electric compartment. A DifShell is an iso-concentration region with respect to "
+ "the ion. Adjoining DifShells exchange flux of this ion, and also keep track of "
+ "changes in concentration due to pumping, buffering and channel currents, by "
+ "talking to the appropriate objects.",
+ };
+ static ZeroSizeDinfo< int > dinfo;
+ //static Dinfo< DifShellBase> dinfo;
+ static Cinfo difShellBaseCinfo(
+ "DifShellBase",
+ Neutral::initCinfo(),
+ difShellBaseFinfos,
+ sizeof( difShellBaseFinfos ) / sizeof( Finfo* ),
+ &dinfo,
+ doc,
+ sizeof( doc ) / sizeof( string ));
+
+ return &difShellBaseCinfo;
+}
+
+//Cinfo *object* corresponding to the class.
+static const Cinfo* difShellBaseCinfo = DifShellBase::initCinfo();
+
+DifShellBase::DifShellBase()
+{ ; }
+
+void DifShellBase::setC( const Eref& e, double C)
+{
+ vSetC(e,C);
+}
+double DifShellBase::getC(const Eref& e) const
+{
+ return vGetC(e);
+}
+
+void DifShellBase::setCeq( const Eref& e, double Ceq )
+{
+ vSetCeq(e,Ceq);
+}
+
+double DifShellBase::getCeq(const Eref& e ) const
+{
+ return vGetCeq(e);
+}
+
+void DifShellBase::setD(const Eref& e, double D )
+{
+ vSetD(e,D);
+}
+
+double DifShellBase::getD(const Eref& e ) const
+{
+ return vGetD(e);
+}
+
+void DifShellBase::setValence(const Eref& e, double valence )
+{
+ vSetValence(e,valence);
+}
+
+double DifShellBase::getValence(const Eref& e ) const
+{
+ return vGetValence(e);
+}
+
+void DifShellBase::setLeak(const Eref& e, double leak )
+{
+ vSetLeak(e,leak);
+}
+
+double DifShellBase::getLeak(const Eref& e ) const
+{
+ return vGetLeak(e);
+}
+
+void DifShellBase::setShapeMode(const Eref& e, unsigned int shapeMode )
+{
+ vSetShapeMode(e,shapeMode);
+}
+
+unsigned int DifShellBase::getShapeMode(const Eref& e) const
+{
+ return vGetShapeMode(e);
+}
+
+void DifShellBase::setLength(const Eref& e, double length )
+{
+ vSetLength(e,length);
+}
+
+double DifShellBase::getLength(const Eref& e ) const
+{
+ return vGetLength(e);
+}
+
+void DifShellBase::setDiameter(const Eref& e, double diameter )
+{
+ vSetDiameter(e,diameter);
+}
+
+double DifShellBase::getDiameter(const Eref& e ) const
+{
+ return vGetDiameter(e);
+}
+
+void DifShellBase::setThickness( const Eref& e, double thickness )
+{
+ vSetThickness(e,thickness);
+}
+
+double DifShellBase::getThickness(const Eref& e) const
+{
+ return vGetThickness(e);
+}
+
+void DifShellBase::setVolume(const Eref& e, double volume )
+{
+ vSetVolume(e,volume);
+}
+
+double DifShellBase::getVolume(const Eref& e ) const
+{
+ return vGetVolume(e);
+}
+
+void DifShellBase::setOuterArea(const Eref& e, double outerArea )
+{
+ vSetOuterArea(e,outerArea);
+}
+
+double DifShellBase::getOuterArea(const Eref& e ) const
+{
+ return vGetOuterArea(e);
+}
+
+void DifShellBase::setInnerArea(const Eref& e, double innerArea )
+{
+ vSetInnerArea(e,innerArea);
+}
+
+double DifShellBase::getInnerArea(const Eref& e) const
+{
+ return vGetInnerArea(e);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Dest functions
+////////////////////////////////////////////////////////////////////////////////
+
+void DifShellBase::reinit( const Eref& e, ProcPtr p )
+{
+ vReinit( e, p );
+}
+
+void DifShellBase::process( const Eref& e, ProcPtr p )
+{
+ vProcess( e, p );
+}
+
+void DifShellBase::buffer(
+ const Eref& e,
+ double kf,
+ double kb,
+ double bFree,
+ double bBound )
+{
+ vBuffer( e, kf, kb, bFree, bBound );
+}
+
+void DifShellBase::fluxFromOut(const Eref& e,
+ double outerC,
+ double outerThickness )
+{
+ vFluxFromOut(e, outerC, outerThickness );
+}
+
+void DifShellBase::fluxFromIn(
+ const Eref& e,
+ double innerC,
+ double innerThickness )
+{
+ vFluxFromIn( e, innerC, innerThickness );
+}
+
+void DifShellBase::influx(const Eref& e,
+ double I )
+{
+ vInflux( e, I );
+}
+
+void DifShellBase::outflux(const Eref& e,
+ double I )
+{
+ vOutflux(e, I );
+}
+
+void DifShellBase::fInflux(const Eref& e,
+ double I,
+ double fraction )
+{
+ vFInflux(e, I, fraction );
+}
+
+void DifShellBase::fOutflux(const Eref& e,
+ double I,
+ double fraction )
+{
+ vFOutflux(e, I, fraction );
+}
+
+void DifShellBase::storeInflux(const Eref& e,
+ double flux )
+{
+ vStoreInflux( e, flux );
+}
+
+void DifShellBase::storeOutflux(const Eref& e,
+ double flux )
+{
+ vStoreOutflux(e, flux );
+}
+
+void DifShellBase::tauPump(const Eref& e,
+ double kP,
+ double Ceq )
+{
+ vTauPump( e, kP, Ceq );
+}
+
+void DifShellBase::eqTauPump(const Eref& e,
+ double kP )
+{
+ vEqTauPump(e, kP );
+}
+
+void DifShellBase::mmPump(const Eref& e,
+ double vMax,
+ double Kd )
+{
+ vMMPump(e, vMax, Kd );
+}
+
+void DifShellBase::hillPump(const Eref& e,
+ double vMax,
+ double Kd,
+ unsigned int hill )
+{
+ vHillPump(e, vMax, Kd, hill );
+}
+void DifShellBase::vSetSolver( const Eref& e, Id hsolve )
+{;}
+
+void DifShellBase::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;
+ unsigned int len = 12;
+ vector< double > data( num * len );
+
+ unsigned int j = 0;
+
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ const DifShellBase* ds =
+ reinterpret_cast< const DifShellBase* >( er.data() );
+ data[j + 0] = ds->getC( er );
+ data[j + 1] = ds->getCeq( er );
+ data[j + 2] = ds->getD( er );
+ data[j + 3] = ds->getValence( er );
+ data[j + 4] = ds->getLeak( er );
+ data[j + 5] = ds->getShapeMode( er );
+ data[j + 6] = ds->getLength( er );
+ data[j + 7] = ds->getDiameter( er );
+ data[j + 8] = ds->getThickness( er );
+ data[j + 9] = ds->getVolume( er );
+ data[j + 10] = ds->getOuterArea( er );
+ data[j + 11] = ds->getInnerArea( er );
+ j += len;
+ }
+ orig->zombieSwap( zClass );
+ j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ DifShellBase* ds =
+ reinterpret_cast< DifShellBase* >( er.data() );
+ ds->vSetSolver(er,hsolve);
+ ds->setC(er, data[j+0]);
+ ds->setCeq(er, data[j + 1]);
+ ds->setD(er, data[j + 2]);
+ ds->setValence(er, data[j + 3]);
+ ds->setLeak(er, data[j + 4]);
+ ds->setShapeMode(er, data[j + 5]);
+ ds->setLength(er, data[j + 6]);
+ ds->setDiameter(er, data[j + 7]);
+ ds->setThickness(er, data[j + 8]);
+ ds->setVolume(er, data[j + 9]);
+ ds->setOuterArea(er, data[j + 10]);
+ ds->setInnerArea(er, data[j + 11]);
+ j += len; //??
+ }
+
+}
diff --git a/biophysics/DifShellBase.h b/biophysics/DifShellBase.h
new file mode 100644
index 00000000..db74c536
--- /dev/null
+++ b/biophysics/DifShellBase.h
@@ -0,0 +1,144 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani 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 _DIFSHELL_BASE_H
+#define _DIFSHELL_BASE_H
+/*This is base class for DifShell*/
+
+class DifShellBase
+{
+ public:
+ DifShellBase();
+ /////////////////////////////////////////////////////////////
+ // Dest functions
+ /////////////////////////////////////////////////////////////
+ void reinit( const Eref & e, ProcPtr p );
+ void process(const Eref & e, ProcPtr p );
+ void buffer( const Eref& e, double kf, double kb, double bFree, double bBound );
+ void fluxFromOut(const Eref& e, double outerC, double outerThickness );
+ void fluxFromIn(const Eref& e, double innerC, double innerThickness );
+ void influx(const Eref& e, double I );
+ void outflux(const Eref& e, double I );
+ void fInflux(const Eref& e, double I, double fraction );
+ void fOutflux(const Eref& e, double I, double fraction );
+ void storeInflux(const Eref& e, double flux );
+ void storeOutflux(const Eref& e, double flux );
+ void tauPump(const Eref& e, double kP, double Ceq );
+ void eqTauPump(const Eref& e, double kP );
+ void mmPump(const Eref& e, double vMax, double Kd );
+ void hillPump(const Eref& e, double vMax, double Kd, unsigned int hill);
+
+ virtual void vReinit(const Eref & e, ProcPtr p ) = 0;
+ virtual void vProcess(const Eref & e, ProcPtr p ) = 0;
+ virtual void vBuffer(const Eref& e, double kf, double kb, double bFree, double bBound ) = 0;
+ virtual void vFluxFromOut(const Eref& e, double outerC, double outerThickness ) = 0;
+ virtual void vFluxFromIn(const Eref& e, double innerC, double innerThickness ) = 0;
+ virtual void vInflux(const Eref& e, double I ) = 0;
+ virtual void vOutflux(const Eref& e, double I ) = 0;
+ virtual void vFInflux(const Eref& e, double I, double fraction ) = 0;
+ virtual void vFOutflux(const Eref& e, double I, double fraction ) = 0;
+ virtual void vStoreInflux(const Eref& e, double flux ) = 0;
+ virtual void vStoreOutflux(const Eref& e, double flux ) = 0;
+ virtual void vTauPump(const Eref& e, double kP, double Ceq ) = 0;
+ virtual void vEqTauPump(const Eref& e, double kP ) = 0;
+ virtual void vMMPump(const Eref& e, double vMax, double Kd ) = 0;
+ virtual void vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill ) = 0;
+ /////////////////////////////////////////////////////////////
+ // Field access functions
+ /////////////////////////////////////////////////////////////
+
+ void setC(const Eref& e,double C);
+ double getC( const Eref& e) const;
+
+ void setCeq(const Eref& e,double Ceq );
+ double getCeq(const Eref& e) const;
+
+ void setD(const Eref& e, double D );
+ double getD(const Eref& e) const;
+
+ void setValence(const Eref& e, double valence );
+ double getValence(const Eref& e) const;
+
+ void setLeak(const Eref& e, double leak );
+ double getLeak(const Eref& e) const;
+
+ void setShapeMode(const Eref& e, unsigned int shapeMode );
+ unsigned int getShapeMode(const Eref& e) const;
+
+ void setLength(const Eref& e, double length );
+ double getLength(const Eref& e) const;
+
+ void setDiameter(const Eref& e, double diameter );
+ double getDiameter(const Eref& e) const;
+
+ void setThickness(const Eref& e, double thickness );
+ double getThickness(const Eref& e) const;
+
+ void setVolume(const Eref& e, double volume );
+ double getVolume(const Eref& e) const;
+
+ void setOuterArea(const Eref& e, double outerArea );
+ double getOuterArea(const Eref& e) const;
+
+ void setInnerArea(const Eref& e, double innerArea );
+ double getInnerArea(const Eref& e) const;
+
+ virtual void vSetC(const Eref& e,double C) = 0;
+ virtual double vGetC( const Eref& e) const = 0;
+
+ virtual void vSetCeq(const Eref& e,double Ceq ) = 0;
+ virtual double vGetCeq(const Eref& e) const = 0;
+
+ virtual void vSetD(const Eref& e, double D ) = 0;
+ virtual double vGetD(const Eref& e) const = 0;
+
+ virtual void vSetValence(const Eref& e, double valence ) = 0;
+ virtual double vGetValence(const Eref& e) const = 0;
+
+ virtual void vSetLeak(const Eref& e, double leak ) = 0;
+ virtual double vGetLeak(const Eref& e) const = 0;
+
+ virtual void vSetShapeMode(const Eref& e, unsigned int shapeMode ) = 0;
+ virtual unsigned int vGetShapeMode(const Eref& e) const = 0;
+
+ virtual void vSetLength(const Eref& e, double length ) = 0;
+ virtual double vGetLength(const Eref& e) const = 0;
+
+ virtual void vSetDiameter(const Eref& e, double diameter ) = 0;
+ virtual double vGetDiameter(const Eref& e) const = 0;
+
+ virtual void vSetThickness(const Eref& e, double thickness ) = 0;
+ virtual double vGetThickness(const Eref& e) const = 0;
+
+ virtual void vSetVolume(const Eref& e, double volume ) = 0;
+ virtual double vGetVolume(const Eref& e) const = 0;
+
+ virtual void vSetOuterArea(const Eref& e, double outerArea ) = 0;
+ virtual double vGetOuterArea(const Eref& e) const = 0;
+
+ virtual void vSetInnerArea(const Eref& e, double innerArea ) = 0;
+ virtual double vGetInnerArea(const Eref& e) const = 0;
+
+
+ virtual void vSetSolver( const Eref& e, Id hsolve );
+ static void zombify( Element* orig, const Cinfo* zClass,
+ Id hsolve );
+ static SrcFinfo1< double >* concentrationOut();
+ static SrcFinfo2< double, double >* innerDifSourceOut();
+ static SrcFinfo2< double, double >* outerDifSourceOut();
+ static const Cinfo * initCinfo();
+
+ private:
+
+
+
+};
+
+#endif //_DIFSHELL_BASE_H
diff --git a/biophysics/HHGate.cpp b/biophysics/HHGate.cpp
index 1f975b66..db9042d2 100644
--- a/biophysics/HHGate.cpp
+++ b/biophysics/HHGate.cpp
@@ -1,11 +1,11 @@
/**********************************************************************
-** 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.
-**********************************************************************/
+ ** 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 "header.h"
#include "ElementValueFinfo.h"
@@ -15,209 +15,209 @@ static const double SINGULARITY = 1.0e-6;
const Cinfo* HHGate::initCinfo()
{
- ///////////////////////////////////////////////////////
- // Field definitions.
- ///////////////////////////////////////////////////////
- static ReadOnlyLookupValueFinfo< HHGate, double, double > A( "A",
- "lookupA: Look up the A gate value from a double. Usually 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 dV fields.",
- &HHGate::lookupA );
- static ReadOnlyLookupValueFinfo< HHGate, double, double > B( "B",
- "lookupB: Look up the B gate value from a double."
- "Note that this looks up the raw tables, which are transformed"
- "from the reference parameters.",
- &HHGate::lookupB );
-
- static ElementValueFinfo< HHGate, vector< double > > alpha( "alpha",
- "Parameters for voltage-dependent rates, alpha:"
- "Set up alpha term using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setAlpha,
- &HHGate::getAlpha
- );
-
- static ElementValueFinfo< HHGate, vector< double > > beta( "beta",
- "Parameters for voltage-dependent rates, beta:"
- "Set up beta term using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setBeta,
- &HHGate::getBeta
- );
-
- static ElementValueFinfo< HHGate, vector< double > > tau( "tau",
- "Parameters for voltage-dependent rates, tau:"
- "Set up tau curve using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))",
- &HHGate::setTau,
- &HHGate::getTau
- );
-
- static ElementValueFinfo< HHGate, vector< double > > mInfinity(
- "mInfinity",
- "Parameters for voltage-dependent rates, mInfinity:"
- "Set up mInfinity curve using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setMinfinity,
- &HHGate::getMinfinity
- );
-
- static ElementValueFinfo< HHGate, double > min( "min",
- "Minimum range for lookup",
- &HHGate::setMin,
- &HHGate::getMin
- );
-
- static ElementValueFinfo< HHGate, double > max( "max",
- "Minimum range for lookup",
- &HHGate::setMax,
- &HHGate::getMax
- );
-
- static ElementValueFinfo< HHGate, unsigned int > divs( "divs",
- "Divisions for lookup. Zero means to use linear interpolation",
- &HHGate::setDivs,
- &HHGate::getDivs
- );
-
- static ElementValueFinfo< HHGate, vector< double > > tableA(
- "tableA",
- "Table of A entries",
- &HHGate::setTableA,
- &HHGate::getTableA
- );
-
- static ElementValueFinfo< HHGate, vector< double > > tableB(
- "tableB",
- "Table of alpha + beta entries",
- &HHGate::setTableB,
- &HHGate::getTableB
- );
-
- static ElementValueFinfo< HHGate, bool > useInterpolation(
- "useInterpolation",
- "Flag: use linear interpolation if true, else direct lookup",
- &HHGate::setUseInterpolation,
- &HHGate::getUseInterpolation
- );
-
- static ElementValueFinfo< HHGate, vector< double > > alphaParms(
- "alphaParms",
- "Set up both gates using 13 parameters, as follows:"
- "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
- "Here AA-AF are Coefficients A to F of the alpha (forward) term"
- "Here BA-BF are Coefficients A to F of the beta (reverse) term"
- "Here xdivs is the number of entries in the table,"
- "xmin and xmax define the range for lookup."
- "Outside this range the returned value will be the low [high]"
- "entry of the table."
- "The equation describing each table is:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setupAlpha,
- &HHGate::getAlphaParms
- );
-
- ///////////////////////////////////////////////////////
- // DestFinfos
- ///////////////////////////////////////////////////////
- static DestFinfo setupAlpha( "setupAlpha",
- "Set up both gates using 13 parameters, as follows:"
- "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
- "Here AA-AF are Coefficients A to F of the alpha (forward) term"
- "Here BA-BF are Coefficients A to F of the beta (reverse) term"
- "Here xdivs is the number of entries in the table,"
- "xmin and xmax define the range for lookup."
- "Outside this range the returned value will be the low [high]"
- "entry of the table."
- "The equation describing each table is:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- new EpFunc1< HHGate, vector< double > >( &HHGate::setupAlpha )
- );
- static DestFinfo setupTau( "setupTau",
- "Identical to setupAlpha, except that the forms specified by"
- "the 13 parameters are for the tau and m-infinity curves rather"
- "than the alpha and beta terms. So the parameters are:"
- "setupTau TA TB TC TD TF MA MB MC MD MF xdivs xmin xmax"
- "As before, the equation describing each curve is:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))",
- new EpFunc1< HHGate, vector< double > >( &HHGate::setupTau )
- );
- static DestFinfo tweakAlpha( "tweakAlpha",
- "Dummy function for backward compatibility. It used to convert"
- "the tables from alpha, beta values to alpha, alpha+beta"
- "because the internal calculations used these forms. Not"
- "needed now, deprecated.",
- new OpFunc0< HHGate >( &HHGate::tweakAlpha )
- );
- static DestFinfo tweakTau( "tweakTau",
- "Dummy function for backward compatibility. It used to convert"
- "the tables from tau, minf values to alpha, alpha+beta"
- "because the internal calculations used these forms. Not"
- "needed now, deprecated.",
- new OpFunc0< HHGate >( &HHGate::tweakTau )
- );
- static DestFinfo setupGate( "setupGate",
- "Sets up one gate at a time using the alpha/beta form."
- "Has 9 parameters, as follows:"
- "setupGate A B C D F xdivs xmin xmax is_beta"
- "This sets up the gate using the equation:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "Deprecated.",
- new EpFunc1< HHGate, vector< double > >( &HHGate::setupGate )
- );
- static Finfo* HHGateFinfos[] =
- {
- &A, // ReadOnlyLookupValue
- &B, // ReadOnlyLookupValue
- &alpha, // ElementValue
- &beta, // ElementValue
- &tau, // ElementValue
- &mInfinity, // ElementValue
- &min, // ElementValue
- &max, // ElementValue
- &divs, // ElementValue
- &tableA, // ElementValue
- &tableB, // ElementValue
- &useInterpolation, // ElementValue
- &alphaParms, // ElementValue
- &setupAlpha, // Dest
- &setupTau, // Dest
- &tweakAlpha, // Dest
- &tweakTau, // Dest
- &setupGate, // Dest
- };
-
- static string doc[] =
- {
- "Name", "HHGate",
- "Author", "Upinder S. Bhalla, 2011, NCBS",
- "Description", "HHGate: 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.",
- };
-
- static Dinfo< HHGate > dinfo;
- static Cinfo HHGateCinfo(
- "HHGate",
- Neutral::initCinfo(),
- HHGateFinfos, sizeof(HHGateFinfos)/sizeof(Finfo *),
- &dinfo,
- doc,
- sizeof(doc)/sizeof(string)
- );
-
- return &HHGateCinfo;
+ ///////////////////////////////////////////////////////
+ // Field definitions.
+ ///////////////////////////////////////////////////////
+ static ReadOnlyLookupValueFinfo< HHGate, double, double > A( "A",
+ "lookupA: Look up the A gate value from a double. Usually 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 dV fields.",
+ &HHGate::lookupA );
+ static ReadOnlyLookupValueFinfo< HHGate, double, double > B( "B",
+ "lookupB: Look up the B gate value from a double."
+ "Note that this looks up the raw tables, which are transformed"
+ "from the reference parameters.",
+ &HHGate::lookupB );
+
+ static ElementValueFinfo< HHGate, vector< double > > alpha( "alpha",
+ "Parameters for voltage-dependent rates, alpha:"
+ "Set up alpha term using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setAlpha,
+ &HHGate::getAlpha
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > beta( "beta",
+ "Parameters for voltage-dependent rates, beta:"
+ "Set up beta term using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setBeta,
+ &HHGate::getBeta
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > tau( "tau",
+ "Parameters for voltage-dependent rates, tau:"
+ "Set up tau curve using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))",
+ &HHGate::setTau,
+ &HHGate::getTau
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > mInfinity(
+ "mInfinity",
+ "Parameters for voltage-dependent rates, mInfinity:"
+ "Set up mInfinity curve using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setMinfinity,
+ &HHGate::getMinfinity
+ );
+
+ static ElementValueFinfo< HHGate, double > min( "min",
+ "Minimum range for lookup",
+ &HHGate::setMin,
+ &HHGate::getMin
+ );
+
+ static ElementValueFinfo< HHGate, double > max( "max",
+ "Minimum range for lookup",
+ &HHGate::setMax,
+ &HHGate::getMax
+ );
+
+ static ElementValueFinfo< HHGate, unsigned int > divs( "divs",
+ "Divisions for lookup. Zero means to use linear interpolation",
+ &HHGate::setDivs,
+ &HHGate::getDivs
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > tableA(
+ "tableA",
+ "Table of A entries",
+ &HHGate::setTableA,
+ &HHGate::getTableA
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > tableB(
+ "tableB",
+ "Table of alpha + beta entries",
+ &HHGate::setTableB,
+ &HHGate::getTableB
+ );
+
+ static ElementValueFinfo< HHGate, bool > useInterpolation(
+ "useInterpolation",
+ "Flag: use linear interpolation if true, else direct lookup",
+ &HHGate::setUseInterpolation,
+ &HHGate::getUseInterpolation
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > alphaParms(
+ "alphaParms",
+ "Set up both gates using 13 parameters, as follows:"
+ "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
+ "Here AA-AF are Coefficients A to F of the alpha (forward) term"
+ "Here BA-BF are Coefficients A to F of the beta (reverse) term"
+ "Here xdivs is the number of entries in the table,"
+ "xmin and xmax define the range for lookup."
+ "Outside this range the returned value will be the low [high]"
+ "entry of the table."
+ "The equation describing each table is:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setupAlpha,
+ &HHGate::getAlphaParms
+ );
+
+ ///////////////////////////////////////////////////////
+ // DestFinfos
+ ///////////////////////////////////////////////////////
+ static DestFinfo setupAlpha( "setupAlpha",
+ "Set up both gates using 13 parameters, as follows:"
+ "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
+ "Here AA-AF are Coefficients A to F of the alpha (forward) term"
+ "Here BA-BF are Coefficients A to F of the beta (reverse) term"
+ "Here xdivs is the number of entries in the table,"
+ "xmin and xmax define the range for lookup."
+ "Outside this range the returned value will be the low [high]"
+ "entry of the table."
+ "The equation describing each table is:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ new EpFunc1< HHGate, vector< double > >( &HHGate::setupAlpha )
+ );
+ static DestFinfo setupTau( "setupTau",
+ "Identical to setupAlpha, except that the forms specified by"
+ "the 13 parameters are for the tau and m-infinity curves rather"
+ "than the alpha and beta terms. So the parameters are:"
+ "setupTau TA TB TC TD TF MA MB MC MD MF xdivs xmin xmax"
+ "As before, the equation describing each curve is:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))",
+ new EpFunc1< HHGate, vector< double > >( &HHGate::setupTau )
+ );
+ static DestFinfo tweakAlpha( "tweakAlpha",
+ "Dummy function for backward compatibility. It used to convert"
+ "the tables from alpha, beta values to alpha, alpha+beta"
+ "because the internal calculations used these forms. Not"
+ "needed now, deprecated.",
+ new OpFunc0< HHGate >( &HHGate::tweakAlpha )
+ );
+ static DestFinfo tweakTau( "tweakTau",
+ "Dummy function for backward compatibility. It used to convert"
+ "the tables from tau, minf values to alpha, alpha+beta"
+ "because the internal calculations used these forms. Not"
+ "needed now, deprecated.",
+ new OpFunc0< HHGate >( &HHGate::tweakTau )
+ );
+ static DestFinfo setupGate( "setupGate",
+ "Sets up one gate at a time using the alpha/beta form."
+ "Has 9 parameters, as follows:"
+ "setupGate A B C D F xdivs xmin xmax is_beta"
+ "This sets up the gate using the equation:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "Deprecated.",
+ new EpFunc1< HHGate, vector< double > >( &HHGate::setupGate )
+ );
+ static Finfo* HHGateFinfos[] =
+ {
+ &A, // ReadOnlyLookupValue
+ &B, // ReadOnlyLookupValue
+ &alpha, // ElementValue
+ &beta, // ElementValue
+ &tau, // ElementValue
+ &mInfinity, // ElementValue
+ &min, // ElementValue
+ &max, // ElementValue
+ &divs, // ElementValue
+ &tableA, // ElementValue
+ &tableB, // ElementValue
+ &useInterpolation, // ElementValue
+ &alphaParms, // ElementValue
+ &setupAlpha, // Dest
+ &setupTau, // Dest
+ &tweakAlpha, // Dest
+ &tweakTau, // Dest
+ &setupGate, // Dest
+ };
+
+ static string doc[] =
+ {
+ "Name", "HHGate",
+ "Author", "Upinder S. Bhalla, 2011, NCBS",
+ "Description", "HHGate: 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.",
+ };
+
+ static Dinfo< HHGate > dinfo;
+ static Cinfo HHGateCinfo(
+ "HHGate",
+ Neutral::initCinfo(),
+ HHGateFinfos, sizeof(HHGateFinfos)/sizeof(Finfo *),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string)
+ );
+
+ return &HHGateCinfo;
}
static const Cinfo* hhGateCinfo = HHGate::initCinfo();
@@ -225,22 +225,22 @@ static const Cinfo* hhGateCinfo = HHGate::initCinfo();
// Core class functions
///////////////////////////////////////////////////
HHGate::HHGate()
- : xmin_(0), xmax_(1), invDx_(1),
- originalChanId_(0),
- originalGateId_(0),
- lookupByInterpolation_(0),
- isDirectTable_(0)
+ : xmin_(0), xmax_(1), invDx_(1),
+ originalChanId_(0),
+ originalGateId_(0),
+ lookupByInterpolation_(0),
+ isDirectTable_(0)
{;}
HHGate::HHGate( Id originalChanId, Id originalGateId )
- :
- A_( 1, 0.0 ),
- B_( 1, 0.0 ),
- xmin_(0), xmax_(1), invDx_(1),
- originalChanId_( originalChanId ),
- originalGateId_( originalGateId ),
- lookupByInterpolation_(0),
- isDirectTable_(0)
+ :
+ A_( 1, 0.0 ),
+ B_( 1, 0.0 ),
+ xmin_(0), xmax_(1), invDx_(1),
+ originalChanId_( originalChanId ),
+ originalGateId_( originalGateId ),
+ lookupByInterpolation_(0),
+ isDirectTable_(0)
{;}
///////////////////////////////////////////////////
@@ -249,272 +249,272 @@ HHGate::HHGate( Id originalChanId, Id originalGateId )
double HHGate::lookupTable( const vector< double >& tab, double v ) const
{
- if ( v <= xmin_ ) return tab[0];
- if ( v >= xmax_ ) return tab.back();
- if ( lookupByInterpolation_ ) {
- unsigned int index =
- static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
- assert( tab.size() > index );
- double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
- return tab[ index ] * ( 1 - frac ) + tab[ index + 1 ] * frac;
- } else {
- return tab[ static_cast< unsigned int >( (v - xmin_) * invDx_ ) ];
- }
+ if ( v <= xmin_ ) return tab[0];
+ if ( v >= xmax_ ) return tab.back();
+ if ( lookupByInterpolation_ ) {
+ unsigned int index =
+ static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
+ assert( tab.size() > index );
+ double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
+ return tab[ index ] * ( 1 - frac ) + tab[ index + 1 ] * frac;
+ } else {
+ return tab[ static_cast< unsigned int >( (v - xmin_) * invDx_ ) ];
+ }
}
double HHGate::lookupA( double v ) const
{
- return lookupTable( A_, v );
+ return lookupTable( A_, v );
}
double HHGate::lookupB( double v ) const
{
- return lookupTable( B_, v );
+ return lookupTable( B_, v );
}
void HHGate::lookupBoth( double v, double* A, double* B ) const
{
- if ( v <= xmin_ ) {
- *A = A_[0];
- *B = B_[0];
- } else if ( v >= xmax_ ) {
- *A = A_.back();
- *B = B_.back();
- } else {
- unsigned int index =
- static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
- assert( A_.size() > index && B_.size() > index );
- if ( lookupByInterpolation_ ) {
- double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
- *A = A_[ index ] * ( 1 - frac ) + A_[ index + 1 ] * frac;
- *B = B_[ index ] * ( 1 - frac ) + B_[ index + 1 ] * frac;
- } else {
- *A = A_[ index ];
- *B = B_[ index ];
- }
- }
+ if ( v <= xmin_ ) {
+ *A = A_[0];
+ *B = B_[0];
+ } else if ( v >= xmax_ ) {
+ *A = A_.back();
+ *B = B_.back();
+ } else {
+ unsigned int index =
+ static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
+ assert( A_.size() > index && B_.size() > index );
+ if ( lookupByInterpolation_ ) {
+ double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
+ *A = A_[ index ] * ( 1 - frac ) + A_[ index + 1 ] * frac;
+ *B = B_[ index ] * ( 1 - frac ) + B_[ index + 1 ] * frac;
+ } else {
+ *A = A_[ index ];
+ *B = B_[ index ];
+ }
+ }
}
vector< double > HHGate::getAlpha( const Eref& e) const
{
- return alpha_;
+ return alpha_;
}
void HHGate::setAlpha( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setAlpha on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "alpha" ) ) {
- alpha_ = val;
- updateTauMinf();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setAlpha on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "alpha" ) ) {
+ alpha_ = val;
+ updateTauMinf();
+ updateTables();
+ }
}
vector< double > HHGate::getBeta( const Eref& e) const
{
- return beta_;
+ return beta_;
}
void HHGate::setBeta( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setBeta on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "beta" ) ) {
- beta_ = val;
- updateTauMinf();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setBeta on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "beta" ) ) {
+ beta_ = val;
+ updateTauMinf();
+ updateTables();
+ }
}
vector< double > HHGate::getTau( const Eref& e) const
{
- return tau_;
+ return tau_;
}
void HHGate::setTau( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setTau on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "tau" ) ) {
- tau_ = val;
- updateAlphaBeta();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setTau on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "tau" ) ) {
+ tau_ = val;
+ updateAlphaBeta();
+ updateTables();
+ }
}
vector< double > HHGate::getMinfinity( const Eref& e) const
{
- return mInfinity_;
+ return mInfinity_;
}
void HHGate::setMinfinity( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setMinfinity on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "mInfinity" ) ) {
- mInfinity_ = val;
- updateAlphaBeta();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setMinfinity on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "mInfinity" ) ) {
+ mInfinity_ = val;
+ updateAlphaBeta();
+ updateTables();
+ }
}
double HHGate::getMin( const Eref& e) const
{
- return xmin_;
+ return xmin_;
}
void HHGate::setMin( const Eref& e, double val )
{
- if ( checkOriginal( e.id(), "min" ) ) {
- xmin_ = val;
- unsigned int xdivs = A_.size() - 1;
- if ( isDirectTable_ && xdivs > 0 ) {
- // Stuff here to stretch out table using interpolation.
- invDx_ = static_cast< double >( xdivs ) / ( xmax_ - val );
- tabFill( A_, xdivs, val, xmax_ );
- tabFill( B_, xdivs, val, xmax_ );
- } else {
- updateTables();
- }
- }
+ if ( checkOriginal( e.id(), "min" ) ) {
+ xmin_ = val;
+ unsigned int xdivs = A_.size() - 1;
+ if ( isDirectTable_ && xdivs > 0 ) {
+ // Stuff here to stretch out table using interpolation.
+ invDx_ = static_cast< double >( xdivs ) / ( xmax_ - val );
+ tabFill( A_, xdivs, val, xmax_ );
+ tabFill( B_, xdivs, val, xmax_ );
+ } else {
+ updateTables();
+ }
+ }
}
double HHGate::getMax( const Eref& e) const
{
- return xmax_;
+ return xmax_;
}
void HHGate::setMax( const Eref& e, double val )
{
- if ( checkOriginal( e.id(), "max" ) ) {
- xmax_ = val;
- unsigned int xdivs = A_.size() - 1;
- if ( isDirectTable_ && xdivs > 0 ) {
- // Set up using direct assignment of table values.
- invDx_ = static_cast< double >( xdivs ) / ( val - xmin_ );
- tabFill( A_, xdivs, xmin_, val );
- tabFill( B_, xdivs, xmin_, val );
- } else {
- // Set up using functional form. here we just recalculate.
- updateTables();
- }
- }
+ if ( checkOriginal( e.id(), "max" ) ) {
+ xmax_ = val;
+ unsigned int xdivs = A_.size() - 1;
+ if ( isDirectTable_ && xdivs > 0 ) {
+ // Set up using direct assignment of table values.
+ invDx_ = static_cast< double >( xdivs ) / ( val - xmin_ );
+ tabFill( A_, xdivs, xmin_, val );
+ tabFill( B_, xdivs, xmin_, val );
+ } else {
+ // Set up using functional form. here we just recalculate.
+ updateTables();
+ }
+ }
}
unsigned int HHGate::getDivs( const Eref& e) const
{
- return A_.size() - 1;
+ return A_.size() - 1;
}
void HHGate::setDivs( const Eref& e, unsigned int val )
{
- if ( checkOriginal( e.id(), "divs" ) ) {
- if ( isDirectTable_ ) {
- invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
- tabFill( A_, val, xmin_, xmax_ );
- tabFill( B_, val, xmin_, xmax_ );
- } else {
- /// Stuff here to redo sizes.
- A_.resize( val + 1 );
- B_.resize( val + 1 );
- invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
- updateTables();
- }
- }
+ if ( checkOriginal( e.id(), "divs" ) ) {
+ if ( isDirectTable_ ) {
+ invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
+ tabFill( A_, val, xmin_, xmax_ );
+ tabFill( B_, val, xmin_, xmax_ );
+ } else {
+ /// Stuff here to redo sizes.
+ A_.resize( val + 1 );
+ B_.resize( val + 1 );
+ invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
+ updateTables();
+ }
+ }
}
vector< double > HHGate::getTableA( const Eref& e) const
{
- return A_;
+ return A_;
}
void HHGate::setTableA( const Eref& e, vector< double > v )
{
- if ( v.size() < 2 ) {
- cout << "Warning: HHGate::setTableA: size must be >= 2 entries on "
- << e.id().path() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "tableA" ) ) {
- isDirectTable_ = 1;
- A_ = v;
- unsigned int xdivs = A_.size() - 1;
- invDx_ = static_cast< double >( xdivs ) / ( xmax_ - xmin_ );
- }
+ if ( v.size() < 2 ) {
+ cout << "Warning: HHGate::setTableA: size must be >= 2 entries on "
+ << e.id().path() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "tableA" ) ) {
+ isDirectTable_ = 1;
+ A_ = v;
+ unsigned int xdivs = A_.size() - 1;
+ invDx_ = static_cast< double >( xdivs ) / ( xmax_ - xmin_ );
+ }
}
vector< double > HHGate::getTableB( const Eref& e) const
{
- return B_;
+ return B_;
}
void HHGate::setTableB( const Eref& e, vector< double > v )
{
- if ( checkOriginal( e.id(), "tableB" ) ) {
- isDirectTable_ = 1;
- if ( A_.size() != v.size() ) {
- cout << "Warning: HHGate::setTableB: size should be same as table A: " << v.size() << " != " << A_.size() << ". Ignoring.\n";
- return;
- }
- B_ = v;
- }
+ if ( checkOriginal( e.id(), "tableB" ) ) {
+ isDirectTable_ = 1;
+ if ( A_.size() != v.size() ) {
+ cout << "Warning: HHGate::setTableB: size should be same as table A: " << v.size() << " != " << A_.size() << ". Ignoring.\n";
+ return;
+ }
+ B_ = v;
+ }
}
bool HHGate::getUseInterpolation( const Eref& e) const
{
- return lookupByInterpolation_;
+ return lookupByInterpolation_;
}
void HHGate::setUseInterpolation( const Eref& e, bool val )
{
- if ( checkOriginal( e.id(), "useInterpolation" ) )
- lookupByInterpolation_ = val;
+ if ( checkOriginal( e.id(), "useInterpolation" ) )
+ lookupByInterpolation_ = val;
}
void HHGate::setupAlpha( const Eref& e,
- vector< double > parms )
-{
- if ( checkOriginal( e.id(), "setupAlpha" ) ) {
- if ( parms.size() != 13 ) {
- cout << "HHGate::setupAlpha: Error: parms.size() != 13\n";
- return;
- }
- setupTables( parms, false );
- alpha_.resize( 5, 0 );
- beta_.resize( 5, 0 );
- for ( unsigned int i = 0; i < 5; ++i )
- alpha_[i] = parms[i];
- for ( unsigned int i = 5; i < 10; ++i )
- beta_[i - 5] = parms[i];
- }
+ vector< double > parms )
+{
+ if ( checkOriginal( e.id(), "setupAlpha" ) ) {
+ if ( parms.size() != 13 ) {
+ cout << "HHGate::setupAlpha: Error: parms.size() != 13\n";
+ return;
+ }
+ setupTables( parms, false );
+ alpha_.resize( 5, 0 );
+ beta_.resize( 5, 0 );
+ for ( unsigned int i = 0; i < 5; ++i )
+ alpha_[i] = parms[i];
+ for ( unsigned int i = 5; i < 10; ++i )
+ beta_[i - 5] = parms[i];
+ }
}
vector< double > HHGate::getAlphaParms( const Eref& e ) const
{
- vector< double > ret = alpha_;
- ret.insert( ret.end(), beta_.begin(), beta_.end() );
- ret.push_back( A_.size() );
- ret.push_back( xmin_ );
- ret.push_back( xmax_ );
+ vector< double > ret = alpha_;
+ ret.insert( ret.end(), beta_.begin(), beta_.end() );
+ ret.push_back( A_.size() );
+ ret.push_back( xmin_ );
+ ret.push_back( xmax_ );
- return ret;
+ return ret;
}
///////////////////////////////////////////////////
@@ -522,25 +522,25 @@ vector< double > HHGate::getAlphaParms( const Eref& e ) const
///////////////////////////////////////////////////
void HHGate::setupTau( const Eref& e,
- vector< double > parms )
+ vector< double > parms )
{
- if ( checkOriginal( e.id(), "setupTau" ) ) {
- if ( parms.size() != 13 ) {
- cout << "HHGate::setupTau: Error: parms.size() != 13\n";
- return;
- }
- setupTables( parms, true );
- }
+ if ( checkOriginal( e.id(), "setupTau" ) ) {
+ if ( parms.size() != 13 ) {
+ cout << "HHGate::setupTau: Error: parms.size() != 13\n";
+ return;
+ }
+ setupTables( parms, true );
+ }
}
void HHGate::tweakAlpha()
{
- ; // Dummy
+ ; // Dummy
}
void HHGate::tweakTau()
{
- ; // Dummy
+ ; // Dummy
}
/**
@@ -550,96 +550,96 @@ void HHGate::tweakTau()
*/
void HHGate::setupTables( const vector< double >& parms, bool doTau )
{
- assert( parms.size() == 13 );
- static const int XDIVS = 10;
- static const int XMIN = 11;
- static const int XMAX = 12;
- if ( parms[XDIVS] < 1 ) return;
- unsigned int xdivs = static_cast< unsigned int >( parms[XDIVS] );
-
- A_.resize( xdivs + 1 );
- B_.resize( xdivs + 1 );
- xmin_ = parms[XMIN];
- xmax_ = parms[XMAX];
- assert( xmax_ > xmin_ );
- invDx_ = xdivs / (xmax_ - xmin_ );
- double dx = ( xmax_ - xmin_ ) / xdivs;
-
- double x = xmin_;
- double prevAentry = 0.0;
- double prevBentry = 0.0;
- double temp;
- double temp2 = 0.0;
- unsigned int i;
-
- for( i = 0; i <= xdivs; i++ ) {
- if ( fabs( parms[4] ) < SINGULARITY ) {
- temp = 0.0;
- A_[i] = temp;
- } else {
- temp2 = parms[2] + exp( ( x + parms[3] ) / parms[4] );
- if ( fabs( temp2 ) < SINGULARITY ) {
- temp2 = parms[2] + exp( ( x + dx/10.0 + parms[3] ) / parms[4] );
- temp = ( parms[0] + parms[1] * (x + dx/10 ) ) / temp2;
-
- temp2 = parms[2] + exp( ( x - dx/10.0 + parms[3] ) / parms[4] );
- temp += ( parms[0] + parms[1] * (x - dx/10 ) ) / temp2;
- temp /= 2.0;
- // cout << "interpolated temp = " << temp << ", prev = " << prevAentry << endl;
-
- // temp = prevAentry;
- A_[i] = temp;
- } else {
- temp = ( parms[0] + parms[1] * x) / temp2;
- A_[i] = temp;
- }
- }
- if ( fabs( parms[9] ) < SINGULARITY ) {
- B_[i] = 0.0;
- } else {
- temp2 = parms[7] + exp( ( x + parms[8] ) / parms[9] );
- if ( fabs( temp2 ) < SINGULARITY ) {
- temp2 = parms[7] + exp( ( x + dx/10.0 + parms[8] ) / parms[9] );
- temp = (parms[5] + parms[6] * (x + dx/10) ) / temp2;
- temp2 = parms[7] + exp( ( x - dx/10.0 + parms[8] ) / parms[9] );
- temp += (parms[5] + parms[6] * (x - dx/10) ) / temp2;
- temp /= 2.0;
- B_[i] = temp;
- // B_[i] = prevBentry;
- } else {
- B_[i] = (parms[5] + parms[6] * x ) / temp2;
- // B_.table_[i] = ( parms[5] + parms[6] * x ) / temp2;
- }
- }
- // There are cleaner ways to do this, but this keeps
- // the relation to the GENESIS version clearer.
- // Note the additional SINGULARITY check, to fix a bug
- // in the earlier code.
- if ( doTau == 0 && fabs( temp2 ) > SINGULARITY )
- B_[i] += temp;
-
- prevAentry = A_[i];
- prevBentry = B_[i];
- x += dx;
- }
-
- prevAentry = 0.0;
- prevBentry = 0.0;
- if ( doTau ) {
- for( i = 0; i <= xdivs; i++ ) {
- temp = A_[i];
- temp2 = B_[i];
- if ( fabs( temp ) < SINGULARITY ) {
- A_[i] = prevAentry;
- B_[i] = prevBentry;
- } else {
- A_[i] = temp2 / temp;
- B_[i] = 1.0 / temp;
- }
- prevAentry = A_[i];
- prevBentry = B_[i];
- }
- }
+ assert( parms.size() == 13 );
+ static const int XDIVS = 10;
+ static const int XMIN = 11;
+ static const int XMAX = 12;
+ if ( parms[XDIVS] < 1 ) return;
+ unsigned int xdivs = static_cast< unsigned int >( parms[XDIVS] );
+
+ A_.resize( xdivs + 1 );
+ B_.resize( xdivs + 1 );
+ xmin_ = parms[XMIN];
+ xmax_ = parms[XMAX];
+ assert( xmax_ > xmin_ );
+ invDx_ = xdivs / (xmax_ - xmin_ );
+ double dx = ( xmax_ - xmin_ ) / xdivs;
+
+ double x = xmin_;
+ double prevAentry = 0.0;
+ double prevBentry = 0.0;
+ double temp;
+ double temp2 = 0.0;
+ unsigned int i;
+
+ for( i = 0; i <= xdivs; i++ ) {
+ if ( fabs( parms[4] ) < SINGULARITY ) {
+ temp = 0.0;
+ A_[i] = temp;
+ } else {
+ temp2 = parms[2] + exp( ( x + parms[3] ) / parms[4] );
+ if ( fabs( temp2 ) < SINGULARITY ) {
+ temp2 = parms[2] + exp( ( x + dx/10.0 + parms[3] ) / parms[4] );
+ temp = ( parms[0] + parms[1] * (x + dx/10 ) ) / temp2;
+
+ temp2 = parms[2] + exp( ( x - dx/10.0 + parms[3] ) / parms[4] );
+ temp += ( parms[0] + parms[1] * (x - dx/10 ) ) / temp2;
+ temp /= 2.0;
+ // cout << "interpolated temp = " << temp << ", prev = " << prevAentry << endl;
+
+ // temp = prevAentry;
+ A_[i] = temp;
+ } else {
+ temp = ( parms[0] + parms[1] * x) / temp2;
+ A_[i] = temp;
+ }
+ }
+ if ( fabs( parms[9] ) < SINGULARITY ) {
+ B_[i] = 0.0;
+ } else {
+ temp2 = parms[7] + exp( ( x + parms[8] ) / parms[9] );
+ if ( fabs( temp2 ) < SINGULARITY ) {
+ temp2 = parms[7] + exp( ( x + dx/10.0 + parms[8] ) / parms[9] );
+ temp = (parms[5] + parms[6] * (x + dx/10) ) / temp2;
+ temp2 = parms[7] + exp( ( x - dx/10.0 + parms[8] ) / parms[9] );
+ temp += (parms[5] + parms[6] * (x - dx/10) ) / temp2;
+ temp /= 2.0;
+ B_[i] = temp;
+ // B_[i] = prevBentry;
+ } else {
+ B_[i] = (parms[5] + parms[6] * x ) / temp2;
+ // B_.table_[i] = ( parms[5] + parms[6] * x ) / temp2;
+ }
+ }
+ // There are cleaner ways to do this, but this keeps
+ // the relation to the GENESIS version clearer.
+ // Note the additional SINGULARITY check, to fix a bug
+ // in the earlier code.
+ if ( doTau == 0 && fabs( temp2 ) > SINGULARITY )
+ B_[i] += temp;
+
+ prevAentry = A_[i];
+ prevBentry = B_[i];
+ x += dx;
+ }
+
+ prevAentry = 0.0;
+ prevBentry = 0.0;
+ if ( doTau ) {
+ for( i = 0; i <= xdivs; i++ ) {
+ temp = A_[i];
+ temp2 = B_[i];
+ if ( fabs( temp ) < SINGULARITY ) {
+ A_[i] = prevAentry;
+ B_[i] = prevBentry;
+ } else {
+ A_[i] = temp2 / temp;
+ B_[i] = 1.0 / temp;
+ }
+ prevAentry = A_[i];
+ prevBentry = B_[i];
+ }
+ }
}
/**
@@ -649,97 +649,97 @@ void HHGate::setupTables( const vector< double >& parms, bool doTau )
*/
void HHGate::tweakTables( bool doTau )
{
- unsigned int i;
- unsigned int size = A_.size();
- assert( size == B_.size() );
- if ( doTau ) {
- for ( i = 0; i < size; i++ ) {
- double temp = A_[ i ];
- double temp2 = B_[ i ];
- if ( fabs( temp ) < SINGULARITY ) {
- if ( temp < 0.0 )
- temp = -SINGULARITY;
- else
- temp = SINGULARITY;
- }
- A_[ i ] = temp2 / temp;
- B_[ i ] = 1.0 / temp;
- }
- } else {
- for ( i = 0; i < size; i++ )
- B_[i] = A_[i] + B_[i];
- }
+ unsigned int i;
+ unsigned int size = A_.size();
+ assert( size == B_.size() );
+ if ( doTau ) {
+ for ( i = 0; i < size; i++ ) {
+ double temp = A_[ i ];
+ double temp2 = B_[ i ];
+ if ( fabs( temp ) < SINGULARITY ) {
+ if ( temp < 0.0 )
+ temp = -SINGULARITY;
+ else
+ temp = SINGULARITY;
+ }
+ A_[ i ] = temp2 / temp;
+ B_[ i ] = 1.0 / temp;
+ }
+ } else {
+ for ( i = 0; i < size; i++ )
+ B_[i] = A_[i] + B_[i];
+ }
}
void HHGate::setupGate( const Eref& e,
- vector< double > parms )
-{
- // The nine arguments are :
- // A B C D F size min max isbeta
- // If size == 0 then we check that the gate has already been allocated.
- // If isbeta is true then we also have to do the conversion to
- // HHGate form of alpha, alpha+beta, assuming that the alpha gate
- // has already been setup. This uses tweakTables.
- // We may need to resize the tables if they don't match here.
- if ( !checkOriginal( e.id(), "setupGate" ) )
- return;
-
- if ( parms.size() != 9 ) {
- cout << "HHGate::setupGate: Error: parms.size() != 9\n";
- return;
- }
-
- double A = parms[0];
- double B = parms[1];
- double C = parms[2];
- double D = parms[3];
- double F = parms[4];
- int size = static_cast< int > (parms[5] );
- double min = parms[6];
- double max = parms[7];
- bool isBeta = static_cast< bool >( parms[8] );
-
- vector< double >& ip = ( isBeta ) ? B_ : A_;
-
- if ( size <= 0 ) { // Look up size, min, max from the interpol
- size = ip.size() - 1;
- if ( size <= 0 ) {
- cout << "Error: setupGate has zero size\n";
- return;
- }
- } else {
- ip.resize( size + 1 );
- }
-
- double dx = ( max - min ) / static_cast< double >( size );
- double x = min + dx / 2.0;
- for ( int i = 0; i <= size; i++ ) {
- if ( fabs ( F ) < SINGULARITY ) {
- ip[i] = 0.0;
- } else {
- double temp2 = C + exp( ( x + D ) / F );
- if ( fabs( temp2 ) < SINGULARITY )
- ip[i] = ip[i-1];
- else
- ip[i] = ( A + B * x ) / temp2;
- }
- }
-
- if ( isBeta ) {
- assert( A_.size() > 0 );
- // Here we ensure that the tables are the same size
- if ( A_.size() != B_.size() ) {
- if ( A_.size() > B_.size() ) {
- // Note that the tabFill expects to allocate the
- // terminating entry, so we put in size - 1.
- tabFill( B_, A_.size() - 1, xmin_, xmax_ );
- } else {
- tabFill( A_, B_.size() - 1, xmin_, xmax_ );
- }
- }
- // Then we do the tweaking to convert to HHChannel form.
- tweakTables( 0 );
- }
+ vector< double > parms )
+{
+ // The nine arguments are :
+ // A B C D F size min max isbeta
+ // If size == 0 then we check that the gate has already been allocated.
+ // If isbeta is true then we also have to do the conversion to
+ // HHGate form of alpha, alpha+beta, assuming that the alpha gate
+ // has already been setup. This uses tweakTables.
+ // We may need to resize the tables if they don't match here.
+ if ( !checkOriginal( e.id(), "setupGate" ) )
+ return;
+
+ if ( parms.size() != 9 ) {
+ cout << "HHGate::setupGate: Error: parms.size() != 9\n";
+ return;
+ }
+
+ double A = parms[0];
+ double B = parms[1];
+ double C = parms[2];
+ double D = parms[3];
+ double F = parms[4];
+ int size = static_cast< int > (parms[5] );
+ double min = parms[6];
+ double max = parms[7];
+ bool isBeta = static_cast< bool >( parms[8] );
+
+ vector< double >& ip = ( isBeta ) ? B_ : A_;
+
+ if ( size <= 0 ) { // Look up size, min, max from the interpol
+ size = ip.size() - 1;
+ if ( size <= 0 ) {
+ cout << "Error: setupGate has zero size\n";
+ return;
+ }
+ } else {
+ ip.resize( size + 1 );
+ }
+
+ double dx = ( max - min ) / static_cast< double >( size );
+ double x = min + dx / 2.0;
+ for ( int i = 0; i <= size; i++ ) {
+ if ( fabs ( F ) < SINGULARITY ) {
+ ip[i] = 0.0;
+ } else {
+ double temp2 = C + exp( ( x + D ) / F );
+ if ( fabs( temp2 ) < SINGULARITY )
+ ip[i] = ip[i-1];
+ else
+ ip[i] = ( A + B * x ) / temp2;
+ }
+ }
+
+ if ( isBeta ) {
+ assert( A_.size() > 0 );
+ // Here we ensure that the tables are the same size
+ if ( A_.size() != B_.size() ) {
+ if ( A_.size() > B_.size() ) {
+ // Note that the tabFill expects to allocate the
+ // terminating entry, so we put in size - 1.
+ tabFill( B_, A_.size() - 1, xmin_, xmax_ );
+ } else {
+ tabFill( A_, B_.size() - 1, xmin_, xmax_ );
+ }
+ }
+ // Then we do the tweaking to convert to HHChannel form.
+ tweakTables( 0 );
+ }
}
///////////////////////////////////////////////////////////////////////
@@ -753,54 +753,54 @@ void HHGate::setupGate( const Eref& e,
* subdivisions that the table represents.
*/
void HHGate::tabFill( vector< double >& table,
- unsigned int newXdivs, double newXmin, double newXmax )
-{
- if ( newXdivs < 3 ) {
- cout << "Error: tabFill: # divs must be >= 3. Not filling table.\n";
- return;
- }
-
- vector< double > old = table;
- double newDx = ( newXmax - newXmin ) / newXdivs;
- table.resize( newXdivs + 1 );
- bool origLookupMode = lookupByInterpolation_;
- lookupByInterpolation_ = 1;
-
- for( unsigned int i = 0; i <= newXdivs; ++i ) {
- table[i] = lookupTable( table, newXmin + i * newDx );
- }
+ unsigned int newXdivs, double newXmin, double newXmax )
+{
+ if ( newXdivs < 3 ) {
+ cout << "Error: tabFill: # divs must be >= 3. Not filling table.\n";
+ return;
+ }
+
+ vector< double > old = table;
+ double newDx = ( newXmax - newXmin ) / newXdivs;
+ table.resize( newXdivs + 1 );
+ bool origLookupMode = lookupByInterpolation_;
+ lookupByInterpolation_ = 1;
+
+ for( unsigned int i = 0; i <= newXdivs; ++i ) {
+ table[i] = lookupTable( table, newXmin + i * newDx );
+ }
- lookupByInterpolation_ = origLookupMode;
+ lookupByInterpolation_ = origLookupMode;
}
bool HHGate::checkOriginal( Id id, const string& field ) const
{
- if ( id == originalGateId_ )
- return 1;
+ if ( id == originalGateId_ )
+ return 1;
- cout << "Warning: HHGate: attempt to set field '" << field << "' on " <<
- id.path() << "\nwhich is not the original Gate element. Ignored.\n";
- return 0;
+ cout << "Warning: HHGate: attempt to set field '" << field << "' on " <<
+ id.path() << "\nwhich is not the original Gate element. Ignored.\n";
+ return 0;
}
bool HHGate::isOriginalChannel( Id id ) const
{
- return ( id == originalChanId_ );
+ return ( id == originalChanId_ );
}
bool HHGate::isOriginalGate( Id id ) const
{
- return ( id == originalGateId_ );
+ return ( id == originalGateId_ );
}
Id HHGate::originalChannelId() const
{
- return originalChanId_;
+ return originalChanId_;
}
Id HHGate::originalGateId() const
{
- return originalGateId_;
+ return originalGateId_;
}
void HHGate::updateAlphaBeta()
@@ -813,13 +813,13 @@ void HHGate::updateTauMinf()
void HHGate::updateTables()
{
- if ( alpha_.size() == 0 || beta_.size() == 0 )
- return;
- vector< double > parms = alpha_;
- parms.insert( parms.end(), beta_.begin(), beta_.end() );
- parms.push_back( A_.size() );
- parms.push_back( xmin_ );
- parms.push_back( xmax_ );
+ if ( alpha_.size() == 0 || beta_.size() == 0 )
+ return;
+ vector< double > parms = alpha_;
+ parms.insert( parms.end(), beta_.begin(), beta_.end() );
+ parms.push_back( A_.size() );
+ parms.push_back( xmin_ );
+ parms.push_back( xmax_ );
- setupTables( parms, 0 );
+ setupTables( parms, 0 );
}
diff --git a/biophysics/MMPump.cpp b/biophysics/MMPump.cpp
new file mode 100644
index 00000000..f9ed90ef
--- /dev/null
+++ b/biophysics/MMPump.cpp
@@ -0,0 +1,141 @@
+#include "header.h"
+#include "MMPump.h"
+#include "ElementValueFinfo.h"
+#include "../utility/numutil.h"
+#include <cmath>
+
+SrcFinfo2< double,double >* MMPump::PumpOut()
+{
+ static SrcFinfo2< double,double > pumpOut( "PumpOut",
+ "Sends out MMPump parameters.");
+ return &pumpOut;
+}
+
+const Cinfo * MMPump::initCinfo()
+{
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< MMPump>(&MMPump::process ) );
+ static DestFinfo reinit( "reinit",
+ "Reinit happens only in stage 0",
+ new ProcOpFunc< MMPump>( &MMPump::reinit ));
+
+ static Finfo* processShared[] = {
+ &process, &reinit
+ };
+
+ static SharedFinfo proc(
+ "proc",
+ "Shared message to receive Process message from scheduler",
+ processShared, sizeof( processShared ) / sizeof( Finfo* ));
+
+
+ ////////////////////////////
+ // Field defs
+ ////////////////////////////
+
+
+
+ static ElementValueFinfo<MMPump, double> Vmax("Vmax",
+ "maximum pump velocity, scaled by mebrane"
+ "surface area. i.e., max ion flux in moles/sec",
+ &MMPump::setVmax,
+ &MMPump::getVmax);
+
+ static ElementValueFinfo<MMPump, double> Kd("Kd",
+ "half-maximal activating concentration in mM",
+ &MMPump::setKd,
+ &MMPump::getKd);
+
+ ////
+ // DestFinfo
+ ////
+ static Finfo * difMMPumpFinfos[] = {
+ //////////////////////////////////////////////////////////////////
+ // Field definitions
+ //////////////////////////////////////////////////////////////////
+
+
+ &Vmax,
+ &Kd,
+ //////////////////////////////////////////////////////////////////
+ // SharedFinfo definitions
+ /////////////////////////////////////////////////////////////////
+
+ &proc,
+ PumpOut(),
+ //////////////////////////////////////////////////////////////////
+ // DestFinfo definitions
+ //////////////////////////////////////////////////////////////////
+
+ };
+
+ static string doc[] = {
+ "Name", "MMPump",
+ "Author", "Subhasis Ray (ported from GENESIS2)",
+ "Description", "Models Michaelis-Menten pump. It is coupled with a DifShell.",
+ };
+ static ZeroSizeDinfo<int> dinfo;
+ static Cinfo MMPumpCinfo(
+ "MMPump",
+ Neutral::initCinfo(),
+ difMMPumpFinfos,
+ sizeof(difMMPumpFinfos)/sizeof(Finfo*),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string));
+
+ return &MMPumpCinfo;
+}
+
+static const Cinfo * MMpumpCinfo = MMPump::initCinfo();
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Class functions
+////////////////////////////////////////////////////////////////////////////////
+
+MMPump::MMPump()
+{ ; }
+
+
+
+
+double MMPump::getVmax(const Eref& e) const
+{
+
+ return Vmax_;
+}
+void MMPump::setVmax(const Eref& e,double value)
+{
+if ( value < 0.0 ) {
+ cerr << "Error: MMPump: Vmax cannot be negative!\n";
+ return;
+ }
+ Vmax_ = value;
+}
+
+
+double MMPump::getKd(const Eref& e) const
+{
+ return Kd_;
+}
+void MMPump::setKd(const Eref& e,double value)
+{
+if ( value < 0.0 ) {
+ cerr << "Error: MMPump: Kd cannot be negative!\n";
+ return;
+ }
+ Kd_ = value;
+}
+
+
+void MMPump::process(const Eref& e, ProcPtr p)
+{
+ PumpOut()->send(e,Vmax_,Kd_);
+}
+
+void MMPump::reinit(const Eref& e, ProcPtr p)
+{
+ return;
+}
diff --git a/biophysics/MMPump.h b/biophysics/MMPump.h
new file mode 100644
index 00000000..911fb1d4
--- /dev/null
+++ b/biophysics/MMPump.h
@@ -0,0 +1,45 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani 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 _MMPUMP_H
+#define _MMPUMP_H
+/*This is base class for MMPump*/
+class MMPump
+{
+public:
+ MMPump();
+ void process( const Eref & e, ProcPtr p);
+ void reinit(const Eref &e, ProcPtr p);
+
+ double getVmax(const Eref& e) const; //Vmax of the pump
+ void setVmax(const Eref& e,double value);
+
+ int getVal(const Eref& e) const; //Valence
+ void setVal(const Eref& e,int value);
+
+ double getKd(const Eref& e) const; // Pump's Kd
+ void setKd(const Eref& e,double value);
+
+
+
+ static SrcFinfo2< double,double >* PumpOut(); //Pump parameters;
+
+ static const Cinfo * initCinfo();
+
+private:
+
+ double Vmax_;
+ double Kd_;
+};
+
+
+
+
+#endif //_MMPUMP_BASE_H
diff --git a/biophysics/Makefile b/biophysics/Makefile
index ae8391ef..6659630c 100644
--- a/biophysics/Makefile
+++ b/biophysics/Makefile
@@ -38,7 +38,10 @@ OBJ = \
NMDAChan.o \
testBiophysics.o \
IzhikevichNrn.o \
+ DifShellBase.o \
DifShell.o \
+ DifBufferBase.o \
+ DifBuffer.o \
Leakage.o \
VectorTable.o \
MarkovRateTable.o \
@@ -50,6 +53,7 @@ OBJ = \
VClamp.o \
+ # MMpump.o \
# LeakyIaF.o \
# GapJunction.o \
# Nernst.o \
@@ -94,7 +98,11 @@ ReadCell.o: CompartmentBase.h Compartment.h SymCompartment.h ReadCell.h ../shell
SwcSegment.o: SwcSegment.h ../utility/Vec.h
ReadSwc.o: CompartmentBase.h Compartment.h SymCompartment.h SwcSegment.h ReadSwc.h ../shell/Shell.h ../utility/Vec.h
IzhikevichNrn.o: IzhikevichNrn.h
-DifShell.o: DifShell.h
+DifShellBase.o: DifShellBase.h
+DifShell.o: DifShellBase.h DifShell.h
+DifBufferBase.o: DifBufferBase.h
+DifBuffer.o: DifBufferBase.h DifBuffer.h
+MMPump.o: MMPump.h
testBiophysics.o: IntFire.h CompartmentBase.h Compartment.h HHChannel.h HHGate.h
VectorTable.o : VectorTable.h
MarkovGslSolver.o : MarkovGslSolver.h
diff --git a/python/moose/chemUtil/graphUtils.py b/python/moose/chemUtil/graphUtils.py
index 1d9ee1d5..5bf05f54 100644
--- a/python/moose/chemUtil/graphUtils.py
+++ b/python/moose/chemUtil/graphUtils.py
@@ -1,10 +1,20 @@
import moose
-import pygraphviz as pgv
-#import networkx as nx
+
+pygraphvizFound_ = True
+try:
+ import pygraphviz as pgv
+except Exception as e:
+ pygraphvizFound_ = False
def autoCoordinates(meshEntry,srcdesConnection):
+ global pygraphvizFound_
positionInfo = {}
+ if not pygraphvizFound_:
+ print( '[warn] python-pygraphviz could not be found.' )
+ print( '\tMOOSE Install pygraphviz to use this feature' )
+ return positionInfo
+
if meshEntry:
#G = nx.Graph()
G = pgv.AGraph()
@@ -67,7 +77,7 @@ def autoCoordinates(meshEntry,srcdesConnection):
ann.x = float(valuelist[0])
ann.y = float(valuelist[1])
- return(positionInfo)
+ return positionInfo
def find_index(value, key):
""" Value.get(key) to avoid expection which would raise if empty value in dictionary for a given key """
diff --git a/python/moose/moose.py b/python/moose/moose.py
index cda22e8f..904d6c32 100644
--- a/python/moose/moose.py
+++ b/python/moose/moose.py
@@ -33,7 +33,13 @@ except Exception as e:
print('\tError was %s' % e)
genesisSupport_ = False
-import chemUtil.add_Delete_ChemicalSolver
+chemUtilSupport_ = True
+try:
+ import chemUtil.add_Delete_ChemicalSolver
+except Exception as e:
+ chemUtilSupport_ = False
+ print( 'Failed to import utility module chemUtil' )
+ print( '\tError was %s' % e )
sequence_types = ['vector<double>',
'vector<int>',
diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp
index e48b6a34..434d2f62 100644
--- a/scheduling/Clock.cpp
+++ b/scheduling/Clock.cpp
@@ -358,6 +358,10 @@ const Cinfo* Clock::initCinfo()
" CaConc 1 50e-6\n"
" CaConcBase 1 50e-6\n"
" DifShell 1 50e-6\n"
+ " DifShellBase 1 50e-6\n"
+ " MMPump 1 50e-6\n"
+ " DifBuffer 1 50e-6\n"
+ " DifBufferBase 1 50e-6\n"
" MgBlock 1 50e-6\n"
" Nernst 1 50e-6\n"
" RandSpike 1 50e-6\n"
@@ -837,6 +841,10 @@ void Clock::buildDefaultTick()
defaultTick_["CaConc"] = 1;
defaultTick_["CaConcBase"] = 1;
defaultTick_["DifShell"] = 1;
+ defaultTick_["DifShellBase"] = 1;
+ defaultTick_["MMPump"] = 1;
+ defaultTick_["DifBuffer"] = 1;
+ defaultTick_["DifBufferBase"] = 1;
defaultTick_["MgBlock"] = 1;
defaultTick_["Nernst"] = 1;
defaultTick_["RandSpike"] = 1;
--
GitLab