A splitting, discontinuous Galerkin solver for the cell-by-cell electroneutral Nernst-Planck framework

TL;DR

The paper addresses the challenge of simulating ionic electrodiffusion in brain tissue with explicit cell-scale geometry by introducing a splitting-based, multiplier-free formulation of the KNP-EMI model and a discontinuous Galerkin discretization. By decoupling NP ion transport from EMI electrostatics and eliminating Lagrange multipliers, the authors obtain a single-dimensional, mass-conserving system that is amenable to standard solvers with scalable AMG preconditioning. The method achieves optimal spatial and temporal convergence, exhibits robustness to discretization and problem size, and demonstrates near-optimal parallel scalability across idealized 2D/3D and morphologically realistic brain geometries. These advances enable high-fidelity, large-scale simulations of brain ion homeostasis and activity-driven ionic shifts in complex geometries, with potential applications in neuroscience and energy storage modeling alike.

Abstract

Mathematical models for excitable tissue with explicit representation of individual cells are highly detailed and can, unlike classical homogenized models, represent complex cellular geometries and local membrane variations. However, these cell-based models are challenging to approximate numerically, partly due to their mixed-dimensional nature with unknowns both in the bulk and at the lower-dimensional cellular membranes. We here develop and evaluate a novel solution strategy for the cell-based KNP-EMI model describing ionic electrodiffusion in and between intra- and extracellular compartments with explicit representation of individual cells. The strategy is based on operator splitting, a multiplier-free formulation of the coupled dynamics across sub-regions, and a discontinuous Galerkin discretization. In addition to desirable theoretical properties, such as local mass conservation, the scheme is practical as it requires no specialized functionality in the finite element assembly and order optimal solvers for the resulting linear systems can be realized with black-box algebraic multigrid preconditioners. Numerical investigations show that the proposed solution strategy is accurate, robust with respect to discretization parameters, and that the parallel scalability of the solver is close to optimal - both for idealized and realistic two and three dimensional geometries.

Figures (7)

• Figure 1: Finite element discretization of the KNP-EMI model by continuous linear and discontinuous linear Lagrange elements. Setup with single intracellular domain $\Omega_i$ (in green) is considered. Triangulation conforms to the interface $\Gamma$. Mesh vertices are shown with black circle markers (A). Intra-/extracellular variables are represented in separate$H^1$-conforming FE spaces over different meshes of $\Omega_i$ and $\Omega_e$. The coupling requires specialized implementation as on $\Gamma$ the degrees of freedom (shown with red and orange markers) interact (B). With DG discretization the problem unknowns are represented in a single FE spaces posed over the global mesh. Coupling on $\Gamma$ requires no special treatment compared to the remaining facets. However, the scheme leads to more degrees of freedom (shown in blue markers, C).
• Figure 1: Robustness of the preconditioner with respect to mesh size for model B (idealized 2D axon). The figure displays the temporal evolution of the membrane potential at point $(25, 1)$$\mu$m (A), the total number of iterations over time during refinement in space for the full KNP-EMI system (B), and an illustration of the 2D domain with ICS in yellow and ECS in purple (C).
• Figure 2: Robustness of the preconditioner with respect to mesh size for model C (idealized 3D axons). The figure displays the temporal evolution of the membrane potential at point $(25, 0.3, 0.4)$$\mu$m (A), the total number of iterations over time during refinement in space for the full KNP-EMI system (B), and an illustration of the 3D domain with ICS in yellow and ECS in purple (C).
• Figure 3: Parallel scalability of assembly and solve times for the EMI and the KNP sub-problems, and for solving the ODEs. Strong scaling data for model C (idealized 3D geometry), with ideal scaling reported (A, log-log plot). Weak scaling data for model B (idealized 2D geometry), with both ideal and linear scaling reported (B, log-log plot). The timings represent total runtime and assembly over 10 timesteps.
• Figure 4: Diffusive and advective contributions of ion transport in model C (idealized 3D axons). The two upper panels display the temporal (A) and spatial (B) evolution of the membrane potential recorded at respectively the point $x=(100, 0.5, 0.5)$$\mu$m and at 20 ms (pink), 50 ms (orange) and 70 ms (blue). The lower panels display the ratio between diffusion and advection for transport of $\rm Na^+$ (C) and $\rm K^+$ (D) where solid lines and dashed lines are respectively ECS and ICS traces of the ratio along the $x$-axis. Recall that the ratio is positive where diffusion dominates, and negative where drift dominates.