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

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.