Axon Neural Simulation Framework

• Axon Neural Simulation Framework is a modular simulation environment that integrates electrodiffusion, membrane kinetics, and axon guidance mechanisms to capture realistic axonal behavior.
• It employs scalable high-performance computing techniques including finite element methods and Newton solvers to resolve detailed biophysical processes across multiple scales.
• The framework supports diverse models such as growth-cone dynamics, stochastic CTRW guidance, and mechanotactic durotaxis, enabling comprehensive studies of neural pathfinding and signal propagation.

The Axon Neural Simulation Framework denotes a reusable and extensible computational environment for simulating axonal dynamics across multiple scales, from electrodiffusion and membrane kinetics to growth-cone guidance, interstitial branching, stochastic trajectory formation, and substrate-mechanics-dependent steering. In the electrodiffusion formulation implemented as the dune-ax1 module atop DUNE, PDELab, ISTL, and MultidomainGrid/Multidomain, the framework couples Poisson–Nernst–Planck (PNP) electrodiffusion and Hodgkin–Huxley (HH) channel kinetics across a realistic axon geometry, resolves Debye layers and diffusion layers, and scales in parallel (Pods, 2019). Related modeling lines add activity-independent axon guidance with receptor regulation (Costa, 2015), interstitial branching controlled by axon guidance molecules (Suleymanov et al., 2013), GPU cue-diffusion guidance in 3D (Mihut et al., 2014), continuous-time random walk descriptions of growth-cone trajectories (Zhizhina et al., 2015), and a three-scale theory of axonal durotaxis (Kassianides et al., 3 Jan 2025). This suggests a framework concept defined less by a single codebase than by a common modular organization of geometry, fields, decision rules, constitutive laws, and numerical solvers.

1. Constituent model families and scope

The framework concept spans several distinct but partially compatible model classes. One line treats the axon as an electrodiffusive cable embedded in extracellular fluid and computes intracellular action potentials, extracellular action potentials, and near-/far-field signals by solving the full PNP system rather than assuming electroneutrality (Pods, 2019). A second line models axon guidance (AG) by representing neurons with soma, axon and growth cone, glial cells acting as guideposts, attraction/repulsion guidance cues, growth cone adaptation, tissue-gradient intersections, axonal transport, changes in the growth cone complexity and a range of responses for each receptor (Costa, 2015). Further lines address branching from the axon shaft (Suleymanov et al., 2013), 3D cue-diffusion guidance on massively parallel processors (Mihut et al., 2014), stochastic growth-cone motion via continuous time random walk (CTRW) (Zhizhina et al., 2015), and durotaxis through a three-scale model that links molecular clutch mechanics to growth-cone traction and axon-scale morphoelasticity (Kassianides et al., 3 Jan 2025).

Model family	Paper	Core content
Electrodiffusion framework	(Pods, 2019)	PNP equations, HH membrane coupling, DUNE implementation
Computational axon guidance	(Costa, 2015)	Receptors, ligands, regulatory network, midline crossing
Interstitial branching	(Suleymanov et al., 2013)	AGM-driven growth, pause-triggered active points, branch initiation
GPU cue diffusion guidance	(Mihut et al., 2014)	Particle-based diffusion surrogate, 3D steering, collision-based synapse formation
CTRW growth model	(Zhizhina et al., 2015)	Waiting-time laws, diffusive and subdiffusive regimes
Multiscale durotaxis	(Kassianides et al., 3 Jan 2025)	Molecular clutch, growth-cone traction, morphoelastic filament

A plausible implication is that the term Axon Neural Simulation Framework is most precise when used for a modular stack in which biophysical transport, extracellular fields, growth-cone sensing, and morphology updates are implemented as separable services rather than as a single monolithic model.

2. Electrodiffusion formulation and membrane biophysics

In the PNP formulation, the primary variables are the electric potential $\phi(\mathbf{x},t)$, ionic concentrations $c_i(\mathbf{x},t)$ for $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$, and the charge density $\rho = F \sum_i z_i c_i$. The governing equations are the Poisson equation

$$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$

and the Nernst–Planck transport equations

$$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$

The model is isothermal, includes no fluid convection and no bulk reactions, and is solved in both intra- and extracellular electrolytes without assuming electroneutrality (Pods, 2019).

Membrane coupling is imposed through HH-type current balance. With transmembrane potential $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$, specific capacitance $C_m$, and injected stimulus current $I_{\mathrm{stim}}$, the membrane relation is

$$C_m \frac{\partial V_m}{\partial t} + \sum_i I_i = I_{\mathrm{stim}},$$

