A computational investigation of preconditioning strategies and iterative methods for finite element based neurostimulation simulations

https://doi.org/10.1016/j.camwa.2020.02.021Get rights and content

Abstract

Computational simulations of transcranial direct current stimulation (tDCS) enable researchers and medical practitioners to investigate this form of neurostimulation with in silico experiments. For these computer-based simulations to be of practical use to the medical community, patient-specific head geometries and finely discretized computational grids must be used. As a result, solving the partial differential equation based mathematical model that governs tDCS can be computationally burdensome. Further, consider the task of identifying optimal electrode configurations and parameters for a particular patient’s condition, head geometry, and therapeutic objectives; it is conceivable that hundreds of tDCS simulations could be executed. It is therefore important and necessary to identify efficient solution methods for medically-based tDCS simulations. To address this requirement, we exhaustively compare the convergence performance of geometric multigrid and the preconditioned conjugate gradient method when solving the linear system of equations generated from a finite element discretization of the tDCS governing equations. Our simulations consist of three commonly used real-world tDCS electrode montages on MRI-derived three-dimensional head models with physiologically-based inhomogeneous tissue conductivities. Simulations are realized on fine computational meshes with resolutions deemed applicable to the medical community, and as a result, our finite element solutions highlight tDCS-specific phenomena such as electric field shunting that contributes to a notable intensification of the stimulation’s electric current dosage. Convergence metrics of each linear system solver are examined, and compared and linked to theoretical estimates. It is shown that the conjugate gradient method achieves superior convergence rates only when preconditioned with an appropriately configured multigrid algorithm. In addition, it is demonstrated that physiological characteristics of tDCS simulations make multigrid as a stand-alone solver highly ineffective, despite the fact that this method is typically effective in solving the tDCS-based Poisson problem. By identifying the solution methods optimal for medically-driven tDCS simulations, our results extend simulation support to high-resolution and high-volume computing applications, and will ultimately help guide tDCS numerical simulations towards becoming an integrated aspect of the patient-specific tDCS treatment protocol.

Introduction

Transcranial direct current stimulation (tDCS) is a non-invasive medical procedure that is designed to strategically stimulate specific areas of the brain. Low current, typically on the order of 1 mA, is delivered via electrodes placed on a patient’s scalp. The electrical current propagates intra-cerebrally, and alters the excitability of proximal neurons. The result is an increase (or decrease) in spontaneous neuron action potential generation. Numerous recent research findings demonstrate the potential of tDCS as a medical intervention. Individuals with Parkinson’s disease have demonstrated improved gait and extremity use [1], as well as enhanced working memory [2]. Research with Alzheimer disease patients has shown improved recognition memory [3], [4]. Chronic pain management can in part be addressed with tDCS [5], as well as facilitating and expediting post-stroke recover [6], [7].

In parallel to these medical findings, researchers in mathematics and computation have been producing simulations of tDCS to enhance the efficiency of tDCS treatments, and to expand our overall comprehension of this neurostimulation technique. Simulations have demonstrated how cerebral current density distribution can be predicted prior to treatment [8], [9], and how tDCS electrode positioning correlates with region stimulation [10], [11]. Other simulations have demonstrated the necessity of current amplitude dosage specificity [12]. Recently, software tools have been developed to facilitate patient-specific transcranial stimulation computer-based simulations by automating the generation of computational grids directly from a patient’s MRI data [13].

The quasi-static form of the Maxwell–Faraday equation with appropriate boundary conditions realistically depicts the current density distribution within the head from tDCS, and the finite element method is commonly used in solving this problem [10]. To most accurately model tDCS phenomena and to provide maximal benefit the medical community, high-resolution computational meshes are a necessity; resulting finite element discretizations therefore yield a large system of linear equations to be solved. The sheer magnitude of these linear systems virtually eliminates the use of direct methods to solve them, and iterative methods are much more appropriate. The effectiveness of the overall tDCS simulation is therefore directly correlated to the efficiency of the chosen iterative scheme. Iterative solver performance has been analyzed for non-neurostimulation EEG applications [14]. However, distinctive characteristics of tDCS simulations for medical applications, including: (i) an externally-delivered electrical stimulation, (ii) non-idealized geometries and boundaries derived from MRI data, (iii) extremely fine computational meshes, (iv) physiology-based, inhomogeneous tissue conductivities, and (v) multiple homogeneous and inhomogeneous boundary conditions, suggest that a comprehensive investigation of iterative solution strategies for tDCS simulations is warranted.

The conjugate gradient (CG) method is an ideal iterative method for symmetric, positive-definite linear systems [15], like those generated from tDCS numerical simulations. However, poor conditioning of the linear system coefficient matrix can result in extremely slow convergence; proper preconditioning of the linear system is needed to increase the solution convergence rate [16]. In this paper, we compare the convergence performance of the CG method with the following preconditioning strategies: none, Jacobi, symmetric Gauss–Seidel (SGS), incomplete LU decomposition (ILU), modified ILU (MILU), relaxed ILU (RILU), and geometric multigrid (MG). We also assess the robustness of MG as a stand-alone solver, as this method has shown to be highly effective in solving the tDCS-like classical Poisson problem [17]. When appropriate, comparisons are made with theoretical estimates of the CG method and also to results from related Poisson equation based simulations.

