From 8af9d74b8f95303c9468bf65d45e1c34ff267988 Mon Sep 17 00:00:00 2001 From: adajel Date: Fri, 23 Aug 2024 11:15:55 +0200 Subject: [PATCH] add code for ODE pre-calibration and full system (ODE/PDE) calibration --- emix-simulations/make_mesh_3D_two_tags.py | 122 +++++++++ .../mm_hh_glial_calibration_ODE.py | 257 ++++++++++++++++++ emix-simulations/run_calibration_two_tags.py | 158 +++++++++++ emix-simulations/run_two_tags.py | 162 +++++++++++ 4 files changed, 699 insertions(+) create mode 100755 emix-simulations/make_mesh_3D_two_tags.py create mode 100644 emix-simulations/mm_hh_glial_calibration_ODE.py create mode 100644 emix-simulations/run_calibration_two_tags.py create mode 100644 emix-simulations/run_two_tags.py diff --git a/emix-simulations/make_mesh_3D_two_tags.py b/emix-simulations/make_mesh_3D_two_tags.py new file mode 100755 index 0000000..fe73831 --- /dev/null +++ b/emix-simulations/make_mesh_3D_two_tags.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 + +""" +This script generates a 3D mesh representing 4 axons. +""" + +from dolfin import * +import sys + +class Boundary(SubDomain): + # define exterior boundary + def inside(self, x, on_boundary): + return on_boundary + +def add_axon(mesh, subdomains, surfaces, a, b, tag): + # define interior domain + in_interior = """ (x[0] >= %g && x[0] <= %g && + x[1] >= %g && x[1] <= %g && + x[2] >= %g && x[2] <= %g) """ % ( + a[0], + b[0], + a[1], + b[1], + a[2], + b[2], + ) + + interior = CompiledSubDomain(in_interior) + + # mark interior and exterior domain + for cell in cells(mesh): + x = cell.midpoint().array() + if (int(interior.inside(x, False))): + subdomains[cell] = tag + assert sum(1 for _ in SubsetIterator(subdomains, 1)) > 0 + + for facet in facets(mesh): + x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()] + side_1 = near(x[0], a[0]) and a[1] <= x[1] <= b[1] and a[2] <= x[2] <= b[2] + side_2 = near(x[0], b[0]) and a[1] <= x[1] <= b[1] and a[2] <= x[2] <= b[2] + side_3 = near(x[1], a[1]) and a[0] <= x[0] <= b[0] and a[2] <= x[2] <= b[2] + side_4 = near(x[1], b[1]) and a[0] <= x[0] <= b[0] and a[2] <= x[2] <= b[2] + side_5 = near(x[2], a[2]) and a[0] <= x[0] <= b[0] and a[1] <= x[1] <= b[1] + side_6 = near(x[2], b[2]) and a[0] <= x[0] <= b[0] and a[1] <= x[1] <= b[1] + if (side_1 or side_2 or side_3 or side_4 or side_5 or side_6): + surfaces[facet] = tag + + return + +import argparse +from pathlib import Path + + +class CustomParser( + argparse.ArgumentDefaultsHelpFormatter, argparse.RawDescriptionHelpFormatter +): ... + +def main(argv=None): + parser = argparse.ArgumentParser(formatter_class=CustomParser) + parser.add_argument( + "-r", + "--resolution", + dest="resolution_factor", + default=0, + type=int, + help="Mesh resolution factor", + ) + parser.add_argument( + "-d", + "--directory", + dest="mesh_dir", + type=Path, + default=Path("meshes/3D_two_tags"), + help="Directory to save the mesh", + ) + args = parser.parse_args(argv) + resolution_factor = args.resolution_factor + out_dir = args.mesh_dir + + l = 1 + nx = l * 20 * 2 ** resolution_factor + ny = 20 * 2 ** resolution_factor + nz = 20 * 2 ** resolution_factor + + # box mesh + mesh = BoxMesh(Point(0, 0.0, 0.0), Point(l * 20, 2.0, 2.0), nx, ny, nz) + subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0) + surfaces = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0) + + a = Point(5, 0.2, 0.2) + b = Point(l * 20 - 5, 0.6, 0.6) + add_axon(mesh, subdomains, surfaces, a, b, 1) + + a = Point(5, 0.8, 0.8) + b = Point(l * 20 - 5, 1.4, 1.4) + add_axon(mesh, subdomains, surfaces, a, b, 2) + + # mark exterior boundary + Boundary().mark(surfaces, 5) + + # convert mesh from um to cm + mesh.coordinates()[:, :] *= 1e-4 + + # save .xml files + mesh_file = File(str((out_dir / f"mesh_{resolution_factor}.xml").absolute())) + mesh_file << mesh + + subdomains_file = File(str((out_dir / f"subdomains_{resolution_factor}.xml").absolute())) + subdomains_file << subdomains + + surfaces_file = File(str((out_dir / f"surfaces_{resolution_factor}.xml").absolute())) + surfaces_file << surfaces + + # save .pvd files + mesh_file = File(str((out_dir / f"mesh_{resolution_factor}.pvd").absolute())) + mesh_file << mesh + + subdomains_file = File(str((out_dir / f"subdomains_{resolution_factor}.pvd").absolute())) + subdomains_file << subdomains + + surfaces_file = File(str((out_dir / f"surfaces_{resolution_factor}.pvd").absolute())) + surfaces_file << surfaces diff --git a/emix-simulations/mm_hh_glial_calibration_ODE.py b/emix-simulations/mm_hh_glial_calibration_ODE.py new file mode 100644 index 0000000..2f3916a --- /dev/null +++ b/emix-simulations/mm_hh_glial_calibration_ODE.py @@ -0,0 +1,257 @@ +# Gotran generated code for the "hodgkin_huxley_squid_axon_model_1952_original" model + +import numpy as np +import math + +def init_state_values(**values): + """ + Initialize state values + """ + + # Init values + phi_M_n_init = -74.38609374461858 + phi_M_g_init = -83.08511451850003 + K_e_init = 3.3236967382613933 + K_n_init = 124.15397583492471 + K_g_init = 102.75563828644862 + Na_e_init = 100.71925900028181 + Na_n_init = 12.838513108606818 + Na_g_init = 12.39731187972181 + n_init = 0.18820202480418377 + m_init = 0.016648440745826054 + h_init = 0.8542015627820454 + + init_values = np.array([m_init, h_init, n_init, phi_M_n_init, phi_M_g_init, \ + K_e_init, K_n_init, K_g_init, \ + Na_e_init, Na_n_init, Na_g_init], dtype=np.float_) + + # State indices and limit checker + state_inds = dict([("m", 0), ("h", 1), ("n", 2), ("V_n", 3), ("V_g", 4), + ("K_e", 5), ("K_n", 6), ("K_g", 7), + ("Na_e", 8), ("Na_n", 9), ("Na_g", 10)]) + + for state_name, value in values.items(): + if state_name not in state_ind: + raise ValueError("{0} is not a state.".format(state_name)) + ind = state_ind[state_name] + + # Assign value + init_values[ind] = value + + return init_values + +def init_parameter_values(**values): + """ + Initialize parameter values + """ + + # Membrane parameters + g_Na_bar = 120 # Na max conductivity (mS/cm**2) + g_K_bar = 36 # K max conductivity (mS/cm**2) + g_leak_Na_n = 0.1 # Na leak conductivity (mS/cm**2) + g_leak_K_n = 0.4 # K leak conductivity (mS/cm**2) + + g_leak_Na_g = 0.1 # Na leak conductivity (mS/cm**2) + g_leak_K_g = 1.7 # K leak conductivity (mS/cm**2) + I_max_g = 50 # max pump strength (muA/cm^2) + + m_K = 2 # threshold ECS K (mol/m^3) + m_Na = 7.7 # threshold ICS Na (mol/m^3) + I_max_n = 44.9 # max pump strength (muA/cm^2) + C_M = 2.0 # Faraday's constant (mC/ mol) + + # Set initial parameter values + init_values = np.array([g_Na_bar, g_K_bar, \ + g_leak_Na_n, g_leak_K_n, \ + g_leak_Na_g, g_leak_K_g, \ + C_M, 0, \ + m_K, m_Na, I_max_n, I_max_g], dtype=np.float_) + + # Parameter indices and limit checker + param_ind = dict([("g_Na_bar", 0), ("g_K_bar", 1), + ("g_leak_Na_n", 2), ("g_leak_K_n", 3), + ("g_leak_Na_g", 4), ("g_leak_K_g", 5), + ("Cm", 6), ("stim_amplitude", 7), + ("m_K", 8), ("m_Na", 9), ("I_max_n", 10), + ("I_max_g", 11)]) + + for param_name, value in values.items(): + if param_name not in param_ind: + raise ValueError("{0} is not a parameter.".format(param_name)) + ind = param_ind[param_name] + + # Assign value + init_values[ind] = value + + return init_values + +def state_indices(*states): + """ + State indices + """ + # State indices and limit checker + state_inds = dict([("m", 0), ("h", 1), ("n", 2), ("V_n", 3), ("V_g", 4), + ("K_e", 5), ("K_n", 6), ("K_g", 7), + ("Na_e", 8), ("Na_n", 9), ("Na_g", 10)]) + + indices = [] + for state in states: + if state not in state_inds: + raise ValueError("Unknown state: '{0}'".format(state)) + indices.append(state_inds[state]) + if len(indices)>1: + return indices + else: + return indices[0] + +def parameter_indices(*params): + """ + Parameter indices + """ + + param_inds = dict([("g_Na_bar", 0), ("g_K_bar", 1), + ("g_leak_Na_n", 2), ("g_leak_K_n", 3), + ("g_leak_Na_g", 4), ("g_leak_K_g", 5), + ("Cm", 6), ("stim_amplitude", 7), + ("m_K", 8), ("m_Na", 9), ("I_max_n", 10), + ("I_max_g", 11)]) + + indices = [] + for param in params: + if param not in param_inds: + raise ValueError("Unknown param: '{0}'".format(param)) + indices.append(param_inds[param]) + if len(indices)>1: + return indices + else: + return indices[0] + +from numbalsoda import lsoda_sig +from numba import njit, cfunc, jit +import numpy as np +import timeit +import math + +@cfunc(lsoda_sig, nopython=True) +def rhs_numba(t, states, values, parameters): + """ + Compute the right hand side of the\ + hodgkin_huxley_squid_axon_model_1952_original ODE + """ + + # Assign states + #assert(len(states)) == 4 + + # Assign parameters + #assert(len(parameters)) == 11 + + # # Init return args + # if values is None: + # values = np.zeros((4,), dtype=np.float_) + # else: + # assert isinstance(values, np.ndarray) and values.shape == (4,) + + # Physical parameters (PDEs) + temperature = 300e3 # temperature (m K) + R = 8.314e3 # Gas Constant (m J/(K mol)) + F = 96485e3 # Faraday's constant (mC/ mol) + + ICS_vol = 3.42e-11/2.0 # ICS volume (cm^3) + ECS_vol = 7.08e-11 # ECS volume (cm^3) + surface = 2.29e-6 # membrane surface (cm²) + + K_g_init = 102.74050220804774 + K_e_init = 3.32597273958481 + + K_e = states[5] + K_n = states[6] + K_g = states[7] + + Na_e = states[8] + Na_n = states[9] + Na_g = states[10] + + E_Na_n = R * temperature / F * np.log(Na_e/Na_n) #4 + E_K_n = R * temperature / F * np.log(K_e/K_n) #5 + + E_Na_g = R * temperature / F * np.log(Na_e/Na_g) #4 + E_K_g = R * temperature / F * np.log(K_e/K_g) #5 + E_K_init = R * temperature / F * np.log(K_e_init/K_g_init) #5 + + alpha_m = 0.1 * (states[3] + 40.0)/(1.0 - math.exp(-(states[3] + 40.0) / 10.0)) + beta_m = 4.0 * math.exp(-(states[3] + 65.0) / 18.0) + + alpha_h = 0.07 * math.exp(-(states[3] + 65.0) / 20.0) + beta_h = 1.0 / (1.0 + math.exp(-(states[3] + 35.0) / 10.0)) + + alpha_n = 0.01 * (states[3] + 55.0)/(1.0 - math.exp(-(states[3] + 55.0) / 10.0)) + beta_n = 0.125 * math.exp(-(states[3] + 65) / 80.0) + + # Expressions for the m gate component + values[0] = (1 - states[0])*alpha_m - states[0]*beta_m + + # Expressions for the h gate component + values[1] = (1 - states[1])*alpha_h - states[1]*beta_h + + # Expressions for the n gate component + values[2] = (1 - states[2])*alpha_n - states[2]*beta_n + + # Expressions for the Membrane component + i_Stim = parameters[7] * np.exp(-np.mod(t, 20.0)/2.0) + + i_pump_n = parameters[10] / ((1 + parameters[8] / K_e) ** 2 \ + * (1 + parameters[9] / Na_n) ** 3) + + i_pump_g = parameters[11] / ((1 + parameters[8] / K_e) ** 2 \ + * (1 + parameters[9] / Na_g) ** 3) + + # set conductance + dphi = states[4] - E_K_g + A = 1 + np.exp(18.4/42.4) # shorthand + B = 1 + np.exp(-(0.1186e3 + E_K_init)/0.0441e3) # shorthand + C = 1 + np.exp((dphi + 0.0185e3)/0.0425e3) # shorthand + D = 1 + np.exp(-(0.1186e3 + states[4])/0.0441e3) # shorthand + g_Kir = np.sqrt(K_e/K_e_init)*(A*B)/(C*D) + + # define and return current + I_Kir = parameters[5]*g_Kir*(states[4] - E_K_g) # umol/(cm^2*ms) + + # Expressions for the Sodium channel component + i_Na_n = (parameters[2] + parameters[0]*states[1]*math.pow(states[0], 3) + i_Stim) * \ + (states[3] - E_Na_n) + 3 * i_pump_n + + # Expressions for the Potassium channel component + i_K_n = (parameters[3] + parameters[1]*math.pow(states[2], 4)) * \ + (states[3] - E_K_n) - 2 * i_pump_n + + # Expressions for the Sodium channel component + i_Na_g = parameters[4] * (states[4] - E_Na_g) + 3 * i_pump_g + + # Expressions for the Potassium channel component + i_K_g = I_Kir - 2 * i_pump_g + + # Expression for phi_M_n + values[3] = (- i_K_n - i_Na_n)/parameters[6] + + # Expression for phi_M_g + values[4] = (- i_K_g - i_Na_g)/parameters[6] + + # Expression for K_e + values[5] = i_K_n * surface / (F * ECS_vol) \ + + i_K_g * surface / (F * ECS_vol) + + # Expression for K_n + values[6] = - i_K_n * surface / (F * ICS_vol) + + # Expression for K_g + values[7] = - i_K_g * surface / (F * ICS_vol) + + # Expression for Na_e + values[8] = i_Na_n * surface / (F * ECS_vol) \ + + i_Na_g * surface / (F * ECS_vol) + + # Expression for Na_n + values[9] = - i_Na_n * surface / (F * ICS_vol) + + # Expression for Na_g + values[10] = - i_Na_g * surface / (F * ICS_vol) diff --git a/emix-simulations/run_calibration_two_tags.py b/emix-simulations/run_calibration_two_tags.py new file mode 100644 index 0000000..a286313 --- /dev/null +++ b/emix-simulations/run_calibration_two_tags.py @@ -0,0 +1,158 @@ +#!/usr/bin/python3 + +import os +import sys +import time + +from dolfin import * +import numpy as np + +from collections import namedtuple + +from knpemidg import Solver +import mm_hh as mm_hh +import mm_glial as mm_glial +import mm_leak as mm_leak + +# Define colors for printing +class bcolors: + HEADER = '\033[95m' + OKBLUE = '\033[94m' + OKCYAN = '\033[96m' + OKGREEN = '\033[92m' + WARNING = '\033[93m' + FAIL = '\033[91m' + ENDC = '\033[0m' + BOLD = '\033[1m' + UNDERLINE = '\033[4m' + +if __name__ == "__main__": + + # Resolution factor of mesh + for resolution in [0]: + + # Time variables PDEs + dt = 0.1 # global time step (ms) + Tstop = 2 # global end time (ms) + t = Constant(0.0) # time constant + + # Time variables ODEs + n_steps_ODE = 20 # number of ODE steps + + # Physical parameters (PDEs) + C_M = 2.0 # capacitance + temperature = 300e3 # temperature (m K) + R = 8.314e3 # Gas Constant (m J/(K mol)) + F = 96485e3 # Faraday's constant (mC/ mol) + + D_Na = Constant(1.33e-8) # diffusion coefficients Na (cm/ms) + D_K = Constant(1.96e-8) # diffusion coefficients K (cm/ms) + D_Cl = Constant(2.03e-8) # diffusion coefficients Cl (cm/ms) + + psi = F / (R * temperature) # shorthand + C_phi = C_M / dt # shorthand + + # Initial values potassium + K_e_init = 3.3236967382613933 + K_n_init = 124.15397583492471 + K_g_init = 102.75563828644862 + + # Initial values sodium + Na_e_init = 100.71925900028181 + Na_n_init = 12.838513108606818 + Na_g_init = 12.39731187972181 + + # Initial values chloride + Cl_e_init = Na_e_init + K_e_init + Cl_n_init = Na_n_init + K_n_init + Cl_g_init = Na_g_init + K_g_init + + # Set physical parameters + params = namedtuple('params', ('dt', 'n_steps_ODE', 'F', 'psi', \ + 'C_phi', 'C_M', 'R', 'temperature'))(dt, \ + n_steps_ODE, F, psi, C_phi, C_M, R, temperature) + + # diffusion coefficients for each sub-domain + D_Na_sub = {0:D_Na, 1:D_Na, 2:D_Na} + D_K_sub = {0:D_K, 1:D_K, 2:D_K} + D_Cl_sub = {0:D_Cl, 1:D_Cl, 2:D_Cl} + + # initial concentrations for each sub-domain + Na_init_sub = {0:Constant(Na_e_init), 1:Constant(Na_n_init), 2:Constant(Na_g_init)} + K_init_sub = {0:Constant(K_e_init), 1:Constant(K_n_init), 2:Constant(K_g_init)} + Cl_init_sub = {0:Constant(Cl_e_init), 1:Constant(Cl_n_init), 2:Constant(Cl_g_init)} + c_init_sub_type = 'constant' + + # Create ions (channel conductivity is set below in the membrane model) + Na = {'c_init_sub':Na_init_sub, 'c_init_sub_type':c_init_sub_type, + 'bdry': Constant((0, 0)), 'z':1.0, 'name':'Na', 'D_sub':D_Na_sub} + K = {'c_init_sub':K_init_sub, 'c_init_sub_type':c_init_sub_type, + 'bdry': Constant((0, 0)), 'z':1.0, 'name':'K', 'D_sub':D_K_sub} + Cl = {'c_init_sub':Cl_init_sub, 'c_init_sub_type':c_init_sub_type, + 'bdry': Constant((0, 0)), 'z':-1.0, 'name':'Cl', 'D_sub':D_Cl_sub} + + # Create ion list. NB! The last ion in list will be eliminated, and + # should be the ion with the smallest diffusion coefficient D_k + ion_list = [K, Cl, Na] + + # Membrane parameters + g_syn_bar = 0.0 # synaptic conductivity (mS/cm**2) + + # Set stimulus ODE + stimulus = {'stim_amplitude': g_syn_bar} + stimulus_locator = lambda x: (x[0] < 7.0e-4) + + stim_params = namedtuple('membrane_params', ('g_syn_bar', \ + 'stimulus', 'stimulus_locator'))(g_syn_bar, \ + stimulus, stimulus_locator) + + # Get mesh, subdomains, surfaces paths + mesh_prefix = "meshes/3D_two_tags/" + mesh_path = mesh_prefix + "mesh_" + str(resolution) + ".xml" + subdomains_path = mesh_prefix + "subdomains_" + str(resolution) + ".xml" + surfaces_path = mesh_prefix + "surfaces_" + str(resolution) + ".xml" + + # Generate mesh if it does not exist + if not os.path.isfile(mesh_path): + from make_mesh_3D_two_tags import main + main(["-r", str(resolution), "-d", mesh_prefix]) + + mesh = Mesh(mesh_path) + subdomains = MeshFunction('size_t', mesh, subdomains_path) + surfaces = MeshFunction('size_t', mesh, surfaces_path) + + # Set solver parameters EMI (True is direct, and False is iterate) + direct_emi = False + rtol_emi = 1E-5 + atol_emi = 1E-40 + threshold_emi = 0.9 + + # Set solver parameters KNP (True is direct, and False is iterate) + direct_knp = False + rtol_knp = 1E-7 + atol_knp = 2E-40 + threshold_knp = 0.75 + + # Set parameters + solver_params = namedtuple('solver_params', ('direct_emi', + 'direct_knp', 'resolution', + 'rtol_emi', 'rtol_knp', + 'atol_emi', 'atol_knp', + 'threshold_emi', 'threshold_knp' + ))(direct_emi, direct_knp, resolution, \ + rtol_emi, rtol_knp, atol_emi, atol_knp, \ + threshold_emi, threshold_knp) + + # File for results + fname = "results/data/calibration_two_tags/" + + # Dictionary with membrane models (key is facet tag, value is ode model) + ode_models = {1: mm_hh, 2:mm_glial} + + # Solve system + S = Solver(params, ion_list) # create solver + S.setup_domain(mesh, subdomains, surfaces) # setup meshes + S.setup_parameters() # setup physical parameters + S.setup_FEM_spaces() # setup function spaces and numerical parameters + S.setup_membrane_model(stim_params, ode_models) # setup membrane model(s) + S.solve_system_active(Tstop, t, solver_params, filename=fname) # solve diff --git a/emix-simulations/run_two_tags.py b/emix-simulations/run_two_tags.py new file mode 100644 index 0000000..72315ab --- /dev/null +++ b/emix-simulations/run_two_tags.py @@ -0,0 +1,162 @@ +#!/usr/bin/python3 + +import os +import sys +import time + +from dolfin import * +import numpy as np + +from collections import namedtuple + +from knpemidg import Solver +import mm_hh as mm_hh +import mm_glial as mm_glial +import mm_leak as mm_leak + +# Define colors for printing +class bcolors: + HEADER = '\033[95m' + OKBLUE = '\033[94m' + OKCYAN = '\033[96m' + OKGREEN = '\033[92m' + WARNING = '\033[93m' + FAIL = '\033[91m' + ENDC = '\033[0m' + BOLD = '\033[1m' + UNDERLINE = '\033[4m' + +if __name__ == "__main__": + + # Time variables PDEs + dt = 0.1 # global time step (ms) + Tstop = 10 # global end time (ms) + t = Constant(0.0) # time constant + + # Time variables ODEs + n_steps_ODE = 20 # number of ODE steps + + # Physical parameters (PDEs) + C_M = 2.0 # capacitance + temperature = 300e3 # temperature (m K) + R = 8.314e3 # Gas Constant (m J/(K mol)) + F = 96485e3 # Faraday's constant (mC/ mol) + + D_Na = Constant(1.33e-8) # diffusion coefficients Na (cm/ms) + D_K = Constant(1.96e-8) # diffusion coefficients K (cm/ms) + D_Cl = Constant(2.03e-8) # diffusion coefficients Cl (cm/ms) + + psi = F / (R * temperature) # shorthand + C_phi = C_M / dt # shorthand + + n = 20 + fname_callibrate = 'results/data/calibration_two_tags/results.h5' + hdf5file = HDF5File(MPI.comm_world, fname_callibrate, "r") + + mesh = Mesh() + subdomains = MeshFunction("size_t", mesh, 2) + surfaces = MeshFunction("size_t", mesh, 1) + hdf5file.read(mesh, '/mesh', False) + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + P1 = FiniteElement('DG', mesh.ufl_cell(), 1) + W = FunctionSpace(mesh, MixedElement(2*[P1])) + V = FunctionSpace(mesh, P1) + + u = Function(W) + v = Function(V) + w = Function(V) + + Na_init = Function(V) + K_init = Function(V) + Cl_init = Function(V) + + hdf5file.read(u, "/concentrations/vector_" + str(n)) + + # K concentrations + assign(K_init, u.sub(0)) + # Cl concentrations + assign(Cl_init, u.sub(1)) + + # Na concentrations + hdf5file.read(v, "/elim_concentration/vector_" + str(n)) + assign(Na_init, v) + + hdf5file.close() + + # Set physical parameters + params = namedtuple('params', ('dt', 'n_steps_ODE', 'F', 'psi', \ + 'C_phi', 'C_M', 'R', 'temperature'))(dt, \ + n_steps_ODE, F, psi, C_phi, C_M, R, temperature) + + # diffusion coefficients for each sub-domain + D_Na_sub = {0:D_Na, 1:D_Na, 2:D_Na} + D_K_sub = {0:D_K, 1:D_K, 2:D_K} + D_Cl_sub = {0:D_Cl, 1:D_Cl, 2:D_Cl} + + # initial concentrations for each ion + Na_init_sub = Na_init + K_init_sub = K_init + Cl_init_sub = Cl_init + c_init_sub_type = 'function' # type - function or constant + + # Create ions (channel conductivity is set below in the membrane model) + Na = {'c_init_sub':Na_init_sub, 'c_init_sub_type':c_init_sub_type, + 'bdry': Constant((0, 0)), 'z':1.0, 'name':'Na', 'D_sub':D_Na_sub} + K = {'c_init_sub':K_init_sub, 'c_init_sub_type':c_init_sub_type, + 'bdry': Constant((0, 0)), 'z':1.0, 'name':'K', 'D_sub':D_K_sub} + Cl = {'c_init_sub':Cl_init_sub, 'c_init_sub_type':c_init_sub_type, + 'bdry': Constant((0, 0)), 'z':-1.0, 'name':'Cl', 'D_sub':D_Cl_sub} + + # Create ion list. NB! The last ion in list will be eliminated, and + # should be the ion with the smallest diffusion coefficient D_k + ion_list = [K, Cl, Na] + + # Membrane parameters + g_syn_bar = 5.0 # synaptic conductivity (mS/cm**2) + #g_syn_bar = 0.0 # synaptic conductivity (mS/cm**2) + + # Set stimulus ODE + stimulus = {'stim_amplitude': g_syn_bar} + stimulus_locator = lambda x: (x[0] < 7.0e-4) + + stim_params = namedtuple('membrane_params', ('g_syn_bar', \ + 'stimulus', 'stimulus_locator'))(g_syn_bar, \ + stimulus, stimulus_locator) + + # Set solver parameters EMI (True is direct, and False is iterate) + direct_emi = False + rtol_emi = 1E-5 + atol_emi = 1E-40 + threshold_emi = 0.9 + + # Set solver parameters KNP (True is direct, and False is iterate) + direct_knp = False + rtol_knp = 1E-7 + atol_knp = 2E-40 + threshold_knp = 0.75 + + # Set parameters + solver_params = namedtuple('solver_params', ('direct_emi', + 'direct_knp', 'resolution', + 'rtol_emi', 'rtol_knp', + 'atol_emi', 'atol_knp', + 'threshold_emi', 'threshold_knp' + ))(direct_emi, direct_knp, 0, \ + rtol_emi, rtol_knp, atol_emi, atol_knp, \ + threshold_emi, threshold_knp) + + # File for results + fname = "results/data/two_tags/" + + # Dictionary with membrane models (key is facet tag, value is ode model) + ode_models = {1: mm_hh, 2:mm_glial} + + # Solve system + S = Solver(params, ion_list) # create solver + S.setup_domain(mesh, subdomains, surfaces) # setup meshes + S.setup_parameters() # setup physical parameters + S.setup_FEM_spaces() # setup function spaces and numerical parameters + S.setup_membrane_model(stim_params, ode_models) # setup membrane model(s) + S.solve_system_active(Tstop, t, solver_params, filename=fname) # solve