From a418fd300a44a7ed9b9f25172691a84211f668b5 Mon Sep 17 00:00:00 2001 From: adajel Date: Thu, 8 Aug 2024 12:01:42 +0200 Subject: [PATCH] add code for calibrating ODE models in EMIx simulations --- emix-simulations/make_figures_calibration.py | 634 +++++++++++++++++++ emix-simulations/make_mesh_3D.py | 117 ++++ emix-simulations/mm_glial.py | 170 +++++ emix-simulations/mm_glial_calibration_ODE.py | 211 ++++++ emix-simulations/mm_hh.py | 174 +++++ emix-simulations/mm_hh_calibration_ODE.py | 219 +++++++ emix-simulations/run_calibration.py | 171 +++++ emix-simulations/run_calibration_ODE.py | 118 ++++ 8 files changed, 1814 insertions(+) create mode 100644 emix-simulations/make_figures_calibration.py create mode 100755 emix-simulations/make_mesh_3D.py create mode 100644 emix-simulations/mm_glial.py create mode 100644 emix-simulations/mm_glial_calibration_ODE.py create mode 100644 emix-simulations/mm_hh.py create mode 100644 emix-simulations/mm_hh_calibration_ODE.py create mode 100644 emix-simulations/run_calibration.py create mode 100644 emix-simulations/run_calibration_ODE.py diff --git a/emix-simulations/make_figures_calibration.py b/emix-simulations/make_figures_calibration.py new file mode 100644 index 0000000..00cab83 --- /dev/null +++ b/emix-simulations/make_figures_calibration.py @@ -0,0 +1,634 @@ +import matplotlib as mpl +from mpl_toolkits.axes_grid1 import make_axes_locatable + +import matplotlib.pyplot as plt +import numpy as np +import os +import sys + +from fenics import * +import string + +from surf_plot.vtk_io import DltWriter +from surf_plot.dlt_embedding import P0surf_to_DLT0_map + +from knpemidg import pcws_constant_project +from knpemidg import interface_normal, plus, minus + +JUMP = lambda f, n: minus(f, n) - plus(f, n) + +# set font & text parameters +font = {'family' : 'serif', + 'weight' : 'bold', + 'size' : 13} + +plt.rc('font', **font) +plt.rc('text', usetex=True) +mpl.rcParams['image.cmap'] = 'jet' + +path = 'results/data/' + +def write_to_pvd(dt, T, fname): + # read data file + hdf5file = HDF5File(MPI.comm_world, fname, "r") + + mesh = Mesh() + subdomains = MeshFunction("size_t", mesh, 2) + surfaces = MeshFunction("size_t", mesh, 1) + hdf5file.read(mesh, '/mesh', False) + mesh.coordinates()[:] *= 1e4 + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + P1 = FiniteElement('CG', mesh.ufl_cell(), 1) + W = FunctionSpace(mesh, MixedElement(2*[P1])) + V = FunctionSpace(mesh, P1) + + u = Function(W) + v = Function(V) + w = Function(V) + + Na = Function(V) + K = Function(V) + Cl = Function(V) + phi = Function(V) + + f_phi = File('results/pvd/pot.pvd') + f_K = File('results/pvd/K.pvd') + f_Na = File('results/pvd/Na.pvd') + f_Cl = File('results/pvd/Cl.pvd') + + for n in range(1, int(T/dt)): + + # read file + hdf5file.read(u, "/concentrations/vector_" + str(n)) + + # K concentrations + assign(K, u.sub(0)) + # Cl concentrations + assign(Cl, u.sub(1)) + + # Na concentrations + hdf5file.read(v, "/elim_concentration/vector_" + str(n)) + assign(Na, v) + + # potential + hdf5file.read(w, "/potential/vector_" + str(n)) + assign(phi, w) + + f_Na << Na, n + f_K << K, n + f_Cl << Cl, n + f_phi << phi, n + + return + + +def get_time_series(dt, T, fname, x, y, z): + # read data file + hdf5file = HDF5File(MPI.comm_world, fname, "r") + + mesh = Mesh() + subdomains = MeshFunction("size_t", mesh, 2) + surfaces = MeshFunction("size_t", mesh, 1) + hdf5file.read(mesh, '/mesh', False) + mesh.coordinates()[:] *= 1e4 + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + P1 = FiniteElement('CG', mesh.ufl_cell(), 1) + W = FunctionSpace(mesh, MixedElement(2*[P1])) + V = FunctionSpace(mesh, P1) + + u = Function(W) + v = Function(V) + w = Function(V) + + f_Na = Function(V) + f_K = Function(V) + f_Cl = Function(V) + f_phi = Function(V) + + Na = [] + K = [] + Cl = [] + phi = [] + + for n in range(1, int(T/dt)): + + # read file + hdf5file.read(u, "/concentrations/vector_" + str(n)) + + # K concentrations + assign(f_K, u.sub(0)) + K.append(f_K(x, y, z)) + + # Cl concentrations + assign(f_Cl, u.sub(1)) + Cl.append(f_Cl(x, y, z)) + + # Na concentrations + hdf5file.read(v, "/elim_concentration/vector_" + str(n)) + assign(f_Na, v) + Na.append(f_Na(x, y, z)) + + # potential + hdf5file.read(w, "/potential/vector_" + str(n)) + assign(f_phi, w) + phi.append(f_phi(x, y, z)) + + return Na, K, Cl, phi + +def get_time_series_membrane(dt, T, fname, x, y, z): + # read data file + hdf5file = HDF5File(MPI.comm_world, fname, "r") + + mesh = Mesh() + subdomains = MeshFunction("size_t", mesh, 2) + surfaces = MeshFunction("size_t", mesh, 1) + hdf5file.read(mesh, '/mesh', False) + mesh.coordinates()[:] *= 1e4 + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + x_min = x - 0.5; x_max = x + 0.1 + y_min = y - 0.5; y_max = y + 0.1 + z_min = z; z_max = z + 0.04 + + # define one facet to 10 for getting membrane potential + for facet in facets(mesh): + x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()] + point_1 = y_min <= x[1] <= y_max \ + and x_min <= x[0] <= x_max \ + and z_min <= x[2] <= z_max + if point_1 and surfaces[facet] == 1: + #if surfaces[facet] == 1: + print(x[0], x[1], x[2]) + surfaces[facet] = 10 + + surfacesfile = File('surfaces_plot.pvd') + surfacesfile << surfaces + + # define function space of piecewise constants on interface gamma for solution to ODEs + Q = FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0) + phi_M = Function(Q) + E_Na = Function(Q) + E_K = Function(Q) + + # interface normal + n_g = interface_normal(subdomains, mesh) + + dS = Measure('dS', domain=mesh, subdomain_data=surfaces) + iface_size = assemble(Constant(1)*dS(10)) + + P1 = FiniteElement('DG', mesh.ufl_cell(), 1) + W = FunctionSpace(mesh, MixedElement(2*[P1])) + V = FunctionSpace(mesh, P1) + + v = Function(W) + w_Na = Function(V) + w_phi = Function(V) + + w_phi_ODE = Function(Q) + f_phi_ODE = Function(Q) + + f_phi = Function(V) + f_Na = Function(V) + f_K = Function(V) + + phi_M_s = [] + phi_M_ODE_s = [] + E_Na_s = [] + E_K_s = [] + + z_Na = 1; z_K = 1; temperature = 300e3; F = 96485e3; R = 8.314e3 + + for n in range(1, int(T/dt)): + + # update membrane potential + hdf5file.read(w_phi, "/potential/vector_" + str(n)) + assign(f_phi, w_phi) + phi_M_step = JUMP(f_phi, n_g) + assign(phi_M, pcws_constant_project(phi_M_step, Q)) + phi_M_s.append(assemble(1.0/iface_size*avg(phi_M)*dS(10))) + + hdf5file.read(w_phi_ODE, "/membranepotential/vector_" + str(n)) + assign(f_phi_ODE, w_phi_ODE) + phi_M_ODE_s.append(assemble(1.0/iface_size*avg(f_phi_ODE)*dS(10))) + + # K concentrations + hdf5file.read(v, "/concentrations/vector_" + str(n)) + assign(f_K, v.sub(0)) + E = R * temperature / (F * z_K) * ln(plus(f_K, n_g) / minus(f_K, n_g)) + assign(E_K, pcws_constant_project(E, Q)) + E_K_ = assemble(1.0/iface_size*avg(E_K)*dS(10)) + E_K_s.append(E_K_) + + # Na concentrations + hdf5file.read(w_Na, "/elim_concentration/vector_" + str(n)) + assign(f_Na, w_Na) + E = R * temperature / (F * z_Na) * ln(plus(f_Na, n_g) / minus(f_Na, n_g)) + assign(E_Na, pcws_constant_project(E, Q)) + E_Na_ = assemble(1.0/iface_size*avg(E_Na)*dS(10)) + E_Na_s.append(E_Na_) + + return phi_M_s, phi_M_ODE_s, E_Na_s, E_K_s + +def plot_3D_concentration(res, T, dt, fname): + + time = np.arange(0, T-dt, dt) + + # original mesh + #x_M = 20 + #y_M = 0.2 + #z_M = 0.63 + + # new mesh + x_M = 8 + y_M = 0.8 + z_M = 1.23 + + # 0.05 um above axon A (ECS) + x_e = x_M; y_e = y_M + 0.1; z_e = z_M + # mid point inside axon A (ICS) + x_i = x_M; y_i = y_M - 0.1; z_i = z_M + + # trace concentrations + phi_M, phi_M_ODE, E_Na, E_K = get_time_series_membrane(dt, T, fname, x_M, y_M, z_M) + + # bulk concentrations + Na_e, K_e, Cl_e, _ = get_time_series(dt, T, fname, x_e, y_e, z_e) + Na_i, K_i, Cl_i, _ = get_time_series(dt, T, fname, x_i, y_i, z_i) + + ################################################################# + # get data axons BC are stimulated + + # Concentration plots + fig = plt.figure(figsize=(12*0.9,12*0.9)) + ax = plt.gca() + + ax1 = fig.add_subplot(3,3,1) + plt.title(r'Na$^+$ concentration (ECS)') + plt.ylabel(r'[Na]$_e$ (mM)') + #plt.ylim([90, 110]) + plt.plot(Na_e, linewidth=3, color='b') + + ax3 = fig.add_subplot(3,3,2) + plt.title(r'K$^+$ concentration (ECS)') + plt.ylabel(r'[K]$_e$ (mM)') + plt.plot(K_e, linewidth=3, color='b') + + ax3 = fig.add_subplot(3,3,3) + plt.title(r'Cl$^-$ concentration (ECS)') + plt.ylabel(r'[Cl]$_e$ (mM)') + plt.plot(Cl_e, linewidth=3, color='b') + + ax2 = fig.add_subplot(3,3,4) + plt.title(r'Na$^+$ concentration (ICS)') + plt.ylabel(r'[Na]$_i$ (mM)') + #plt.ylim([99, 101]) + plt.plot(Na_i,linewidth=3, color='r') + + ax2 = fig.add_subplot(3,3,5) + plt.title(r'K$^+$ concentration (ICS)') + plt.ylabel(r'[K]$_i$ (mM)') + #plt.ylim([4, 4.2]) + plt.plot(K_i,linewidth=3, color='r') + + ax2 = fig.add_subplot(3,3,6) + plt.title(r'Cl$^-$ concentration (ICS)') + plt.ylabel(r'[Cl]$_i$ (mM)') + #plt.ylim([103, 105]) + plt.plot(Cl_i,linewidth=3, color='r') + + ax5 = fig.add_subplot(3,3,7) + plt.title(r'Membrane potential') + plt.ylabel(r'$\phi_M$ (mV)') + #plt.ylim([-71.3, -71.2]) + plt.xlabel(r'time (ms)') + plt.plot(phi_M, linewidth=3) + + ax6 = fig.add_subplot(3,3,8) + plt.title(r'Na$^+$ reversal potential') + plt.ylabel(r'E$_Na$ (mV)') + plt.xlabel(r'time (ms)') + plt.plot(E_K, linewidth=3) + plt.plot(E_Na, linewidth=3) + + ax7 = fig.add_subplot(3,3,9) + plt.title(r'Membrane potential ODE') + #plt.ylim([-71.3, -71.2]) + plt.plot(phi_M_ODE, linewidth=3) + + print("membrane potential", phi_M[-1]) + print("membrane potential", phi_M_ODE[-1]) + #print("Na_e", Na_e[-1]) + #print("Na_i", Na_i[-1]) + #print("K_i", K_i[-1]) + #print("K_e", K_e[-1]) + #print("Cl_i", Cl_i[-1]) + #print("Cl_e", Cl_e[-1]) + + # make pretty + ax.axis('off') + plt.tight_layout() + + # save figure to file + plt.savefig('results/figures/callibrate.svg', format='svg') + + f_phi_M = open('phi_M_3D.txt', "w") + for p in phi_M: + f_phi_M.write("%.10f \n" % p*1000) + f_phi_M.close() + + return + +def plot_surface(fname, T, dt): + + # read data file + hdf5file = HDF5File(MPI.comm_world, fname, "r") + + mesh = Mesh() + subdomains = MeshFunction("size_t", mesh, 2) + surfaces = MeshFunction("size_t", mesh, 1) + hdf5file.read(mesh, '/mesh', False) + mesh.coordinates()[:] *= 1e7 + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + ds = Measure('ds', domain=mesh, subdomain_data=surfaces) + dS = Measure('dS', domain=mesh, subdomain_data=surfaces) + + surface_tags = (1, 2) + + surface_mesh, L, cell2Ldofs = P0surf_to_DLT0_map(surfaces, tags=surface_tags) + + P1 = FiniteElement('DG', mesh.ufl_cell(), 1) + V = FunctionSpace(mesh, P1) + + u = Function(V) + uh = Function(V) + + # read file + n = 0 + hdf5file.read(u, "/potential/vector_" + str(n)) + + # K concentrations + assign(uh, u) + + # To get the restriction as we want we build a P0(mesh) tagging function + dx = Measure('dx', domain=mesh, subdomain_data=subdomains) + + V = FunctionSpace(mesh, 'DG', 0) + v = TestFunction(V) + + subdomain_tag = Function(V) + hK = CellVolume(mesh) + # Project coefs to P0 based on subdomain + assemble((1/hK)*inner(v, Constant(1))*dx(1) + (1/hK)*inner(v, Constant(2))*dx(2), subdomain_tag.vector()) + # With the cell taggging function we will pick sides based on the marker value + small_side = lambda u, K=subdomain_tag: conditional(gt(K('+'), K('-')), u('+'), u('-')) + big_side = lambda u, K=subdomain_tag: conditional(le(K('+'), K('-')), u('+'), u('-')) + + # So finally for the projection of P0 data + dl = TestFunction(L) + fK = avg(FacetArea(mesh)) + # To get the quantity on the surface + vh = Function(L) + + fK = FacetArea(mesh) + dS0 = dS.reconstruct(metadata={'quadrature_degree': 0}) + ds0 = ds.reconstruct(metadata={'quadrature_degree': 0}) + + #for truth, side in enumerate((big_side, small_side), 1): + #print(side) + # NOTE: that DLT is "continuous" on the facet, avg is just to shout up FFC + #assemble(sum((1.0e3*1/avg(fK))*inner(small_side(uh) - big_side(uh), avg(dl))*dS0(tag) + #for tag in surface_tags), tensor=vh.vector()) + #as_backend_type(vh.vector()).update_ghost_values() + L_form = sum((1/avg(fK))*inner(small_side(uh) - big_side(uh), avg(dl))*dS0(tag) + for tag in surface_tags) + + assemble(L_form, tensor=vh.vector()) + as_backend_type(vh.vector()).update_ghost_values() + + values = vh.vector().get_local()[cell2Ldofs] + cell_data = {'uh': values} + + dlt = DltWriter('testing/bar', mesh, surface_mesh) + with dlt as output: + output.write(cell_data, t=0) + + # Some more time steps + for n in range(1, int(T/dt)): + # read file + hdf5file.read(u, "/potential/vector_" + str(n)) + # K concentrations + assign(uh, u) + + assemble(L_form, tensor=vh.vector()) + as_backend_type(vh.vector()).update_ghost_values() + + # Now we build the data for P0 function on the mesh + values[:] = vh.vector().get_local()[cell2Ldofs] + print(np.linalg.norm(values)) + output.write(cell_data, t=n) + + # Check the data + #true = truth*np.ones_like(values) + #assert np.linalg.norm(true - values) < 1E-12 + + return + +def plot_surface_time(fname, T, dt): + + # read data file + hdf5file = HDF5File(MPI.comm_world, fname, "r") + + mesh = Mesh() + subdomains = MeshFunction("size_t", mesh, 2) + surfaces = MeshFunction("size_t", mesh, 1) + hdf5file.read(mesh, '/mesh', False) + mesh.coordinates()[:] *= 1e7 + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + ds = Measure('ds', domain=mesh, subdomain_data=surfaces) + dS = Measure('dS', domain=mesh, subdomain_data=surfaces) + + surface_tags = (1, 2) + + surface_mesh, L, cell2Ldofs = P0surf_to_DLT0_map(surfaces, tags=surface_tags) + + P1 = FiniteElement('DG', mesh.ufl_cell(), 1) + V = FunctionSpace(mesh, P1) + + u = Function(V) + uh = Function(V) + + # read file + hdf5file.read(u, "/potential/vector_" + str(0)) + # K concentrations + assign(uh, u) + + # To get the restriction as we want we build a P0(mesh) tagging function + dx = Measure('dx', domain=mesh, subdomain_data=subdomains) + + V = FunctionSpace(mesh, 'DG', 0) + v = TestFunction(V) + + subdomain_tag = Function(V) + hK = CellVolume(mesh) + # Project coefs to P0 based on subdomain + assemble((1/hK)*inner(v, Constant(1))*dx(1) + (1/hK)*inner(v, Constant(2))*dx(2), subdomain_tag.vector()) + # With the cell taggging function we will pick sides based on the marker value + small_side = lambda u, K=subdomain_tag: conditional(gt(K('+'), K('-')), u('+'), u('-')) + big_side = lambda u, K=subdomain_tag: conditional(le(K('+'), K('-')), u('+'), u('-')) + + # So finally for the projection of P0 data + dl = TestFunction(L) + fK = avg(FacetArea(mesh)) + # To get the quantity on the surface + vh = Function(L) + + fK = FacetArea(mesh) + dS0 = dS.reconstruct(metadata={'quadrature_degree': 0}) + ds0 = ds.reconstruct(metadata={'quadrature_degree': 0}) + + #for truth, side in enumerate((big_side, small_side), 1): + #print(side) + # NOTE: that DLT is "continuous" on the facet, avg is just to shout up FFC + L_form = sum((1/avg(fK))*inner(small_side(uh) - big_side(uh), avg(dl))*dS0(tag) + for tag in surface_tags) + + assemble(L_form, tensor=vh.vector()) + as_backend_type(vh.vector()).update_ghost_values() + + # At this point vh has the quantitiy of interest. Iguess it typically is + # some restriction, see `example_P0.py` for how to get them + # For plotting we always want grab subset of DLT dofs + values = vh.vector().get_local()[cell2Ldofs] + + # And dump - note that we can have several quantities (e.g. neg_ih) in the + # same file + cell_data = {'uh': values, 'neg_uh': values} + + dlt = DltWriter(f'testing/cux_world{mesh.mpi_comm().size}', mesh, surface_mesh) + with dlt as output: + output.write(cell_data, t=0) + + # Some more time steps + for n in range(1, int(T/dt)): + # read file + hdf5file.read(u, "/potential/vector_" + str(n)) + # potential + assign(uh, u) + + assemble(L_form, tensor=vh.vector()) + as_backend_type(vh.vector()).update_ghost_values() + + # Now we build the data for P0 function on the mesh + values[:] = vh.vector().get_local()[cell2Ldofs] + print(np.linalg.norm(values)) + output.write(cell_data, t=n) + return + +def get_velocity(fname, T, dt): + # read data file + hdf5file = HDF5File(MPI.comm_world, fname, "r") + + mesh = Mesh() + subdomains = MeshFunction("size_t", mesh, 2) + surfaces = MeshFunction("size_t", mesh, 1) + hdf5file.read(mesh, '/mesh', False) + mesh.coordinates()[:] *= 1e6 + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + #x_min = 0.1; x_max = x + 0.1 + #y_min = 0.1; y_max = y + 0.1 + #z_min = 0.1; z_max = z + 0.1 + + # define one facet to 10 for getting membrane potential + for facet in facets(mesh): + x = [facet.midpoint().x(), facet.midpoint().y(), facet.midpoint().z()] + # mark top point + if surfaces[facet] == 2 and (x[1] > 539): + print(x[0], x[1], x[2]) + surfaces[facet] = 10 + # mark bottom point + if surfaces[facet] == 2 and (-0.004 < x[1] < 0.0045): + print(x[0], x[1], x[2]) + surfaces[facet] = 20 + + # define function space of piecewise constants on interface gamma for solution to ODEs + Q = FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0) + phi_M = Function(Q) + + # interface normal + n_g = interface_normal(subdomains, mesh) + + dS = Measure('dS', domain=mesh, subdomain_data=surfaces) + iface_size_10 = assemble(Constant(1)*dS(10)) + iface_size_20 = assemble(Constant(1)*dS(20)) + + P1 = FiniteElement('DG', mesh.ufl_cell(), 1) + V = FunctionSpace(mesh, P1) + + w_phi = Function(V) + f_phi = Function(V) + + phi_M_s_1 = [] + phi_M_s_2 = [] + + time_point_1 = 0 + time_point_2 = 0 + + for n in range(1, int(T/dt)): + + # potential + hdf5file.read(w_phi, "/potential/vector_" + str(n)) + assign(f_phi, w_phi) + + # update membrane potential + phi_M_step = JUMP(f_phi, n_g) + assign(phi_M, pcws_constant_project(phi_M_step, Q)) + + if 1.0e3*assemble(1.0/iface_size_10*avg(phi_M)*dS(10)) > 0 and time_point_1 == 0: + time_point_1 = n + + if 1.0e3*assemble(1.0/iface_size_10*avg(phi_M)*dS(20)) > 0 and time_point_2 == 0: + time_point_2 = n + + print(time_point_1) + print(time_point_2) + + # if membrane potential has reach 0 in both points then break + if (time_point_1 > 0) and (time_point_2 > 0): + break + + delta_t = (time_point_2 - time_point_1)*1.0e-4 # s + delta_x = 539e-6 # m + + print("velocity (m/s)", delta_x/delta_t) + + return phi_M_s_1, phi_M_s_2 + +# create directory for figures +if not os.path.isdir('results/figures'): + os.mkdir('results/figures') + +# create figures +res_3D = '0' # mesh resolution for 3D axon bundle +dt = 0.1 +T = 10 + +#fname = 'results/data/EMIX/results.h5' +#plot_surface(fname, T, dt) +#plot_surface_time(fname, T, dt) + +fname = 'results/data/calibration/results.h5' +plot_3D_concentration(res_3D, T, dt, fname) + +#write_to_pvd(dt, T, fname) +#get_velocity(fname, T, dt) diff --git a/emix-simulations/make_mesh_3D.py b/emix-simulations/make_mesh_3D.py new file mode 100755 index 0000000..86a2df1 --- /dev/null +++ b/emix-simulations/make_mesh_3D.py @@ -0,0 +1,117 @@ +#!/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): + # 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] = 1 + 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] + surfaces[facet] += side_1 or side_2 or side_3 or side_4 or side_5 or side_6 + + 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"), + 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 * 105 * 2 ** resolution_factor + ny = 10 * 2 ** resolution_factor + nz = 10 * 2 ** resolution_factor + + # box mesh + mesh = BoxMesh(Point(0, 0.0, 0.0), Point(l * 105, 1.0, 1.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 * 105 - 5, 0.8, 0.8) + add_axon(mesh, subdomains, surfaces, a, b) + + # 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_glial.py b/emix-simulations/mm_glial.py new file mode 100644 index 0000000..056080d --- /dev/null +++ b/emix-simulations/mm_glial.py @@ -0,0 +1,170 @@ +# 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_init = -81.73542331243227 + + init_values = np.array([phi_M_init], dtype=np.float_) + + # State indices and limit checker + state_ind = dict([("V", 0)]) + + 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 = 0 # Na max conductivity (mS/cm**2) + g_K_bar = 0 # K max conductivity (mS/cm**2) + g_leak_Na = 0.1 # Na leak conductivity (mS/cm**2) + g_leak_K = 1.7 # K leak conductivity (mS/cm**2) + + m_K = 2.0 # threshold ECS K (mol/m^3) - yao 2011 + m_Na = 7.7 # threshold ICS Na (mol/m^3) - yao 2011 + I_max = 43 # max pump strength (muA/cm^2) + + K_i_init = 101.45093275680335 + K_e_init = 3.369320617474813 + + # Set initial parameter values + init_values = np.array([g_Na_bar, g_K_bar, \ + g_leak_Na, g_leak_K, \ + 0, 0, 0, 0, \ + 0, 0, 0, 0, 0, + m_K, m_Na, I_max, K_e_init, K_i_init], dtype=np.float_) + + # Parameter indices and limit checker + param_ind = dict([("g_Na_bar", 0), ("g_K_bar", 1), \ + ("g_leak_Na", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), \ + ("Cm", 6), ("stim_amplitude", 7), \ + ("I_ch_Na", 8), ("I_ch_K", 9), ("I_ch_Cl", 10), \ + ("K_e", 11), ("Na_i", 12), \ + ("m_K", 13), ("m_Na", 14), ("I_max", 15), \ + ("K_e_init", 16), ("K_i_init", 17)]) + + 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_inds = dict([("V", 0)]) + + 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", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), ("Cm", 6), \ + ("stim_amplitude", 7), ("I_ch_Na", 8), \ + ("I_ch_K", 9), ("I_ch_Cl", 10), ("K_e", 11), \ + ("Na_i", 12), ("m_K", 13), ("m_Na", 14), ("I_max", 15), \ + ("K_e_init", 16), ("K_i_init", 17)]) + + 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,) + + i_pump = parameters[15] / ((1 + parameters[13] / parameters[11]) ** 2 \ + * (1 + parameters[14] / parameters[12]) ** 3) + + # Physical parameters (PDEs) + temperature = 300e3 # temperature (m K) + R = 8.314e3 # Gas Constant (m J/(K mol)) + F = 96485e3 # Faraday's constant (mC/ mol) + + # set conductance + E_K_init = R * temperature / F * np.log(parameters[16]/parameters[17]) #5 + dphi = states[0] - parameters[5] + 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[0])/0.0441e3) # shorthand + g_Kir = np.sqrt(parameters[11]/parameters[16])*(A*B)/(C*D) + + # define and return current + i_Kir = parameters[3]*g_Kir*(states[0] - parameters[5]) # umol/(cm^2*ms) + + # Expressions for the Sodium channel component + i_Na = parameters[2] * (states[3] - parameters[4]) + 3 * i_pump + + # Expressions for the Potassium channel component + i_K = i_Kir - 2 * i_pump + + # set I_ch_Na + parameters[8] = i_Na + # set I_ch_K + parameters[9] = i_K + # set I_ch_Cl + parameters[10] = 0.0 + + # update membrane potential + values[0] = (- i_K - i_Na)/parameters[6] diff --git a/emix-simulations/mm_glial_calibration_ODE.py b/emix-simulations/mm_glial_calibration_ODE.py new file mode 100644 index 0000000..587abc2 --- /dev/null +++ b/emix-simulations/mm_glial_calibration_ODE.py @@ -0,0 +1,211 @@ +# 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_init = -81.73542331243227 + + # Initial values + #Na_i_init = 15.189 # Intracellular Na concentration + #Na_e_init = 99.959 # extracellular Na concentration + #K_i_init = 99.959 # intracellular K concentration + #K_e_init = 4.0900 # extracellular K concentration + + K_i_init = 101.45093275680335 + Na_i_init = 13.701595802750164 + K_e_init = 3.369320617474813 + Na_e_init = 100.67749185800223 + + init_values = np.array([phi_M_init, K_e_init, K_i_init, \ + Na_e_init, Na_i_init], dtype=np.float_) + + # State indices and limit checker + state_ind = dict([("V", 0), ("K_e", 1), ("K_i", 2), ("Na_e", 3), ("Na_i", 4)]) + + 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 + """ + + K_i_init = 99.959 # intracellular K concentration + K_e_init = 4.0900 # extracellular K concentration + + # Membrane parameters + g_Na_bar = 0 # Na max conductivity (mS/cm**2) + g_K_bar = 0 # K max conductivity (mS/cm**2) + g_leak_Na = 0.1 # Na leak conductivity (mS/cm**2) + g_leak_K = 1.7 # K leak conductivity (mS/cm**2) + + m_K = 2.0 # threshold ECS K (mol/m^3) - yao 2011 + m_Na = 7.7 # threshold ICS Na (mol/m^3) - yao 2011 + I_max = 43 # 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, g_leak_K, \ + 0, 0, C_M, 0, \ + 0, 0, 0, 0, 0, + m_K, m_Na, I_max, K_e_init, K_i_init], dtype=np.float_) + + # Parameter indices and limit checker + param_ind = dict([("g_Na_bar", 0), ("g_K_bar", 1), \ + ("g_leak_Na", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), \ + ("Cm", 6), ("stim_amplitude", 7), \ + ("I_ch_Na", 8), ("I_ch_K", 9), ("I_ch_Cl", 10), \ + ("K_e", 11), ("Na_i", 12), \ + ("m_K", 13), ("m_Na", 14), ("I_max", 15), + ("K_e_init", 16), ("K_i_init", 17)]) + + 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_inds = dict([("V", 0), ("K_e", 1), ("K_i", 2), ("Na_e", 3), ("Na_i", 4)]) + + 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", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), ("Cm", 6), \ + ("stim_amplitude", 7), ("I_ch_Na", 8), \ + ("I_ch_K", 9), ("I_ch_Cl", 10), ("K_e", 11), \ + ("Na_i", 12), ("m_K", 13), ("m_Na", 14), ("I_max", 15),\ + ("K_e_init", 16), ("K_i_init", 17)]) + + 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 # ICS volume (cm^3) + ECS_vol = 7.08e-11 # ECS volume (cm^3) + surface = 2.29e-6 # membrane surface (cm²) + + K_e = states[1] #11 + K_i = states[2] + Na_e = states[3] + Na_i = states[4] #12 + + K_e_init = parameters[16] + K_i_init = parameters[17] + + E_Na = R * temperature / F * np.log(Na_e/Na_i) #4 + E_K = R * temperature / F * np.log(K_e/K_i) #5 + E_K_init = R * temperature / F * np.log(K_e_init/K_i_init) #5 + + i_pump = parameters[15] / ((1 + parameters[13] / K_e) ** 2 \ + * (1 + parameters[14] / Na_i) ** 3) + + # set conductance + dphi = states[0] - E_K + 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[0])/0.0441e3) # shorthand + g_Kir = np.sqrt(K_e/K_e_init)*(A*B)/(C*D) + + # define and return current + I_Kir = parameters[3]*g_Kir*(states[0] - E_K) # umol/(cm^2*ms) + + # Expressions for the Sodium channel component + i_Na = parameters[2] * (states[0] - E_Na) + 3 * i_pump + + # Expressions for the Potassium channel component + i_K = I_Kir - 2 * i_pump + + # set I_ch_Na + parameters[8] = i_Na + # set I_ch_K + parameters[9] = i_K + # set I_ch_Cl + parameters[10] = 0.0 + + # Expression for phi_M + values[0] = (- i_K - i_Na)/parameters[6] + + # Expression for K_e + values[1] = i_K * surface / (F * ECS_vol) + + # Expression for K_i + values[2] = - i_K * surface / (F * ICS_vol) + + # Expression for Na_e + values[3] = i_Na * surface / (F * ECS_vol) + + # Expression for Na_i + values[4] = - i_Na * surface / (F * ICS_vol) diff --git a/emix-simulations/mm_hh.py b/emix-simulations/mm_hh.py new file mode 100644 index 0000000..2bb5b65 --- /dev/null +++ b/emix-simulations/mm_hh.py @@ -0,0 +1,174 @@ +# 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 + n_init = 0.19059822295360918 + m_init = 0.01708076504334425 + h_init = 0.8504559822908828 + + phi_M_init = -74.18510727096373 + + init_values = np.array([m_init, h_init, n_init, phi_M_init], dtype=np.float_) + + # State indices and limit checker + state_ind = dict([("m", 0), ("h", 1), ("n", 2), ("V", 3)]) + + 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 = 0.1 # Na leak conductivity (mS/cm**2) + g_leak_K = 0.4 # K leak conductivity (mS/cm**2) + + m_K = 2 # threshold ECS K (mol/m^3) + m_Na = 7.7 # threshold ICS Na (mol/m^3) + I_max = 44.9 # max pump strength (muA/cm^2) + + # Set initial parameter values + init_values = np.array([g_Na_bar, g_K_bar, \ + g_leak_Na, g_leak_K, \ + 0, 0, 0, 0, \ + 0, 0, 0, 0, 0, + m_K, m_Na, I_max], dtype=np.float_) + + # Parameter indices and limit checker + param_ind = dict([("g_Na_bar", 0), ("g_K_bar", 1), \ + ("g_leak_Na", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), \ + ("Cm", 6), ("stim_amplitude", 7), \ + ("I_ch_Na", 8), ("I_ch_K", 9), ("I_ch_Cl", 10), \ + ("K_e", 11), ("Na_i", 12), \ + ("m_K", 13), ("m_Na", 14), ("I_max", 15)]) + + 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_inds = dict([("m", 0), ("h", 1), ("n", 2), ("V", 3)]) + + 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", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), ("Cm", 6), \ + ("stim_amplitude", 7), ("I_ch_Na", 8), \ + ("I_ch_K", 9), ("I_ch_Cl", 10), ("K_e", 11), \ + ("Na_i", 12), ("m_K", 13), ("m_Na", 14), ("I_max", 15)]) + + 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,) + + 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 = parameters[15] / ((1 + parameters[13] / parameters[11]) ** 2 \ + * (1 + parameters[14] / parameters[12]) ** 3) + + # Expressions for the Sodium channel component + i_Na = (parameters[2] + parameters[0]*states[1]*math.pow(states[0], 3) + i_Stim) * \ + (states[3] - parameters[4]) + 3 * i_pump + + # Expressions for the Potassium channel component + i_K = (parameters[3] + parameters[1]*math.pow(states[2], 4)) * \ + (states[3] - parameters[5]) - 2 * i_pump + + # set I_ch_Na + parameters[8] = i_Na + # set I_ch_K + parameters[9] = i_K + # set I_ch_Cl + parameters[10] = 0.0 + + values[3] = (- i_K - i_Na)/parameters[6] diff --git a/emix-simulations/mm_hh_calibration_ODE.py b/emix-simulations/mm_hh_calibration_ODE.py new file mode 100644 index 0000000..6b85f2c --- /dev/null +++ b/emix-simulations/mm_hh_calibration_ODE.py @@ -0,0 +1,219 @@ +# 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 + #n_init = 0.2276174382146544 # gating variable n + #m_init = 0.02472911223827967 # gating variable m + #h_init = 0.7867106718874558 # gating variable h + + n_init = 0.19059822295360918 + m_init = 0.01708076504334425 + h_init = 0.8504559822908828 + + phi_M_init = -74.18510727096373 + + # Initial values + K_i_init = 124.25309301681166 + Na_i_init = 12.73795855170082 + + K_e_init = 3.369320617474813 + Na_e_init = 100.67749185800223 + + init_values = np.array([m_init, h_init, n_init, phi_M_init, \ + K_e_init, K_i_init, Na_e_init, Na_i_init], dtype=np.float_) + + # State indices and limit checker + state_ind = dict([("m", 0), ("h", 1), ("n", 2), ("V", 3), + ("K_e", 4), ("K_i", 5), ("Na_e", 6), ("Na_i", 7)]) + + 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 = 0.1 # Na leak conductivity (mS/cm**2) + g_leak_K = 0.4 # K leak conductivity (mS/cm**2) + + m_K = 2 # threshold ECS K (mol/m^3) + m_Na = 7.7 # threshold ICS Na (mol/m^3) + I_max = 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, g_leak_K, \ + 0, 0, C_M, 0, \ + 0, 0, 0, 0, 0, + m_K, m_Na, I_max], dtype=np.float_) + + # Parameter indices and limit checker + param_ind = dict([("g_Na_bar", 0), ("g_K_bar", 1), \ + ("g_leak_Na", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), \ + ("Cm", 6), ("stim_amplitude", 7), \ + ("I_ch_Na", 8), ("I_ch_K", 9), ("I_ch_Cl", 10), \ + ("K_e", 11), ("Na_i", 12), \ + ("m_K", 13), ("m_Na", 14), ("I_max", 15)]) + + 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_inds = dict([("m", 0), ("h", 1), ("n", 2), ("V", 3), + ("K_e", 4), ("K_i", 5), ("Na_e", 6), ("Na_i", 7)]) + + 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", 2), ("g_leak_K", 3), \ + ("E_Na", 4), ("E_K", 5), ("Cm", 6), \ + ("stim_amplitude", 7), ("I_ch_Na", 8), \ + ("I_ch_K", 9), ("I_ch_Cl", 10), ("K_e", 11), \ + ("Na_i", 12), ("m_K", 13), ("m_Na", 14), ("I_max", 15)]) + + 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 # ICS volume (cm^3) + ECS_vol = 7.08e-11 # ECS volume (cm^3) + surface = 2.29e-6 # membrane surface (cm²) + + K_e = states[4] #11 + K_i = states[5] + Na_e = states[6] + Na_i = states[7] #12 + E_Na = R * temperature / F * np.log(Na_e/Na_i) #4 + E_K = R * temperature / F * np.log(K_e/K_i) #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 = parameters[15] / ((1 + parameters[13] / K_e) ** 2 \ + * (1 + parameters[14] / Na_i) ** 3) + + # Expressions for the Sodium channel component + i_Na = (parameters[2] + parameters[0]*states[1]*math.pow(states[0], 3) + i_Stim) * \ + (states[3] - E_Na) + 3 * i_pump + + # Expressions for the Potassium channel component + i_K = (parameters[3] + parameters[1]*math.pow(states[2], 4)) * \ + (states[3] - E_K) - 2 * i_pump + + # set I_ch_Na + parameters[8] = i_Na + # set I_ch_K + parameters[9] = i_K + # set I_ch_Cl + parameters[10] = 0.0 + + # Expression for phi_M + values[3] = (- i_K - i_Na)/parameters[6] + + # Expression for K_e + values[4] = i_K * surface / (F * ECS_vol) + + # Expression for K_i + values[5] = - i_K * surface / (F * ICS_vol) + + # Expression for Na_e + values[6] = i_Na * surface / (F * ECS_vol) + + # Expression for Na_i + values[7] = - i_Na * surface / (F * ICS_vol) diff --git a/emix-simulations/run_calibration.py b/emix-simulations/run_calibration.py new file mode 100644 index 0000000..71a3773 --- /dev/null +++ b/emix-simulations/run_calibration.py @@ -0,0 +1,171 @@ +#!/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 + +# 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 = 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 + + # neuron + # Initial values + K_i_init = 124.25309301681166 + Na_i_init = 12.73795855170082 + K_e_init = 3.369320617474813 + Na_e_init = 100.67749185800223 + Cl_e_init = Na_e_init + K_e_init + Cl_i_init = Na_i_init + K_i_init + + phi_M_init = Constant(-74.18510727096373) + + """ + + # Glial + # Init values + K_i_init = 101.45093275680335 + Na_i_init = 13.701595802750164 + K_e_init = 3.369320617474813 + Na_e_init = 100.67749185800223 + Cl_e_init = Na_e_init + K_e_init + Cl_i_init = Na_i_init + K_i_init + + phi_M_init = Constant(-81.73542331243227) + """ + + # Set physical parameters + params = namedtuple('params', ('dt', 'n_steps_ODE', 'F', 'psi', \ + 'C_phi', 'C_M', 'R', 'temperature', 'phi_M_init'))(dt, \ + n_steps_ODE, F, psi, C_phi, C_M, R, temperature, phi_M_init) + + # diffusion coefficients for each sub-domain + D_Na_sub = {1:D_Na, 0:D_Na} + D_K_sub = {1:D_K, 0:D_K} + D_Cl_sub = {1:D_Cl, 0:D_Cl} + + # initial concentrations for each sub-domain + Na_init_sub = {1:Constant(Na_i_init), 0:Constant(Na_e_init)} + K_init_sub = {1:Constant(K_i_init), 0:Constant(K_e_init)} + Cl_init_sub = {1:Constant(Cl_i_init), 0:Constant(Cl_e_init)} + + # Create ions (channel conductivity is set below in the membrane model) + Na = {'c_init_sub':Na_init_sub, 'bdry': Constant((0, 0)), + 'z':1.0, 'name':'Na', 'D_sub':D_Na_sub} + K = {'c_init_sub':K_init_sub, 'bdry': Constant((0, 0)), + 'z':1.0, 'name':'K', 'D_sub':D_K_sub} + Cl = {'c_init_sub':Cl_init_sub, '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] < 20.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/" + mesh_prefix = "meshes/3D_v_II/" + 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 import main + from make_mesh_3D_v_II 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 = True + 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 = True + 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/" + + # Dictionary with membrane models (key is facet tag, value is ode model) + ode_models = {1: mm_hh} + #ode_models = {1: 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_calibration_ODE.py b/emix-simulations/run_calibration_ODE.py new file mode 100644 index 0000000..b12b84a --- /dev/null +++ b/emix-simulations/run_calibration_ODE.py @@ -0,0 +1,118 @@ +import os +import dolfin as df +import numpy as np +import matplotlib.pyplot as plt +import mm_hh_calibration_ODE as ode +#import mm_glial_calibration_ODE as ode +from knpemidg.membrane import MembraneModel +from collections import namedtuple + +mesh = df.UnitSquareMesh(2, 2) +V = df.FunctionSpace(mesh, 'Discontinuous Lagrange Trace', 0) + +facet_f = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0) +tag = 0 + +#stimulus = None + +g_syn_bar = 5 +#g_syn_bar = 0 +stimulus = {'stim_amplitude': g_syn_bar} + +membrane = MembraneModel(ode, facet_f=facet_f, tag=tag, V=V) + +V_index = ode.state_indices('V') +K_e_index = ode.state_indices('K_e') +K_i_index = ode.state_indices('K_i') +Na_e_index = ode.state_indices('Na_e') +Na_i_index = ode.state_indices('Na_i') + +#n_index = ode.state_indices('n') +#m_index = ode.state_indices('m') +#h_index = ode.state_indices('h') + +potential_history = [] +K_e_history = [] +K_i_history = [] +Na_e_history = [] +Na_i_history = [] + +#n_history = [] +#m_history = [] +#h_history = [] + +#for _ in range(50000): +for _ in range(500): +#for _ in range(50000): + membrane.step_lsoda(dt=0.1, stimulus=stimulus) + + potential_history.append(1*membrane.states[:, V_index]) + K_e_history.append(1*membrane.states[:, K_e_index]) + K_i_history.append(1*membrane.states[:, K_i_index]) + Na_e_history.append(1*membrane.states[:, Na_e_index]) + Na_i_history.append(1*membrane.states[:, Na_i_index]) + + #n_history.append(1*membrane.states[:, n_index]) + #m_history.append(1*membrane.states[:, m_index]) + #h_history.append(1*membrane.states[:, h_index]) + +potential_history = np.array(potential_history) +K_e_history = np.array(K_e_history) +K_i_history = np.array(K_i_history) +Na_e_history = np.array(Na_e_history) +Na_i_history = np.array(Na_i_history) + +#n_history = np.array(n_history) +#m_history = np.array(m_history) +#h_history = np.array(h_history) + +print("V", potential_history[-1, 2]) +print("K_e", K_e_history[-1, 2]) +print("K_i", K_i_history[-1, 2]) +print("Na_e", Na_e_history[-1, 2]) +print("Na_i", Na_i_history[-1, 2]) + +#print("n", n_history[-1, 2]) +#print("m", m_history[-1, 2]) +#print("h", h_history[-1, 2]) + +fig, ax = plt.subplots(5, 1, sharex=True) +ax[0].plot(potential_history[:, 2]) +ax[1].plot(K_e_history[:, 2]) +ax[2].plot(K_i_history[:, 2]) +ax[3].plot(Na_e_history[:, 2]) +ax[4].plot(Na_i_history[:, 2]) +ax[0].set_title("V") +ax[1].set_title("K_e") +ax[2].set_title("K_i") +ax[3].set_title("Na_e") +ax[4].set_title("Na_i") +plt.tight_layout() +fig.savefig("ode.png") +#plt.show() +plt.close() + +fig = plt.figure() +plt.plot(potential_history[:, 2]) +#plt.ylim([-72, -70]) +plt.tight_layout() +fig.savefig("phiM.png") + +fig = plt.figure() +plt.plot(K_e_history[:, 2]) +#plt.ylim([4, 4.1]) +plt.tight_layout() +fig.savefig("K_e.png") + +fig = plt.figure() +plt.plot(Na_e_history[:, 2]) +#plt.ylim([99, 100]) +plt.tight_layout() +fig.savefig("Na_e.png") + +# TODO: +# - consider a test where we have dy/dt = A(x)y with y(t=0) = y0 +# - after stepping u should be fine +# - add forcing: dy/dt = A(x)y + f(t) with y(t=0) = y0 +# - things are currently quite slow -> multiprocessing? +# - rely on cbc.beat?