diff --git a/include/mechanisms/hh_verb.hpp b/include/mechanisms/hh_verb.hpp
index c3999107115f7d21c1c6ac4a70c616796a016e2d..8b67a69d56a0f61b2db8baaa37d03fc34d7e9687 100644
--- a/include/mechanisms/hh_verb.hpp
+++ b/include/mechanisms/hh_verb.hpp
@@ -73,11 +73,10 @@ public:
htau = data_(14*field_size, 15*size());
// set initial values for variables and parameters
- std::fill(gnabar.data(), gnabar.data()+size(), 1.2);
- std::fill(gl.data(), gl.data()+size(), 0.0029999999999999997);
- std::fill(gkbar.data(), gkbar.data()+size(), 0.35999999999999997);
- std::fill(el.data(), el.data()+size(), -54.299999999999997);
-
+ std::fill(gnabar.data(), gnabar.data()+size(), 0.12);
+ std::fill(gkbar.data(), gkbar.data()+size(), 0.036);
+ std::fill(gl.data(), gl.data()+size(), 0.0003);
+ std::fill(el.data(), el.data()+size(), -54.3);
}
using base::size;
@@ -138,23 +137,24 @@ public:
m[i_] = minf[i_]+(m[i_]-minf[i_])*exp( -dt/mtau[i_]);
h[i_] = hinf[i_]+(h[i_]-hinf[i_])*exp( -dt/htau[i_]);
n[i_] = ninf[i_]+(n[i_]-ninf[i_])*exp( -dt/ntau[i_]);
+ printf("m h n : %18.14f %18.14f %18.14f\n", m[i_], h[i_], n[i_]);
}
}
void rates(const int i_, value_type v) {
value_type ll3_, ll1_, ll0_, alpha, beta, ll2_, sum, q10;
q10 = std::pow( 3, (celsius- 6.2999999999999998)/ 10);
- //ll2_ = -v+ 40;
+ printf("q10 %18.14f\n", q10);
ll2_ = -(v+ 40);
ll0_ = ll2_/(exp(ll2_/ 10)- 1);
alpha = 0.10000000000000001*ll0_;
beta = 4*exp( -(v+ 65)/ 18);
- std::cout << "v " << v << "\n";
+ //std::cout << "v " << v << " dt " << dt << "\n";
//std::cout << "m (alpha, beta) : " << alpha << " " << beta << "\n";
sum = alpha+beta;
- mtau[i_] = 1/q10*sum;
+ mtau[i_] = 1/(q10*sum);
minf[i_] = alpha/sum;
//////////////////////////////////////////////////////////
@@ -165,7 +165,7 @@ public:
//std::cout << "h (alpha, beta) : " << alpha << " " << beta << "\n";
sum = alpha+beta;
- htau[i_] = 1/q10*sum;
+ htau[i_] = 1/(q10*sum);
hinf[i_] = alpha/sum;
//////////////////////////////////////////////////////////
@@ -178,8 +178,11 @@ public:
//std::cout << "n (alpha, beta) : " << alpha << " " << beta << "\n";
+ //printf("m_inf, h_inf, n_inf : %18.16f %18.16f %18.16f\n",
+ //minf[i_], hinf[i_], ninf[i_]);
+
sum = alpha+beta;
- ntau[i_] = 1/q10*sum;
+ ntau[i_] = 1/(q10*sum); // TODO : make modparser insert parenthesis
ninf[i_] = alpha/sum;
}
diff --git a/mechanisms/mod/hh.mod b/mechanisms/mod/hh.mod
index c96280fd490ce70fbfbeb2a66f64fc42761187a5..8b23f44a4dc3e27cdbb2d78a161b9504c8651073 100644
--- a/mechanisms/mod/hh.mod
+++ b/mechanisms/mod/hh.mod
@@ -29,7 +29,6 @@ NEURON {
}
PARAMETER {
- : neuron uses S/cm2, i have to scale by a factor of 10
gnabar = .12 (S/cm2) : <0,1e9>
gkbar = .036 (S/cm2) : <0,1e9>
gl = .0003 (S/cm2) : <0,1e9>
diff --git a/src/fvm.hpp b/src/fvm.hpp
index d680edc8d290ffce11a2b000d3edcf287713d8ed..7abf6ed4b2c465122f8b30e0b5c97c1d9a135df0 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -400,9 +400,9 @@ void fvm_cell<T, I>::setup_matrix(T dt)
// . . .
// l[i] . . d[i]
//
- d(all) = cv_areas_;
+ d(all) = 1.0;
for(auto i=1u; i<d.size(); ++i) {
- auto a = dt * face_alpha_[i];
+ auto a = dt * face_alpha_[i]; // TODO get this right
d[i] += a;
l[i] = -a;
@@ -413,9 +413,10 @@ void fvm_cell<T, I>::setup_matrix(T dt)
}
// the RHS of the linear system is
- // sigma_i * (V[i] - dt/cm*(im - ie))
+ // V[i] - dt/cm*(im - ie)
+ auto factor = 10.*dt;
for(auto i=0u; i<d.size(); ++i) {
- rhs[i] = cv_areas_[i] * (voltage_[i] - dt/cv_capacitance_[i]*current_[i]);
+ rhs[i] = voltage_[i] - factor/cv_capacitance_[i]*current_[i];
}
}
@@ -441,21 +442,24 @@ void fvm_cell<T, I>::advance(T dt)
m->nrn_current();
}
- //if(t_>=10.) {
- current_[0] -= 0.1;
- //}
+ if(t_>=5. && t_<8.) {
+ current_[0] -= 0.01;
+ }
+
//std::cout << "t " << t_ << " current " << current_;
// set matrix diagonals and rhs
setup_matrix(dt);
- //std::cout << " rhs " << matrix_.rhs() << " d " << matrix_.d();
+ //printf("rhs %18.14f d %18.14f\n", matrix_.rhs()[0], matrix_.d()[0]);
// solve the linear system
matrix_.solve();
voltage_(all) = matrix_.rhs();
+ //printf("v solve %18.14f\n", voltage_[0]);
+
//std::cout << " v " << voltage_ << "\n";
// update states
@@ -464,6 +468,7 @@ void fvm_cell<T, I>::advance(T dt)
}
t_ += dt;
+ //std::cout << "******************\n";
}
} // namespace fvm
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 9ed84900c97aaca267a02c6aed398fe13d6cae79..e4002736c32b1075c3a2e6e01bedd27a4923e1d8 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -3,7 +3,7 @@
#include "../src/cell.hpp"
#include "../src/fvm.hpp"
-// based on hh/Neuron/steps_B.py
+// based on hh/Neuron/steps_A.py
TEST(run, single_compartment)
{
using namespace nest::mc;
@@ -15,6 +15,7 @@ TEST(run, single_compartment)
// Soma with diameter 18.8um -> 1.88e-3 cm and HH channel
auto soma = cell.add_soma(1.88e-3/2.0);
+ soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
soma->add_mechanism(hh_parameters());
std::cout << soma->mechanism("membrane");
@@ -22,23 +23,21 @@ TEST(run, single_compartment)
// make the lowered finite volume cell
fvm::fvm_cell<double, int> model(cell);
- std::cout << "CV areas " << model.cv_areas() << "\n";
+ //std::cout << "CV areas " << model.cv_areas() << "\n";
- std::cout << "-----------------------------\n";
+ //std::cout << "-----------------------------\n";
// set initial conditions
using memory::all;
model.voltage()(all) = -65.;
model.initialize(); // have to do this _after_ initial conditions are set
- std::cout << "-----------------------------\n";
+ //std::cout << "-----------------------------\n";
// run the simulation
- //auto dt = 0.02 / 1000.; // convert ms to s
- //auto tfinal = 100./1000.;
auto dt = 0.02; // ms
- auto tfinal = 5.0; // ms
//auto tfinal = 0.10; // ms
+ auto tfinal = 10.; // ms
int nt = tfinal/dt;
std::vector<double> result;
result.push_back(model.voltage()[0]);