Skip to content

Commit

Permalink
add code for ODE pre-calibration and full system (ODE/PDE) calibration
Browse files Browse the repository at this point in the history
  • Loading branch information
adajel committed Aug 23, 2024
1 parent 18f7832 commit 8af9d74
Show file tree
Hide file tree
Showing 4 changed files with 699 additions and 0 deletions.
122 changes: 122 additions & 0 deletions emix-simulations/make_mesh_3D_two_tags.py
View file
Original file line number Diff line number Diff line change
@@ -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
257 changes: 257 additions & 0 deletions emix-simulations/mm_hh_glial_calibration_ODE.py
View file
Original file line number Diff line number Diff line change
@@ -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)
Loading

0 comments on commit 8af9d74

Please sign in to comment.