PaperlayerView exact revision on arXiv

Source-linked AI summary

Second-Order IMEX Time-Stepping Methods for Efficient Bidomain Multi-Electrode Array Simulations of Cardiac Stem Cell Monolayers

Sofia Tonali, Sofia Botti, Luca Franco Pavarino

arXiv:2609.12650v1math.NAcs.NA

TL;DR

Bidomain–MEA simulations must represent stiff ionic dynamics and the electrode–tissue interface, while higher-order IMEX methods in this detailed setting remain to be assessed. The paper evaluates second-order IMEX–RK schemes with Strang-based splitting and finds substantially lower errors than the first-order approach at comparable computational cost. The authors identify the lack of rigorous convergence and stability analysis for the fully coupled nonlinear MEA discretization as an open limitation.

  • Problem

    Detailed MEA simulations of hiPSC-CM monolayers require accurate and efficient treatment of stiff Bidomain dynamics and the electrode–tissue interface, but higher-order IMEX schemes in this setting remain underexplored.

    Source: §1 · p. 2

  • Method

    The study combines second-order Strang-based operator splitting with IMEX Runge–Kutta schemes and evaluates them in detailed MEA simulations.

    Source: §1 · p. 2, §3.1 · p. 5, §3.2 · p. 6

  • Results

    Up to 2–3 orders of magnitude accuracy gains were achieved by higher-order schemes over the first-order approach for a fixed time step size in most test cases.

    Source: §5 · p. 20

  • Takeaways & Limitations

    Second-order IMEX–RK schemes achieve substantially lower errors for comparable computational cost across the investigated 0D and 2D settings, variables, and error measures.

    Source: §5 · pp. 19–20

    Rigorous convergence, stability, and error estimates for the proposed IMEX operator-splitting discretization of the fully coupled nonlinear MEA model remain unavailable.

    Source: Limitations and Future Work · p. 20

Abstract

from arXiv · show

Multi-electrode Arrays (MEAs) enable tissue-level electrophysiological studies of human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) by recording extracellular field potentials. Recent computational models have improved MEA simulations by coupling detailed electrode descriptions with the Bidomain framework, a parabolic-elliptic system of nonlinear PDEs coupled with a stiff ionic model. The standard numerical strategy relies on operator splitting techniques that decouple the PDE and ODE components. In most implementations, the ionic subsystem is treated explicitly while the diffusive operator is handled implicitly, resulting in a first-order implicit-explicit (IMEX) time discretization. Although computationally convenient, this approach limits temporal accuracy and may reduce efficiency in large-scale simulations. In this work, we investigate higher-order IMEX Runge-Kutta schemes within a Strang-based operator splitting, specifically tailored for the MEA model of hiPSC-CMs monolayers. First- and second-order schemes are compared in terms of computational cost and global error against a high-fidelity reference solution obtained with a very small time step. For a fixed computational cost, the second-order schemes achieve errors that are 2-3 orders of magnitude smaller than those of the first-order method. These results demonstrate that higher-order IMEX integration significantly improves the accuracy-to-cost ratio of MEA simulations, providing a practical and reliable approach for large-scale electrophysiological studies.

1 Introduction

hiPSC-CM tissue simulations couple Bidomain reaction–diffusion equations with stiff ionic ODEs, while MEA recordings require models that capture both tissue dynamics and the electrode–tissue interface. This study tests whether second-order IMEX Runge–Kutta schemes within Strang-based splitting improve efficiency over the first-order approach for detailed MEA simulations.

