Scalable approximation and solvers for ionic electrodiffusion in cellular geometries

TL;DR

The paper tackles scalable simulation of ionic electrodiffusion in highly detailed brain tissue by formulating the KNP-EMI model with explicit intracellular and extracellular compartments and membrane coupling. It introduces a monolithic, variational finite element discretization with an implicit-explicit time scheme and a block-diagonal AMG-preconditioned GMRES solver, enabling robust performance up to $10^8$ degrees of freedom on large HPC resources. Key contributions include a detailed solver strategy, a four-model benchmark suite spanning idealized and image-based geometries, and comprehensive numerical experiments demonstrating near-optimal strong and weak scalability. The work enables high-fidelity, large-scale in-silico investigations of ion dynamics in brain tissue and provides open-source tooling for future benchmarks and extensions.

Abstract

The activity and dynamics of excitable cells are fundamentally regulated and moderated by extracellular and intracellular ion concentrations and their electric potentials. The increasing availability of dense reconstructions of excitable tissue at extreme geometric detail pose a new and clear scientific computing challenge for computational modelling of ion dynamics and transport. In this paper, we design, develop and evaluate a scalable numerical algorithm for solving the time-dependent and nonlinear KNP-EMI equations describing ionic electrodiffusion for excitable cells with an explicit geometric representation of intracellular and extracellular compartments and interior interfaces. We also introduce and specify a set of model scenarios of increasing complexity suitable for benchmarking. Our solution strategy is based on an implicit-explicit discretization and linearization in time, a mixed finite element discretization of ion concentrations and electric potentials in intracellular and extracellular domains, and an algebraic multigrid-based, inexact block-diagonal preconditioner for GMRES. Numerical experiments with up to $10^8$ unknowns per time step and up to 256 cores demonstrate that this solution strategy is robust and scalable with respect to the problem size, time discretization and number of cores.

Figures (7)

• Figure 1: Model A. (a) Sample 2D geometry ($N_x = 4$). (b) Sample 3D geometry ($N_x = 20$). (c) Membrane current stimulus (acting on $\Gamma$) $g_{\rm stim}^{\rm Na}$ versus time $t$. (d) Evolution of ion concentrations over time. Extracellular and intracellular quantities are sampled respectively at $\bm{x}_e=(0.15,0.15)\,\mu\mathrm{m}\in\Omega_e$ and $\bm{x}_i=(0.5,0.5)\,\mu\mathrm{m}\in\Omega_i$. (f) Evolution of membrane potential, which is homogeneous in $\Gamma$, for different ionic models: the default Hodgkin-Huxley (HH); the Kir-Na/K model; a leak model obtained by setting $f\mathrm{_{Kir}}=j_\mathrm{pump}=0$ in the Kir-Na/K model.
• Figure 1: Number of iterations (left) and runtimes (right) over $N_t=300$ time steps for Model A with $d = 2$ and $(N_x, p) = (512, 1)$ ($N=1\,056\,772$), using 32 MPI processes.
• Figure 2: Model B. (a) Astrocyte geometry with corresponding tetrahedral mesh. (b) Ionic source terms over time, imposed for $x<0$. (c),(d) Intracellular potassium [K]$_i$ (mM) at two different time steps. Potassium is injected in the extracellular space for $x<0$. At $t=2$ ms the extracellular potassium [K]$_e$ is increased by approximately 10 mM w.r.t. the initial condition [K]$^0_e=4$ mM. As the system evolves, potassium is electro-diffusing in both $\Omega_i$ and $\Omega_e$ through $\Gamma$. (e)--(g) Time evolution of ionic concentrations $[k]_r(\bm{x}_\star,t)$ (mM) with $\bm{x}_\star\in\Gamma$ labeled by $\star$ in panel (c).
• Figure 2: Parallel performance: strong scaling for Models A, B, and C. Plots show runtimes for iterative solves, finite element assembly and total runtime versus number of MPI processes (cores). For all models, we use $p = 1$ and $N_t = 20$ time steps. Average AMG($P_0$)-GMRES iterations over $n=1, \ldots, N_t$ are indicated for each data point for Models A ($d=3$), B, and C.
• Figure 3: Model C. (a) Dendritic segment (red), spines heads $\Gamma_{\rm head}$ (green), and glial cells (blue). A portion of the meshed extracellular domain is shown as long as a particular of few meshed dendritic spines. (b),(c) [K]$_i$ (mM) for two different time steps, given a membrane stimulus located in the spines heads. (d)--(f) Concentrations $[k]_r(\bm{x}_\star,t)$ (mM) for $k\in\{\text{K$^+$, Na$^+$}\}$ and potential $\phi_M(\bm{x}_\star,t)$ (mV) with $\bm{x}_\star\in\Gamma_{\rm head}$. The location of $\bm{x}_\star\in\Gamma_{\rm head}$ is labeled in panels (b) and (c) with $\star$.

Theorems & Definitions (1)

• Remark 4.1