diff --git a/moose-core/biophysics/DifBufferBase.cpp b/moose-core/biophysics/DifBufferBase.cpp
index 5c6c43a20dd7034673fa8b0e10fe31fde7a5a839..1be36d79e2b1aeae6fbe98fe0cae9109eaaa554f 100644
--- a/moose-core/biophysics/DifBufferBase.cpp
+++ b/moose-core/biophysics/DifBufferBase.cpp
@@ -397,65 +397,3 @@ void DifBufferBase:: fluxFromIn(const Eref& e,double innerC, double innerThickne
{
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/moose-core/biophysics/DifBufferBase.h b/moose-core/biophysics/DifBufferBase.h
index e1bc23b616caa153d9cde8affe5b8808fb0fa1d4..00618d8c5349863cebd15dae2a803739c80a8aed 100644
--- a/moose-core/biophysics/DifBufferBase.h
+++ b/moose-core/biophysics/DifBufferBase.h
@@ -119,9 +119,7 @@ public:
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:
diff --git a/moose-core/biophysics/DifShellBase.cpp b/moose-core/biophysics/DifShellBase.cpp
index e6d1e5338b9242c03fcf6a7c49bc475faeabbea9..8577a767c4d9dbbd6fef4ba549e30a7c4159e530 100644
--- a/moose-core/biophysics/DifShellBase.cpp
+++ b/moose-core/biophysics/DifShellBase.cpp
@@ -453,61 +453,3 @@ void DifShellBase::hillPump(const Eref& e,
{
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/moose-core/biophysics/DifShellBase.h b/moose-core/biophysics/DifShellBase.h
index db74c53696c1d133d4bfbf96432da1e6a46b4ccb..57b208105bdf0fb5dcdb2e6ce4d7e0d8dc0d146c 100644
--- a/moose-core/biophysics/DifShellBase.h
+++ b/moose-core/biophysics/DifShellBase.h
@@ -127,9 +127,7 @@ class DifShellBase
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();
diff --git a/moose-core/hsolve/HSolve.h b/moose-core/hsolve/HSolve.h
index 1b26b856c91f2d373de8521c0ce2e362ae6ef023..3db7d2ad006ebeec6f92e81e2ace903db1e0e568 100644
--- a/moose-core/hsolve/HSolve.h
+++ b/moose-core/hsolve/HSolve.h
@@ -87,7 +87,7 @@ public:
//~ const vector< Id >& getCompartments() const;
void addGkEk( Id id, double v1, double v2 );
-
+ void addConc( Id id, double conc );
/// Interface to channels
//~ const vector< Id >& getHHChannels() const;
void setPowers(
diff --git a/moose-core/hsolve/HSolveActive.cpp b/moose-core/hsolve/HSolveActive.cpp
index 004392eb480711a594386b9f89f4c9df5e56e600..2e80688339e3388b8683678a3860e185b8f4e7c6 100644
--- a/moose-core/hsolve/HSolveActive.cpp
+++ b/moose-core/hsolve/HSolveActive.cpp
@@ -18,6 +18,7 @@
#include "../biophysics/CompartmentBase.h"
#include "../biophysics/Compartment.h"
#include "../biophysics/CaConcBase.h"
+#include "../biophysics/ChanBase.h"
#include "ZombieCaConc.h"
using namespace moose;
//~ #include "ZombieCompartment.h"
@@ -227,9 +228,12 @@ void HSolveActive::advanceChannels( double dt )
vector< LookupColumn >::iterator icolumn = column_.begin();
vector< LookupRow >::iterator icarowcompt;
vector< LookupRow* >::iterator icarow = caRow_.begin();
-
+ vector< double >::iterator iextca = externalCalcium_.begin();
+
LookupRow vRow;
+ LookupRow dRow;
double C1, C2;
+
for ( iv = V_.begin(); iv != V_.end(); ++iv )
{
vTable_.row( *iv, vRow );
@@ -238,6 +242,7 @@ void HSolveActive::advanceChannels( double dt )
for ( ; ica < caBoundary; ++ica )
{
caTable_.row( *ica, *icarowcompt );
+
++icarowcompt;
}
@@ -252,6 +257,9 @@ void HSolveActive::advanceChannels( double dt )
chanBoundary = ichan + *ichannelcount;
for ( ; ichan < chanBoundary; ++ichan )
{
+
+ caTable_.row( *iextca, dRow );
+
if ( ichan->Xpower_ > 0.0 )
{
vTable_.lookup( *icolumn, vRow, C1, C2 );
@@ -279,21 +287,29 @@ void HSolveActive::advanceChannels( double dt )
{
double temp = 1.0 + dt / 2.0 * C2;
*istate = ( *istate * ( 2.0 - temp ) + dt * C1 ) / temp;
- }
-
+
+}
++icolumn, ++istate;
}
if ( ichan->Zpower_ > 0.0 )
{
LookupRow* caRow = *icarow;
+
if ( caRow )
{
caTable_.lookup( *icolumn, *caRow, C1, C2 );
+
}
- else
+ else if (*iextca >0)
+
+ {
+ caTable_.lookup( *icolumn, dRow, C1, C2 );
+ }
+ else
{
- vTable_.lookup( *icolumn, vRow, C1, C2 );
+ vTable_.lookup( *icolumn, vRow, C1, C2 );
+
}
//~ *istate = *istate * C1 + C2;
@@ -307,7 +323,9 @@ void HSolveActive::advanceChannels( double dt )
}
++icolumn, ++istate, ++icarow;
+
}
+ ++iextca;
}
++ichannelcount, ++icacount;
@@ -339,12 +357,24 @@ void HSolveActive::sendValues( ProcPtr info )
vector< unsigned int >::iterator i;
for ( i = outVm_.begin(); i != outVm_.end(); ++i )
- moose::Compartment::VmOut()->send(
+ Compartment::VmOut()->send(
//~ ZombieCompartment::VmOut()->send(
compartmentId_[ *i ].eref(),
V_[ *i ]
);
+
+ for ( i = outIk_.begin(); i != outIk_.end(); ++i ){
+
+ unsigned int comptIndex = chan2compt_[ *i ];
+
+ assert( comptIndex < V_.size() );
+
+ ChanBase::IkOut()->send(channelId_[*i].eref(),
+ (current_[ *i ].Ek - V_[ comptIndex ]) * current_[ *i ].Gk);
+
+ }
+
for ( i = outCa_.begin(); i != outCa_.end(); ++i )
//~ CaConc::concOut()->send(
CaConcBase::concOut()->send(
diff --git a/moose-core/hsolve/HSolveActive.h b/moose-core/hsolve/HSolveActive.h
index eff1e59220a1bd37cede6bd7e1b816dc52e3dc5f..b17741436257b425dd8309e3d400850ad054ebcd 100644
--- a/moose-core/hsolve/HSolveActive.h
+++ b/moose-core/hsolve/HSolveActive.h
@@ -119,6 +119,7 @@ protected:
vector< double > externalCurrent_; ///< External currents from
///< channels that HSolve
///< cannot internalize.
+ vector< double > externalCalcium_; /// calcium from difshells
vector< Id > caConcId_; ///< Used for localIndex-ing.
vector< Id > channelId_; ///< Used for localIndex-ing.
vector< Id > gateId_; ///< Used for localIndex-ing.
@@ -131,6 +132,7 @@ protected:
* Tells you which compartments have external calcium-dependent
* channels so that you can send out Calcium concentrations in only
* those compartments. */
+ vector< unsigned int > outIk_;
private:
/**
diff --git a/moose-core/hsolve/HSolveActiveSetup.cpp b/moose-core/hsolve/HSolveActiveSetup.cpp
index b7d649958b3496388bc1558499115ccb3910545a..e9a2720788d74a968afd8c772df042d66f77e415 100644
--- a/moose-core/hsolve/HSolveActiveSetup.cpp
+++ b/moose-core/hsolve/HSolveActiveSetup.cpp
@@ -88,23 +88,32 @@ void HSolveActive::reinitChannels()
vector< LookupColumn >::iterator icolumn = column_.begin();
vector< LookupRow >::iterator icarowcompt;
vector< LookupRow* >::iterator icarow = caRow_.begin();
-
+ vector< double>::iterator iextca = externalCalcium_.begin();
+
LookupRow vRow;
+ LookupRow dRow;
double C1, C2;
for ( iv = V_.begin(); iv != V_.end(); ++iv )
{
+
vTable_.row( *iv, vRow );
icarowcompt = caRowCompt_.begin();
+
caBoundary = ica + *icacount;
for ( ; ica < caBoundary; ++ica )
{
caTable_.row( *ica, *icarowcompt );
+
++icarowcompt;
+
}
chanBoundary = ichan + *ichannelcount;
for ( ; ichan < chanBoundary; ++ichan )
{
+
+ caTable_.row( *iextca, dRow );
+
if ( ichan->Xpower_ > 0.0 )
{
vTable_.lookup( *icolumn, vRow, C1, C2 );
@@ -126,10 +135,13 @@ void HSolveActive::reinitChannels()
if ( ichan->Zpower_ > 0.0 )
{
LookupRow* caRow = *icarow;
+
if ( caRow )
{
caTable_.lookup( *icolumn, *caRow, C1, C2 );
}
+ else if (*iextca>0)
+ caTable_.lookup( *icolumn, dRow, C1, C2 );
else
{
vTable_.lookup( *icolumn, vRow, C1, C2 );
@@ -139,6 +151,8 @@ void HSolveActive::reinitChannels()
++icolumn, ++istate, ++icarow;
}
+
+ ++iextca;
}
++ichannelcount, ++icacount;
@@ -241,16 +255,18 @@ void HSolveActive::readCalcium()
vector< Id > caConcId;
vector< int > caTargetIndex;
map< Id, int > caConcIndex;
- int nTarget, nDepend;
+ int nTarget, nDepend = 0;
vector< Id >::iterator iconc;
caCount_.resize( nCompt_ );
unsigned int ichan = 0;
+
for ( unsigned int ic = 0; ic < nCompt_; ++ic )
{
+
unsigned int chanBoundary = ichan + channelCount_[ ic ];
unsigned int nCa = caConc_.size();
-
+
for ( ; ichan < chanBoundary; ++ichan )
{
caConcId.clear();
@@ -261,16 +277,21 @@ void HSolveActive::readCalcium()
caTargetIndex.push_back( -1 );
nDepend = HSolveUtils::caDepend( channelId_[ ichan ], caConcId );
- if ( nDepend == 0 )
+
+ if ( nDepend == 0)
// Channel does not depend on calcium.
- caDependIndex_.push_back( -1 );
+ caDependIndex_.push_back( -1 );
+
+
+ externalCalcium_.push_back(0);
+
for ( iconc = caConcId.begin(); iconc != caConcId.end(); ++iconc )
if ( caConcIndex.find( *iconc ) == caConcIndex.end() )
{
caConcIndex[ *iconc ] = caCount_[ ic ];
++caCount_[ ic ];
-
+
Ca = Field< double >::get( *iconc, "Ca" );
CaBasal = Field< double >::get( *iconc, "CaBasal" );
tau = Field< double >::get( *iconc, "tau" );
@@ -293,9 +314,13 @@ void HSolveActive::readCalcium()
caTargetIndex.push_back( caConcIndex[ caConcId.front() ] + nCa );
if ( nDepend != 0 )
caDependIndex_.push_back( caConcIndex[ caConcId.back() ] );
+
+
}
}
-
+
+
+
caTarget_.resize( channel_.size() );
ca_.resize( caConc_.size() );
caActivation_.resize( caConc_.size() );
@@ -307,6 +332,8 @@ void HSolveActive::readCalcium()
else
caTarget_[ ichan ] = &caActivation_[ caTargetIndex[ ichan ] ];
}
+
+
}
void HSolveActive::createLookupTables()
@@ -418,8 +445,8 @@ void HSolveActive::createLookupTables()
// *ia = ( 2.0 - dt_ * b ) / ( 2.0 + dt_ * b );
// *ib = dt_ * a / ( 1.0 + dt_ * b / 2.0 );
- // *ia = dt_ * a;
- // *ib = 1.0 + dt_ * b / 2.0;
+ // *ia = dt_ * a;
+ // *ib = 1.0 + dt_ * b / 2.0;
++ia, ++ib;
}
@@ -435,6 +462,7 @@ void HSolveActive::createLookupTables()
//~ for ( int igrid = 0; igrid <= vDiv_; ++igrid )
//~ grid[ igrid ] = vMin_ + igrid * dv;
//~ }
+
for ( unsigned int ig = 0; ig < vGate.size(); ++ig )
{
@@ -467,7 +495,7 @@ void HSolveActive::createLookupTables()
for ( unsigned int ig = 0; ig < gateId_.size(); ++ig )
{
unsigned int species = gateSpecies[ gateId_[ ig ] ];
-
+
LookupColumn column;
if ( gCaDepend_[ ig ] )
caTable_.column( species, column );
@@ -480,23 +508,27 @@ void HSolveActive::createLookupTables()
///////////////////!!!!!!!!!!
unsigned int maxN = *( max_element( caCount_.begin(), caCount_.end() ) );
caRowCompt_.resize( maxN );
+
for ( unsigned int ichan = 0; ichan < channel_.size(); ++ichan )
{
if ( channel_[ ichan ].Zpower_ > 0.0 )
{
int index = caDependIndex_[ ichan ];
+
if ( index == -1 )
caRow_.push_back( 0 );
else
caRow_.push_back( &caRowCompt_[ index ] );
}
}
+
+
}
/**
* Reads in SynChans and SpikeGens.
- *
- * Unlike Compartments, HHChannels, etc., neither of these are zombified.
+ *
+* Unlike Compartments, HHChannels, etc., neither of these are zombified.
* In other words, their fields are not managed by HSolve, and their "process"
* functions are invoked to do their calculations. For SynChans, the process
* calls are made by their respective clocks, and hence the process message is
@@ -615,6 +647,30 @@ void HSolveActive::manageOutgoingMessages()
if ( nTargets )
outCa_.push_back( ica );
}
+
+ filter.clear();
+ filter.push_back( "CaConc" );
+
+ for ( unsigned int ik = 0; ik < channelId_.size(); ++ik )
+ {
+ targets.clear();
+
+ int nTargets = HSolveUtils::targets(
+ channelId_[ ik ],
+ "IkOut",
+ targets,
+ filter,
+ false // include = false. That is, use filter to exclude.
+ );
+
+ if ( nTargets )
+ outIk_.push_back( ik );
+
+
+ }
+
+
+
}
void HSolveActive::cleanup()
diff --git a/moose-core/hsolve/HSolveInterface.cpp b/moose-core/hsolve/HSolveInterface.cpp
index 4499b0062de94670e5c82753359ed21cb8fc07dc..30b39db19f52a81009dea6a827e28f283a798697 100644
--- a/moose-core/hsolve/HSolveInterface.cpp
+++ b/moose-core/hsolve/HSolveInterface.cpp
@@ -243,6 +243,14 @@ void HSolve::addGkEk( Id id, double Gk, double Ek )
externalCurrent_[ 2 * index + 1 ] += Gk * Ek;
}
+void HSolve::addConc( Id id, double conc )
+{
+ unsigned int index = localIndex( id );
+ assert( index + 1 < externalCalcium_.size() );
+ externalCalcium_[ index ] = conc;
+}
+
+
void HSolve::setPowers(
Id id,
double Xpower,
@@ -280,7 +288,7 @@ void HSolve::setHHChannelGbar( Id id, double value )
unsigned int index = localIndex( id );
assert( index < channel_.size() );
channel_[ index ].Gbar_ = value;
- // cout << "HSolve::setHHChannelGbar( " << id.path() << ", " << value << " ), index = " << index << endl;
+
}
double HSolve::getEk( Id id ) const
diff --git a/moose-core/hsolve/Makefile b/moose-core/hsolve/Makefile
index da4e33a3a96c04e163d07ea2c2c454f08f9787cd..48ad098f6c7c1e99a9442d35371f90707f6fd933 100644
--- a/moose-core/hsolve/Makefile
+++ b/moose-core/hsolve/Makefile
@@ -35,7 +35,7 @@ HSolveStruct.o: HSolveStruct.h
HinesMatrix.o: HinesMatrix.h TestHSolve.h
HSolvePassive.o: HSolvePassive.h HinesMatrix.h HSolveStruct.h HSolveUtils.h TestHSolve.h ../biophysics/Compartment.h
RateLookup.o: RateLookup.h
-HSolveActive.o: HSolveActive.h RateLookup.h HSolvePassive.h HinesMatrix.h HSolveStruct.h
+HSolveActive.o: HSolveActive.h RateLookup.h HSolvePassive.h HinesMatrix.h HSolveStruct.h ../biophysics/ChanBase.h
HSolveActiveSetup.o: HSolveActive.h RateLookup.h HSolvePassive.h HinesMatrix.h HSolveStruct.h HSolveUtils.h ../biophysics/HHChannelBase.h ../biophysics/HHChannel.h ../biophysics/ChanBase.h ../biophysics/ChanCommon.h ../biophysics/HHGate.h ../biophysics/CaConc.h
HSolveInterface.o: HSolve.h HSolveActive.h RateLookup.h HSolvePassive.h HinesMatrix.h HSolveStruct.h
HSolve.o: ../biophysics/Compartment.h ZombieCompartment.h ../biophysics/CaConc.h ZombieCaConc.h ../biophysics/HHGate.h ../biophysics/ChanBase.h ../biophysics/ChanCommon.h ../biophysics/HHChannelBase.h ../biophysics/HHChannel.h ZombieHHChannel.h HSolve.h HSolveActive.h RateLookup.h HSolvePassive.h HinesMatrix.h HSolveStruct.h ../basecode/ElementValueFinfo.h
diff --git a/moose-core/hsolve/ZombieHHChannel.cpp b/moose-core/hsolve/ZombieHHChannel.cpp
index d72cd6e3d8ff76ac74e7087a849ac3599bdae4a7..7876f9e991e6109f95f7bba8ea65e47c3185d907 100644
--- a/moose-core/hsolve/ZombieHHChannel.cpp
+++ b/moose-core/hsolve/ZombieHHChannel.cpp
@@ -177,7 +177,9 @@ void ZombieHHChannel::vReinit( const Eref& er, ProcPtr info )
void ZombieHHChannel::vHandleConc( const Eref& e, double conc )
{
- ;// cout << "Warning: ZombieHHChannel::vHandleConc\n";
+
+ hsolve_->addConc(e.id(),conc);
+
}
void ZombieHHChannel::vCreateGate(const Eref& e, string name)
diff --git a/moose-core/tests/python/chan_proto.py b/moose-core/tests/python/chan_proto.py
index ae37024f5de98402d8df87822af6258aaeef7b54..52dceee92812e7e032c0e9dbbe811d0bf27cc368 100644
--- a/moose-core/tests/python/chan_proto.py
+++ b/moose-core/tests/python/chan_proto.py
@@ -12,36 +12,37 @@ If no inactivation, just send in empty Yparam array.
from __future__ import print_function, division
import moose
import numpy as np
-from util import NamedList
-
-SSTauChannelParams = NamedList('SSTauChannelParams', '''
- Arate
- A_B
- A_C
- Avhalf
- Avslope
- taumin
- tauVdep
- tauPow
- tauVhalf
- tauVslope''')
-
-AlphaBetaChannelParams = NamedList('AlphaBetaChannelParams', '''
- A_rate
- A_B
- A_C
- Avhalf
- A_vslope
- B_rate
- B_B
- B_C
- Bvhalf
- B_vslope''')
-
-ZChannelParams = NamedList('ZChannelParams', 'Kd power tau')
-
-
-ChannelSettings = NamedList('ChannelSettings', 'Xpow Ypow Zpow Erev name')
+from collections import namedtuple
+
+ChannelSettings = namedtuple('''ChannelSettings''', '''
+Xpow
+Ypow
+Zpow
+Erev
+name''')
+ZChannelParams = namedtuple('''ZChannelParams''',
+ '''Kd
+ power
+ tau''')
+
+ChannelParams = namedtuple('''ChannelParams''', '''channel
+X
+Y
+Z
+''')
+
+GateParams = namedtuple('''GateParams''', '''
+A_rate
+A_B
+A_C
+Avhalf
+A_vslope
+B_rate
+B_B
+B_C
+Bvhalf
+B_vslope''')
+
def interpolate_values_in_table(tabA,V_0,l=40):
import param_chan
@@ -92,13 +93,13 @@ def chan_proto(chanpath, params):
chan.Xpower = params.channel.Xpow
if params.channel.Xpow > 0:
xGate = moose.HHGate(chan.path + '/gateX')
- xGate.setupAlpha(params.X + [param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX])
+ xGate.setupAlpha(params.X + (param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX))
xGate = fix_singularities(params.X, xGate)
chan.Ypower = params.channel.Ypow
if params.channel.Ypow > 0:
yGate = moose.HHGate(chan.path + '/gateY')
- yGate.setupAlpha(params.Y + [param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX])
+ yGate.setupAlpha(params.Y + (param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX))
yGate = fix_singularities(params.Y, yGate)
if params.channel.Zpow > 0:
chan.Zpower = params.channel.Zpow
@@ -116,25 +117,3 @@ def chan_proto(chanpath, params):
return chan
-TypicalOneDalpha = NamedList('TypicalOneDalpha',
- '''channel X Y Z=[] calciumPermeable=False calciumPermeable2=False''')
-
-_FUNCTIONS = {
- TypicalOneDalpha: chan_proto,
-
-}
-
-#*params... passes the set of values not as a list but as individuals
-def make_channel(chanpath, params):
- func = _FUNCTIONS[params.__class__]
- return func(chanpath, params)
-
-def chanlib():
- import param_chan
- if not moose.exists('/library'):
- moose.Neutral('/library')
- #Adding all the channels to the library. *list removes list elements from the list,
- #so they can be used as function arguments
- chan = [make_channel('/library/'+key, value) for key, value in param_chan.ChanDict.items()]
-
- return chan
diff --git a/moose-core/tests/python/param_chan.py b/moose-core/tests/python/param_chan.py
index 0db3f697226b9e30aa0a766d5e4b41ad8255cad8..a6531f97093036c77b541d7cd4e3737123157888 100644
--- a/moose-core/tests/python/param_chan.py
+++ b/moose-core/tests/python/param_chan.py
@@ -1,25 +1,6 @@
# -*- coding: utf-8 -*-
+from chan_proto import (ChannelSettings, ZChannelParams,ChannelParams,GateParams)
-from util import NamedDict
-from chan_proto import (
- SSTauChannelParams,
- AlphaBetaChannelParams,
- ZChannelParams,
- ChannelSettings,
- TypicalOneDalpha,)
-
-#chanDictSP.py
-#contains all gating parameters and reversal potentials
-# Gate equations have the form:
-#
-# y(x) = (rate + B * x) / (C + exp((x + vhalf) / vslope))
-#
-# OR
-# y(x) = tau_min + tau_vdep / (1 + exp((x + vhalf) / vslope))
-#
-# where x is membrane voltage and y is the rate constant
-#KDr params used by Sriram, RE paper1, Krp params used by RE paper 2
-#Parameters for Ca channels may need to be shifted - see Dorman model
krev=-87e-3
narev=50e-3
carev=140e-3 #assumes CaExt=2 mM and CaIn=50e-3
@@ -32,49 +13,16 @@ CAMIN=0.01e-3 #10 nM
CAMAX=40e-3 #40 uM, might want to go up to 100 uM with spines
CADIVS=3000 #10 nM steps
-#mtau: Ogata fig 5, no qfactor accounted in mtau, 1.2 will improve spike shape
-#activation minf fits Ogata 1990 figure 3C (which is cubed root)
-#inactivation hinf fits Ogata 1990 figure 6B
-#htau fits the main -50 through -10 slope of Ogata figure 9 (log tau), but a qfact of 2 is already taken into account.
-
-
-CaTparam = ChannelSettings(Xpow=3, Ypow=1, Zpow=1, Erev=carev, name='CaT')
-qfactCaT = 1
+CaL12param = ChannelSettings(Xpow=1, Ypow=0, Zpow=0, Erev=0.048, name='CaL12')
+qfactCaL = 1
+CaL12_X_params = GateParams(A_rate=-880,A_B= -220e3,A_C = -1.0,Avhalf= 4.0003e-3, A_vslope=-7.5e-3,B_rate=-284,B_B=71e3,B_C=-1.0,Bvhalf=-4.0003e-3,B_vslope=5e-3)
-#Original inactivation ws too slow compared to activation, made closder the alpha1G
-CaT_X_params = AlphaBetaChannelParams(A_rate = 5342.4*qfactCaT,
- A_B = 2100.*qfactCaT,
- A_C = 11.9,
- Avhalf = 1e-3,
- A_vslope = -12e-3,
- B_rate = 289.7*qfactCaT,
- B_B = 0.*qfactCaT,
- B_C = 1.,
- Bvhalf = 0.0969,
- B_vslope = 0.0141)
-#CaT_Y_params = CaT_X_params
-CaT_Y_params = AlphaBetaChannelParams(A_rate = 0*qfactCaT,
- A_B = -74.*qfactCaT,
- A_C = 1.,
- Avhalf = 9e-2,
- A_vslope = 5e-3,
- B_rate = 15*qfactCaT,
- B_B = -1.5*qfactCaT,
- B_C = 1.5,
- Bvhalf = 0.05,
- B_vslope = -15e-3)
+SKparam = ChannelSettings(Xpow=0, Ypow=0, Zpow=1, Erev=-.087, name='SKCa')
-#These CDI params can be used with every channel, make ZpowCDI=2
-#If ZpowCDI=0 the CDI will not be used, power=-4 is to transform
-#(Ca/Kd)^pow/(1+(Ca/Kd)^pow) to 1/(1+(ca/Kd)^-pow)
-CDI_Z_params = ZChannelParams(Kd = 0.12e-3,
- power = -4,
- tau = 142e-3)
+SK_Z_params = ZChannelParams(Kd = 0.57e-3,
+ power = 5.2,
+ tau = 4.9e-3)
-#Dictionary of "standard" channels, to create channels using a loop
-#NaF doesn't fit since it uses different prototype form
-#will need separate dictionary for BK
-ChanDict = NamedDict(
- 'ChannelParams',
- CaT = TypicalOneDalpha(CaTparam, CaT_X_params, CaT_Y_params, CDI_Z_params, calciumPermeable=True))
+Cal = ChannelParams(channel=CaL12param,X=CaL12_X_params,Y=[],Z=[])
+SK = ChannelParams(channel=SKparam,X=[],Y=[],Z=SK_Z_params)
diff --git a/moose-core/tests/python/soma.p b/moose-core/tests/python/soma.p
new file mode 100644
index 0000000000000000000000000000000000000000..01e87c41384c074a3534b7a895de3bff77a3f98a
--- /dev/null
+++ b/moose-core/tests/python/soma.p
@@ -0,0 +1,15 @@
+
+*relative
+*cartesian
+*asymmetric
+
+*set_global RA 4.0
+//change Cm to account for no spines - make 3x higher?
+*set_global CM 0.03
+*set_global RM 2.8
+*set_global EREST_ACT -80e-3
+*set_global ELEAK -50e-3
+
+*start_cell
+soma none 16.000 0 0 16.000
+
diff --git a/moose-core/tests/python/test_difshells.py b/moose-core/tests/python/test_difshells.py
index 669f46a891717c87cac7892988f86b652af8af05..9f4f35e99488066a3a3e4f292effbc4576e51606 100644
--- a/moose-core/tests/python/test_difshells.py
+++ b/moose-core/tests/python/test_difshells.py
@@ -2,7 +2,7 @@ import moose
import numpy as np
import matplotlib.pyplot as plt
import chan_proto
-
+import param_chan
from params import *
@@ -47,13 +47,13 @@ def add_difbuffer_to_dishell(comp, i, j, shell_thickness, shell_radius):
return buf
-def add_difshells_and_buffers(comp):
+def add_difshells_and_buffers(comp,difshell_no,difbuff_no):
if difshell_no < 1:
return [], []
- difshell = []
- shell_thickness = dend.diameter / difshell_no / 2.
+ difshell = []
+ shell_thickness = comp.diameter / difshell_no / 2.
difbuffer = []
for i in range(difshell_no):
shell_radius = comp.diameter / 2 - i * shell_thickness
@@ -84,7 +84,7 @@ def add_difshells_and_buffers(comp):
return difshell, difbuffer
-def addOneChan(chanpath, gbar, comp):
+def addOneChan(chanpath, gbar,comp):
SA = np.pi * comp.length * comp.diameter
proto = moose.element('/library/' + chanpath)
@@ -97,7 +97,7 @@ def addOneChan(chanpath, gbar, comp):
if __name__ == '__main__':
-
+ lib = moose.Neutral('/library')
for tick in range(0, 7):
moose.setClock(tick, dt)
moose.setClock(8, 0.005) # set output clock
@@ -120,19 +120,21 @@ if __name__ == '__main__':
pulse.width[0] = 500e-3
pulse.level[0] = inject
pulse.delay[1] = 1e9
-
+
+ chan = chan_proto.chan_proto('/library/CaL12',param_chan.Cal)
+
m = moose.connect(pulse, 'output', dend, 'injectMsg')
moose.connect(vmtab, 'requestOut', dend, 'getVm')
- chan_proto.chanlib()
- chan = addOneChan('CaT', gbar, dend)
+
+ chan = addOneChan('CaL12', gbar, dend)
moose.connect(gktab, 'requestOut', chan, 'getGk')
moose.connect(iktab, 'requestOut', chan, 'getIk')
diftab = []
buftab = []
- difs, difb = add_difshells_and_buffers(dend)
+ difs, difb = add_difshells_and_buffers(dend,difshell_no,difbuff_no)
if pumps:
pump = moose.MMPump('/model/dend/pump')
pump.Vmax = kcat
diff --git a/moose-core/tests/python/test_hsolve_externalCalcium.py b/moose-core/tests/python/test_hsolve_externalCalcium.py
new file mode 100644
index 0000000000000000000000000000000000000000..1a831b607abc3d146f1d0a081bc138410192f1b0
--- /dev/null
+++ b/moose-core/tests/python/test_hsolve_externalCalcium.py
@@ -0,0 +1,69 @@
+import sys
+import numpy as np
+import time
+import moose
+import matplotlib.pyplot as plt
+import test_difshells as td
+import chan_proto
+import param_chan
+print('Using moose from %s' % moose.__file__ )
+
+difshell_no = 3
+difbuf_no = 0
+
+
+p_file = "soma.p"
+cond = {'CaL12':30*0.35e-5,
+ 'SK':0.5*0.35e-6}
+
+
+if __name__ =='__main__':
+ for tick in range(0, 7):
+ moose.setClock(tick,10e-6)
+ moose.setClock(8, 0.005)
+
+ lib = moose.Neutral('/library')
+ model = moose.loadModel(p_file,'neuron')
+ pulse = moose.PulseGen('/neuron/pulse')
+ inject = 100e-10
+ chan_proto.chan_proto('/library/SK',param_chan.SK)
+ chan_proto.chan_proto('/library/CaL12',param_chan.Cal)
+ pulse.delay[0] = 8.
+ pulse.width[0] = 500e-12
+ pulse.level[0] = inject
+ pulse.delay[1] = 1e9
+
+ for comp in moose.wildcardFind('/neuron/#[TYPE=Compartment]'):
+ new_comp = moose.element(comp)
+ new_comp.initVm = -.08
+ difs, difb = td.add_difshells_and_buffers(new_comp,difshell_no,difbuf_no)
+ for name in cond:
+ chan = td.addOneChan(name,cond[name],new_comp)
+ if 'Ca' in name:
+ moose.connect(chan, "IkOut", difs[0], "influx")
+
+ if 'SK' in name:
+ moose.connect(difs[0], 'concentrationOut', chan, 'concen')
+
+
+
+ data = moose.Neutral('/data')
+ vmtab = moose.Table('/data/Vm')
+ shelltab = moose.Table('/data/Ca')
+ caltab = moose.Table('/data/CaL_Gk')
+ sktab = moose.Table('/data/SK_Gk')
+ moose.connect(vmtab, 'requestOut',moose.element('/neuron/soma') , 'getVm')
+ moose.connect(shelltab, 'requestOut', difs[0], 'getC')
+ moose.connect(caltab,'requestOut',moose.element('/neuron/soma/CaL12') ,'getGk')
+ moose.connect(sktab,'requestOut', moose.element('/neuron/soma/SK'),'getGk')
+
+ hsolve = moose.HSolve('/neuron/hsolve')
+ hsolve.dt = 10e-6
+ hsolve.target = ('/neuron/soma')
+ t_stop = 10.
+ moose.reinit()
+ moose.start(t_stop)
+ print len(np.where(sktab.vector<1e-19)[0]), len(np.where(shelltab.vector>50e-6)[0])
+
+
+