1 Introduction

  • hiPSC-CMs express main cardiac ion channels and exhibit electrophysiological responses comparable to adult human cells, supporting their use in drug safety and disease studies.Applications include both in-vitro and in-silico studies.Source: §1 · pp. 1–2
  • Bidomain tissue simulations couple reaction–diffusion equations with stiff ionic-current ODEs, creating a multiscale discretization challenge.Numerical methods must balance stability, accuracy, and computational efficiency.Source: §1 · p. 2
  • Finite elements and operator splitting decouple nonlinear PDEs into simpler subproblems, with IMEX schemes treating diffusion implicitly and ionic dynamics explicitly.Prior cardiac studies have explored higher-order splitting and IMEX schemes for improved accuracy without substantial cost increases.Source: §1 · p. 2
  • MEA recordings of extracellular field potentials motivate models that capture both intracellular dynamics and the electrode–tissue interface.The Bidomain formulation is extended to incorporate a detailed electrode description.Source: §1 · p. 2
  • Second-order IMEX–RK methods with Strang-based splitting are evaluated against the first-order approach in detailed MEA simulations using the stiff Paci 2020 hiPSC-CM model.The model includes 23 ionic currents and 3 ionic concentrations.Source: §1 · p. 2
  • The paper first isolates temporal integration error in a zero-dimensional ionic model, then evaluates the schemes in the full two-dimensional MEA setting.The study also describes its reference solution and error-analysis framework.Source: §1 · p. 2

2 Mathematical Model

The model couples Bidomain tissue dynamics to recording-electrode behavior to represent a single-well MEA and its extracellular field potentials. It uses Paci 2020 ionic dynamics, simplified isotropic conductivities, and prescribed stimulation and boundary conditions.

2 Mathematical Model

  • The framework couples Bidomain equations to a recording-electrode model and represents a 60–6 well MEA with nine electrodes per well; simulations use a single well.The framework is described as applicable to other MEA types.Source: §2 · pp. 2–3

The MEA model

  • The two-dimensional Bidomain model describes intracellular and extracellular electrical activity in a hiPSC-CM monolayer connected to a MEA.The full coupled system is presented in parabolic-elliptic form.Source: The MEA model · p. 3

Initial and Boundary Conditions

  • The model initializes transmembrane potential, gating variables, and concentrations, with no-flux intracellular-potential boundaries and grounded or no-flux extracellular-potential boundaries.The extracellular-potential boundary is grounded on one portion and has no flux on the remaining portion.Source: Initial and Boundary Conditions · pp. 3–4

Ionic model: the Paci 2020 model

  • The Paci 2020 model represents ventricular-like hiPSC-CMs and spontaneous automaticity with a stiff, nonlinear ionic system.It contains 23 gating variables and three concentrations: intracellular sodium, intracellular calcium, and sarcoplasmic-reticulum calcium.Source: Ionic model: the Paci 2020 model · p. 4

Tissue conductivity

  • Because directional conductivities are difficult to characterize in hiPSC-CM monolayers, the study assumes an isotropic tissue medium.The difficulty is linked to variation in cell orientation, gap-junction distribution, and maturation state.Source: Tissue conductivity · p. 5

2.1 Applied Current

  • A 150 µA/cm2 stimulus is applied to a 3600 µm2 square at the domain’s bottom-left corner for 0.51 ms, with a smoothed transition to avoid sharp temporal discontinuities.The smoothing uses a hyperbolic tangent function.Source: §2.1 · p. 5

3 Numerical Methods

The numerical method combines bilinear finite elements for spatial discretization with a Strang-based splitting framework and IMEX Runge–Kutta time integration. Electrode potentials are computed from the extracellular potential over each electrode.

3.1 MEA spatial discretization

  • The MEA system is spatially discretized with piecewise bilinear finite elements on a uniform 15 µm mesh, yielding mass, stiffness, and electrode-specific local mass matrices.The matrices are constructed using a Galerkin procedure.Source: §3.1 · p. 5
  • At each time step, the average extracellular potential on each electrode is computed by interpolating nodal values and integrating over the electrode with a composite two-dimensional trapezoidal rule.For a uniform grid, the quadrature weights differ for interior, edge, and corner points.Source: §3.1 · pp. 5–6
  • Figure 2 presents transmembrane-potential snapshots at 20, 40, and 300 ms and electrode-9 traces of v, ue, and U9.The figure aligns the temporal traces with the spatial snapshots.Source: Figure 2 · p. 6, §3.1 · p. 6