with ionic current

$c_i(\mathbf{x},t)$0

and Nernst reversal potential

$c_i(\mathbf{x},t)$1

The flux boundary condition on the membrane interface maps current to ionic flux,

$c_i(\mathbf{x},t)$2

The displacement current is included through

$c_i(\mathbf{x},t)$3

with parallel-plate approximation

$c_i(\mathbf{x},t)$4

A defining feature of the framework is the explicit treatment of screening. The Debye length

$c_i(\mathbf{x},t)$5

sets the screening scale, and the model reports a Debye-layer coupling relation

$c_i(\mathbf{x},t)$6

with $c_i(\mathbf{x},t)$7 for the parameters used (Pods, 2019). This is the basis for the reported Debye-layer “AP echo”, in which the extracellular potential mirrors the intracellular action potential with reduced amplitude.

3. Geometry, discretization, and high-performance implementation

The baseline geometry is a straight cylindrical axon embedded in extracellular fluid. By rotational symmetry, the 3D problem is reduced to a 2D axisymmetric domain with axial coordinate $c_i(\mathbf{x},t)$8 and radial coordinate $c_i(\mathbf{x},t)$9. For the unmyelinated axon baseline, the computational domain is $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$0, the axon radius is $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$1, the membrane thickness is $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$2, and the bulk Debye length is approximately $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$3 (Pods, 2019). The extracellular upper boundary uses Dirichlet conditions $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$4 and $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$5, while other external boundaries use Neumann conditions.

The ionic and dielectric parameters are specified explicitly. The concentrations are

$$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$6

$$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$7

with diffusion coefficients

$$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$8

The framework uses continuous Galerkin FEM with bilinear quadrilateral elements $$i \in \{ \mathrm{Na^+}, \mathrm{K^+}, \mathrm{Cl^-} \}$$9 on a tensor grid. Radial refinement reaches $\rho = F \sum_i z_i c_i$0 to resolve $\rho = F \sum_i z_i c_i$1, while axial spacing is coarse away from active regions and graded near nodes in the myelinated case (Pods, 2019).

Time integration is implicit Euler (first order) with adaptive step size. The reported bounds are $\rho = F \sum_i z_i c_i$2 and $\rho = F \sum_i z_i c_i$3, capped at $\rho = F \sum_i z_i c_i$4 during the action potential. The nonlinear system is solved by a fully coupled Newton method with line search; HH gates are integrated implicitly per time step and held fixed during Newton linearization. Linear algebra uses ISTL BiCGStab with ILU(0) for the unmyelinated case and GMRes with AMG preconditioner for myelinated and large-grid simulations (Pods, 2019).

Implementation is organized as a parallel scientific-computing stack: Dune for grid and geometry, PDELab for grid operator and Newton, ISTL for solvers, MultidomainGrid and Multidomain for subdomain coupling, and LocalFunctions, with application module dune-ax1. Domain decomposition is in the $\rho = F \sum_i z_i c_i$5-direction only, with the membrane and adjacent Debye layers remaining on a single rank (Pods, 2019). The reported unmyelinated baseline has approximately 73,124 DOFs, sequential runtime of approximately 19 h to $\rho = F \sum_i z_i c_i$6, and parallel speedup of approximately 7.6× with 10 ranks. The reported myelinated 48-node axon has approximately 5,250,448 DOFs, and dense output can produce full dataset sizes on disk up to approximately 419 GB (Pods, 2019).

A distinct GPU-oriented implementation strategy appears in the 3D cue-diffusion guidance model, where OpenGL/GLSL is used for rendering and some parallel computations, while CUDA kernels implement calculateHash, radixSortPairs, identifyCellStart, updateParticles, updateAxons, and optional bfsTraversal (Mihut et al., 2014). There, cue fields are represented as particle systems rather than PDE solutions, neighbor queries inspect up to 27 cells, and the simulation exports connectivity and positional data for external analysis. This suggests two numerically different but architecturally compatible back ends: a field-resolving PDE path and a particle-based diffusion surrogate.

4. Guidance, branching, and mechanotactic modules

In the computational axon-guidance model, the essential objects are neurons, with soma, axon and growth cone, and glial cells acting as guideposts. The mechanisms include attraction/repulsion guidance cues, growth cone adaptation, tissue-gradient intersections, axonal transport, changes in the growth cone complexity and a range of responses for each receptor (Costa, 2015). Growth-cone guidance is defined as a function that maps receptor activation by ligands into a repulsive or attractive force, and this force is converted into a turning angle using spherical coordinates. The model also includes a regulatory network between the receptors and the intracellular proteins, intended to generate more complex and realistic behaviors, while ligand diffusion through the extracellular environment is modeled with linear or exponential functions (Costa, 2015).

