diff --git a/diffusion/ConcChanInfo.h b/diffusion/ConcChanInfo.h index 9aeb984b5d..160915fe87 100644 --- a/diffusion/ConcChanInfo.h +++ b/diffusion/ConcChanInfo.h @@ -14,22 +14,26 @@ */ class ConcChanInfo { -public: - ConcChanInfo(); - ConcChanInfo( unsigned int my, unsigned int other, - unsigned int chan, double perm, bool isSwapped ) - : myPool( my ), otherPool( other ), chanPool( chan ), - swapped( isSwapped ), - permeability( perm ) - {;} + public: + ConcChanInfo(); + ConcChanInfo( unsigned int my, unsigned int other, + unsigned int chan, double perm, bool isSwapped, bool _isLocal ) + : myPool( my ), otherPool( other ), chanPool( chan ), + swapped( isSwapped ), + isLocal( _isLocal ), + permeability( perm ) + {;} - unsigned int myPool; - unsigned int otherPool; - unsigned int chanPool; - bool swapped; // Flag records whether myPool/otherPool are - // swapped wrt inPool/outPool. System expects - // that inPool should be in same compt as chanPool. - double permeability; + unsigned int myPool; /// This is not Id, it is internal PoolIndex + /// Obtained using convertIdToPoolIndex. + unsigned int otherPool; /// Internal PoolIndex + unsigned int chanPool; /// Internal PoolIndex + bool swapped; // Flag records whether myPool/otherPool are + // swapped wrt inPool/outPool. System expects + // that inPool should be in same compt as chanPool. + bool isLocal; // Flag, records odd cases where the channel is + // entirely local to one compartment. + double permeability; }; #endif // _CONC_CHAN_INFO_H diff --git a/diffusion/Dsolve.cpp b/diffusion/Dsolve.cpp index d65b528acc..bfbd2d387c 100644 --- a/diffusion/Dsolve.cpp +++ b/diffusion/Dsolve.cpp @@ -412,6 +412,8 @@ void Dsolve::calcJnChan( const DiffJunction& jn, Dsolve* other, double dt ) for ( unsigned int i = 0; i < jn.myChannels.size(); ++i ) { + if ( channels_[ jn.myChannels[i] ].isLocal ) + continue; ConcChanInfo& myChan = channels_[ jn.myChannels[i] ]; DiffPoolVec& myDv = pools_[ myChan.myPool ]; DiffPoolVec& otherDv = other->pools_[ myChan.otherPool ]; @@ -426,7 +428,7 @@ void Dsolve::calcJnChan( const DiffJunction& jn, Dsolve* other, double dt ) double chanN = chanDv.getN( j->first ); // Stick in a conversion factor for the myN and otherN into // concentrations. Note that SI is millimolar. - double perm = myChan.permeability * chanN * 1000.0 / NA; + double perm = myChan.permeability * chanN / NA; myN = integ( myN, perm * myN/j->firstVol, perm * otherN/j->secondVol, dt ); otherN += lastN - myN; // Mass consv @@ -446,6 +448,8 @@ void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt ) { for ( unsigned int i = 0; i < jn.otherChannels.size(); ++i ) { + if ( other->channels_[ jn.otherChannels[i] ].isLocal ) + continue; ConcChanInfo& otherChan = other->channels_[ jn.otherChannels[i] ]; // This is the DiffPoolVec for the pools on the other Dsolve, // the one with the channel. @@ -465,7 +469,7 @@ void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt ) double chanN = chanDv.getN( j->second ); // Stick in a conversion factor for the myN and otherN into // concentrations. Note that SI is millimolar. - double perm = otherChan.permeability * chanN * 1000.0 / NA; + double perm = otherChan.permeability * chanN / NA; myN = integ( myN, perm * myN/j->firstVol, perm * otherN/j->secondVol, dt ); otherN += lastN - myN; // Mass consv @@ -480,6 +484,46 @@ void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt ) } } +void Dsolve::calcLocalChan( double dt ) +{ + // There may be some chans which connect within the compartment. + // This is odd biologically, but happens for modelling convenience. + // The channels have a flag indicating if this is so. It is a simple + // calculation, but lacks numerical accuracy as it is outside the + // solver framework. + + ChemCompt* cc = reinterpret_cast< ChemCompt* >( compartment_.eref().data() ); + const vector< double >& vols = cc->vGetVoxelVolume(); + for (auto ch = channels_.begin(); ch != channels_.end(); ++ch ) { + if ( ch->isLocal ) { + DiffPoolVec& myDv = pools_[ ch->myPool ]; + DiffPoolVec& otherDv = pools_[ ch->otherPool ]; + DiffPoolVec& chanDv = pools_[ ch->chanPool ]; + assert( myDv.getNumVoxels() == numVoxels_ ); + for (unsigned int i = 0; i < numVoxels_; ++i ) { + double myN = myDv.getN( i ); + double lastN = myN; + double otherN = otherDv.getN( i ); + double chanN = chanDv.getN( i ); + double vol = vols[i]; + // Stick in a conversion factor for the myN and otherN into + // concentrations. Note that SI is millimolar. + double perm = ch->permeability * chanN / NA; + myN = integ( myN, perm * myN/vol, perm * otherN/vol, dt ); + otherN += lastN - myN; // Mass consv + if ( otherN < 0.0 ) // Avoid negatives + { + myN += otherN; + otherN = 0.0; + } + myDv.setN( i, myN ); + otherDv.setN( i, otherN ); + } + } + } +} + + /** * Computes flux through a junction between diffusion solvers. * Most used at junctions on spines and PSDs, but can also be used @@ -525,6 +569,7 @@ void Dsolve::reinit( const Eref& e, ProcPtr p ) void Dsolve::updateJunctions( double dt ) { + calcLocalChan( dt ); for (auto i = junctions_.begin(); i != junctions_.end(); ++i ) calcJunction( *i, dt ); } @@ -614,7 +659,8 @@ void Dsolve::fillConcChans( const vector< ObjId >& chans ) ret.clear(); unsigned int outPoolValue = outPool.id.value(); bool swapped = false; - if ( !( inPool.bad() or chanPool.bad() ) ) + bool isLocal = false; + if ( !( inPool.bad() || chanPool.bad() ) ) { unsigned int inPoolIndex = convertIdToPoolIndex( inPool.id ); unsigned int chanPoolIndex = convertIdToPoolIndex(chanPool.id); @@ -623,14 +669,22 @@ void Dsolve::fillConcChans( const vector< ObjId >& chans ) inPoolIndex = convertIdToPoolIndex( outPool.id ); outPoolValue = inPool.id.value(); swapped = true; - } + } else { + unsigned int outPoolIndex = convertIdToPoolIndex( outPool.id ); + // edge case where both in and out are on same compartment. + // Not much sense in bio, but a common convenience hack. + if ( outPoolIndex != ~0U ) { + outPoolValue = outPoolIndex; + isLocal = true; + } + + } if ( ( inPoolIndex != ~0U) && (chanPoolIndex != ~0U ) ) { ConcChanInfo cci( inPoolIndex, outPoolValue, chanPoolIndex, Field< double >::get( *i, "permeability" ), - ////// Fix it below //////// - swapped + swapped, isLocal ); channels_.push_back( cci ); } @@ -972,27 +1026,29 @@ void Dsolve::mapChansBetweenDsolves( DiffJunction& jn, Id self, Id other) unsigned int outIndex; for ( unsigned int i = 0; i < ch.size(); ++i ) { - unsigned int chanIndex = ch[i].chanPool; - outIndex = otherSolve->convertIdToPoolIndex( ch[i].otherPool ); - if ( (outIndex != ~0U) && (chanIndex != ~0U ) ) - { + if ( ch[i].isLocal ) { jn.myChannels.push_back(i); - ch[i].otherPool = outIndex; // replace the Id with the index. - ch[i].chanPool = chanIndex; //chanIndex may be on either Dsolve + } else { + outIndex = otherSolve->convertIdToPoolIndex( ch[i].otherPool ); + if (outIndex != ~0U) { + ch[i].otherPool = outIndex; // replace Id with the index. + jn.myChannels.push_back(i); + } } } // Now set up the other Dsolve. vector< ConcChanInfo >& ch2 = otherSolve->channels_; for ( unsigned int i = 0; i < ch2.size(); ++i ) { - unsigned int chanIndex = ch2[i].chanPool; - outIndex = selfSolve->convertIdToPoolIndex( ch2[i].otherPool ); - if ( (outIndex != ~0U) && (chanIndex != ~0U) ) - { + if ( ch2[i].isLocal ) { jn.otherChannels.push_back(i); - ch2[i].otherPool = outIndex; // replace the Id with the index - ch2[i].chanPool = chanIndex; //chanIndex may be on either Dsolve - } + } else { + outIndex = selfSolve->convertIdToPoolIndex( ch2[i].otherPool ); + if ( outIndex != ~0U ) { + ch2[i].otherPool = outIndex; // replace Id with the index + jn.otherChannels.push_back(i); + } + } } } @@ -1063,15 +1119,18 @@ void Dsolve::setNumAllVoxels( unsigned int num ) unsigned int Dsolve::convertIdToPoolIndex( const Id id ) const { + bool verbose = false; unsigned int i = id.value() - poolMapStart_; if ( i < poolMap_.size() ) { return poolMap_[i]; } - cout << "Warning: Dsolve::convertIdToPoollndex: Id out of range, (" << + if ( verbose ) { + cout << "Warning: Dsolve::convertIdToPoollndex: Id out of range, (" << poolMapStart_ << ", " << id << ", " << id.path() << ", " << poolMap_.size() + poolMapStart_ << "\n"; - return 0; + } + return ~0U; } unsigned int Dsolve::convertIdToPoolIndex( const Eref& e ) const @@ -1083,7 +1142,7 @@ void Dsolve::setN( const Eref& e, double v ) { unsigned int pid = convertIdToPoolIndex( e ); // Ignore silently, as this may be a valid pid for the ksolve to use. - if ( pid >= pools_.size() ) + if ( pid == ~0U || pid >= pools_.size() ) return; unsigned int vox = e.dataIndex(); if ( vox < numVoxels_ ) @@ -1098,7 +1157,7 @@ void Dsolve::setN( const Eref& e, double v ) double Dsolve::getN( const Eref& e ) const { unsigned int pid = convertIdToPoolIndex( e ); - if ( pid >= pools_.size() ) return 0.0; // ignore silently + if ( pid == ~0U || pid >= pools_.size() ) return 0.0; //ignore silently unsigned int vox = e.dataIndex(); if ( vox < numVoxels_ ) { @@ -1112,7 +1171,7 @@ double Dsolve::getN( const Eref& e ) const void Dsolve::setNinit( const Eref& e, double v ) { unsigned int pid = convertIdToPoolIndex( e ); - if ( pid >= pools_.size() ) // Ignore silently + if ( pid == ~0U || pid >= pools_.size() ) // Ignore silently return; unsigned int vox = e.dataIndex(); if ( vox < numVoxels_ ) @@ -1127,7 +1186,7 @@ void Dsolve::setNinit( const Eref& e, double v ) double Dsolve::getNinit( const Eref& e ) const { unsigned int pid = convertIdToPoolIndex( e ); - if ( pid >= pools_.size() ) return 0.0; // ignore silently + if ( pid == ~0U || pid >= pools_.size() ) return 0.0; //ignore silently unsigned int vox = e.dataIndex(); if ( vox < numVoxels_ ) { @@ -1141,36 +1200,52 @@ double Dsolve::getNinit( const Eref& e ) const void Dsolve::setDiffConst( const Eref& e, double v ) { unsigned int pid = convertIdToPoolIndex( e ); - if ( pid >= pools_.size() ) // Ignore silently, out of range. + if ( pid == ~0U || pid >= pools_.size() ) // Ignore silently, out of range. return; - pools_[ convertIdToPoolIndex( e ) ].setDiffConst( v ); + pools_[ pid ].setDiffConst( v ); } double Dsolve::getDiffConst( const Eref& e ) const { unsigned int pid = convertIdToPoolIndex( e ); - if ( pid >= pools_.size() ) // Ignore silently, out of range. + if ( pid == ~0U || pid >= pools_.size() ) // Ignore silently, out of range. return 0.0; - return pools_[ convertIdToPoolIndex( e ) ].getDiffConst(); + return pools_[ pid ].getDiffConst(); } void Dsolve::setMotorConst( const Eref& e, double v ) { unsigned int pid = convertIdToPoolIndex( e ); - if ( pid >= pools_.size() ) // Ignore silently, out of range. + if ( pid == ~0U || pid >= pools_.size() ) // Ignore silently, out of range. return; - pools_[ convertIdToPoolIndex( e ) ].setMotorConst( v ); + pools_[ pid ].setMotorConst( v ); } -void Dsolve::setNumPools( unsigned int numPoolSpecies ) +void Dsolve::setNumVarTotPools( unsigned int var, unsigned int tot ) { // Decompose numPoolSpecies here, assigning some to each node. - numTotPools_ = numPoolSpecies; - numLocalPools_ = numPoolSpecies; + numTotPools_ = tot; + numLocalPools_ = var; poolStartIndex_ = 0; - pools_.resize( numLocalPools_ ); - for ( unsigned int i = 0 ; i < numLocalPools_; ++i ) + pools_.resize( numTotPools_ ); + for ( unsigned int i = 0 ; i < numTotPools_; ++i ) + { + pools_[i].setNumVoxels( numVoxels_ ); + // pools_[i].setId( reversePoolMap_[i] ); + // pools_[i].setParent( me ); + } +} + +void Dsolve::setNumPools( unsigned int numVarPoolSpecies ) +{ + // Decompose numPoolSpecies here, assigning some to each node. + numTotPools_ = numVarPoolSpecies; + numLocalPools_ = numVarPoolSpecies; + poolStartIndex_ = 0; + + pools_.resize( numTotPools_ ); + for ( unsigned int i = 0 ; i < numTotPools_; ++i ) { pools_[i].setNumVoxels( numVoxels_ ); // pools_[i].setId( reversePoolMap_[i] ); diff --git a/diffusion/Dsolve.h b/diffusion/Dsolve.h index dc013d3cea..e6b911ff67 100644 --- a/diffusion/Dsolve.h +++ b/diffusion/Dsolve.h @@ -136,6 +136,8 @@ class Dsolve: public ZombiePoolInterface void setNumPools( unsigned int num ); unsigned int getNumPools() const; + void setNumVarTotPools( unsigned int var, unsigned int tot ); + unsigned int getNumLocalVoxels() const; VoxelPoolsBase* pools( unsigned int i ); double volume( unsigned int i ) const; @@ -185,6 +187,7 @@ class Dsolve: public ZombiePoolInterface void calcJnChan( const DiffJunction& jn, Dsolve* other, double dt ); void calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt ); + void calcLocalChan( double dt ); void fillConcChans( const vector< ObjId >& chans ); /** diff --git a/hsolve/HSolveInterface.cpp b/hsolve/HSolveInterface.cpp index 455c8cb7b9..05cd95c21e 100644 --- a/hsolve/HSolveInterface.cpp +++ b/hsolve/HSolveInterface.cpp @@ -89,8 +89,10 @@ void HSolve::setCm( Id id, double value ) { unsigned int index = localIndex( id ); assert( index < tree_.size() ); - // Also update data structures used for calculations. tree_[ index ].Cm = value; + // Also update data structures used for calculations. + assert( tree_.size() == compartment_.size() ); + compartment_[index].CmByDt = 2.0 * value / dt_; } double HSolve::getEm( Id id ) const @@ -104,8 +106,10 @@ void HSolve::setEm( Id id, double value ) { unsigned int index = localIndex( id ); assert( index < tree_.size() ); - // Also update data structures used for calculations. tree_[ index ].Em = value; + // Also update data structures used for calculations. + assert( tree_.size() == compartment_.size() ); + compartment_[index].EmByRm = value / tree_[index].Rm; } double HSolve::getRm( Id id ) const @@ -119,8 +123,10 @@ void HSolve::setRm( Id id, double value ) { unsigned int index = localIndex( id ); assert( index < tree_.size() ); - // Also update data structures used for calculations. tree_[ index ].Rm = value; + // Also update data structures used for calculations. + assert( tree_.size() == compartment_.size() ); + compartment_[index].EmByRm = tree_[index].Em / value; } double HSolve::getRa( Id id ) const diff --git a/kinetics/ReadKkit.cpp b/kinetics/ReadKkit.cpp index a33986ac76..366fe7809b 100644 --- a/kinetics/ReadKkit.cpp +++ b/kinetics/ReadKkit.cpp @@ -1192,7 +1192,8 @@ Id ReadKkit::buildChan( const vector< string >& args ) // double permeability = atof( args[ chanMap_["perm"] ].c_str() ); Id chan = shell_->doCreate( "ConcChan", pa, tail, 1 ); - Field< double >::set( chan, "permeability", permeability ); + // Convert from perm in uM in GENESIS, to mM for MOOSE. + Field< double >::set( chan, "permeability", permeability *1000.0 ); assert( chan != Id() ); string chanPath = clean.substr( 10 ); chanIds_[ chanPath ] = chan; diff --git a/ksolve/Gsolve.cpp b/ksolve/Gsolve.cpp index 4b31deec50..22c630b31b 100644 --- a/ksolve/Gsolve.cpp +++ b/ksolve/Gsolve.cpp @@ -1018,6 +1018,10 @@ void Gsolve::setNumPools( unsigned int numPoolSpecies ) } } +void Gsolve::setNumVarTotPools( unsigned int var, unsigned int tot ) { + setNumPools( tot ); +} + unsigned int Gsolve::getNumPools() const { if ( pools_.size() > 0 ) diff --git a/ksolve/Gsolve.h b/ksolve/Gsolve.h index 594b560751..71b74da791 100644 --- a/ksolve/Gsolve.h +++ b/ksolve/Gsolve.h @@ -46,6 +46,7 @@ class Gsolve: public ZombiePoolInterface */ void setNumPools( unsigned int num ); /// Inherited. unsigned int getNumPools() const; /// Inherited. + void setNumVarTotPools( unsigned int var, unsigned int tot );//Inherited VoxelPoolsBase* pools( unsigned int i ); /// Inherited. double volume( unsigned int i ) const; diff --git a/ksolve/Ksolve.cpp b/ksolve/Ksolve.cpp index e81b4d8f7a..f26c84c18f 100644 --- a/ksolve/Ksolve.cpp +++ b/ksolve/Ksolve.cpp @@ -726,6 +726,15 @@ void Ksolve::setNumPools( unsigned int numPoolSpecies ) } } +void Ksolve::setNumVarTotPools( unsigned int var, unsigned int tot ) +{ + unsigned int numVoxels = pools_.size(); + for ( unsigned int i = 0 ; i < numVoxels; ++i ) + { + pools_[i].resizeArrays( tot ); + } +} + unsigned int Ksolve::getNumPools() const { if ( pools_.size() > 0 ) diff --git a/ksolve/Ksolve.h b/ksolve/Ksolve.h index 89b87cd76b..fb810094f3 100644 --- a/ksolve/Ksolve.h +++ b/ksolve/Ksolve.h @@ -107,6 +107,7 @@ class Ksolve: public ZombiePoolInterface */ void setNumPools( unsigned int num ); unsigned int getNumPools() const; + void setNumVarTotPools( unsigned int var, unsigned int tot ); VoxelPoolsBase* pools( unsigned int i ); double volume( unsigned int i ) const; diff --git a/ksolve/Stoich.cpp b/ksolve/Stoich.cpp index 504f8ec10b..e23583c0bb 100644 --- a/ksolve/Stoich.cpp +++ b/ksolve/Stoich.cpp @@ -947,8 +947,8 @@ void Stoich::resizeArrays() N_.setSize( totNumPools, totNumRates ); if ( kinterface_ ) kinterface_->setNumPools( totNumPools ); - if ( dinterface_ ) // Only set up var pools managed locally. - dinterface_->setNumPools( varPoolVec_.size() ); + if ( dinterface_ ) // Need to set both the numVar and numTot + dinterface_->setNumVarTotPools( varPoolVec_.size(), totNumPools ); } /// Calculate sizes of all arrays, and allocate them. diff --git a/ksolve/ZombiePoolInterface.h b/ksolve/ZombiePoolInterface.h index 70d4995475..8cae5d5644 100644 --- a/ksolve/ZombiePoolInterface.h +++ b/ksolve/ZombiePoolInterface.h @@ -46,6 +46,9 @@ class ZombiePoolInterface /// gets number of pools (species) handled by system. virtual unsigned int getNumPools() const = 0; + /// Specifies number of var pools and total pools including buffered . + virtual void setNumVarTotPools( unsigned int var, unsigned int tot )=0; + /// Assign number of voxels (size of pools_ vector ) virtual void setNumAllVoxels( unsigned int numVoxels ) = 0; /// Number of voxels here. pools_.size() == getNumLocalVoxels diff --git a/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py index c1cb25ec96..5deea13e9b 100644 --- a/python/rdesigneur/rdesigneur.py +++ b/python/rdesigneur/rdesigneur.py @@ -460,6 +460,9 @@ def _buildVclampOnCompt( self, dendCompts, spineCompts, stimInfo ): stimObj = [] for i in dendCompts + spineCompts: vclamp = make_vclamp( name = 'vclamp', parent = i.path ) + + # Assume SI units. Scale by Cm to get reasonable gain. + vclamp.gain = i.Cm * 1e4 moose.connect( i, 'VmOut', vclamp, 'sensedIn' ) moose.connect( vclamp, 'currentOut', i, 'injectMsg' ) stimObj.append( vclamp ) @@ -1267,6 +1270,7 @@ def _assignComptNamesFromKkit( self ): #print( "SortedClist= {}".format( sortedComptList[i] )) if sortedComptList[i].name != sortedNames[i]: sortedComptList[i].name = sortedNames[i] + #print( "SortedClist= {} {}".format( sortedNames[i], sortedComptList[i].volume )) return sortedComptList diff --git a/python/rdesigneur/rdesigneurProtos.py b/python/rdesigneur/rdesigneurProtos.py index 991bf4e8f3..b71ce8e4d3 100644 --- a/python/rdesigneur/rdesigneurProtos.py +++ b/python/rdesigneur/rdesigneurProtos.py @@ -240,7 +240,6 @@ def make_Ca( name ): yA[i] = 5.0 * math.exp( -50 * (x - EREST_ACT) ) else: yA[i] = 5.0 - #yB[i] = 6.0 - yA[i] yB[i] = 5.0 x += dx ygate.tableA = yA diff --git a/tests/python/test_Xchan1.py b/tests/python/test_Xchan1.py index d1eb35114c..0eda4dfad1 100644 --- a/tests/python/test_Xchan1.py +++ b/tests/python/test_Xchan1.py @@ -63,7 +63,7 @@ def makeModel(): s.concInit = 0.001 chanPool.nInit = 1000.0 # Flux (mM/s) = permeability * N * (conc_out - conc_in ) - chan.permeability = 0.1 * chanPool.volume * 6.022e23 *0.001 / chanPool.nInit + chan.permeability = 0.1 * chanPool.volume * 6.022e23 / chanPool.nInit ##################################################################### fixXreacs.fixXreacs( '/model' )