3.2 MEA Operator Splitting Framework

The MEA time integrator uses second-order Strang-based operator splitting to separate reaction and diffusion. It advances the ionic system for half a time step, solves diffusion for a full step, then advances the ionic system for the remaining half step.

3.2 MEA Operator Splitting Framework

  • The second-order Strang-based splitting advances the ionic system over a half step, solves diffusion over the full step, then advances the ionic system over the second half step.The first ionic solve provides intermediate values for the diffusion solve, which produces an updated intermediate potential before the final ionic solve.Source: §3.2 · pp. 6–7

Step 4 — Extracellular and electrodes update

The MEA time integrator uses IMEX-RK for reaction-diffusion steps and Crank–Nicolson for the coupled electrode-extracellular update.

Step 4 — Extracellular and electrodes update

3.3 MEA time discretization

  • Reaction-diffusion steps use IMEX-RK, while the coupled electrode equations use Crank–Nicolson.Source: §3.3 · p. 7
  • The hybrid time integration is designed to preserve the splitting framework’s second-order accuracy.Source: §3.3 · p. 7

IMEX–RK time discretization

The reaction-diffusion system uses an IMEX Runge–Kutta scheme, with nonlinear ionic dynamics treated explicitly and linear diffusion implicitly. The coupled electrode-extracellular system is advanced with Crank–Nicolson.

IMEX–RK time discretization

  • Nonlinear ionic dynamics are integrated explicitly, while the linear diffusion operator is handled implicitly.Source: IMEX–RK time discretization · p. 7
  • An s-stage IMEX-RK method uses distinct explicit and implicit contributions, with stage solutions at t_n + c_i∆t and a final approximation at t_n+1.Source: IMEX–RK time discretization · p. 8
  • The method’s explicit and implicit coefficient matrices define the treatments of the two contributions, with internal nodes and weight vectors also specifying the scheme.Source: IMEX–RK time discretization · p. 8

Crank-Nicolson method for electrode dynamics

3.4 Error computation and order of convergence

The convergence study assesses the proposed framework using a reference solution because analytical solutions for the coupled MEA system are unavailable. The reference is computed with a fine time step to make its temporal discretization error negligible relative to measured errors.

3.4 Error computation and order of convergence

  • The convergence study evaluates the numerical implementation using a reference solution, error metrics, and experimental order of convergence.Source: §3.4 · p. 9

Reference solution

  • Because analytical solutions for the coupled MEA system are unavailable, the reference solution is computed using the second-order IMEX-SSP2(2,2,2) scheme.Source: Reference solution · p. 10
  • A very fine time step makes the reference solution’s temporal discretization error negligible relative to errors measured in the convergence study.Source: Reference solution · p. 10

Accuracy evaluation

Accuracy is assessed with a global mixed root mean square error and discrete L2 errors evaluated at selected spatial nodes or time instants.

Accuracy evaluation

  • The mixed root mean square error is used as a global error measure for the MEA quantities over the fully discrete space-time set.Source: Accuracy evaluation · p. 10
  • The test and reference solutions’ quantities include v, u_e, and electrode field potentials U^k_FP at each space-time point.Source: Accuracy evaluation · p. 10
  • Discrete temporal L2 error is evaluated at a selected spatial node, while spatial L2 error is evaluated at a selected time instant.Source: Accuracy evaluation · p. 10

Order of Convergence

The framework’s overall convergence order is limited by the lower of the IMEX integrator and operator-splitting orders. The study estimates experimental order by refining the time step and comparing errors across levels.

Order of Convergence

  • The overall temporal order is min(pIMEX, psplit); with second-order Strang splitting, pIMEX ≥ 2 gives expected second-order accuracy.The experimental order is estimated by repeatedly halving the time step and comparing consecutive-level errors.Source: Order of Convergence · pp. 10–11

