Merge branch 'gap-junctions' of https://github.com/eth-cscs/arbor into gap-junctions

example/gap_junctions/gap_junctions.cpp

@@ -45,9 +45,9 @@ void write_trace_json(const std::vector<arb::trace_data<double>>& trace);
 // Generate a cell.
 arb::mc_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params, bool stim);
 
-class ring_recipe: public arb::recipe {
+class gj_recipe: public arb::recipe {
 public:
-    ring_recipe(unsigned num_cells, cell_parameters params, unsigned min_delay):
+    gj_recipe(unsigned num_cells, cell_parameters params, unsigned min_delay):
         num_cells_(num_cells),
         cell_params_(params),
         min_delay_(min_delay)
@@ -61,8 +61,8 @@ public:
             }
         }
         for (unsigned i = 0; i < num_cells - 1; i++) {
-            cells[i].add_gap_junction(i, {0, 1}, i+1, {0,1}, 0.0003430471145874217);
-            cells[i+1].add_gap_junction(i+1, {0, 1}, i, {0,1}, 0.0003430471145874217);
+            cells[i].add_gap_junction(i, {0, 1}, i+1, {0,1}, 0.05551822537423883);
+            cells[i+1].add_gap_junction(i+1, {0, 1}, i, {0,1}, 0.05551822537423883);
         }
     }
 
@@ -207,7 +207,7 @@ int main(int argc, char** argv) {
         };
 
         // Create an instance of our recipe.
