New Approximate Matrix Exponential Solver Speeds Up Ion Channel Simulations in Brain Models

Simulating the human brain at scale has always forced computational neuroscientists into a painful trade-off: either use simplified neuron models that run fast but miss crucial biology, or use biologically realistic conductance-based models that crawl through simulation time. A new study published in Frontiers in Computational Neuroscience introduces a method that could dramatically ease this dilemma. The technique, called the approximate matrix exponential (AME) method, makes it possible to simulate biologically detailed Markov kinetic models of ion channels an order of magnitude faster and more accurately than established numerical approaches, achieving real-time performance on a state-of-the-art GPU for tens to hundreds of thousands of model instances simultaneously.

Ion channels are the protein macromolecules that conduct electrically charged ions across neuronal membranes, and their conductance is controlled by their conformational state. Markov models of molecular kinetics sit at the very foundation of neuroscience, describing the physical components that implement the brain’s computations. They were first applied to neuroscience to describe the ion channels responsible for generating action potentials, and because brains are composed of trillions of repeated components, computational neuroscience has what the study’s author calls a voracious appetite for these models. Yet despite their biological fidelity, Markov kinetic models have remained underutilized because they are notoriously difficult to simulate efficiently.

The mathematical core of the problem is the master equation, which describes how a probability distribution over the channel’s conformational states evolves over time. This initial value problem has an analytic solution expressed through the matrix exponential function, yielding what is known as the propagator matrix. The propagator matrix advances the state of the model over each time step, but computing it exactly is prohibitively expensive for large-scale simulations. Direct evaluation of the Taylor series expansion converges slowly and is numerically unstable, while better methods such as scaling and squaring with Padé approximants still involve matrix-matrix multiplications whose time complexity scales cubically with the number of model states.

Established numerical solvers fall into two camps, each with serious drawbacks. Explicit integration methods, such as the forward Euler and Runge–Kutta–Chebyshev approaches, evaluate the derivative and run the state forward in time. They are simple and efficient, with quadratic time complexity, but they can be numerically unstable, allowing small errors to accumulate and compound exponentially into sudden catastrophic failure. Implicit methods, including the backward Euler and Crank–Nicolson schemes, evaluate the derivative at the end of the time step and offer accurate, stable integration, but they require solving a system of equations via Gaussian elimination, giving them cubic time complexity. The NEURON simulator, the state-of-the-art tool for conductance-based neuronal modeling, uses implicit methods by default.

The AME method takes a fundamentally different approach. Instead of computing the propagator matrix from scratch at every time step, it constructs a piecewise polynomial approximation of the propagator matrix function. The input space is partitioned into a regular grid of non-overlapping regions, each containing a polynomial approximation fitted to sampled data using linear least squares regression. Within each region, the input vector is transformed into a normalized range, and because the function being approximated returns a full matrix, the method builds a separate approximation for every cell of the propagator matrix, with all approximations sharing the same partitions and polynomial form for computational efficiency.

What makes the method particularly elegant is how it measures and controls its own accuracy. The integration error is defined as the maximum absolute error in any component of the state vector after integrating over one unit of time, across all possible initial states and inputs. Exploiting the Chapman–Kolmogorov equation, the method squares propagator matrices to extend sample integration periods from one time step to a full unit of time. Because Markov models are linear and probability is a conserved quantity, the maximum error occurs when the entire probability mass sits in a single state, which means the error of the standard basis vectors bounds the error of all possible state vectors. A greedy algorithm then adjusts the number of input partitions until the residual integration error falls below a user-specified error target, and a hill-climbing algorithm searches the space of polynomial forms to maximize throughput.

The method also handles a subtlety specific to neurobiology: chemical concentrations such as neurotransmitter levels span many orders of magnitude, which is difficult to approximate with fixed-size piecewise regions. The AME method preconditions concentration inputs with a base-2 logarithmic transform and a small positive offset, determined experimentally by testing 50 offsets in a geometric progression and selecting the one that minimizes the root-mean-square residual error. The base-2 logarithm is preferred because efficient algorithms and even dedicated hardware instructions, such as the lg2 instruction in Nvidia PTX assembly, can compute it quickly.

To validate the approach, the author assembled a benchmark of 99 Markov models drawn from four independent sources, covering AMPA and NMDA receptors, all nine subtypes of voltage-gated sodium channels, and 88 models of voltage-gated potassium channels built from seven kinetic schemes. The models were obtained from ModelDB and GitHub, and some numerically unstable models were detected and rejected by the tooling. Accuracy comparisons using the analytic solution as ground truth showed that the backward Euler method is only first-order accurate, with error directly proportional to the time step, whereas the AME method controls its accuracy independently of the time step, achieving an error target of 10 to the power of minus 3 for all state and input vectors tested.

The practical consequences were demonstrated in simulations of action potentials, the electrical signals that broadcast information throughout the brain, where timing errors quickly cascade. A single-compartment neuron containing eight ion channels from the benchmark was stimulated to produce trains of action potentials. With the backward Euler method, action potentials were delayed in proportion to the time step, while simulations using the exact matrix exponential and the AME method delayed action potentials only in proportion to the square of the time step, limited by NEURON’s second-order cable solver. Varying the AME error target showed convergence by 10 to the power of minus 3, while the loosest target of 10 to the power of minus 1 caused simulations to fail to produce action potentials altogether. In speed benchmarks with 1,000 single-compartment neurons, the AME method consistently ranked fastest, followed by backward Euler and then the exact matrix exponential.

The accompanying open-source program, called ame and available under the MIT license on the Python Package Index, integrates the method into the NEURON simulator, and for comparison the exact analytic solution was also added to NEURON as a new solver option. On an NVIDIA B200 GPU, the method achieved integration speeds as low as 0.059 nanoseconds per step per instance for a four-state potassium channel model, corresponding to real-time simulation of roughly 420,000 instances. The method does have limitations: it is restricted to models with one or two input variables, its time complexity scales quadratically with the number of model states, and building the approximation can take up to an hour per model on state-of-the-art hardware. Nevertheless, because most voltage-gated ion channels genuinely have only one input and most phenomena of interest can be modeled with a dozen states or less, the method’s trade-offs are well suited to its purpose. By decoupling accuracy from the time step and delivering GPU-friendly, branch-free execution, the AME method paves the way for large-scale neuronal simulations to adopt far more realistic models of the ion channels that generate the brain’s electrical activity.

Subject of Research: A fast and accurate numerical method for simulating Markov kinetic models of ion channels in large-scale neuronal simulations

Article Title: Approximate matrix exponential: a fast and accurate solver for Markov kinetic models of ion channels

Article References: McDougall, D. A. (2026). Approximate matrix exponential: a fast and accurate solver for Markov kinetic models of ion channels. Frontiers in Computational Neuroscience, 20, Article 1901593. https://doi.org/10.3389/fncom.2026.1901593

Image Credits: AI Generated

DOI: 10.3389/fncom.2026.1901593

Keywords: ion channels, Markov kinetic models, matrix exponential, NEURON simulator, GPU computing, numerical integration, computational neuroscience, propagator matrix, action potentials, conductance-based models, piecewise polynomial approximation, neuronal simulation