diff --git a/moose-core/biophysics/DifBuffer.cpp b/moose-core/biophysics/DifBuffer.cpp index 0d695f807d8043500445553f10bf62a57a129311..5b460fd1b906efbe4d4eca39465508d2f345aa21 100644 --- a/moose-core/biophysics/DifBuffer.cpp +++ b/moose-core/biophysics/DifBuffer.cpp @@ -89,6 +89,8 @@ DifBuffer::DifBuffer() : Bf_(0), bFree_(0), bBound_(0), + prevFree_(0), + prevBound_(0), bTot_(0), kf_(0), kb_(0), @@ -137,6 +139,9 @@ void DifBuffer::vSetBFree(const Eref& e,double value) } bFree_ = value; bBound_ = bTot_ - bFree_; + prevFree_= bFree_; + prevBound_ = bBound_; + } double DifBuffer::vGetBBound(const Eref& e) const @@ -156,6 +161,8 @@ void DifBuffer::vSetBBound(const Eref& e,double value) } bBound_ = value; bFree_ = bTot_ - bBound_; + prevFree_= bFree_; + prevBound_ = bBound_; } @@ -172,6 +179,7 @@ void DifBuffer::vSetBTot(const Eref& e,double value) } bTot_ = value; bFree_ = bTot_; + bBound_ = 0; } @@ -340,6 +348,7 @@ void DifBuffer::vBuffer(const Eref& e, double C ) { activation_ = C; + //cout<<"Buffer"<< C<<" "; } @@ -360,14 +369,19 @@ void DifBuffer::vProcess( const Eref & e, ProcPtr p ) * then compute their incoming fluxes. */ - innerDifSourceOut()->send( e, bFree_, thickness_ ); - outerDifSourceOut()->send( e, bFree_, thickness_ ); + innerDifSourceOut()->send( e, prevFree_, thickness_ ); + outerDifSourceOut()->send( e, prevFree_, thickness_ ); + reactionOut()->send(e,kf_,kb_,bFree_,bBound_); + Af_ += kb_ * bBound_; Bf_ += kf_ * activation_; + bFree_ = integrate(bFree_,p->dt,Af_,Bf_); bBound_ = bTot_ - bFree_; - - reactionOut()->send(e,kf_,kb_,bFree_,bBound_); + 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) @@ -380,10 +394,11 @@ void DifBuffer::vProcess( 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; + + const double rOut = diameter_/2.; + + const double rIn = rOut - thickness_; + switch ( shapeMode_ ) { /* @@ -391,15 +406,15 @@ void DifBuffer::vReinit( const Eref& e, ProcPtr p ) */ case 0: if ( length_ == 0.0 ) { // Spherical shell - volume_ = ( M_PI / 6.0 ) * ( dOut * dOut * dOut - dIn * dIn * dIn ); - outerArea_ = M_PI * dOut * dOut; - innerArea_ = M_PI * dIn * dIn; + 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_ / 4.0 ) * ( dOut * dOut - dIn * dIn ); - outerArea_ = M_PI * dOut * length_; - innerArea_ = M_PI * dIn * length_; + volume_ = ( M_PI * length_ ) * ( rOut * rOut - rIn * rIn ); + outerArea_ = 2*M_PI * rOut * length_; + innerArea_ = 2*M_PI * rIn * length_; } - + break; /* @@ -422,6 +437,7 @@ void DifBuffer::vReinit( const Eref& e, ProcPtr p ) default: assert( 0 ); } + } void DifBuffer::vFluxFromIn(const Eref& e,double innerC, double innerThickness) diff --git a/moose-core/biophysics/DifBuffer.h b/moose-core/biophysics/DifBuffer.h index fe3f7772a0ced00af6b64a7d915308a73928493c..60f574458385487a732b3653d568e99fbba1a420 100644 --- a/moose-core/biophysics/DifBuffer.h +++ b/moose-core/biophysics/DifBuffer.h @@ -74,6 +74,8 @@ class DifBuffer: public DifBufferBase{ 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 diff --git a/moose-core/biophysics/DifShell.cpp b/moose-core/biophysics/DifShell.cpp index c0e3040b16843bd2afdbf0b0ce9f18c8fcf378c1..d4e37d085579d810d6d5fe784fbff94e673c98d0 100644 --- a/moose-core/biophysics/DifShell.cpp +++ b/moose-core/biophysics/DifShell.cpp @@ -58,6 +58,7 @@ DifShell::DifShell() : Cmultiplier_(0.0), C_( 0.0 ), Ceq_( 0.0 ), + prevC_(0.0), D_( 0.0 ), valence_( 0.0 ), leak_( 0.0 ), @@ -84,6 +85,7 @@ void DifShell::vSetC( const Eref& e, double C) } C_ = C; + prevC_ = C_; } double DifShell::vGetC(const Eref& e) const { @@ -270,6 +272,7 @@ double DifShell::integrate( double state, double dt, double A, double B ) double x = exp( -B * dt ); return state * x + ( A / B ) * ( 1 - x ); } + return state + A * dt ; } @@ -277,7 +280,7 @@ void DifShell::vReinit( const Eref& e, ProcPtr p ) { dCbyDt_ = leak_; Cmultiplier_ = 0; - C_ = Ceq_; + const double rOut = diameter_/2.; const double rIn = rOut - thickness_; @@ -329,19 +332,23 @@ void DifShell::vProcess( const Eref & e, ProcPtr p ) * then compute their incoming fluxes. */ - innerDifSourceOut()->send( e, C_, thickness_ ); - outerDifSourceOut()->send( e, C_, thickness_ ); + innerDifSourceOut()->send( e, prevC_, thickness_ ); + outerDifSourceOut()->send( e, prevC_, thickness_ ); + concentrationOut()->send( e, C_ ); 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. */ - concentrationOut()->send( e, C_ ); + prevC_ = C_; + + //cout<<"Shell "<< C_<<" "; dCbyDt_ = leak_; Cmultiplier_ = 0; + } void DifShell::vBuffer(const Eref& e, double kf, diff --git a/moose-core/biophysics/DifShell.h b/moose-core/biophysics/DifShell.h index 6bef08039138e5dd5bfeae0be1a4a80d4c0f64de..42cf643d69ca04e9af18f9cd5da8ea1729dcdf83 100644 --- a/moose-core/biophysics/DifShell.h +++ b/moose-core/biophysics/DifShell.h @@ -84,6 +84,7 @@ class DifShell: public DifShellBase{ double Cmultiplier_; double C_; double Ceq_; + double prevC_; double D_; double valence_; double leak_;