From 256c5421b863ba3123e7bb4e9bdb22de7a17b6ab Mon Sep 17 00:00:00 2001
From: Ben Cumming <louncharf@gmail.com>
Date: Tue, 14 Nov 2017 15:02:17 +0100
Subject: [PATCH] Use Nernst equation to calculate reversal potentials (#387)
Replace constant values for ion species reversal potentials taken from HH model with values that are updated on each time step via the Nernst equation.
* Implement Nernst equation in multicore and gpu back-ends.
* Extend interface of `ion` type to proved a method `update_reversal_potential` which calls the back-end nernst routine.
* Add `valency` and default concentration values to the ion species.
* Add interface for resetting ion species state (for restarting simulations).
Fixes #376
---
src/CMakeLists.txt | 1 +
src/backends/gpu/fvm.hpp | 12 +++++
src/backends/gpu/kernels/nernst.cu | 39 ++++++++++++++++
src/backends/gpu/nernst.hpp | 19 ++++++++
src/backends/multicore/fvm.hpp | 17 +++++++
src/constants.hpp | 32 +++++++++++++
src/fvm_multicell.hpp | 62 +++++++++++++++++--------
src/ion.hpp | 41 +++++++++++++---
validation/ref/neuron/nrn_validation.py | 20 ++++++++
validation/ref/numeric/HHChannels.jl | 20 ++++++--
10 files changed, 234 insertions(+), 29 deletions(-)
create mode 100644 src/backends/gpu/kernels/nernst.cu
create mode 100644 src/backends/gpu/nernst.hpp
create mode 100644 src/constants.hpp
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 6712d047..9fc7f84a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -30,6 +30,7 @@ set(CUDA_SOURCES
backends/gpu/multi_event_stream.cu
backends/gpu/kernels/assemble_matrix.cu
backends/gpu/kernels/interleave.cu
+ backends/gpu/kernels/nernst.cu
backends/gpu/kernels/solve_matrix.cu
backends/gpu/kernels/stim_current.cu
backends/gpu/kernels/take_samples.cu
diff --git a/src/backends/gpu/fvm.hpp b/src/backends/gpu/fvm.hpp
index a0aec621..4dd50b4e 100644
--- a/src/backends/gpu/fvm.hpp
+++ b/src/backends/gpu/fvm.hpp
@@ -13,6 +13,7 @@
#include "kernels/take_samples.hpp"
#include "matrix_state_interleaved.hpp"
#include "multi_event_stream.hpp"
+#include "nernst.hpp"
#include "stimulus.hpp"
#include "threshold_watcher.hpp"
#include "time_ops.hpp"
@@ -135,6 +136,17 @@ struct backend {
arb::gpu::take_samples(s, time.data(), sample_time.data(), sample_value.data());
}
+ // Calculate the reversal potential eX (mV) using Nernst equation
+ // eX = RT/zF * ln(Xo/Xi)
+ // R: universal gas constant 8.3144598 J.K-1.mol-1
+ // T: temperature in Kelvin
+ // z: valency of species (K, Na: +1) (Ca: +2)
+ // F: Faraday's constant 96485.33289 C.mol-1
+ // Xo/Xi: ratio of out/in concentrations
+ static void nernst(int valency, value_type temperature, const_view Xo, const_view Xi, view eX) {
+ arb::gpu::nernst(eX.size(), valency, temperature, Xo.data(), Xi.data(), eX.data());
+ }
+
private:
using maker_type = mechanism_ptr (*)(size_type, const_iview, const_view, const_view, const_view, view, view, array&&, iarray&&);
static std::map<std::string, maker_type> mech_map_;
diff --git a/src/backends/gpu/kernels/nernst.cu b/src/backends/gpu/kernels/nernst.cu
new file mode 100644
index 00000000..5e84d63b
--- /dev/null
+++ b/src/backends/gpu/kernels/nernst.cu
@@ -0,0 +1,39 @@
+#include <cstdint>
+
+#include <constants.hpp>
+
+#include "../nernst.hpp"
+#include "detail.hpp"
+
+namespace arb {
+namespace gpu {
+
+namespace kernels {
+ template <typename T>
+ __global__ void nernst(std::size_t n, int valency, T temperature, const T* Xo, const T* Xi, T* eX) {
+ auto i = threadIdx.x+blockIdx.x*blockDim.x;
+
+ // factor 1e3 to scale from V -> mV
+ constexpr T RF = 1e3*constant::gas_constant/constant::faraday;
+ T factor = RF*temperature/valency;
+ if (i<n) {
+ eX[i] = factor*std::log(Xo[i]/Xi[i]);
+ }
+ }
+} // namespace kernels
+
+void nernst(std::size_t n,
+ int valency,
+ fvm_value_type temperature,
+ const fvm_value_type* Xo,
+ const fvm_value_type* Xi,
+ fvm_value_type* eX)
+{
+ constexpr int block_dim = 128;
+ const int grid_dim = impl::block_count(n, block_dim);
+ kernels::nernst<<<grid_dim, block_dim>>>
+ (n, valency, temperature, Xo, Xi, eX);
+}
+
+} // namespace gpu
+} // namespace arb
diff --git a/src/backends/gpu/nernst.hpp b/src/backends/gpu/nernst.hpp
new file mode 100644
index 00000000..2597dada
--- /dev/null
+++ b/src/backends/gpu/nernst.hpp
@@ -0,0 +1,19 @@
+#pragma once
+
+#include <cstdint>
+
+#include <backends/fvm_types.hpp>
+
+namespace arb {
+namespace gpu {
+
+// prototype for nernst equation cacluation
+void nernst(std::size_t n, int valency,
+ fvm_value_type temperature,
+ const fvm_value_type* Xo,
+ const fvm_value_type* Xi,
+ fvm_value_type* eX);
+
+} // namespace gpu
+} // namespace arb
+
diff --git a/src/backends/multicore/fvm.hpp b/src/backends/multicore/fvm.hpp
index af90096b..a6f758a1 100644
--- a/src/backends/multicore/fvm.hpp
+++ b/src/backends/multicore/fvm.hpp
@@ -6,6 +6,7 @@
#include <backends/event.hpp>
#include <backends/fvm_types.hpp>
#include <common_types.hpp>
+#include <constants.hpp>
#include <event_queue.hpp>
#include <mechanism.hpp>
#include <memory/memory.hpp>
@@ -148,6 +149,22 @@ struct backend {
}
}
+ // Calculate the reversal potential eX (mV) using Nernst equation
+ // eX = RT/zF * ln(Xo/Xi)
+ // R: universal gas constant 8.3144598 J.K-1.mol-1
+ // T: temperature in Kelvin
+ // z: valency of species (K, Na: +1) (Ca: +2)
+ // F: Faraday's constant 96485.33289 C.mol-1
+ // Xo/Xi: ratio of out/in concentrations
+ static void nernst(int valency, value_type temperature, const_view Xo, const_view Xi, view eX) {
+ // factor 1e3 to scale from V -> mV
+ constexpr value_type RF = 1e3*constant::gas_constant/constant::faraday;
+ value_type factor = RF*temperature/valency;
+ for (std::size_t i=0; i<Xi.size(); ++i) {
+ eX[i] = factor*std::log(Xo[i]/Xi[i]);
+ }
+ }
+
private:
using maker_type = mechanism_ptr (*)(value_type, const_iview, const_view, const_view, const_view, view, view, array&&, iarray&&);
static std::map<std::string, maker_type> mech_map_;
diff --git a/src/constants.hpp b/src/constants.hpp
new file mode 100644
index 00000000..02da7014
--- /dev/null
+++ b/src/constants.hpp
@@ -0,0 +1,32 @@
+#pragma once
+
+namespace arb {
+namespace constant {
+
+// TODO: handle change of constants over time. Take, for example, the values
+// for the universal gas constant (R) and Faraday's constant as given by NIST:
+//
+// R F
+// 1973 8.31441 96486.95
+// 1986 8.314510 96485.309
+// 1998 8.314472 96485.3415
+// 2002 8.314472 96485.3383
+// 2010 8.3144621 96485.3365
+// 2014 8.3144598 96485.33289
+
+// TODO: use value templates from C++14
+
+// Universal gas constant (R)
+// https://physics.nist.gov/cgi-bin/cuu/Value?r
+constexpr double gas_constant = 8.3144598; // J.°K^-1.mol^-1
+
+// Faraday's constant (F)
+// https://physics.nist.gov/cgi-bin/cuu/Value?f
+constexpr double faraday = 96485.33289; // C.mol^-1
+
+// Temperature used in original Hodgkin-Huxley paper
+// doi:10.1113/jphysiol.1952.sp004764
+constexpr double hh_squid_temp = 6.3+273.15; // °K
+
+} // namespace arb
+} // namespace arb
diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp
index a58961c2..a36ee18b 100644
--- a/src/fvm_multicell.hpp
+++ b/src/fvm_multicell.hpp
@@ -12,6 +12,7 @@
#include <backends/fvm_types.hpp>
#include <cell.hpp>
#include <compartment.hpp>
+#include <constants.hpp>
#include <event_queue.hpp>
#include <ion.hpp>
#include <math.hpp>
@@ -1053,24 +1054,28 @@ void fvm_multicell<Backend>::initialize(
}
}
- // FIXME: Hard code parameters for now.
- // Take defaults for reversal potential of sodium and potassium from
- // the default values in Neuron.
- // Neuron's defaults are defined in the file
- // nrn/src/nrnoc/membdef.h
- constexpr value_type DEF_vrest = -65.0; // same name as #define in Neuron
+ // Note: NEURON defined default values for reversal potential as follows,
+ // with units mV:
+ //
+ // const auto DEF_vrest = -65.0
+ // ena = 115.0 + DEF_vrest
+ // ek = -12.0 + DEF_vrest
+ // eca = 12.5*std::log(2.0/5e-5)
+ //
+ // Whereas we use the Nernst equation to calculate reversal potentials at
+ // the start of each time step.
- memory::fill(ion_na().reversal_potential(), 115+DEF_vrest); // mV
- memory::fill(ion_na().internal_concentration(), 10.0); // mM
- memory::fill(ion_na().external_concentration(), 140.0); // mM
+ ion_na().default_internal_concentration = 10;
+ ion_na().default_external_concentration =140;
+ ion_na().valency = 1;
- memory::fill(ion_k().reversal_potential(), -12.0+DEF_vrest);// mV
- memory::fill(ion_k().internal_concentration(), 54.4); // mM
- memory::fill(ion_k().external_concentration(), 2.5); // mM
+ ion_k().default_internal_concentration =54.4;
+ ion_k().default_external_concentration = 2.5;
+ ion_k().valency = 1;
- memory::fill(ion_ca().reversal_potential(), 12.5*std::log(2.0/5e-5));// mV
- memory::fill(ion_ca().internal_concentration(), 5e-5); // mM
- memory::fill(ion_ca().external_concentration(), 2.0); // mM
+ ion_ca().default_internal_concentration =5e-5;
+ ion_ca().default_external_concentration = 2.0;
+ ion_ca().valency = 2;
// initialize mechanism and voltage state
reset();
@@ -1083,13 +1088,25 @@ void fvm_multicell<Backend>::reset() {
set_time_global(0);
set_time_to_global(0);
+ // Update ion species:
+ // - clear currents
+ // - reset concentrations to defaults
+ // - recalculate reversal potentials
+ for (auto& i: ions_) {
+ i.second.reset();
+ }
+
for (auto& m : mechanisms_) {
- // TODO : the parameters have to be set before the nrn_init
- // for now use a dummy value of dt.
m->set_params();
m->nrn_init();
}
+ // Update reversal potential to account for changes to concentrations made
+ // by calls to nrn_init() in mechansisms.
+ for (auto& i: ions_) {
+ i.second.nernst_reversal_potential(constant::hh_squid_temp); // TODO: use temperature specfied in model
+ }
+
// Reset state of the threshold watcher.
// NOTE: this has to come after the voltage_ values have been reinitialized,
// because these values are used by the watchers to set their initial state.
@@ -1110,11 +1127,18 @@ template <typename Backend>
void fvm_multicell<Backend>::step_integration() {
EXPECTS(!integration_complete());
+ // mark pending events for delivery
+ events_.mark_until_after(time_);
+
PE("current");
memory::fill(current_, 0.);
- // mark pending events for delivery
- events_.mark_until_after(time_);
+ // clear currents and recalculate reversal potentials for all ion channels
+ for (auto& i: ions_) {
+ auto& ion = i.second;
+ memory::fill(ion.current(), 0.);
+ ion.nernst_reversal_potential(constant::hh_squid_temp); // TODO: use temperature specfied in model
+ }
// deliver pending events and update current contributions from mechanisms
for (auto& m: mechanisms_) {
diff --git a/src/ion.hpp b/src/ion.hpp
index d719daa8..99530bdf 100644
--- a/src/ion.hpp
+++ b/src/ion.hpp
@@ -1,6 +1,7 @@
#pragma once
#include <array>
+#include <constants.hpp>
#include <memory/memory.hpp>
#include <util/indirect.hpp>
@@ -62,7 +63,10 @@ public :
iX_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()},
eX_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()},
Xi_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()},
- Xo_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()}
+ Xo_{idx.size(), std::numeric_limits<value_type>::quiet_NaN()},
+ valency(0),
+ default_internal_concentration(0),
+ default_external_concentration(0)
{}
std::size_t memory() const {
@@ -74,16 +78,35 @@ public :
view current() {
return iX_;
}
+
view reversal_potential() {
return eX_;
}
+
view internal_concentration() {
return Xi_;
}
+
view external_concentration() {
return Xo_;
}
+ void reset() {
+ // The Nernst equation uses the assumption of nonzero concentrations:
+ EXPECTS(default_internal_concentration> value_type(0));
+ EXPECTS(default_external_concentration> value_type(0));
+ memory::fill(iX_, 0);
+ memory::fill(Xi_, default_internal_concentration);
+ memory::fill(Xo_, default_external_concentration);
+ nernst_reversal_potential(constant::hh_squid_temp); // TODO: use temperature specfied in model
+ }
+
+ /// Calculate the reversal potential for all compartments using Nernst equation
+ /// temperature is in degrees Kelvin
+ void nernst_reversal_potential(value_type temperature) {
+ backend::nernst(valency, temperature, Xo_, Xi_, eX_);
+ }
+
const_iview node_index() const {
return node_index_;
}
@@ -92,13 +115,17 @@ public :
return node_index_.size();
}
-private :
-
+private:
iarray node_index_;
- array iX_;
- array eX_;
- array Xi_;
- array Xo_;
+ array iX_; // (nA) current
+ array eX_; // (mV) reversal potential
+ array Xi_; // (mM) internal concentration
+ array Xo_; // (mM) external concentration
+
+public:
+ int valency; // valency of ionic species
+ value_type default_internal_concentration; // (mM) default internal concentration
+ value_type default_external_concentration; // (mM) default external concentration
};
} // namespace arb
diff --git a/validation/ref/neuron/nrn_validation.py b/validation/ref/neuron/nrn_validation.py
index 961b2771..2e1cf7f3 100644
--- a/validation/ref/neuron/nrn_validation.py
+++ b/validation/ref/neuron/nrn_validation.py
@@ -129,6 +129,26 @@ class VModel:
soma.gl_hh = p['gl_hh']
soma.el_hh = p['el_hh']
+ # For reversal potentials we use those computed using
+ # the Nernst equation with the following values:
+ # R 8.3144598
+ # F 96485.33289
+ # nao 140 mM
+ # nai 10 mM
+ # ko 2.5 mM
+ # ki 64.4 nM
+ # We don't use the default values for ena and ek taken
+ # from the HH paper:
+ # ena = 115.0mV + -65.0mV,
+ # ek = -12.0mV + -65.0mV,
+ soma.ena = 63.55148117386
+ soma.ek = -74.17164678272
+
+ # This is how we would get NEURON to use Nernst equation, when they
+ # correct the Nernst equation implementation.
+ #h.ion_style('k_ion', 3, 2, 1, 1, 1)
+ #h.ion_style('na_ion', 3, 2, 1, 1, 1)
+
self.soma = soma
def add_dendrite(self, name, geom, to=None, **kw):
diff --git a/validation/ref/numeric/HHChannels.jl b/validation/ref/numeric/HHChannels.jl
index bf32c8e8..910a402c 100644
--- a/validation/ref/numeric/HHChannels.jl
+++ b/validation/ref/numeric/HHChannels.jl
@@ -19,12 +19,26 @@ immutable HHParam
# constructor with default values, corresponding
# to a resting potential of -65 mV and temperature 6.3 °C
HHParam(;
- #c_m = 0.01nF*m^-2,
+ # default values from HH paper
+
+ # For reversal potentials we use those computed using
+ # the Nernst equation with the following values:
+ # R 8.3144598
+ # F 96485.33289
+ # nao 140 mM
+ # nai 10 mM
+ # ko 2.5 mM
+ # ki 64.4 nM
+ # We don't use the default values for ena and ek taken
+ # from the HH paper:
+ # ena = 115.0mV + -65.0mV,
+ # ek = -12.0mV + -65.0mV,
+ ena = 63.55148117386mV,
+ ek = -74.17164678272mV,
+
c_m = 0.01F*m^-2,
gnabar = .12S*cm^-2,
- ena = 115.0mV + -65.0mV,
gkbar = .036S*cm^-2,
- ek = -12.0mV + -65.0mV,
gl = .0003S*cm^-2,
el = -54.3mV,
q10 = 1
--
GitLab