-        ring_recipe recipe(params.num_cells, params.cell, params.min_delay);
+        gj_recipe recipe(params.num_cells, params.cell, params.min_delay);
 
         for(unsigned i = 0; i < recipe.num_cells() - 1; i++){
             std::cout << arb::util::any_cast<arb::mc_cell>(recipe.get_cell_description(i)).gap_junctions().size() << std::endl;
@@ -280,7 +280,7 @@ int main(int argc, char** argv) {
         std::cout << report;
     }
     catch (std::exception& e) {
-        std::cerr << "exception caught in ring miniapp:\n" << e.what() << "\n";
+        std::cerr << "exception caught in gap junction miniapp:\n" << e.what() << "\n";
         return 1;
     }
 
@@ -292,7 +292,7 @@ void write_trace_json(const std::vector<arb::trace_data<double>>& trace) {
         std::string path = "./voltages" + std::to_string(i) + ".json";
 
         nlohmann::json json;
-        json["name"] = "ring demo";
+        json["name"] = "gj demo";
         json["units"] = "mV";
         json["cell"] = "0.0";
         json["probe"] = "0";

example/gap_junctions/parameters.hpp

@@ -24,8 +24,8 @@ struct cell_parameters {
     std::array<double,2> lengths = {200, 20};       //  Length of branch in μm.
 };
 
-struct ring_params {
-    ring_params() = default;
+struct gj_params {
+    gj_params() = default;
 
     std::string name = "default";
     unsigned num_cells = 2;
@@ -34,10 +34,10 @@ struct ring_params {
     cell_parameters cell;
 };
 
-ring_params read_options(int argc, char** argv) {
+gj_params read_options(int argc, char** argv) {
     using sup::param_from_json;
 
-    ring_params params;
+    gj_params params;
     if (argc<2) {
         std::cout << "Using default parameters.\n";
         return params;

validation/ref/numeric/HHChannels_multipleSomas.jl

+module HHChannels
+
+export Stim, run_hh
+
+using Sundials
+using Unitful
+using Unitful.DefaultSymbols
+
+struct HHParam
+    c_m       # membrane spacific capacitance
+    gnabar    # Na channel cross-membrane conductivity
+    gkbar     # K channel cross-membrane conductivity
+    gl        # Leak conductivity
+    ena       # Na channel reversal potential
+    ek        # K channel reversal potential
+    el        # Leak reversal potential
+    q10       # temperature dependent rate coefficient
+              # (= 3^((T-T₀)/10K) with T₀ = 6.3 °C)
+
+    # constructor with default values, corresponding
+    # to a resting potential of -65 mV and temperature 6.3 °C
+    HHParam(;
+        # 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,
+        gkbar  = .036S*cm^-2,
+        gl     = .0003S*cm^-2,
+        el     = -54.3mV,
+        q10    = 1
+    ) = new(c_m, gnabar, gkbar, gl, ena, ek, el, q10)
+end
+
+struct Stim
+    t0        # start time of stimulus
+    t1        # stop time of stimulus
+    i_e       # stimulus current density
+
+    Stim() = new(0s, 0s, 0A/m^2)
+    Stim(t0, t1, i_e) = new(t0, t1, i_e)
+end
+
+scale(quantity, unit) = uconvert(NoUnits, quantity/unit)
+
+vtrap(x,y) = x/(exp(x/y) - 1.0)
+
+# "m" sodium activation system
+function m_lims(v, q10)
+    alpha = .1mV^-1 * vtrap(-(v+40mV),10mV)
+    beta =  4 * exp(-(v+65mV)/18mV)
+    sum = alpha + beta
+    mtau = 1ms / (q10*sum)
+    minf = alpha/sum
+    return mtau, minf
+end
+
+# "h" sodium inactivation system
+function h_lims(v, q10)
+    alpha = 0.07*exp(-(v+65mV)/20mV)
+    beta = 1 / (exp(-(v+35mV)/10mV) + 1)
+    sum = alpha + beta
+    htau = 1ms / (q10*sum)
+    hinf = alpha/sum
+    return htau, hinf
+end
+
+# "n" potassium activation system
+function n_lims(v, q10)
+    alpha = .01mV^-1 * vtrap(-(v+55mV),10mV)
+    beta = .125*exp(-(v+65mV)/80mV)
+    sum = alpha + beta
+    ntau = 1ms / (q10*sum)
+    ninf = alpha/sum
+    return ntau, ninf
+end
+
+# Choose initial conditions for the system such that the gating variables
+# are at steady state for the user-specified voltage v
+function initial_conditions(v, q10)
+    mtau, minf = m_lims(v, q10)
+    htau, hinf = h_lims(v, q10)
+    ntau, ninf = n_lims(v, q10)
+
+    return (v, minf, hinf, ninf)
+end
+
+# Get two neighboring somas in a ring of N somas
+function neighbor(m, N)
+    if m == 1
+        l = m + (N-1)
+        r = m + 1
+    elseif m == N
+        l = m-1
+        r = m-(N-1)
+    else
+        l = m-1
+        r = m+1
+    end
+    return l,r
+end
+
+# Given time t and state (v, m, h, n),
+# return (vdot, mdot, hdot, ndot)
+function f(t, state, ggap; p=HHParam(), stim=Stim())
+
+    # state of somas
+    len     = length(state)
+    n_cells = Int64(len/4)
+
+    result = []
+
+    for i=1:n_cells
+        v = state[(i-1)*4 + 1]
+        m = state[(i-1)*4 + 2]
+        h = state[(i-1)*4 + 3]
+        n = state[(i-1)*4 + 4]
+
+        # calculate current density due to ion channels
+        gna = p.gnabar*m*m*m*h
+        gk  = p.gkbar*n*n*n*n
+        ina = gna*(v - p.ena)
+        ik  = gk*(v - p.ek)
+        il  = p.gl*(v - p.el)
+        
+        # calculate gap junction currents of ring
+        l,r = neighbor(i, n_cells)
+        vl  = state[(l-1)*4 + 1]
+        vr  = state[(r-1)*4 + 1]
+        igapl = ggap*(v - vl)
+        igapr = ggap*(v - vr)
+
+        if (n_cells == 2 && i == 1)
+            igapt = igapr
+        elseif (n_cells == 2 && i == 2)
+            igapt = igapl
+        else
+            igapt = igapl + igapr
+        end
+
+        itot = ik + ina + il + igapt
+
+        # calculate current density due to stimulus
+        if t>=stim.t0 && t<stim.t1 && i==1
+            itot -= stim.i_e
+        end
+        
+        # calculate the voltage dependent rates for the gating variables
+        mtau, minf = m_lims(v, p.q10)
+        htau, hinf = h_lims(v, p.q10)
+        ntau, ninf = n_lims(v, p.q10)
+
+        push!(result, -itot/p.c_m, (minf-m)/mtau, (hinf-h)/htau, (ninf-n)/ntau)
+
+    end
+    return result
+end
+
+function run_hh(t_end, ggap, n_cells; v0=-65mV, stim=Stim(), param=HHParam(), sample_dt=0.01ms)
+    v_scale = 1V
+    t_scale = 1s
+
+    # initialize
+    y0 = Float64[]
+    v, m, h, n = initial_conditions(v0, param.q10)
+    for i=1:n_cells
+       push!(y0, v/v_scale, m, h, n)
+    end
+
+    fbis(t, y, ydot) = begin
+        state = []
+        
+        for i=1:n_cells
+           push!(state, y[(i-1)*4 + 1]*v_scale, y[(i-1)*4 + 2], y[(i-1)*4 + 3], y[(i-1)*4 + 4])
+        end
+
+        fdot = f(t*t_scale, state, ggap, stim=stim, p=param)
+        
+        for i=1:n_cells
+           ydot[(i-1)*4 + 1] = fdot[(i-1)*4 + 1]*t_scale/v_scale
+           ydot[(i-1)*4 + 2] = fdot[(i-1)*4 + 2]*t_scale
+           ydot[(i-1)*4 + 3] = fdot[(i-1)*4 + 3]*t_scale
+           ydot[(i-1)*4 + 4] = fdot[(i-1)*4 + 4]*t_scale
+        end
+
+        return Sundials.CV_SUCCESS
+    end
+
+    # Ideally would run with vector absolute tolerance to account for v_scale,
+    # but this would prevent us using the nice cvode wrapper.
+
+    samples = collect(0s: sample_dt: t_end)
+    res = Sundials.cvode(fbis, y0, scale.(samples, t_scale), abstol=1e-6, reltol=5e-10)
+
+    result = []
+    for i=1:n_cells
+      push!(result, res[:,(i-1)*4+1]*v_scale)
+    end
+
+    return samples, result
+
+end
+
+end # module HHChannels

validation/ref/numeric/numeric_multipleSomas.jl

+#!/usr/bin/env julia
+
+include("HHChannels_multipleSomas.jl")
+
+using JSON
+using Unitful
+using Unitful.DefaultSymbols
+using ArgParse
+using DataStructures
+using Main.HHChannels
+
+scale(quantity, unit) = uconvert(NoUnits, quantity/unit)
+
+function parse_commandline()
+    s = ArgParseSettings()
+
+    @add_arg_table s begin
+        "--n_cells", "-n"
+            help = "number of cells"
+            arg_type = Int
+            default = 3
+        "--diameter", "-d"
+            help = "diameter of cells in µm"
+            arg_type = Float64
+            default = 18.8
+        "--ggap", "-g"
+            help = "conductivity of gap junctions in S cm^-2"
+            arg_type = Float64
+            default = 5e-5
+        "--stim_t0"
+            help = "start time of stimulus in ms"
+            arg_type = Float64
+            default = 20.0
+        "--stim_t1"
+            help = "end time of stimulus in ms"
+            arg_type = Float64
+            default = 22.0
+        "--stim", "-s"
+            help = "strength of stimulus in nA"
+            arg_type = Float64
+            default = 0.1
+        "--t_end", "-t"
+            help = "simulation time in ms"
+            arg_type = Float64
+            default = 100.0
+        "--dt"
+            help = "time step in ms"
+            arg_type = Float64
+            default = 0.025
+        "--output", "-o"
+            help = "JSON filename"
+            default = "output.txt"
+    end
+    return parse_args(s)
+end
+
+parsed_args = parse_commandline()
+
+n_cells = parsed_args["n_cells"]
+radius = parsed_args["diameter"]/2*µm
+area = 4*pi*radius^2
+
+ggap = parsed_args["ggap"]*S*cm^-2
+
+ts0 = parsed_args["stim_t0"]*ms
+ts1 = parsed_args["stim_t1"]*ms
+St  = parsed_args["stim"]/area*nA
+
+stim = Stim(ts0, ts1, St)
+
+t_end = parsed_args["t_end"]*ms
+dt = parsed_args["dt"]*ms
+
+ts, vs = run_hh(t_end, ggap, n_cells, stim=stim, sample_dt=dt)
+
+trace = Dict(
+    :name => "membrane voltage",
+    :sim => "numeric",
+    :model => "soma",
+    :units => "mV",
+    :data => OrderedDict{}(
+        :time => scale.(ts, 1ms),
+        Symbol("soma1.mid") => scale.(vs[1], 1mV)
+    )
+)
+
+for i=2:n_cells
+    trace[:data][Symbol("soma.mid"*string(i))] = scale.(vs[i], 1mV)
+end
+
+fname = parsed_args["output"]
+io = open(fname, "w");
+println(io, JSON.json([trace]))