Solving a multitraces formulation problem coupled with ordinary differential equations

Code considerations:

• We try to use numpy arrays for faster computations whenever it is possible (lists are very slow).

Import of the numpy library:

import numpy as np

Here we solve an example based in the problem solved in e2.mtf_time_coupled_Kavian.ipynb but the with FitzHugh-Nagumo and different parameters.

One sphere

We assume that the center of the sphere is in the origin.

Parameters set up.

First, we set up the parameters of the problem:

Geometry and physical parameters

Radius of the sphere. It has to be a float number:

radius = 1.

The variable sigmas is a numpy array with the medium conductivities in the order \(\sigma_0\), \(\sigma_1\). The numbers have to be floats:

sigma_e = 15. # extracellular conductivity
sigma_i =  1.5  # intracellular conductivity
sigmas = np.asarray([sigma_e, sigma_i])

Capacitance of the membrane (needs to be a float):

c_m = 1. # capacitance

Discretization parameters

Maximum degree of the spherical harmonics to be used for discretize (\(L\)). It has to be an integer:

big_l = 17

Maximum degree of the spherical harmonics to be used for quadratures (\(L_c\)). It has to be an integer:

big_l_c = big_l * 3

Parameters for the time scheme used:

initial_time = 0.
final_time = 10.
number_steps = 2**7

Initial conditions

Set up of initial conditions, for this model we need two arrays of length (big_l + 1)**2. Each array has the coefficients of the spherical harmonic expansion of each initial condition. In the following we give a couple of examples for building one of the arrays:

# Initial condition equal to zero
v_0_zero = np.zeros(((big_l+1)**2))
# Intial condition equal to a cte value
cte = 4.
v_0_cte_value = np.zeros(((big_l+1)**2))
v_0_cte_value[0] = cte * 2* np.sqrt(np.pi)

In this example we use initial conditions equal to zero. We need to set up both arrays in one:

initial_conditions = np.zeros((2, (big_l+1)**2))

Set up of the external excitation

Set up the contributions to the right hand side made by \(\phi_e\). For this, we set up a function only dependent on time, to account for changes in time. In the module biosspheres.formulations.mtf.timecouplings.righthands there are implementations on how to do this. We show some of the predetermined, and how to make a new one:

Pre-defined:

# Import of the module where pre-defined right hand sides are.
import biosspheres.formulations.mtf.timecouplings.righthands as righthands
# phi_e cte in time and space
cte=5.
phi_cte_space_time = righthands.phi_part_of_b_cte_space_and_time(
    big_l, n=1, radii=np.asarray([radius]), cte=cte
)
# phi_e point source in space, cte in time
phi_point_source_space_cte_time = righthands.phi_part_of_b_point_source_space_and_cte_time(
  big_l, n=1, radii=np.asarray([radius]), center_positions=[np.zeros((3))], p0 = 25.*np.ones((3)),
  sigma_e = sigmas[0], pii = np.asarray([sigmas[1]/sigmas[0]]), amplitude = 1.
)
# phi_e = cte * z_coordinate
phi_point_source_linear = righthands.phi_part_of_b_linear_space_z_cte_time(
  big_l, n=1, radii=np.asarray([radius]), center_positions=[np.zeros((3))],cte=cte, pii = np.asarray([sigmas[1]/sigmas[0]]))

/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Custom made:

# Space part, cte * z_coordinate
import biosspheres.formulations.mtf.mtf as mtf
import biosspheres.formulations.mtf.righthands as mtfrighthand

# Parameters for the function
n=1
radii=np.asarray([radius])
cte = 5.
ps =[np.zeros((3))]
pii = np.asarray([sigmas[1]/sigmas[0]])
x_dia, x_dia_inv = mtf.x_diagonal_with_its_inv(
        n, big_l, radii, pii, azimuthal=False
    )
# Pre-defined space part
space_phi = mtfrighthand.b_vector_n_spheres_mtf_linear_function_z(
    big_l, n, ps, cte, radii, x_dia)

