A parallel and efficient algorithm for multicompartment neuronal modelling

doi:10.1016/j.neucom.2005.12.040

Neurocomputing

Volume 69, Issues 10–12, June 2006, Pages 1035-1038

https://doi.org/10.1016/j.neucom.2005.12.040 Get rights and content

Abstract

An important metric for simulation algorithms used in compartment modelling is computation efficiency. One algorithmic achievement in efficiency is the Hines method, which substantially reduces the computation cost of solving a system of linear equations arising in each time step of implicit time integration. However, the Hines method does not work for circuits containing gap-junction loops. In this paper, we propose an algorithm with which the Hines method extends its applicability to loop-containing circuits with efficiency almost the same as that of the Hines method for loop-free circuits. Furthermore, our algorithm has good parallelism, promising effective utilization of parallel computing power for large-scale simulations.

Introduction

The dynamics of neuron membrane potential is governed by a partial differential equation of the parabolic type. To quantitatively understand the spatiotemporal properties of the membrane potential, one needs to solve the governing equation. Compartment modelling is an approach with which the governing equation is solved numerically. With the compartmental modelling, the neuronal membrane is divided into a finite number of interconnected anatomical compartments which are short enough so that the potential within each compartment is almost the same. On this set of compartments, the governing equation is spatially discretized, and a temporal integration method with time step size $Δ t$ is then used so that starting with the initial condition at time $t = 0$ , the solution at time $t = (k + 1) Δ t$ is computed using the solution or solutions at time $t = k Δ t$ or earlier. The algorithmically simplest temporal integration method is the forward Euler method, which has a very low computation cost for each time step and also has high parallelism. The major limitation of the forward Euler method is its poor stability that limits the time step to a very small size, which, quoting Pearlmutter and Zador [5], “very rapidly becomes the limiting factor for large numbers of compartments”. Implicit temporal integration methods like the backward Euler and the Crank–Nicolson have good stability and accuracy, but require solving an equation of the form $(I + c Δ t A) V^{k} = r^{k}$ for each time step, where $I$ is the identity matrix, $V^{k}$ , an n-entry vector with n being the number of compartments, holds the values of the potential in all compartments at the $k$ th time step, and $A$ is an $n \times n$ sparse matrix. On an unbranched cable, the matrix A is tridiagonal, which makes Eq. (1) solvable with a low cost of $O (n)$ operations for a cable with n compartments. On a branched tree structure, A is not a tridiagonal matrix, but can form a matrix of some special property when the compartments are ordered in a certain way discovered by Hines [2], a property which enables an efficient solving in an $O (n)$ computation cost. However as pointed out by Hines himself, the Hines method does not work for circuits containing gap-junction loops. In this paper, for Eq. (1) that arises in implicit temporal integration, we propose a solution method which extends the applicability of the Hines method to loop-containing circuits with an efficiency almost the same as that of the Hines method for loop-free circuits. In addition, the solution method also creates parallelism, promising an opportunity for the numerical algorithm to effectively utilizing parallel computing power for large modelling problems.

Section snippets

The one-iteration Schwarz–Hines algorithm

To achieve high efficiency for simulations on loop-containing circuits, we propose to use a non-iterative domain decomposition technique (see [6], [7], [8]) with which the entire circuit is partitioned into several loop-free sub-circuits so that the Hines method is applicable on each sub-circuits. The domain decomposition technique we employ is the one-iteration Schwarz domain decomposition method Kuznetsov [3], [4] and Chen and Lazarov [1] introduced for the elliptic equations obtained from

Numerical experiments

To test the proposed solution algorithm, we constructed two passive neurons with a resting potential of $- 60 mV$ , connected by gap-junctions as in Fig. 1. Each neuron has two dendritic branches originating from the soma. The soma has a diameter of $30 μ m$ and a length of $30 μ m$ . Each dendritic branch has a length of $600 μ m$ from soma to the tip of dendrite, and tapers from $3 μ m$ in diameter at the end connecting soma to $2.4 μ m$ at $300 μ m$ from soma where the dendrite branches into two children of length $300 μ m$

Acknowledgement

The author is grateful to Dr. Frances Skinner for the information on how to model gap junctions in computer programs.