While both stand-alone MG and the CG method preconditioned with MG have perform adequately in solving the linear system of equations arising from discretizations of Poisson equation based models [17], the unique conglomerate of boundary conditions, complex geometries, and inhomogeneous tissue conductivities of tDCS presents a unique challenge for these methods, and results drawn from analyses of idealized scenarios are not necessarily applicable. In addition, to identify optimal iterative solver convergence rates, the (numerous) parameters of MG should be assessed in an environment that most accurately emulates the characteristics of the physical system and medical treatment in a non-simplified and non-idealized fashion. Further, if in fact MG does exhibit advantageous solution speeds, how superior are these convergence results? How do MG parameters affect solution convergence? Do different electrode arrangements impact MG solution performance? Do the computational advantages of multigrid outweigh the need for sophisticated, and potentially challenging, MG programming and software packages? The answers to these questions lie within the unique characteristics of medically-driven tDCS simulations and are the focus of this paper.

The convergence performance of each iterative scheme is assessed on simulations of three commonly used tDCS electrode configurations. Geometries are three-dimensional and derived from patient MRI scans. Five inhomogeneous tissue conductivities are used for the skin, skull, cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM), and three separate mesh refinements are examined. In addition to iterative solver metrics, scalp surface electric potentials and intra-cerebral electric current densities resulting from the simulated neurostimulation treatments are presented and discussed. These results showcase tDCS application-specific tendencies, such as highly non-linear electric fields, inhomogeneous electric current density distributions, and substantial shunting of the electric field on the scalp surface, each of which is quite visible given the fine mesh resolutions utilized in this paper.

In the materials and methods section, we present an overview of the numerics involved with our tDCS simulations including a brief derivation of the governing equations, the resulting weak formulation, and an overview of preconditioning and the preconditioners used in this paper. Then, we specify our computational tools, electrode montages simulated, and a detailed description of each numerical experiment. Results are presented next, and include our finite element solutions of the surface voltages and electric fields, as well as a detailed analysis of the convergence performance of each iterative method. We conclude with closing remarks and a discussion of future work.

To our knowledge, this manuscript is the first comprehensive comparison and assessment of preconditioners and iterative methods for tDCS simulations. Our results demonstrate the benefits of incorporating finely-tuned MG preconditioning to produce efficient, medically applicable simulations. We hope that these findings help promote tDCS simulations towards becoming a valued component of the tDCS treatment process.

Section snippets

tDCS simulation numerics

In this section we present details of the tDCS simulation numerics that we utilize in this paper. We first derive the governing equation and boundary conditions using the quasi-static approximation of the Maxwell–Faraday equation. We then present the weak formulation, followed by an overview of the preconditioning strategies used in this paper. We conclude with a brief description of MG, which is applicable to MG as both a preconditioner and stand-alone iterative method.

Results and discussion

Finite element solutions of the electrical potentials and electric current densities resulting from the simulated tDCS administrations are first presented. Then, convergence performance results of the linear system solvers are presented and discussed.

Conclusions

We have presented a comprehensive comparison of iterative solver strategies for tDCS simulations. Using several grid refinements and tDCS montage arrangements, the CG method for solving the finite element discretized system of equations was evaluated with a range of preconditioners, as well as MG as a standalone iterative solver. Of all methods examined, the CG method preconditioned with an appropriately configured multigrid algorithm produced superior convergence rates.

We found that different

CRediT authorship contribution statement

E.T. Dougherty: Conceptualization, Methodology, Formal analysis, Software, Validation, Visualization. J.C. Turner: Methodology, Formal analysis. F. Vogel: Methodology, Software, Resources.

Acknowledgments

We would like to thank Tobias Wittner of inuTech GmbH for his assistance with Diffpack and Multigrid software support. In addition, E. Dougherty is grateful for the support of the Roger Williams University Foundation to Promote Scholarship and Teaching, USA in providing resources that greatly facilitated this project.

References (43)

  • SchlaugG. et al.

    Transcranial direct current stimulation in stroke recovery

    Arch. Neurol.

    (2008)
  • HesseS. et al.

    Combined transcranial direct current stimulation and robot-assisted arm training in subacute stroke patients: an exploratory, randomized multicenter trial

    Neurorehabil. Neural Repair

    (2011)
  • NeulingT. et al.

    Finite-element model predicts current density distribution for clinical applications of tDCS and tACS

    Front. Psychiatry

    (2012)
  • GascaF. et al.

    Finite element simulation of transcranial current stimulation in realistic rat head model

  • DattaA. et al.

    Validation of finite element model of transcranial electrical stimulation using scalp potentials: implications for clinical dose

    J. Neural Eng.

    (2013)
  • KesslerS.K. et al.

    Dosage considerations for transcranial direct current stimulation in children: a computational modeling study

    PLoS ONE

    (2013)
  • WindhoffMirko et al.

    Electric field calculations in brain stimulation based on finite elements: An optimized processing pipeline for the generation and usage of accurate individual head models

    Hum. Brain Mapp.

    (2013)
  • van der VorstH.A.
  • LangtangenH.P.
  • MardalK.A. et al.

    Software tools for multigrid methods

  • PlonseyRobert et al.

    Considerations of quasi-stationarity in electrophysiological systems

    Bull. Math. Biophys.

    (1967)
  • Cited by (0)

    View full text