# Time part, step function, custom made
def time_function_step(time):
  if time >= 5. or time < 0:
    return 0.
  else:
    return 1.

# Time part, sine function, custom made
parameter_1 = 1.
parameter_2 = 0.
def time_function_sine(time):
  return np.sin(time*parameter_1 + parameter_2)

custom_phi_step = righthands.phi_part_of_b_separable_in_space_time(space_phi, time_function_step)

custom_phi_sine = righthands.phi_part_of_b_separable_in_space_time(space_phi, time_function_sine)


# Which we are using:
b_phi_part_time_function = custom_phi_sine

Parameters for the coupled currents and ordinary differential equations:

For a non-linear current following FitzHugh-Nagumo dynamics:

parameter_theta = 0.8 # must be a float, adimensional parameter
parameter_b = 0.7     # must be a float, adimensional parameter
parameter_a = 0.3     # must be a float, adimensional parameter

Set up of internals

Set up of the time steps:

import biosspheres.formulations.mtf.timecouplings.solvertemplates as solve
tau, times, medium_times = solve.tau_times_medium(
    initial_time, final_time, number_steps
)

Set up of the coupled problem

Set up of the non linear current. Different models for the current can be used, there are models ready to use in the module biosspheres.miscella.forcouplings.currents. See the internals of the module for making custom made currents.

import biosspheres.miscella.forcouplings.currents as currents

FitzHugh-Nagumo current:

i_current_fitzhughnagumo = currents.i_fitzhughnagumo_one_sphere_function_creation(
    big_l,
  big_l_c,
  radius
)

Set up of the semi-implicit scheme for the “gating” variables (or similar)

import biosspheres.miscella.forcouplings.oderesolutions as oderesolutions
ode_next_sol_fitzhugh = (
    oderesolutions.semi_implicit_fitzhughnagumo_one_sphere_next_step(
            tau, parameter_theta, parameter_b, parameter_a
    )
)

Model to be used:

i_current = i_current_fitzhughnagumo
ode_next_sol = ode_next_sol_fitzhugh

Solving

First, we set up internal parameters used for requirements of numpy for threading:

import os
number = "7" # Maximum of threads
os.environ["OMP_NUM_THREADS"] = number
os.environ["OPENBLAS_NUM_THREADS"] = number
os.environ["MKL_NUM_THREADS"] = number
os.environ["VECLIB_MAXIMUM_THREADS"] = number
os.environ["NUMEXPR_NUM_THREADS"] = number

We implemented a direct solver (LU decomposition) that takes as parameters what we previously defined, including the current and the coupled ODE system. The complete list can be seen in the following:

solutions = solve.mtf_time_coupling_one_sphere(
    big_l,
    radius,
    sigmas,
    c_m,
    initial_time,
    final_time,
    number_steps,
    initial_conditions,
    b_phi_part_time_function,
    i_current,
    ode_next_sol,
)

Notice that the functions is defined in the module biosspheres.formulations.mtf.timecouplings.solvertemplates shortened in solve

Saving the results

The following saves all solutions, i.e. traces, transmembrane potential, and recovery variables, in numpy format. We have commented the line for the saving.

filename="example.npy"
directory=""
# np.save(directory+filename, solutions)

To load the solutions:

# solutions = np.load(directory+filename)

To save in txt format:

filename="example.txt"
directory=""
# np.savetxt(directory+filename, solutions)

To load the solutions:

# solutions = np.loadtxt(directory+filename)

Plotting the results

For plotting the evolution of time of the solutions evaluated in the north pole of the sphere:

import matplotlib.pyplot as plt
# Auxiliary arrays for plotting
import biosspheres.utils.auxindexes as auxindexes
num = (big_l + 1)**2
eles, l_square_plus_l, el_square_plus_el_divided_by_two = (
        auxindexes.eles_combination(big_l)
    )
