From 05a77ff454269d82914eeb64728aedaced0e1bd4 Mon Sep 17 00:00:00 2001 From: adajel Date: Fri, 23 Aug 2024 15:07:55 +0200 Subject: [PATCH] add code for running simulation on emix mesh --- emix-simulations/make_figures.py | 725 +++++++++++++++++++ emix-simulations/run_EMIx_simulation.py | 230 ++++++ emix-simulations/run_calibration_two_tags.py | 3 +- 3 files changed, 956 insertions(+), 2 deletions(-) create mode 100644 emix-simulations/make_figures.py create mode 100644 emix-simulations/run_EMIx_simulation.py diff --git a/emix-simulations/make_figures.py b/emix-simulations/make_figures.py new file mode 100644 index 0000000..66f0945 --- /dev/null +++ b/emix-simulations/make_figures.py @@ -0,0 +1,725 @@ +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()[:] *= 1e-7 + 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 = 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()[:] *= 1e-7 + 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) + hdf5file.read(subdomains, '/subdomains') + hdf5file.read(surfaces, '/surfaces') + + x_min = x - 0.5e-5; x_max = x + 0.1e-5 + y_min = y - 0.5e-4; y_max = y + 0.1e-4 + z_min = z; z_max = z + 0.04e-5 + + # 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 or surfaces[facet] == 2): + print("mark 10:", 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_s, E_Na_s, E_K_s + +def plot_3D_concentration_neuron(res, T, dt, fname): + + time = np.arange(0, T-dt, dt) + + # new mesh + x_M = 4.674154869123957e-4 + y_M = 2.208589275930137e-4 + z_M = 1.5693481199926018e-4 + + # 0.05 um above axon A (ECS) + x_e = x_M; y_e = y_M + 0.1e-4; z_e = z_M + # mid point inside axon A (ICS) + x_i = x_M; y_i = y_M - 0.1e-4; 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.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.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.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.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.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.plot(phi_M_ODE, linewidth=3) + + print("membrane potential", phi_M[0], phi_M[-1]) + print("membrane potential", phi_M_ODE[0], phi_M_ODE[-1]) + + print("Na_e", Na_e[0], Na_e[-1]) + print("Na_n", Na_i[0], Na_i[-1]) + + print("K_e", K_e[0], K_e[-1]) + print("K_n", K_i[0], K_i[-1]) + + print("Cl_e", Cl_e[0], Cl_e[-1]) + print("Cl_n", Cl_i[0], Cl_i[-1]) + + # make pretty + ax.axis('off') + plt.tight_layout() + + # save figure to file + plt.savefig('results/figures/callibrate_cell_1.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_3D_concentration_glial(res, T, dt, fname): + + time = np.arange(0, T-dt, dt) + + x_M = 4.813723702751168e-4 + y_M = 0.4791020242777646e-4 + z_M = 0.6595906191667912e-4 + + # 0.05 um above axon A (ECS) + x_e = x_M; y_e = y_M; z_e = z_M + 0.1e-4 + # mid point inside axon A (ICS) + x_i = x_M; y_i = y_M; z_i = z_M - 0.1e-4 + + # 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.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.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.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.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.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.plot(phi_M_ODE, linewidth=3) + + print("membrane potential", phi_M[0], phi_M[-1]) + print("membrane potential", phi_M_ODE[0], phi_M_ODE[-1]) + + print("Na_e", Na_e[0], Na_e[-1]) + print("Na_n", Na_i[0], Na_i[-1]) + + print("K_e", K_e[0], K_e[-1]) + print("K_n", K_i[0], K_i[-1]) + + print("Cl_e", Cl_e[0], Cl_e[-1]) + print("Cl_n", Cl_i[0], Cl_i[-1]) + + # make pretty + ax.axis('off') + plt.tight_layout() + + # save figure to file + plt.savefig('results/figures/callibrate_cell_2.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()[:] *= 1e-7 + 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()[:] *= 1e-7 + 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()[:] *= 1e-7 + 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 = 2 + +fname = 'results/data/EMIx/results.h5' + +#plot_surface(fname, T, dt) +plot_3D_concentration_neuron(res_3D, T, dt, fname) +plot_3D_concentration_glial(res_3D, T, dt, fname) +#write_to_pvd(dt, T, fname) + +#get_velocity(fname, T, dt) +#plot_surface_time(fname, T, dt) diff --git a/emix-simulations/run_EMIx_simulation.py b/emix-simulations/run_EMIx_simulation.py new file mode 100644 index 0000000..38a9055 --- /dev/null +++ b/emix-simulations/run_EMIx_simulation.py @@ -0,0 +1,230 @@ +#!/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 + resolution = 0 + + # Time variables (PDEs) + dt = 0.1 # global time step (s) + Tstop = 2 # global end time (s) + t = Constant(0.0) # time constant + + # Time variables (ODEs) + n_steps_ODE = 25 # number of ODE steps + + # Physical parameters + C_M = 2.0 # capacitance + temperature = 300e3 # temperature (mK) + F = 96485e3 # Faraday's constant (mC/mol) + R = 8.314e3 # Gas Constant (mJ/(K*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_g_init), \ + 2:Constant(Na_n_init)} + K_init_sub = {0:Constant(K_e_init), 1:Constant(K_g_init), \ + 2:Constant(K_n_init)} + Cl_init_sub = {0:Constant(Cl_e_init), 1:Constant(Cl_g_init), \ + 2:Constant(Cl_n_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 # synaptic conductivity (mS/cm**2) + + # Set stimulus ODE + stimulus = {'stim_amplitude': g_syn_bar} + stimulus_locator = lambda x: (x[0] < 2.0e-6) + + stim_params = namedtuple('membrane_params', ('g_syn_bar', \ + 'stimulus', 'stimulus_locator'))(g_syn_bar, \ + stimulus, stimulus_locator) + + # Get mesh, subdomains, surfaces paths + mesh_path = 'meshes/emix_meshes/volume_ncells_5_size_5000/' + + mesh = Mesh() + infile = XDMFFile(mesh_path + 'mesh.xdmf') + infile.read(mesh) + cdim = mesh.topology().dim() + subdomains = MeshFunction("size_t", mesh, cdim) + infile.read(subdomains, "label") + infile.close() + + print(np.unique(subdomains.array())) + + # Remark subdomains + for cell in cells(mesh): + if subdomains[cell] == 1: + subdomains[cell] = 0 + elif subdomains[cell] == 2: + subdomains[cell] = 2 + elif subdomains[cell] == 3: + subdomains[cell] = 2 + elif subdomains[cell] == 4: + subdomains[cell] = 1 + elif subdomains[cell] == 5: + subdomains[cell] = 1 + elif subdomains[cell] == 6: + subdomains[cell] = 1 + + print(np.unique(subdomains.array())) + + # get all local labels + File("meshes/emix_meshes/subdomains.pvd") << subdomains + + infile = XDMFFile(mesh_path + 'tags.xdmf') + surfaces = MeshFunction("size_t", mesh, cdim - 1) + infile.read(surfaces, "boundaries") + infile.close() + + print(np.unique(surfaces.array())) + + unique, counts = np.unique(surfaces.array(), return_counts=True) + + #print(dict(zip(unique, counts))) + + # Remark facets + for facet in facets(mesh): + if surfaces[facet] > 10: + surfaces[facet] = 10 + elif surfaces[facet] == 2: + surfaces[facet] = 2 + elif surfaces[facet] == 3: + surfaces[facet] = 2 + elif surfaces[facet] == 4: + surfaces[facet] = 1 + elif surfaces[facet] == 5: + surfaces[facet] = 1 + elif surfaces[facet] == 6: + surfaces[facet] = 1 + + print(np.unique(surfaces.array())) + + File("meshes/emix_meshes/surfaces.pvd") << surfaces + + # convert mesh from nm to cm + mesh.coordinates()[:,:] *= 1e-7 + + """ + + # Get mesh, subdomains, surfaces paths + mesh_prefix = "meshes/3D/" + 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 + 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 = None + + # Set solver parameters KNP (True is direct, and False is iterate) + direct_knp = False + rtol_knp = 1E-7 + atol_knp = 1E-40 + threshold_knp = None + + # 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/EMIx/" + + # Dictionary with membrane models (key is facet tag, value is ode model) + ode_models = {1: mm_glial, 2: mm_hh} + + # 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_two_tags.py b/emix-simulations/run_calibration_two_tags.py index a286313..d8ee89f 100644 --- a/emix-simulations/run_calibration_two_tags.py +++ b/emix-simulations/run_calibration_two_tags.py @@ -12,7 +12,6 @@ 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: @@ -33,7 +32,7 @@ class bcolors: # Time variables PDEs dt = 0.1 # global time step (ms) - Tstop = 2 # global end time (ms) + Tstop = 10 # global end time (ms) t = Constant(0.0) # time constant # Time variables ODEs