Yu Zhuang is an Assistant Professor of Computer Science at the Texas Tech University. He graduated from Zhejiang University (BS in Applied Mathematics) and Louisiana State University (MS's in Math and System Science, PhDs in Math and Computer Science). His research interests include developing fast, accurate, stable, and/or parallel simulation algorithms and software for multicompartment modelling.

References (8)

H. Chen et al.
Domain splitting algorithm for mixed finite element approximations to parabolic problems
East–West J. Numer. Math.
(1996)
M. Hines
Efficient computation of branched nerve equations
Internat. J. Biomed. Computing
(1984)
Y.A. Kuznetsov
New algorithms for approximate realization of implicit difference schemes
Soviet J. Numer. Anal. Math. Modell.
(1988)
Y.A. Kuznetsov
Domain decomposition methods for unsteady convection–diffusion problems

There are more references available in the full text version of this article.

Cited by (9)

A multiscale numerical method for the heterogeneous cable equation
2012, Neurocomputing
Citation Excerpt :
Although the ever increasing capabilities of computers facilitate such endeavors, a big chunk of the advances is certainly due to better modeling and computational techniques. Several facets of multiple scales and neuroscience pervade the literature, from general discussions on modeling [11,25,38], to topics related to numerical and computational aspects [2,3,8,22,27–30,40–42,45,49]. Powerful mathematical tools from asymptotic analysis, like homogenization, matching asymptotics, singular perturbation, and dimension reduction techniques are commonly employed, along with various degrees of neurological based arguments [4,5,7,17,34,36,39,37,44,46].
Several interesting problems in neuroscience are of multiscale type, i.e. possess different temporal and spatial scales that cannot be disregarded. Such characteristics impose severe burden to numerical simulations since the need to resolve small scale features pushes the computational costs to unreasonable levels. Classical numerical methods that do not resolve the small scales suffer from spurious oscillations and lack of precision.
This paper presents an innovative numerical method of multiscale type that ameliorates these maladies. As an example we consider the case of a cable equation modeling heterogeneous dendrites. Our method is not only easy to parallelize, but it is also nodally exact, i.e., it matches the values of the exact solution at every node of the discretization mesh, for a class of problems.
To show the validity of our scheme under different physiological regimes, we describe how the model behaves whenever the dendrites are thin or long, or the longitudinal conductance is small. We also consider the case of a large number of synapses and of large or low membrane conductance.
An alternating explicit-implicit domain decomposition method for the parallel solution of parabolic equations
2007, Journal of Computational and Applied Mathematics
Explicit–implicit domain decomposition (EIDD) is a class of globally non-iterative, non-overlapping domain decomposition methods for the numerical solution of parabolic problems on parallel computers, which are highly efficient both computationally and communicationally for each time step. In this paper an alternating EIDD method is proposed which is algorithmically simple, efficient for each time step, highly parallel, and satisfies a stability condition that imposes no additional restriction to the time step restriction imposed by the consistency condition, which guarantees a convergence of order $O (Δ {th}^{- 1} \sqrt{N_{B} / N}) + O (h^{2})$ in an $H^{1}$ -type norm, where $N_{B}$ and N, respectively, denote the number of gridpoints on the interface boundaries B and the number of gridpoints on the entire discrete domain.
On the Hamiltonian cycle mapping onto 3-D torus interconnection network based on base-b reflected gray codes
2007, Applied Mathematics and Computation
In this paper we present a 3D large volume simulation decomposition and mapping technique onto a 3D torus interconnection network, based on base-b reflected gray codes. Such simulation leads to fast execution plus an improvement of total execution time. The minimization of the execution time of the simulation is due to minimization of packet routing in the 3D interconnection torus network. Applying the algorithm on real 3D space plasma simulation of the Aurora region of the Earth’s Ionosphere was observed to be almost ideal speedup. The efficiency of the implementations was shown to be constant close to one, as the size of the 3D simulation increases.
Numerical methods and analysis of multiscale problems
2017, SpringerBriefs in Mathematics
Predictor-Corrector Domain Decomposition Algorithm for Parabolic Problems on Graphs
2012, Mathematical Modelling and Analysis
Efficient parabolic solvers scalable across multi-architectural levels
2012, Proceedings of the 2012 10th IEEE International Symposium on Parallel and Distributed Processing with Applications, ISPA 2012