Interstitial branching is treated differently in the AGM-based branching model. There, growth is governed by a dynamic field of axon guidance molecules (AGMs) satisfying a reaction–diffusion equation with degradation and point sources proportional to neuronal activity (Suleymanov et al., 2013). Growth-cone motion follows the AGM gradient,

$\rho = F \sum_i z_i c_i$7

while branch formation depends on pauses, concentration thresholds, shaft length, and a stochastic activity-dependent decision. Active points are recorded when the tip speed is below $\rho = F \sum_i z_i c_i$8, the AGM concentration lies in $\rho = F \sum_i z_i c_i$9, and shaft length exceeds $$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$0; after a delay $$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$1, a branch is created with probability

$$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$2

Candidate branch directions are sampled on a sphere of radius $$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$3 and filtered by an angular window such as $$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$4, $$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$5 (Suleymanov et al., 2013).

The GPU growth-and-guidance model abstracts receptor and intracellular dynamics more strongly. Axons elongate from an initial random direction, local gradients are detected through nearby cue particles, and steering is updated with fixed angular increment $$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$6 rad per step (Mihut et al., 2014). The paper reports two steering modes: a signaling SNR-inspired decision using a Gaussian decision variable and a resultant-vector steering rule in which nearby cue particle velocities and weights are aggregated into a local resultant. Connection formation occurs by sphere–sphere collision or threshold proximity, and growth stops on connection or when leaving the simulation bounds (Mihut et al., 2014).

Mechanotactic guidance enters through the durotaxis model, which is explicitly multiscale. At the molecular scale, clutch traction is generated by a motor–clutch framework with stiffness-dependent attachment/detachment kinetics and an effective substrate compliance. The traction force per filament is

$$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$7

and in steady regime it is written

$$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$8

where $$-\nabla \cdot (\varepsilon \nabla \phi) = \rho = F \sum_i z_i c_i,$$9 is dimensionless stiffness and $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$0 has a unique global maximum at $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$1 with $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$2 (Kassianides et al., 3 Jan 2025). At the growth-cone scale, the net traction is

$$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$3

and at the axon scale the shaft is modeled as a morphoelastic filament growing on an adhesive substrate (Kassianides et al., 3 Jan 2025). The sign of $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$4 determines whether guidance is positive durotaxis or negative durotaxis.

A stochastic alternative to deterministic steering is provided by the CTRW model. There the growth-cone trajectory is a piecewise-constant process with random waiting times $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$5 and jumps $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$6, and the propagator satisfies a generalized master equation with memory kernel $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$7 (Zhizhina et al., 2015). The key distinction is between finite-mean waiting times, which yield diffusive scaling, and heavy-tailed waiting times, which yield subdiffusive scaling. This provides a compact framework for representing differences between normal and mutant axons through waiting-time statistics rather than through explicit receptor networks.

5. Signal phenotypes, mutant scenarios, and validation logic

The electrodiffusion framework produces several signal regimes that are not available in a volume-conductor approximation. In the unmyelinated case, the reported action-potential propagation velocity is approximately $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$8 with baseline HH parameters (Pods, 2019). Near the membrane, the extracellular potential exhibits a characteristic triphasic EAP with P1 (positive), N1 (negative), P2 (small positive), S (slow negative curvature), P3 (rear positive), while in the Debye layer there is an AP echo satisfying $$\frac{\partial c_i}{\partial t} + \nabla \cdot \mathbf{J}_i = 0,\qquad \mathbf{J}_i = -D_i\,\nabla c_i - \frac{D_i z_i F}{R T}\,c_i\,\nabla \phi.$$9 with $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$0 and exponential decay over approximately $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$1 (Pods, 2019). The paper attributes the major near-field deviations from volume-conductor theory to ionic concentration contributions and the diffusion layer, whereas in the far-field the PNP extracellular action potential agrees well with the line-source approximation (LSA) after fitting resistivity $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$2, corresponding to $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$3 (Pods, 2019).

The myelinated formulation preserves the same electrodiffusive logic but redistributes active and passive regions. The reported setup uses nodes of Ranvier of width $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$4, internodes of width $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$5, and channel densities at nodes increased 10× to $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$6 and $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$7 (Pods, 2019). The reported propagation is approximately $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$8 faster than unmyelinated under the chosen parameters, with large ionic contributions at nodes and purely capacitive responses at myelin (Pods, 2019).

