From 7537d1c7e7ce3ccd343b1c9146c446d86d73866a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Asia=20J=C4=99drzejewska-Szmek?= <asia@in.waw.pl> Date: Fri, 20 Jan 2017 15:49:30 -0500 Subject: [PATCH] Rewrite DifShell. Inherit from DifShellBase Most of Cinfo code is moved to DifShellBase, so it got deleted from DifShell. Method's names are changed, because they are now inherited from DifShellBase. --- moose-core/biophysics/CMakeLists.txt | 4 +- moose-core/biophysics/DifBuffer.cpp | 441 ++----------------- moose-core/biophysics/DifBuffer.h | 77 +--- moose-core/biophysics/DifBufferBase.cpp | 435 +++++++++++++++++++ moose-core/biophysics/DifBufferBase.h | 108 +++++ moose-core/biophysics/DifShell.cpp | 547 ++---------------------- moose-core/biophysics/DifShell.h | 171 ++------ moose-core/biophysics/Makefile | 8 +- 8 files changed, 688 insertions(+), 1103 deletions(-) create mode 100644 moose-core/biophysics/DifBufferBase.cpp create mode 100644 moose-core/biophysics/DifBufferBase.h diff --git a/moose-core/biophysics/CMakeLists.txt b/moose-core/biophysics/CMakeLists.txt index 57b3ab45..133a5556 100644 --- a/moose-core/biophysics/CMakeLists.txt +++ b/moose-core/biophysics/CMakeLists.txt @@ -33,8 +33,10 @@ set(BIOPHYSICS_SRCS SynChan.cpp NMDAChan.cpp testBiophysics.cpp - IzhikevichNrn.cpp + IzhikevichNrn.cpp + DifShellBase.cpp DifShell.cpp + DifBufferBase.cpp DifBuffer.cpp Leakage.cpp VectorTable.cpp diff --git a/moose-core/biophysics/DifBuffer.cpp b/moose-core/biophysics/DifBuffer.cpp index c07bfa37..1d6ccf9d 100644 --- a/moose-core/biophysics/DifBuffer.cpp +++ b/moose-core/biophysics/DifBuffer.cpp @@ -25,7 +25,7 @@ // 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 @@ -50,227 +50,11 @@ #include "ElementValueFinfo.h" #include "../utility/numutil.h" #include <cmath> +const double DifBuffer::EPSILON = 1.0e-10; -SrcFinfo4< double, double, double, double >* DifBuffer::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 >* DifBuffer::innerDifSourceOut(){ - static SrcFinfo2< double, double > sourceOut("innerDifSourceOut", - "Sends out source information."); - return &sourceOut; -} - -SrcFinfo2< double, double >* DifBuffer::outerDifSourceOut(){ - static SrcFinfo2< double, double > sourceOut("outerDifSourceOut", - "Sends out source information."); - return &sourceOut; -} const Cinfo * DifBuffer::initCinfo() { - static DestFinfo process( "process", - "Handles process call", - new ProcOpFunc< DifBuffer >( &DifBuffer::process) ); - static DestFinfo reinit( "reinit", - "Reinit happens only in stage 0", - new ProcOpFunc< DifBuffer >( &DifBuffer::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 DifBuffer " - "computations in 2 stages, much as in the Compartment object. " - "In the first stage we send out the concentration value to other DifBuffers and Buffer elements. We also", - processShared, - sizeof( processShared ) / sizeof( Finfo* )); - - /*static DestFinfo process1( "process", - "Handle process call", - new ProcOpFunc< DifBuffer >( &DifBuffer::process_1 ) ); - static DestFinfo reinit1( "reinit", - "Reinit happens only in stage 0", - new ProcOpFunc< DifBuffer >( &DifBuffer::reinit_1 ) ); - - static Finfo* processShared_1[] = { - &process1, &reinit1 - }; - - static SharedFinfo process_1( "process_1", - "Second process call", - processShared_1, - sizeof( processShared_1 ) / sizeof( Finfo* ) ); - */ - //// Diffusion related shared messages - - - static DestFinfo concentration("concentration", - "Receives concentration (from DifShell).", - new EpFunc1<DifBuffer, double>(&DifBuffer::buffer)); - static Finfo* bufferShared[] = { - &concentration, DifBuffer::reactionOut() - }; - static SharedFinfo buffer( "buffer", - "This is a shared message with DifShell. " - "During stage 0:\n " - " - DifBuffer 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< DifBuffer, double, double > ( &DifBuffer::fluxFromOut )); - - static Finfo* innerDifShared[] = { - &fluxFromOut, - DifBuffer::innerDifSourceOut() - - }; - - static SharedFinfo innerDif( "innerDif", - "This shared message (and the next) is between DifBuffers: 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< DifBuffer, double, double> ( &DifBuffer::fluxFromIn) ); - - static Finfo* outerDifShared[] = { - &fluxFromIn, - DifBuffer::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<DifBuffer, double> activation("activation", - "Ion concentration from incoming conc message.", - &DifBuffer::setActivation, - &DifBuffer::getActivation); - static ElementValueFinfo<DifBuffer, double> kf("kf", - "Forward rate constant of buffer molecules 1/mM/s (?)", - &DifBuffer::setKf, - &DifBuffer::getKf); - static ElementValueFinfo<DifBuffer, double> kb("kb", - "Backward rate constant of buffer molecules. 1/s", - &DifBuffer::setKb, - &DifBuffer::getKb); - static ElementValueFinfo<DifBuffer, double> D("D", - "Diffusion constant of buffer molecules. m^2/s", - &DifBuffer::setD, - &DifBuffer::getD); - static ReadOnlyElementValueFinfo<DifBuffer, double> bFree("bFree", - "Free buffer concentration", - &DifBuffer::getBFree); - static ReadOnlyElementValueFinfo<DifBuffer, double> bBound("bBound", - "Bound buffer concentration", - &DifBuffer::getBBound); - /* static ReadOnlyElementValueFinfo<DifBuffer, double> prevFree("prevFree", - "Free buffer concentration in previous step", - &DifBuffer::getPrevFree); - static ReadOnlyElementValueFinfo<DifBuffer, double> prevBound("prevBound", - "Bound buffer concentration in previous step", - &DifBuffer::getPrevBound);*/ - static ElementValueFinfo<DifBuffer, double> bTot("bTot", - "Total buffer concentration.", - &DifBuffer::setBTot, - &DifBuffer::getBTot); - static ElementValueFinfo<DifBuffer, double> length("length", - "Length of shell", - &DifBuffer::setLength, - &DifBuffer::getLength); - static ElementValueFinfo<DifBuffer, double> diameter("diameter", - "Diameter of shell", - &DifBuffer::setDiameter, - &DifBuffer::getDiameter); - static ElementValueFinfo<DifBuffer, int> shapeMode("shapeMode", - "shape of the shell: SHELL=0, SLICE=SLAB=1, USERDEF=3", - &DifBuffer::setShapeMode, - &DifBuffer::getShapeMode); - - static ElementValueFinfo<DifBuffer, double> thickness("thickness", - "Thickness of shell", - &DifBuffer::setThickness, - &DifBuffer::getThickness); - - static ElementValueFinfo<DifBuffer, double> innerArea("innerArea", - "Inner area of shell", - &DifBuffer::setInnerArea, - &DifBuffer::getInnerArea); - static ElementValueFinfo<DifBuffer, double> outerArea("outerArea", - "Outer area of shell", - &DifBuffer::setOuterArea, - &DifBuffer::getOuterArea); - static ElementValueFinfo< DifBuffer, double> volume( "volume", "", - &DifBuffer::setVolume, - &DifBuffer::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", "DifBuffer", @@ -282,8 +66,8 @@ const Cinfo * DifBuffer::initCinfo() static Cinfo difBufferCinfo( "DifBuffer", Neutral::initCinfo(), - difBufferFinfos, - sizeof(difBufferFinfos)/sizeof(Finfo*), + 0, + 0, &dinfo, doc, sizeof(doc)/sizeof(string)); @@ -299,17 +83,6 @@ static const Cinfo * difBufferCinfo = DifBuffer::initCinfo(); //////////////////////////////////////////////////////////////////////////////// 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), diameter_(0), length_(0), @@ -322,12 +95,12 @@ DifBuffer::DifBuffer() : //////////////////////////////////////////////////////////////////////////////// // Field access functions //////////////////////////////////////////////////////////////////////////////// -double DifBuffer::getActivation(const Eref& e) const +double DifBuffer::vGetActivation(const Eref& e) const { return activation_; } -void DifBuffer::setActivation(const Eref& e,double value) +void DifBuffer::vSetActivation(const Eref& e,double value) { if ( value < 0.0 ) { cerr << "Error: DifBuffer: Activation cannot be negative!\n"; @@ -337,47 +110,38 @@ void DifBuffer::setActivation(const Eref& e,double value) } -double DifBuffer::getBFree(const Eref& e) const +double DifBuffer::vGetBFree(const Eref& e) const { return bFree_; } -double DifBuffer::getBBound(const Eref& e) const +double DifBuffer::vGetBBound(const Eref& e) const { return bBound_; } -/*double DifBuffer::getPrevFree(const Eref& e) const -{ - return prevFree_; -} - -double DifBuffer::getPrevBound(const Eref& e) const -{ - return prevBound_; -} -*/ -double DifBuffer::getBTot(const Eref& e) const +double DifBuffer::vGetBTot(const Eref& e) const { return bTot_; } -void DifBuffer::setBTot(const Eref& e,double value) +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_; } -double DifBuffer::getKf(const Eref& e) const +double DifBuffer::vGetKf(const Eref& e) const { return kf_; } -void DifBuffer::setKf(const Eref& e,double value) +void DifBuffer::vSetKf(const Eref& e,double value) { if ( value < 0.0 ) { cerr << "Error: DifBuffer: Kf cannot be negative!\n"; @@ -387,12 +151,12 @@ void DifBuffer::setKf(const Eref& e,double value) } -double DifBuffer::getKb(const Eref& e) const +double DifBuffer::vGetKb(const Eref& e) const { return kb_; } -void DifBuffer::setKb(const Eref& e,double value) +void DifBuffer::vSetKb(const Eref& e,double value) { if ( value < 0.0 ) { cerr << "Error: DifBuffer: Kb cannot be negative!\n"; @@ -401,12 +165,12 @@ void DifBuffer::setKb(const Eref& e,double value) kb_ = value; } -double DifBuffer::getD(const Eref& e) const +double DifBuffer::vGetD(const Eref& e) const { return D_; } -void DifBuffer::setD(const Eref& e,double value) +void DifBuffer::vSetD(const Eref& e,double value) { if ( value < 0.0 ) { @@ -416,152 +180,24 @@ void DifBuffer::setD(const Eref& e,double value) D_ = value; } -int DifBuffer::getShapeMode(const Eref& e) const -{ - return shapeMode_; -} - -void DifBuffer::setShapeMode(const Eref& e,int value) -{ - if ( value != 0 && value !=1 && value != 3 ) { - cerr << "Error: DifBuffer: Shape mode can only be 0, 1 or 3"; - return; - } - shapeMode_ = value; -} - -double DifBuffer::getLength(const Eref& e) const -{ - return length_; -} -void DifBuffer::setLength(const Eref& e,double value) -{ - if ( value < 0.0) { - cerr << "Error: DifBuffer: Length cannot be negative!\n"; - return; - } - length_ = value; -} - -double DifBuffer::getDiameter(const Eref& e) const -{ - return diameter_; -} - -void DifBuffer::setDiameter(const Eref& e,double value) -{ - if ( value < 0.0) { - cerr << "Error: DifBuffer: Diameter cannot be negative!\n"; - return; - } - diameter_ = value; -} - -double DifBuffer::getThickness(const Eref& e) const -{ - return thickness_; -} - -void DifBuffer::setThickness(const Eref& e,double value) -{ - if ( value < 0.0) { - cerr << "Error: DifBuffer: Thickness cannot be negative!\n"; - return; - } - thickness_ = value; -} - -void DifBuffer::setVolume(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::getVolume(const Eref& e ) const -{ - return volume_; -} - -void DifBuffer::setOuterArea(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::getOuterArea(const Eref& e ) const -{ - return outerArea_; -} - -void DifBuffer::setInnerArea(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::getInnerArea(const Eref& e) const -{ - return innerArea_; -} -// Dest function for conenctration received from DifShell. This -// function updates buffer concentration and sends back immediately -// using call back. -// TODO: complete this ... where do we get dt from? ProcPtr p->dt -//Check DifShell - -//In Genesis the following are messages: - - -void DifBuffer::reinit( const Eref& e, ProcPtr p ) -{ - localReinit_0( e, p ); -} - -void DifBuffer::process( const Eref& e, ProcPtr p ) -{ - localProcess_0( e, p ); -} - -/*void DifBuffer::process_1(const Eref& e, ProcPtr p ) +void DifBuffer::vBuffer(const Eref& e, + double C ) { - localProcess_1( e, p ); + activation_ = C; } -void DifBuffer::reinit_1(const Eref& e, ProcPtr p ) +double DifBuffer::integrate( double state, double dt, double A, double B ) { - ; -} -*/ -void DifBuffer::buffer(const Eref& e, - double C ) -{ - activation_ = C; + if ( B > EPSILON ) { + double x = exp( -B * dt ); + return state * x + ( A / B ) * ( 1 - x ); + } + return state + A * dt ; } -void DifBuffer::localProcess_0( const Eref & e, ProcPtr p ) +void DifBuffer::vProcess( const Eref & e, ProcPtr p ) { /** * Send ion concentration and thickness to adjacent DifShells. They will @@ -572,7 +208,9 @@ void DifBuffer::localProcess_0( const Eref & e, ProcPtr p ) outerDifSourceOut()->send( e, bFree_, thickness_ ); Af_ += kb_ * bBound_; Bf_ += kf_ * activation_; - bFree_ += (Af_ - bFree_ * Bf_) * p->dt; + bFree_ = integrate(bFree_,p->dt,Af_,Bf_); + bBound_ = bTot_ - bFree_; + reactionOut()->send(e,kf_,kb_,bFree_,bBound_); /** * Send ion concentration to ion buffers. They will send back information on @@ -583,12 +221,13 @@ void DifBuffer::localProcess_0( const Eref & e, ProcPtr p ) } -void DifBuffer::localReinit_0( const Eref& e, ProcPtr p ) +void DifBuffer::vReinit( const Eref& e, ProcPtr p ) { const double dOut = diameter_; const double dIn = diameter_ - thickness_; - + bFree_ = bTot_; + bBound_ = 0; switch ( shapeMode_ ) { /* @@ -628,23 +267,15 @@ void DifBuffer::localReinit_0( const Eref& e, ProcPtr p ) assert( 0 ); } } -/* -void DifBuffer::localProcess_1(const Eref& e, ProcPtr p ) -{ - //activation_ = conc; - Af_ += kb_ * bBound_; - Bf_ += kf_ * activation_; - bFree_ += (Af_ - bFree_ * Bf_) * p->dt; -} -*/ -void DifBuffer::fluxFromIn(const Eref& e,double innerC, double innerThickness) + +void DifBuffer::vFluxFromIn(const Eref& e,double innerC, double innerThickness) { double dif = 2 * D_ * innerArea_ / ((thickness_ + innerThickness) * volume_); Af_ += dif * innerC; Bf_ += dif; } -void DifBuffer::fluxFromOut(const Eref& e,double outerC, double outerThickness) +void DifBuffer::vFluxFromOut(const Eref& e,double outerC, double outerThickness) { double dif = 2 * D_ * outerArea_ / ((thickness_ + outerThickness) * volume_); Af_ += dif * outerC; diff --git a/moose-core/biophysics/DifBuffer.h b/moose-core/biophysics/DifBuffer.h index 500b4326..57307491 100644 --- a/moose-core/biophysics/DifBuffer.h +++ b/moose-core/biophysics/DifBuffer.h @@ -52,69 +52,40 @@ class DifBuffer 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 getActivation(const Eref& e) const; - void setActivation(const Eref& e,double value); + double vGetActivation(const Eref& e) const; + void vSetActivation(const Eref& e,double value); - double getBFree(const Eref& e) const; - double getBBound(const Eref& e) const; + double vGetBFree(const Eref& e) const; + double vGetBBound(const Eref& e) const; - //double getPrevFree(const Eref& e) const; // Bfree at previous time step - //double getPrevBound(const Eref& e) const; // Bbound at previous time step + double vGetBTot(const Eref& e) const; // total buffer concentration in mM (free + bound) + void vSetBTot(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 vGetKf(const Eref& e) const; // forward rate constant in 1/(mM*sec) + void vSetKf(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); - - int getShapeMode(const Eref& e) const; // Set to one of the predefined global - void setShapeMode(const Eref& e,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); + double vGetKb(const Eref& e) const; // backward rate constant in 1/sec + void vSetKb(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); // + double vGetD(const Eref& e) const; // diffusion constant of buffer molecules, m^2/sec + void vSetD(const Eref& e,double value); //void concentration(); - 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 ); - static SrcFinfo4< double, double, double, double >* reactionOut(); - static SrcFinfo2< double, double >* innerDifSourceOut(); - static SrcFinfo2< double, double >* outerDifSourceOut(); + static const Cinfo * initCinfo(); private: - void localReinit_0( const Eref & e, ProcPtr p ); - void localProcess_0( const Eref & e, ProcPtr p ); + + double integrate(double state, double dt, double A, double B ); double activation_; //ion concentration from incoming CONCEN message double Af_; @@ -127,13 +98,7 @@ class DifBuffer double kf_; //forward rate constant double kb_; //backward rate constant double D_; //diffusion constant - int shapeMode_; //SHELL=0, SLICE=SLAB=1, USERDEF=3 - double diameter_; //shell diameter - double length_; //shell length - double thickness_; //shell thickness - double outerArea_; //area of upper (outer) shell surface - double innerArea_; //area of lower (inner) shell surface - double volume_; + static const double EPSILON; }; #endif // _DifBuffer_h diff --git a/moose-core/biophysics/DifBufferBase.cpp b/moose-core/biophysics/DifBufferBase.cpp new file mode 100644 index 00000000..3b894c35 --- /dev/null +++ b/moose-core/biophysics/DifBufferBase.cpp @@ -0,0 +1,435 @@ +#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 ReadOnlyElementValueFinfo<DifBufferBase, double> bFree("bFree", + "Free buffer concentration", + &DifBufferBase::getBFree); + static ReadOnlyElementValueFinfo<DifBufferBase, double> bBound("bBound", + "Bound buffer concentration", + &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, 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 Dinfo<DifBufferBase> 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() : + shapeMode_(0), + diameter_(0), + length_(0), + thickness_(0), + outerArea_(0), + innerArea_(0), + volume_(0) +{} + + +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); +} + +double DifBufferBase::getBBound(const Eref& e) const +{ + return vGetBBound(e); +} + +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); +} + +int DifBufferBase::getShapeMode(const Eref& e) const +{ + return shapeMode_; +} + +void DifBufferBase::setShapeMode(const Eref& e,int value) +{ + if ( value != 0 && value !=1 && value != 3 ) { + cerr << "Error: DifBuffer: Shape mode can only be 0, 1 or 3"; + return; + } + shapeMode_ = value; +} + +double DifBufferBase::getLength(const Eref& e) const +{ + return length_; +} + +void DifBufferBase::setLength(const Eref& e,double value) +{ + if ( value < 0.0) { + cerr << "Error: DifBuffer: Length cannot be negative!\n"; + return; + } + length_ = value; +} + +double DifBufferBase::getDiameter(const Eref& e) const +{ + return diameter_; +} + +void DifBufferBase::setDiameter(const Eref& e,double value) +{ + if ( value < 0.0) { + cerr << "Error: DifBuffer: Diameter cannot be negative!\n"; + return; + } + diameter_ = value; +} + +double DifBufferBase::getThickness(const Eref& e) const +{ + return thickness_; +} + +void DifBufferBase::setThickness(const Eref& e,double value) +{ + if ( value < 0.0) { + cerr << "Error: DifBuffer: Thickness cannot be negative!\n"; + return; + } + thickness_ = value; +} + +void DifBufferBase::setVolume(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 DifBufferBase::getVolume(const Eref& e ) const +{ + return volume_; +} + +void DifBufferBase::setOuterArea(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 DifBufferBase::getOuterArea(const Eref& e ) const +{ + return outerArea_; +} + +void DifBufferBase::setInnerArea(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 DifBufferBase::getInnerArea(const Eref& e) const +{ + return innerArea_; +} + +} + +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); +} diff --git a/moose-core/biophysics/DifBufferBase.h b/moose-core/biophysics/DifBufferBase.h new file mode 100644 index 00000000..e169ec29 --- /dev/null +++ b/moose-core/biophysics/DifBufferBase.h @@ -0,0 +1,108 @@ +/********************************************************************** + ** 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: + + 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; + double getBBound(const Eref& e) const; + + 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); + + + int getShapeMode(const Eref& e) const; + void setShapeMode(const Eref& e,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 double vGetBBound(const Eref& e) const = 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; + + static SrcFinfo4< double, double, double, double >* reactionOut(); + static SrcFinfo2< double, double >* innerDifSourceOut(); + static SrcFinfo2< double, double >* outerDifSourceOut(); + static const Cinfo * initCinfo(); +private: + + unsigned int shapeMode_; + double length_; + double diameter_; + double thickness_; + double volume_; + double outerArea_; + double innerArea_; + +}; + + + + +#endif //_DIFBUFFER_BASE_H diff --git a/moose-core/biophysics/DifShell.cpp b/moose-core/biophysics/DifShell.cpp index bbb8fefa..b029f3a1 100644 --- a/moose-core/biophysics/DifShell.cpp +++ b/moose-core/biophysics/DifShell.cpp @@ -10,234 +10,15 @@ #include "header.h" #include "DifShell.h" -#include "ElementValueFinfo.h" -#include "../utility/numutil.h" -#include <cmath> +#include "DifShellBase.h" + const double DifShell::EPSILON = 1.0e-10; const double DifShell::F = 96485.3415; /* C / mol like in genesis */ -SrcFinfo1< double >* DifShell::concentrationOut() -{ - static SrcFinfo1< double > concentrationOut("concentrationOut", - "Sends out concentration"); - return &concentrationOut; -} - -SrcFinfo2< double, double >* DifShell::innerDifSourceOut(){ - static SrcFinfo2< double, double > sourceOut("innerDifSourceOut", - "Sends out source information."); - return &sourceOut; -} - -SrcFinfo2< double, double >* DifShell::outerDifSourceOut(){ - static SrcFinfo2< double, double > sourceOut("outerDifSourceOut", - "Sends out source information."); - return &sourceOut; -} - - const Cinfo* DifShell::initCinfo() { - static DestFinfo process( "process", - "Handles process call", - new ProcOpFunc< DifShell >( &DifShell::process ) ); - static DestFinfo reinit( "reinit", - "Reinit happens only in stage 0", - new ProcOpFunc< DifShell >( &DifShell::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 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, - sizeof( processShared ) / sizeof( Finfo* )); - - // static DestFinfo process1( "process", - // "Handle process call", - // new ProcOpFunc< DifShell >( &DifShell::process1 ) ); - // static DestFinfo reinit1( "reinit", - // "Reinit happens only in stage 0", - // new ProcOpFunc< DifShell >( &DifShell::reinit1) - // ); - // 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 EpFunc4< DifShell, double, double, double, double >( &DifShell::buffer )); - - static Finfo* bufferShared[] = { - DifShell::concentrationOut(), &reaction - }; - 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 EpFunc2< DifShell, double, double > ( &DifShell::fluxFromOut )); - - static Finfo* innerDifShared[] = { - &fluxFromOut, - DifShell::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< DifShell, double, double> ( &DifShell::fluxFromIn ) ); - - static Finfo* outerDifShared[] = { - &fluxFromIn, - DifShell::outerDifSourceOut(), - }; - - static SharedFinfo outerDif( "outerDif", - "Using this message, an outer shell sends to, and receives from its inner shell." , - outerDifShared, - sizeof( outerDifShared ) / sizeof( Finfo* )); - //Should be changed it to ElementValueFinfo, to correctly set initial - //conditions (AJS) - - static ElementValueFinfo< DifShell, double> C( "C", - "Concentration C",// is computed by the DifShell", - &DifShell::setC, - &DifShell::getC); - static ElementValueFinfo< DifShell, double> Ceq( "Ceq", "", - &DifShell::setCeq, - &DifShell::getCeq); - static ElementValueFinfo< DifShell, double> D( "D", "", - &DifShell::setD, - &DifShell::getD); - static ElementValueFinfo< DifShell, double> valence( "valence", "", - &DifShell::setValence, - &DifShell::getValence); - static ElementValueFinfo< DifShell, double> leak( "leak", "", - &DifShell::setLeak, - &DifShell::getLeak); - static ElementValueFinfo< DifShell, unsigned int> shapeMode( "shapeMode", "", - &DifShell::setShapeMode, - &DifShell::getShapeMode); - static ElementValueFinfo< DifShell, double> length( "length", "", - &DifShell::setLength, - &DifShell::getLength); - static ElementValueFinfo< DifShell, double> diameter( "diameter", "", - &DifShell::setDiameter, - &DifShell::getDiameter ); - static ElementValueFinfo< DifShell, double> thickness( "thickness", "", - &DifShell::setThickness, - &DifShell::getThickness ); - static ElementValueFinfo< DifShell, double> volume( "volume", "", - &DifShell::setVolume, - &DifShell::getVolume ); - static ElementValueFinfo< DifShell, double> outerArea( "outerArea", "", - &DifShell::setOuterArea, - &DifShell::getOuterArea); - static ElementValueFinfo< DifShell, double> innerArea( "innerArea", "", - &DifShell::setInnerArea, - &DifShell::getInnerArea ); - - - static DestFinfo influx( "influx", "", - new EpFunc1< DifShell, double > (&DifShell::influx )); - static DestFinfo outflux( "outflux", "", - new EpFunc1< DifShell, double >( &DifShell::influx )); - static DestFinfo fInflux( "fInflux", "", - new EpFunc2< DifShell, double, double >( &DifShell::fInflux )); - static DestFinfo fOutflux( "fOutflux", "", - new EpFunc2< DifShell, double, double> (&DifShell::fOutflux )); - static DestFinfo storeInflux( "storeInflux", "", - new EpFunc1< DifShell, double >( &DifShell::storeInflux ) ); - static DestFinfo storeOutflux( "storeOutflux", "", - new EpFunc1< DifShell, double > ( &DifShell::storeOutflux ) ); - static DestFinfo tauPump( "tauPump","", - new EpFunc2< DifShell, double, double >( &DifShell::tauPump ) ); - static DestFinfo eqTauPump( "eqTauPump", "", - new EpFunc1< DifShell, double >( &DifShell::eqTauPump ) ); - static DestFinfo mmPump( "mmPump", "", - new EpFunc2< DifShell, double, double >( &DifShell::mmPump ) ); - static DestFinfo hillPump( "hillPump", "", - new EpFunc3< 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 - ////////////////////////////////////////////////////////////////// - &proc, - &buffer, - // concentrationOut(), - //innerDifSourceOut(), - //outerDifSourceOut(), - &innerDif, - &outerDif, - ////////////////////////////////////////////////////////////////// - // DestFinfo definitions - ////////////////////////////////////////////////////////////////// - &influx, - &outflux, - &fInflux, - &fOutflux, - &storeInflux, - &storeOutflux, - &tauPump, - &eqTauPump, - &mmPump, - &hillPump, - }; - static string doc[] = { "Name", "DifShell", @@ -252,8 +33,8 @@ const Cinfo* DifShell::initCinfo() static Cinfo difShellCinfo( "DifShell", Neutral::initCinfo(), - difShellFinfos, - sizeof( difShellFinfos ) / sizeof( Finfo* ), + 0, + 0, &dinfo, doc, sizeof( doc ) / sizeof( string )); @@ -278,14 +59,7 @@ DifShell::DifShell() : 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 ) + leak_( 0.0 ) { ; } //////////////////////////////////////////////////////////////////////////////// @@ -294,7 +68,7 @@ DifShell::DifShell() : /// C is a read-only field -void DifShell::setC( const Eref& e, double C) +void DifShell::vSetC( const Eref& e, double C) { if ( C < 0.0 ) { cerr << "Error: DifShell: C cannot be negative!\n"; @@ -303,12 +77,12 @@ void DifShell::setC( const Eref& e, double C) C_ = C; } -double DifShell::getC(const Eref& e) const +double DifShell::vGetC(const Eref& e) const { return C_; } -void DifShell::setCeq( const Eref& e, double Ceq ) +void DifShell::vSetCeq( const Eref& e, double Ceq ) { if ( Ceq < 0.0 ) { cerr << "Error: DifShell: Ceq cannot be negative!\n"; @@ -318,12 +92,12 @@ void DifShell::setCeq( const Eref& e, double Ceq ) Ceq_ = Ceq; } -double DifShell::getCeq(const Eref& e ) const +double DifShell::vGetCeq(const Eref& e ) const { return Ceq_; } -void DifShell::setD(const Eref& e, double D ) +void DifShell::vSetD(const Eref& e, double D ) { if ( D < 0.0 ) { cerr << "Error: DifShell: D cannot be negative!\n"; @@ -333,12 +107,12 @@ void DifShell::setD(const Eref& e, double D ) D_ = D; } -double DifShell::getD(const Eref& e ) const +double DifShell::vGetD(const Eref& e ) const { return D_; } -void DifShell::setValence(const Eref& e, double valence ) +void DifShell::vSetValence(const Eref& e, double valence ) { if ( valence < 0.0 ) { cerr << "Error: DifShell: valence cannot be negative!\n"; @@ -348,249 +122,21 @@ void DifShell::setValence(const Eref& e, double valence ) valence_ = valence; } -double DifShell::getValence(const Eref& e ) const +double DifShell::vGetValence(const Eref& e ) const { return valence_; } -void DifShell::setLeak(const Eref& e, double leak ) +void DifShell::vSetLeak(const Eref& e, double leak ) { leak_ = leak; } -double DifShell::getLeak(const Eref& e ) const +double DifShell::vGetLeak(const Eref& e ) const { return leak_; } -void DifShell::setShapeMode(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; -} - -unsigned int DifShell::getShapeMode(const Eref& e) const -{ - return shapeMode_; -} - -void DifShell::setLength(const Eref& e, double length ) -{ - if ( length < 0.0 ) { - cerr << "Error: DifShell: length cannot be negative!\n"; - return; - } - - length_ = length; -} - -double DifShell::getLength(const Eref& e ) const -{ - return length_; -} - -void DifShell::setDiameter(const Eref& e, double diameter ) -{ - if ( diameter < 0.0 ) { - cerr << "Error: DifShell: diameter cannot be negative!\n"; - return; - } - - diameter_ = diameter; -} - -double DifShell::getDiameter(const Eref& e ) const -{ - return diameter_; -} - -void DifShell::setThickness( const Eref& e, double thickness ) -{ - if ( thickness < 0.0 ) { - cerr << "Error: DifShell: thickness cannot be negative!\n"; - return; - } - - thickness_ = thickness; -} - -double DifShell::getThickness(const Eref& e) const -{ - return thickness_; -} - -void DifShell::setVolume(const Eref& e, double volume ) -{ - 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; - } - - volume_ = volume; -} - -double DifShell::getVolume(const Eref& e ) const -{ - return volume_; -} - -void DifShell::setOuterArea(const Eref& e, double outerArea ) -{ - 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; - } - - outerArea_ = outerArea; -} - -double DifShell::getOuterArea(const Eref& e ) const -{ - return outerArea_; -} - -void DifShell::setInnerArea(const Eref& e, double innerArea ) -{ - 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; - } - - innerArea_ = innerArea; -} - -double DifShell::getInnerArea(const Eref& e) const -{ - return innerArea_; -} - -//////////////////////////////////////////////////////////////////////////////// -// Dest functions -//////////////////////////////////////////////////////////////////////////////// - -void DifShell::reinit( const Eref& e, ProcPtr p ) -{ - localReinit_0( e, p ); -} - -void DifShell::process( 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( - const Eref& e, - double kf, - double kb, - double bFree, - double bBound ) -{ - localBuffer( e, kf, kb, bFree, bBound ); -} - -void DifShell::fluxFromOut(const Eref& e, - double outerC, - double outerThickness ) -{ - localFluxFromOut(e, outerC, outerThickness ); -} - -void DifShell::fluxFromIn( - const Eref& e, - double innerC, - double innerThickness ) -{ - localFluxFromIn( e, innerC, innerThickness ); -} - -void DifShell::influx(const Eref& e, - double I ) -{ - localInflux( e, I ); -} - -void DifShell::outflux(const Eref& e, - double I ) -{ - localOutflux(e, I ); -} - -void DifShell::fInflux(const Eref& e, - double I, - double fraction ) -{ - localFInflux(e, I, fraction ); -} - -void DifShell::fOutflux(const Eref& e, - double I, - double fraction ) -{ - localFOutflux(e, I, fraction ); -} - -void DifShell::storeInflux(const Eref& e, - double flux ) -{ - localStoreInflux( e, flux ); -} - -void DifShell::storeOutflux(const Eref& e, - double flux ) -{ - localStoreOutflux(e, flux ); -} - -void DifShell::tauPump(const Eref& e, - double kP, - double Ceq ) -{ - localTauPump( e, kP, Ceq ); -} - -void DifShell::eqTauPump(const Eref& e, - double kP ) -{ - localEqTauPump(e, kP ); -} - -void DifShell::mmPump(const Eref& e, - double vMax, - double Kd ) -{ - localMMPump(e, vMax, Kd ); -} - -void DifShell::hillPump(const Eref& e, - double vMax, - double Kd, - unsigned int hill ) -{ - localHillPump(e, vMax, Kd, hill ); -} - //////////////////////////////////////////////////////////////////////////////// // Local dest functions //////////////////////////////////////////////////////////////////////////////// @@ -603,7 +149,7 @@ double DifShell::integrate( double state, double dt, double A, double B ) return state + A * dt ; } -void DifShell::localReinit_0( const Eref& e, ProcPtr p ) +void DifShell::vReinit( const Eref& e, ProcPtr p ) { dCbyDt_ = leak_; Cmultiplier_ = 0; @@ -651,7 +197,7 @@ void DifShell::localReinit_0( const Eref& e, ProcPtr p ) } } -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 @@ -660,9 +206,8 @@ void DifShell::localProcess_0( const Eref & e, ProcPtr p ) innerDifSourceOut()->send( e, C_, thickness_ ); outerDifSourceOut()->send( e, C_, thickness_ ); - 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) @@ -671,46 +216,41 @@ void DifShell::localProcess_0( const Eref & e, ProcPtr p ) */ concentrationOut()->send( e, C_ ); dCbyDt_ = leak_; + Cmultiplier_ = 0; } -/* - void DifShell::localProcess_1( const Eref& e, ProcPtr p ) - { - C_ += dCbyDt_ * p->dt; - dCbyDt_ = leak_; - } -*/ -void DifShell::localBuffer(const Eref& e, +void DifShell::vBuffer(const Eref& e, double kf, double kb, double bFree, double bBound ) { dCbyDt_ += kb * bBound; - Cmultiplier_ -= kf * bFree; + Cmultiplier_ += kf * bFree; } -void DifShell::localFluxFromOut(const Eref& e, double outerC, double outerThickness ) +void DifShell::vFluxFromOut(const Eref& e, double outerC, double outerThickness ) { - double dx = ( outerThickness + thickness_ ) / 2.0; + 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_ += ( D_ / volume_ ) * ( outerArea_ / dx ) * outerC; - Cmultiplier_ -= ( D_ / volume_ ) * ( outerArea_ / dx ) ; + dCbyDt_ += diff * outerC; + Cmultiplier_ += diff ; } -void DifShell::localFluxFromIn(const Eref& e, double innerC, double innerThickness ) +void DifShell::vFluxFromIn(const Eref& e, double innerC, double innerThickness ) { //influx from inner shell - double dx = ( innerThickness + thickness_ ) / 2.0; - - dCbyDt_ += ( D_ / volume_ ) * ( innerArea_ / dx ) * innerC ; - Cmultiplier_ -= ( D_ / volume_ ) * ( innerArea_ / dx ) ; + //double dx = ( innerThickness + thickness_ ) / 2.0; + double diff = 2.*( D_ / volume_ ) * ( innerArea_ / (innerThickness + thickness_) ); + + dCbyDt_ += diff * innerC ; + Cmultiplier_ += diff ; } -void DifShell::localInflux(const Eref& e, double I ) +void DifShell::vInflux(const Eref& e, double I ) { /** * I: Amperes @@ -723,54 +263,53 @@ void DifShell::localInflux(const Eref& e, double I ) /** * Same as influx, except subtracting. */ -void DifShell::localOutflux(const Eref& e, double I ) +void DifShell::vOutflux(const Eref& e, double I ) { dCbyDt_ -= I / ( F * valence_ * volume_ ); } -void DifShell::localFInflux(const Eref& e, double I, double fraction ) +void DifShell::vFInflux(const Eref& e, double I, double fraction ) { dCbyDt_ += fraction * I / ( F * valence_ * volume_ ); } -void DifShell::localFOutflux(const Eref& e, double I, double fraction ) +void DifShell::vFOutflux(const Eref& e, double I, double fraction ) { dCbyDt_ -= fraction * I / ( F * valence_ * volume_ ); } -void DifShell::localStoreInflux(const Eref& e, double flux ) +void DifShell::vStoreInflux(const Eref& e, double flux ) { dCbyDt_ += flux / volume_; } -void DifShell::localStoreOutflux(const Eref& e, double flux ) +void DifShell::vStoreOutflux(const Eref& e, double flux ) { dCbyDt_ -= flux / volume_; } -void DifShell::localTauPump(const Eref& e, 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(const Eref& e, double kP ) +void DifShell::vEqTauPump(const Eref& e, double kP ) { - //dCbyDt_ += -kP * ( C_ - Ceq_ ); dCbyDt_ += kP*Ceq_; Cmultiplier_ -= kP; } -void DifShell::localMMPump(const Eref& e, double vMax, double Kd ) +void DifShell::vMMPump(const Eref& e, double vMax, double Kd ) { Cmultiplier_ += -( vMax / volume_ ) / ( C_ + Kd ) ; } -void DifShell::localHillPump(const Eref& e, double vMax, double Kd, unsigned int hill ) +void DifShell::vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill ) { //This needs fixing double ch; diff --git a/moose-core/biophysics/DifShell.h b/moose-core/biophysics/DifShell.h index 35396bc1..3bc4d07c 100644 --- a/moose-core/biophysics/DifShell.h +++ b/moose-core/biophysics/DifShell.h @@ -8,149 +8,57 @@ ** See the file COPYING.LIB for the full notice. **********************************************************************/ -#ifndef _DifShell_h -#define _DifShell_h +#ifndef _DIFSHELL_H +#define _DIFSHELL_H -class DifShell +class DifShell: public DifShellBase { public: DifShell(); - + ///////////////////////////////////////////////////////////// - // Field access functions + // 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 ); - 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; - ///////////////////////////////////////////////////////////// - // Dest functions + // Field access functions ///////////////////////////////////////////////////////////// - void reinit( const Eref & e, ProcPtr p ); - - void process(const Eref & e, ProcPtr p ); - /* - void process_1(const Eref & e, ProcPtr p ); - - void reinit_1(const Eref & e, ProcPtr p ); // dummyFunc - */ - 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 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 tauPump( - const Eref& e, - double kP, - double Ceq ); + void vSetD(const Eref& e, double D ); + double vGetD(const Eref& e) const; - void eqTauPump( - const Eref& e, - double kP ); + void vSetValence(const Eref& e, double valence ); + double vGetValence(const Eref& e) const; - void mmPump( - const Eref& e, - double vMax, - double Kd ); + void vSetLeak(const Eref& e, double leak ); + double vGetLeak(const Eref& e) const; - void hillPump( - const Eref& e, - double vMax, - double Kd, - unsigned int hill ); - static SrcFinfo1< double >* concentrationOut(); - static SrcFinfo2< double, double >* innerDifSourceOut(); - static SrcFinfo2< double, double >* outerDifSourceOut(); 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(const Eref& e, double kf, double kb, double bFree, double bBound ); - void localFluxFromOut( const Eref& e,double outerC, double outerThickness ); - void localFluxFromIn(const Eref& e, double innerC, double innerThickness ); - void localInflux(const Eref& e, double I ); - void localOutflux( const Eref& e,double I ); - void localFInflux(const Eref& e, double I, double fraction ); - void localFOutflux(const Eref& e, double I, double fraction ); - void localStoreInflux(const Eref& e, double flux ); - void localStoreOutflux(const Eref& e, double flux ); - void localTauPump(const Eref& e, double kP, double Ceq ); - void localEqTauPump(const Eref& e, double kP ); - void localMMPump(const Eref& e, double vMax, double Kd ); - void localHillPump(const Eref& e, double vMax, double Kd, unsigned int hill ); + double integrate( double state, double dt, double A, double B ); double dCbyDt_; @@ -160,13 +68,6 @@ class DifShell 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) @@ -174,4 +75,4 @@ class DifShell }; -#endif // _DifShell_h +#endif // _DIFSHELL_H diff --git a/moose-core/biophysics/Makefile b/moose-core/biophysics/Makefile index afe96aa0..bc95b4d2 100644 --- a/moose-core/biophysics/Makefile +++ b/moose-core/biophysics/Makefile @@ -38,7 +38,9 @@ OBJ = \ NMDAChan.o \ testBiophysics.o \ IzhikevichNrn.o \ + DifShellBase.o \ DifShell.o \ + DifBufferBase.o \ DifBuffer.o \ Leakage.o \ VectorTable.o \ @@ -95,8 +97,10 @@ 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 -DifBuffer.o: DifBuffer.h +DifShellBase.o: DifShellBase.h +DifShell.o: DifShellBase.h DifShell.h +DifBufferBase.o: DifBufferBase.h +DifBuffer.o: DifBufferBase.h DifBuffer.h testBiophysics.o: IntFire.h CompartmentBase.h Compartment.h HHChannel.h HHGate.h VectorTable.o : VectorTable.h MarkovGslSolver.o : MarkovGslSolver.h -- GitLab