little_partial_spherical = np.sqrt((2 * eles + 1))
plt.figure()
plt.plot(
    times,
    np.sum(
            solutions[:, (4 * num + l_square_plus_l[0 : big_l + 1])]
            * little_partial_spherical[0 : big_l + 1],
            axis=1,
    )
    / (2.0 * np.sqrt(np.pi)),
    label="Numerical solution",
    marker=".",
)
plt.title("v")

Text(0.5, 1.0, 'v')

plt.figure()
plt.plot(
    times,
    np.sum(
            solutions[:, (5 * num + l_square_plus_l[0 : big_l + 1])]
            * little_partial_spherical[0 : big_l + 1],
            axis=1,
    )
    / (2.0 * np.sqrt(np.pi)),
    label="Numerical solution",
    marker=".",
)
plt.title("Additional variable")

Text(0.5, 1.0, 'Additional variable')

Looking the value of a specific coefficient:

num = (big_l + 1)**2
eles, l_square_plus_l, el_square_plus_el_divided_by_two = (
        auxindexes.eles_combination(big_l)
    )
plt.figure()
l = 10 # Coefficient
plt.plot(
    times,solutions[:, (4 * num + l_square_plus_l[l])],
    label="Numerical solution",
    marker=".",
)
plt.title("Coeficiente l = "+ str(l) + ", m = 0")

Text(0.5, 1.0, 'Coeficiente l = 10, m = 0')

For plotting the transmembrane potential on the surface of the sphere at a given time:

# Auxiliary arrays
import biosspheres.quadratures.sphere as quadratures
(
    quantity_theta_points,
    quantity_phi_points,
    weights,
    pre_vector,
    spherical_harmonics,
) = quadratures.gauss_legendre_trapezoidal_real_sh_mapping_2d(
    big_l, big_l_c
)

time_index = number_steps // 2  # half of the simulation
vector = pre_vector * radius
surface_v = np.zeros((quantity_theta_points, quantity_phi_points))
surface_v[:, :] = np.sum(
    solutions[time_index, 4 * num:5*num, np.newaxis, np.newaxis]
    * spherical_harmonics, axis=0)
surface_v_max = np.max(np.abs(surface_v))
surface_v = (surface_v + surface_v_max) / (2. * surface_v_max)

from matplotlib import cm, colors
fig = plt.figure()
ax_1 = fig.add_subplot(111, projection='3d')
ax_1.plot_surface(
    vector[0, :, :],
    vector[1, :, :],
    vector[2, :, :],
    rstride=1, cstride=1,
    facecolors=cm.coolwarm(surface_v[:, :]))
ax_1.set_xlabel('$x$')
ax_1.set_ylabel('$y$')
ax_1.set_zlabel('$z$')
ax_1.set_aspect('equal')
fig.colorbar(cm.ScalarMappable(
    norm=colors.CenteredNorm(halfrange=surface_v_max), cmap=cm.coolwarm),
    ax=ax_1)

<matplotlib.colorbar.Colorbar at 0x770cc8514a10>

That can be made a video too.

For plotting a plane cut without \(\phi_e\):

import biosspheres.laplace.drawing as draw
from matplotlib import cm, colors
cut = 1
# Size of the window to be plotted
horizontal = 5.
vertical = 5.
# For the resolution of the plot
inter_horizontal = 100
inter_vertical = 100
# For plotting without the external function we make the following function
def zero(x: np.ndarray) -> float:
    return 0.

x1, y1, data = draw.draw_cut_representation_formula_one_sphere_azimuthal_symmetry(
    cut, np.zeros((3)), horizontal, vertical, inter_horizontal, inter_vertical, solutions[time_index, 0 : 4 * num], radius, big_l, zero)
plt.figure()
plt.imshow(data, origin='lower',
           extent=[-horizontal / 2, horizontal / 2,
                   -vertical / 2, vertical / 2],
           norm=colors.CenteredNorm()
           )
plt.xlabel('$y$')
plt.ylabel('$z$')
plt.title('$x = 0.$')
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x770cc8567350>