Implementation details

4 Numerical Results and accuracy study

The 0D and 2D studies assess convergence and computational efficiency for the proposed schemes. In the 2D MEA model, second-order schemes reach target accuracy faster despite higher per-step cost, while profiling identifies the ionic stages as the main computational burden.

4.1 0D Paci 2020 Ionic Model

  • The 0D tests cover 5 s of spontaneous activity and 1 Hz pacing; H(2,2,2) and SSP2(2,2,2) share Heun’s explicit tableau and convergence behavior.The paced stimulus is represented as an instantaneous membrane-potential depolarization.Source: §4.1 · pp. 11–12, Table 2 · p. 11

Error and Order of Convergence

  • Across paced and unpaced 0D runs, errors decrease as ∆t is reduced, with the tested methods exhibiting their expected theoretical convergence orders.Temporal L2 errors are assessed for v and Cai; Figure 5 shows convergence and efficiency measures for paced runs, with unpaced results reported as comparable.Source: Error and Order of Convergence · p. 12, Table 3 · p. 13, Figure 5 · p. 14

Comparative Performance Analysis

4.2 2D MEA results

  • The 2D study considers paced simulations only and limits analysis to a single paced action potential of approximately 600 ms.Unpaced cells would rapidly synchronize and depolarize nearly simultaneously, preventing propagation wavefronts.Source: §4.2 · p. 12

APs and space–time propagation

  • The 2D model reproduces a depolarization wavefront across tissue and corresponding extracellular signals at the MEA electrodes.Figure 2 links snapshots of transmembrane potential v to time traces of v, ue, and U 9 FP at the ninth electrode.Source: APs and space–time propagation · p. 13

Error and order of convergence

Comparative performance analysis

OS CPU time analysis

  • Second-order schemes approximately double total CPU time; ionic-model Steps 1 and 3 dominate, with an initial workload spike during rapid depolarization.After the transient, cost stabilizes during plateau and repolarization; SP(1,1,1) has the lowest per-step cost, while the second-order methods have nearly identical profiles.Source: OS CPU time analysis · pp. 18–19, Figure 7 · p. 18, Figure 8 · p. 19

5 Discussion and conclusions

The discussion finds that the Bidomain–electrode framework reproduces expected field-potential morphology, while numerical tests support expected convergence orders and an efficiency advantage for second-order IMEX schemes. The authors note that tested time orders stop at second order and rigorous analysis of the fully coupled nonlinear model remains open.

5 Discussion and conclusions

  • Simulated field-potential traces display expected morphology, supporting a consistent description of MEA dynamics through Bidomain–electrode coupling.The framework models hiPSC-CM monolayers recorded through MEA devices.Source: §5 · p. 19
  • Convergence analysis in the 0D ionic model and full 2D MEA configuration confirmed the theoretical orders across quantities and error metrics.The investigated quantities were v, ue, and UFP.Source: §5 · p. 19
  • Second-order schemes reach a given accuracy with coarser time steps, reducing overall computational effort despite higher per-step cost.The first-order scheme is less expensive per time step but requires smaller steps for the same error tolerance.Source: §5 · p. 19
  • Second-order schemes achieved accuracy gains of up to 2–3 orders of magnitude over first order at fixed time steps in most test cases.The advantage was observed in both 0D and 2D settings and remained robust across variables and error measures.Source: §5 · pp. 19–20

Limitations and Future Work

  • The analysis covers only first- and second-order time discretizations within the adopted Strang-based operator-splitting framework.Third- or higher-order schemes are identified as a possible extension.Source: Limitations and Future Work · p. 20
  • Rigorous convergence, stability, and error estimates for the fully coupled nonlinear MEA discretization remain open mathematical questions.The numerical experiments confirmed expected convergence orders, but a rigorous theoretical analysis is still missing.Source: Limitations and Future Work · p. 20
Loading 2609.12650v1…