Validation is distributed across the model families rather than concentrated in one criterion. The PNP model reports analytical 1D PNP transient tests, spatial second-order convergence for $V_m = \phi^{\mathrm{CY}}-\phi^{\mathrm{ES}}$9 FEM, temporal first-order convergence for implicit Euler, and far-field consistency with LSA (Pods, 2019). The activity-independent guidance model reports the first computational model and a new theoretical model of the midline crossing of Drosophila axons that focus all the decision points, and states that the model allows describing to a great extent the behaviors reported in the literature for normal, comm mutant, and robo mutant pathfinding scenarios (Costa, 2015). It also suggests the existence of an inhibitory link between the DCC receptor and the Comm protein that is Netrin-mediated or mediated by a third unknown signal (Costa, 2015). In the branching model, simulated branch counts and activity dependence are compared qualitatively and quantitatively with pharmacological observations such as TTX and APV/DNQX manipulations, and the reported network morphologies are described as similar to neural network (Suleymanov et al., 2013). In the durotaxis model, the predictions include attractive zones of preferential stiffness in the substrate domain, reflective and refractive behaviour across interface between regions of different stiffness, and improved success rates in an in vivo-like Xenopus optic tract scenario when the optimal stiffness enters the accessible range (Kassianides et al., 3 Jan 2025).

The CTRW formulation provides a different type of validation target. Its central proposition is that for finite mean waiting time the lateral mean-squared displacement obeys

$C_m$0

whereas for a heavy-tailed waiting-time law $C_m$1 with $C_m$2,

$C_m$3

(Zhizhina et al., 2015). This suggests a compact statistical criterion for separating diffusive and subdiffusive axonal families.

6. Limitations, extensibility, and open problems

Each model family also states its restrictions. The PNP electrodiffusion framework assumes no fluid flow/convection, homogeneous diffusivities and permittivities (except membrane/myelin), no fixed membrane charges, no stochastic channel noise, and cytosol modeled as effective continuum (Pods, 2019). Cylindrical symmetry restricts multi-axon geometries to coaxial configurations, and the far-field Dirichlet boundary introduces a small absolute offset error in $C_m$4 that decays with domain size (Pods, 2019). The GPU cue-diffusion model uses a particle-based diffusion surrogate rather than a PDE solver, abstracts receptor dynamics, and does not implement branching in the reported baseline (Mihut et al., 2014). The branching model omits explicit inhibitory cues, substrate adhesion anisotropy, contact guidance, fasciculation, competition, pruning, and detailed synaptic plasticity (Suleymanov et al., 2013). The durotaxis model uses a 2D, small-deflection beam approximation, deterministic mean-field clutch dynamics, homogeneous adhesion density, and a simplified linear elastic substrate (Kassianides et al., 3 Jan 2025). The CTRW model is deliberately coarse-grained, encoding biochemical reaction times through waiting-time laws rather than explicit intracellular mechanisms (Zhizhina et al., 2015).

At the same time, the literature repeatedly frames these systems as modular and extendable. The PNP implementation is described as readily reusable and extensible as an “Axon Neural Simulation Framework” for high-resolution, parallel PNP simulations (Pods, 2019). The suggested extensions include full 3D PNP or asymptotic electroneutral reductions with Debye-layer boundary conditions, glial buffering, surface charges, heterogeneous extracellular tortuosity, frequency-dependent conductivities, validation against experimental juxtacellular signals, and improved HPC scaling through tailored multigrid preconditioners, block-structured solvers, and in-situ analysis/output (Pods, 2019). The durotaxis work proposes extension to 3D substrate fields and rod mechanics, stochastic clutch simulations, coupling to chemotaxis and electrotaxis, adaptive growth-cone geometry, and rate-and-state friction (Kassianides et al., 3 Jan 2025). The GPU work proposes anisotropic diffusion, receptor dynamics and adaptation, branching, and feedback from growing axons on cues (Mihut et al., 2014).

Taken together, these studies define the Axon Neural Simulation Framework as a layered research program rather than a single canonical model. Its core characteristics are explicit field representations, axon- or growth-cone-level update laws, parameterized constitutive relations, and scalable numerical infrastructure. The common unresolved issue is not whether such a framework can reproduce individual phenomena, but how far its components can be unified without losing the distinct explanatory power of electrodiffusion, chemical guidance, branching dynamics, stochastic waiting-time effects, and substrate mechanics.