Abstract
Computational simulations of transcranial electrical stimulation (TES) are commonly utilized by the neurostimulation community, and while vastly different TES application areas can be investigated, the mathematical equations and physiological characteristics that govern this research are identical. The goal of this work was to develop a robust software framework for TES that efficiently supports the spectrum of computational simulations routinely utilized by the TES community and in addition easily extends to support alternative neurostimulation research objectives. Using well-established object-oriented software engineering techniques, we have designed a software framework based upon the physical and computational aspects of TES. The framework’s versatility is demonstrated with a set of diverse neurostimulation simulations that (i) reinforce the importance of using anisotropic tissue conductivities, (ii) demonstrate the enhanced precision of high-definition stimulation electrodes, and (iii) highlight the benefits of utilizing multigrid solution algorithms. Our approaches result in a framework that facilitates rapid prototyping of real-world, customized TES administrations and supports virtually any clinical, biomedical, or computational aspect of this treatment. Software reuse and maintainability are optimized, and in addition, the same code can be effortlessly augmented to provide support for alternative neurostimulation research endeavors.
1. Introduction
Transcranial electrical stimulation (TES) is a collection of noninvasive neurostimulation techniques that strategically modulate activity in regions of the brain with low magnitude electric current delivered through electrodes positioned on the scalp surface. Forms of TES include the commonly used transcranial direct current stimulation (tDCS), as well as transcranial alternating current stimulation (tACS) [1]. Recently, the use of numerous smaller sized electrodes, termed high-definition tDCS (HD-tDCS), has emerged as a form of TES that enhances electrical current focality [2, 3]. Clinical and biomedical research continue to demonstrate the capabilities of TES as a medical treatment. For example, Alzheimer and Parkinson’s disease patients have demonstrated increased memory abilities [4–6]. In addition, TES has shown to alleviate symptoms of psychiatric disorders including depression [7, 8] and schizophrenia [9–11].
The efficacy and comprehension of TES have been enhanced with mathematical modeling and computational simulation. In particular, simulations can compute the current density distribution for a given patient and TES apparatus configuration [1, 12–15] and have demonstrated the importance of modulating treatment stimulation dosage [16]. Other simulations have illustrated the importance of incorporating anisotropic tissue conductivity data [17, 18], and numerical studies aid in identifying the computational solution methods most efficient in simulating TES mathematical models [19].
Object-oriented design is a software design approach that defines objects, which are simply software entities that encapsulate data and functionality and the relationships among objects. Features of object-oriented design can be utilized to maximize code reusability, simplify software maintenance, and promote application versatility [20]. The prevalence of object-oriented design in scientific and biomedical applications began in the 1990s, and its use has dramatically increased due to its advantages over procedural software implementations [21–26]. In particular, mathematical and physical attributes of a model can be naturally represented with objects [27, 28].
A common approach for simulating partial differential equation (PDE) based mathematical models is to use prebuilt simulation software programs. These “black-boxes” can simplify model implementation; however, there are limitations with this approach. A researcher is often confined to the numerical algorithms and programming controls offered by the simulation program; it can be very difficult, or perhaps not possible, to incorporate numerics and solution techniques not supported by the software. In addition, integrating a prebuilt software application with external data sources and applications can be very challenging, which obstructs the use of its simulations within larger software solutions.
An alternative strategy is to create custom software that utilizes numerical application programming interfaces (APIs). While this approach generally takes longer to implement, it allows the use of problem-specific algorithms and programming logic that facilitate accurate and expedient simulation results [28]. Coupling this philosophy with object-oriented software engineering techniques can produce a final product that is versatile, expandable, and computationally efficient. Interactions with external systems can be seamlessly integrated, and the ability to compile the software to executable machine code simplifies its deployment to alternative hardware platforms, for example, medical imaging machines.
In this paper, we present an object-oriented software framework for finite element based TES simulations. Multiple features of object-oriented design and programming are utilized to create a modular software architecture that encapsulates medical, mathematical, and computational attributes of TES. The result is a program that maximizes code reuse and TES simulation versatility. We demonstrate this versatility with several simulations, each with a distinct TES research focus. These simulations utilize MRI-derived three-dimensional head models, with physiologically based tissue conductivities and real-world electrode configurations. Finally, we show how the same software can be easily extended to address alternative neurostimulation research areas, such as deep brain stimulation (DBS). In addition, all components of the framework are available to the community as a supplement to this paper.
2. Materials and Methods
2.1. Governing Equations
The electric current density within the head and brain from neurostimulation can be modeled by the Poisson equation; namely, , where is the electric potential, is the tissue conductivity tensor, and is a given electrical source term. For isotropic mediums, can be represented as a scalar, which varies for different tissue types.
Electric current delivered by TES anode electrode(s) is given by the nonhomogeneous Neumann boundary condition , where represents the stimulation current density and is the outward boundary normal vector. Cathode electrode(s) are represented by the homogeneous Dirichlet condition . All other points on the skin surface are insulated by the surrounding air; thus .
Our governing equations are as follows:where is the head volume, and represent the areas on the scalp covered by anode and cathode electrodes, respectively, and is the remaining portion of the head surface. Note that, for TES simulations, there is no source term within the volume, and so (1a). Alternatively, DBS is realized with an appropriate definition of and the homogeneous Neumann boundary condition (1d) applied to the entire boundary; namely, [29].
The associated weak formulation (see Appendix) is to find given such thatwhere
2.2. Framework Design
The fundamental task in the design of the framework was to create objects based on attributes of TES and TES simulations. We refer to [30] for a detailed explanation of object-oriented design and provide just a brief overview of the key aspects utilized within our framework.
Objects encapsulate data and functions that operate on the data. A class provides the description of an object type by defining variable and function names, and an object is more formally viewed as a specific instance of a class. Inheritance is a class relationship and provides the ability to extend the functionality of a class with a so-called subclass. A subclass can contain data and functionality from a parent class and can define its own as well. Inheritance is a major object-oriented design concept that exploits code reuse, since all subclasses can reuse data constructs and functions from a parent class. Polymorphism enables a subclass to give specific functionality to an abstract function defined by a parent class. This powerful technique allows parent classes to encapsulate generic and broad ideas by delegating specific function implementations to subclasses [30].
We required that electrode configuration, stimulation parameters, tissue conductivity information, computational domain, and numerical solution methods be specified via an input file. This permits customized TES simulations without the need to rewrite and recompile code. In addition, by mathematically retaining the source term, , in the weak formulation (2) and allowing this function to be defined by a user, different neurostimulation modalities, for example, DBS, can be simulated.
We choose to base our framework design on Diffpack [31], which is a numerical API library for solving PDEs. It is based on the C++ programming language [32] and possesses a vast collection of PDE solution algorithms [27, 28]. The C++ language incorporates well-established software engineering practices [33] and is compiled to machine code, resulting in fast execution speeds and portability to alternative hardware platforms. Despite our use of a specific numerical library and programming language, the object-oriented design and implementation strategies presented in this paper apply to any programming language and API package that supports an object-oriented approach.
Figure 1 displays the main software components in the TES simulation framework. Each box represents a class, and each arrow and diamond tipped line represents inheritance and aggregation, respectively. FEM is a class in the Diffpack library that offers fundamental finite element method data structures and algorithms. Class TES, which is the cornerstone of the framework, inherits FEM functionality to solve simulations given by system (1a)–(1d). Class TES_MG extends TES functionality so that multigrid (MG) algorithms can also be utilized in solving this system.
For TES simulations to be fully realized, tissue conductivity and electrode information are needed. These ideas have been encapsulated with respective classes, and TES possesses instances of each. Since conductivity data can come from a variety of sources with differing storage formats, the Conductivity class merely defines abstract function names that are viewed as common to all conductivity data sources. Then, specific functionality for specific conductivity data sources is implemented in subclasses, for example, Matlab_Data_Conductivity, via polymorphism. This approach allows different conductivity data sources to be incorporated into the framework without needing to modify code in the TES and Conductivity classes. The boundary conditions are also managed by TES, using information from class Electrode, which maintains TES electrode location and size information. In addition, by retaining the source term in the weak formulation (2) and encapsulating it as a class, namely, Source, assignments to enable the simulation of alternative types of neurostimulation, such as DBS.
2.3. Framework Implementation
In this section, key software implementation aspects of the framework and related Diffpack concepts are described. For a complete guide to Diffpack, see [27].
2.3.1. Tissue Conductivity
Class Conductivity defines general function names but delegates the implementation of these functions to its subclasses (Code 1). MatlabConductivity inherits Conductivity and provides specific implementations of the loadConductivities and getConductivity functions based on a Matlab [34] binary data file format produced by the SimNIBS software package [35].
Two additional data members are defined by MatlabConductivity, namely, a MatlabHelper object, mh, which manages a runtime interface with the MatlabEngine [34], and a Matlab matrix, ct. These two members are needed by the MatlabConductivity subclass, not the parent Conductivity class, and are therefore included in only the subclass. The MatlabConductivity implementation of the loadConductivities function simply creates a Matlab command to load the anisotropic conductivity data source, executes this command within the MatlabEngine via the mh object, and then stores the anisotropic conductivity tensor data into the ct matrix which are then accessible throughout a TES simulation via the getConductivity function (Code 1).
The Conductivity-MatlabConductivity relationship demonstrates the advantages of object-oriented inheritance and polymorphism. Conductivity is used to define function names common to all conductivity data sources and serves as the bridge between these repositories and the framework. Polymorphism permits the MatlabConductivity subclass to implement specific loadConductivities and getConductivity functionality for its particular data format. Code within the Conductivity class and all other framework classes does not dependent on these subclass implementations. Thus, alternative conductivity data sources can be easily incorporated into the framework as a subclass of Conductivity, in an identical fashion, requiring no modification to any other framework component.
2.3.2. Source Term
Class Source is a software abstraction of the source term (1a) and like Conductivity defines basic functionality to be implemented in subclasses (Code 2). Source inherits the Diffpack class FieldFunc which allows a scalar function to be defined over the domain. Different subclass implementations of the valuePt function can be used to model different neurostimulation modalities. For TES, valuePt in class TES_Source simply returns zero since in this case (Code 2).
Main TES Class. Class TES is the main class in the framework. Code 3 presents the key elements of this class, which contains numerous objects, functions, and Diffpack concepts. Class TES inherits the Diffpack FEM class, giving it access to finite element data structures and functionality. Handle objects, which are pointers in Diffpack that include memory management features, are defined for the computational grid, grid, and electric potential and current density solution results, u and currDensity, respectively.
Several other class members are needed to implement TES simulations. Variables to store boundary condition values are included, as well as the anodes and cathodes vectors to store electrode information. The vector isoSigma and matrix sigma store isotropic and anisotropic conductivity data, and the choice of conductivity representation is determined by a user via the isotropic boolean variable. The mc data member is a reference to the Conductivity class, providing an interface to conductivity data. Note that mc is type Conductivity and not MatlabConductivity. The mc object can therefore use the loadConductivities and getConductivity functions of any Conductivity subclass, including MatlabConductivity. This design approach allows new Conductivity subclasses to be defined and utilized without the need to modify code in class TES.
There is also a reference to a Source object, for example, TES_Source, as well as a field function for nonconstant stimulation currents, for example, tACS. Finally, TES inherits functions from FEM to perform finite element computations, including fillEssBC to set Dirichlet boundary conditions, calcElmMatVec to assemble the finite element linear system coefficient matrix and load vector, and the integrands and integrands4side functions to evaluate the weak formulation volume and boundary integrals, respectively.
2.3.3. Weak Formulation Integration
The weak formulation volume integrals (2) are computed in the TESintegrands function (Code 4). The Diffpack linear system assembler automatically calls this function as it iterates over the finite elements in the computational grid. Conductivity values for the current element are attained with a call to the updateConductivityTensors function, which determines if isotropic or anisotropic conductivities are desired and then populates the sigma matrix accordingly.
The source term over the current element is retrieved via a call to the valuePt function in a Source class. Then, the finite element basis functions for this element are iterated over, and the weak formulation is computed. Note that, in this implementation, sigma is assumed to be diagonal; however, this can be generalized with the incorporation of an additional loop. Finally, the coefficient matrix, A, and load vector, b, are updated. The boundary integral in the weak formulation is computed similarly in the integrands4side function, and its contribution is incorporated into b.
2.3.4. Multigrid
Multigrid is a linear system solution algorithm that utilizes multiple computational grid resolutions to achieve fast solution convergence [36]. By performing just a few iterations on each grid and then changing between finer and coarser grids, large portions of the error are efficiently removed [37]. In addition, as a preconditioner to the commonly used conjugate gradient (CG) method, MG is highly effective at solving linear systems that result from finite element based TES simulations [38]. Two commonly used MG cycles are the V-cycle and the W-cycle, and identifying the optimal cycle type for a given PDE system, in addition to other MG parameters including the number of grid levels, can be challenging [38].
As a subclass of TES, very few new data members and functions are needed in TES_MG (Code 5), since the entirety of TES is efficiently reused. The new members in TES_MG are a reference to the Diffpack multigrid toolbox, mgtools, and an integer representing the number of MG grid levels, no_of_grids. Finally, just one new function, mgFillEssBC, is needed to set the essential boundary condition (1b) on all grid levels.
2.4. Computational Tools
Numerical simulations were performed on a three-dimensional mesh (Figure 2) generated from human MRI images by the SimNIBS software package [35]. The mesh contains the skin, skull, cerebral spinal fluid (CSF), gray matter (GM), and white matter (WM) tissues of the head. Gmsh [39] enabled mesh visualization, identification of electrode coordinates, and grid conversion to a form supported by Diffpack. Electric potential and current density results were exported from Diffpack and visualized with ParaView [40] and gnuplot [41]. Anisotropic conductivity data for the GM and WM regions are provided by SimNIBS in a Matlab binary data file; these data are accessed by the MatlabEngine [34] via the MatlabHelper object of the MatlabConductivity class (Code 1).
(a) Head surface
(b) Brain region
(c) View of sagittal cross section
2.5. Computational Simulations
Multiple simulations were performed with the TES framework. Simulations were selected to demonstrate the framework’s versatility and ability to target medical, biophysical, and computational research objectives. Finally, to illustrate how alternative forms of neurostimulation can be simulated with the same software, we show a trivial extension to the framework that enables support for DBS applications.
Simulation 1 (comparison of isotropic and anisotropic conductivity data). Previous research suggests that incorporating anisotropic, rather than isotropic, brain tissue conductivity data is important for most accurately modeling TES electrical current distribution [17, 42]. To evaluate the impact that these two conductivity representations have on simulation results, TES simulations were performed with isotropic and anisotropic data. The three-dimensional domain shown in Figure 2 was used with approximately 2.8 million linear tetrahedra finite elements.
Anode and cathode electrodes were positioned at C3 and C4 [43], respectively, each with a surface area of approximately 16 cm2, and the anode electric current magnitude was set to 1.0 mA. This montage has been used to stimulate the motor cortex ipsilateral to the C3 anode electrode [44, 45]. First, isotropic electrical conductivities were assigned to different tissues: skin = 0.465, skull = 0.010, CSF = 1.654, GM = 0.276, and WM = 0.126, each with units (S/m) [46]. Then, anisotropic conductivity data were used via the MatlabConductivity class as previously described (see Code 1).
Simulation 2 (comparison of tDCS and HD-tDCS electrode montages). High-definition TES electrodes have demonstrated a greater ability to focus its electrical current on a targeted brain region than traditional tDCS, which uses two larger electrodes [2, 3]. Using the same computational grid as Simulation 1, the current densities produced by tDCS and HD-tDCS were compared.
For tDCS, the anode was positioned over C3 and the cathode over the contralateral supraorbital, each with a surface area of approximately 16 cm2. In comparison, high-definition electrodes, each circular with a 12 mm diameter, were positioned according to a 4 × 1 configuration; a single anode was positioned over C3, and four cathodes were placed approximately 5 cm radially from the anode, forming a square. Both of these montages are known to target the motor cortex region ipsilateral to the anode electrode [15, 47–50]. Anode stimulation strength for each simulation was once again 1.0 mA, and both utilized anisotropic conductivities.
Simulation 3 (comparison of finite element linear system solvers). Finite element based simulations of TES require the solution of large systems of linear equations, which can become a computational bottleneck for simulations performed on very fine meshes. Therefore, the effectiveness of a TES simulation is directly related to the efficiency of the chosen linear solver. The CG method is ideal for solving these linear systems and appropriate preconditioning can rapidly accelerate numerical results [36, 38].
The numerical efficiency of the preconditioned CG method was evaluated with TES simulations on the three-dimensional volume mesh (Figure 2), with approximately 29 million linear tetrahedra finite elements and roughly 5.1 million unknowns. The anode was positioned at CZ, with a stimulation strength of 1.0 mA, and the cathode at OZ [12, 51]. Simulations were performed with the CG method preconditioned with symmetric successive overrelaxation (SSOR), relaxed incomplete LU decomposition (RILU), and multigrid. The RILU relaxation parameter was set to 0.5 [27]. The MG preconditioner was simulated with a V-cycle with both two and three grid levels, and a W-cycle with three grid levels. The relative residual convergence tolerance was set to for all numerical experiments.
Simulation 4 (impact of conductivity representation on DBS simulation results). The importance of incorporating anisotropic conductivities in TES simulations suggests that accurately simulating other neurostimulation modalities may depend on this conductivity representation as well. In this numerical experiment, we compared DBS simulation results using both isotropic and anisotropic conductivities. A single electrode was positioned in the subthalamic nucleus (STN), the most commonly targeted location for DBS [52]. The electrode is a simplified version of the Medtronic (Model 3387) electrode used in humans [53]. The lower 5.0 mm of the electrode was modeled. The anode is positioned 1.0 mm from the electrode tip and the cathode is separated by 0.5 mm from the anode. Both the anode and cathode are 0.5 mm in height, and the overall electrode diameter was set to 0.75 mm (Figure 3). The anode and cathode metal contacts were modeled as conductors by setting the electrical conductivities of these regions to (S/m), and the remaining electrode shaft was modeled as an insulator with conductivity equal to (S/m) [54]. Because DBS electric current is proximal to the electrode [29, 55], a 10.0 mm × 4.0 mm × 4.0 mm subset of tissue around the electrode was considered as the computational domain.
To simulate DBS with the existing TES software framework, no existing code required modification. Rather, the only addition was a new subclass of class Source that defines (1a) for DBS; in total, less than ten lines of code were added. Code 6 displays this new subclass definition, Source_DBS, and highlights of its valuePt function implementation, which in this simplified scenario simply injects 1.0 mA of current from the anode contact into the surrounding tissue. In practice, DBS stimulations are time-varying pulses, and differing stimulation frequencies impact GM and WM tissue conductivities [56]. However, when examining electric field dispersion around a DBS electrode as in this simulation, a constant-valued source term in (1a) can be utilized [57].
3. Results and Discussion
3.1. Simulation 1
Electric current density results on the surface of the brain tissue are displayed in Figure 4. Viewing perspective is from above the head with the nasion facing up. As expected, the largest electric current density values are attained near the anode and cathode locations. Despite an overall similar pattern of electric current distribution, the simulated current densities, particularly in the motor cortex ipsilateral to the anode electrode, are significantly different between the isotropic and anisotropic simulations. For example, the maximal anisotropic current density value (0.434 (A/m2)) is approximately 18.2% higher than in the isotropic case (0.367 (A/m2)). Similar discrepancies are observed throughout the top surface of the brain, including a substantial impact on the paracentral lobule contralateral to the anode. In addition, these discrepancies extend to the GM and WM interiors. These results reinforce the importance of using anisotropic conductivities to most accurately model TES administrations and in addition showcase the capabilities of the object-oriented framework to support both isotropic and anisotropic TES simulations.
(a) Isotropic conductivities
(b) Anisotropic conductivities
3.2. Simulation 2
Electric potential and electric current density results are displayed in Figure 5. The tDCS montage results in a noticeably larger range of electric potential values (Figures 5(a) and 5(b)). However, the focus of the tDCS current density on the targeted region, in this case the motor cortex under C3, is much less than the HD-tDCS electrode montage (Figures 5(c) and 5(d)). Specifically, the brain tissue outside of the square formed by the HD-tDCS cathodes is virtually unstimulated, whereas tDCS results in a much greater electric current dispersion.
(a) tDCS electric potential
(b) HD-tDCS electric potential
(c) tDCS electric current density
(d) HD-tDCS electric current density
One limitation of this HD-tDCS electrode arrangement is the lower concentration of electric current that reaches the brain tissue. Specifically, the maximal electric current density in the brain from HD-tDCS is just approximately 23% of that achieved with traditional tDCS. The culprit for this effect is the shunting of the electric current around the poorly conducting skull; the HD-tDCS montage used in this simulation yields a minuscule amount of current that penetrates the skull to reach the CSF and brain tissues (Figure 6). Potential remedies for this behavior are a greater anode stimulation or positioning the cathodes at greater distances from the anode [3].
This simulation demonstrates the framework’s ability to support varying TES electrode montages, including high-definition electrode configurations. Results of this simulation corroborate that HD-tDCS offers greater electric current focality; however, its net stimulation of brain tissue is potentially much lower than tDCS.
3.3. Simulation 3
Electric potential and current density results on the head surface are displayed in Figure 7. Viewing perspective is from behind the head, and as anticipated, maximal and minimal electric potential and current density values occur at the anode and cathode, respectively. The shunting of the electric current around the skull, due to its low electrical conductivity, results in the segregated current density pattern around the anode and cathode centers (Figure 7(b)). This phenomenon has been previously observed [1]; however, it is highlighted by the very fine mesh resolution used in this simulation.
(a) Electric potential
(b) Current density
Figure 8 displays convergence history curves and performance metrics for each CG preconditioning strategy. With no preconditioning, the CG method will eventually convergence but will accomplish this extremely slowly, requiring more than 6 hours to solve the linear system. The SSOR and RILU preconditioners demonstrate comparable performances, both in solution time and number of iterations. Multigrid preconditioning with two grid levels has far fewer iterations with 166 than both SSOR and RILU and yet has a greater run time. This observation can be explained by the fact that a single iteration of MG has embedded iterations on the different mesh refinements [37]. However, three-grid-level MG noticeably accelerates convergence rates, with both the V- and W-cycles outperforming the other preconditioners. Three-grid-level MG with a W-cycle pattern has slightly fewer iterations than its corresponding V-cycle, yet this V-cycle preconditioner is approximately 11.3% faster.
(a) Convergence history
(b) Numerical iterations and linear system solve time. Boldface values indicate best convergence results
The ability of the framework to support numerically oriented TES research is presented in this simulation example. The results indicate that the CG method combined with an appropriately configured MG preconditioner is highly efficient in solving the linear systems produced in TES computational simulations.
3.4. Simulation 4
Figure 9 displays the electric current densities produced by the DBS electrode using isotropic (Figure 9(a)) and anisotropic (Figure 9(b)) conductivities. In both cases, the majority of the current density is proximal to the electrode, indicating that accurate placement of electrodes in DBS procedures is paramount [56]. While the maximal current densities of the isotropic and anisotropic scenarios are similar, other differences between them are observable. First, the current density around the electrode is nonsymmetric in the anisotropic case, as is expected with directionally dependent conductivity data. In addition, anisotropic conductivities result in a greater electrical current intensity adjacent to the electrode and more dispersion into the neighbouring tissue. This is also observed in the electrical field lines, where anisotropic conductivities produce a more intense and further-reaching current.
(a) Isotropic
(b) Anisotropic
With a straightforward extension, DBS was simulated with the object-oriented TES framework. The entirety of the framework is reused in the DBS simulation, which demonstrates how object-oriented design can produce efficient and scalable software implementations. This particular example further reinforces the importance that anisotropic conductivities play in accurately modeling neurostimulation.
While this DBS scenario utilizes a constant source term to assess electrode electric field dispersion [57], alternative applications may demand the use of time-dependent pulses, such as those that examine temporal neuronal responses to nonconstant electric current administrations. In these cases, different stimulation frequencies have an impact on the electrical impedance; specifically, WM and GM conductivities are frequency dependent where an increase in stimulation frequency yields an increase in electrical conductivity.
Applications such as these are supported by the software framework. First, since isotropic tissue conductivity values are specified within the configuration input file, and not hard-coded in the software, alternative DBS pulses can be accurately simulated by modifying these values. For anisotropic data, modifications to support frequency-dependent conductivities can be made within the Conductivity classes. For example, the SimNIBS software package volume normalizes the WM and GM anisotropic conductivities such that the mean conductivity value of each tensor is maintained to the associated tissue’s isotropic value [35]. Therefore, for different DBS frequencies, tensors could be updated to be normalized to different isotropic values within the framework’s Conductivity components. To implement the time-varying source itself, the existing definition of the valuePt function in class Source contains an argument for time, namely, dpreal t. Thus, implementations of valuePt in the Source_DBS class can be based on time, and system (1a)–(1d) can then be iteratively solved for with time-based source values.
4. Conclusions
Computational simulations of neurostimulation are a valuable tool that enable researchers to investigate this form of brain therapy in silico, and as simulations become more refined, their utility to medical and biomedical research grows. While prebuilt simulation software programs can simplify and expedite TES model implementation, they possess application and portability limitations. In addition, recreating custom software as applications and research objectives change is inefficient, error-prone, and tedious to maintain.
Since the mathematics and physiology that govern all forms of neurostimulation are the same, a single, well-designed software code can support a versatile range of neurostimulation simulations. In this paper, we have presented one such software framework, described its design and implementation, and demonstrated its abilities to support diverse neurostimulation research objectives. The cornerstone of the design stage was to encapsulate general neurostimulation concepts into modular software objects. In doing so, a multitude of TES simulation areas are supported, and in addition, alternative forms of neurostimulation, for example, DBS, can be simulated with the same software.
Simulation results show the importance of using anisotropic conductivities in both TES and DBS. This is especially important for simulations utilized in selecting patient-specific neurostimulation parameters. In addition, results demonstrate the ability of HD-tDCS to focus the electrical current on a specific brain region; however, this capability must be balanced with a potentially low concentration of electric current actually reaching the target. The object-oriented TES framework is an ideal tool for this scenario, enabling different permutations of HD-tDCS electrode montages and stimulation strengths to be conveniently investigated to identify an optimal configuration for a particular patient’s data and therapeutic objectives.
Results also illustrate that appropriately configured multigrid preconditioning can achieve superior convergence rates. It is conceivable that hundreds of simulations could be run to identify an optimal electrode montage and neurostimulation parameters for a particular patient. Hence, efficiently solving the linear system of equations resulting from TES finite element simulations is crucial, and MG can be used in this capacity to greatly decrease simulation run times. Finally, we demonstrated how a minor addition to the object-oriented TES framework enables simulations of DBS. The DBS simulation example that was performed utilizes a constant stimulation source term; however, the software framework description and simulation ensemble presented in this paper motivate how the framework could be used to address DBS research related to electrode placement and parameter values.
In future work, we plan to extend the framework to support anisotropic transcranial magnetic stimulation (TMS), in a similar fashion as was demonstrated for DBS. In addition, we plan to utilize the framework to more thoroughly investigate correlations between HD-tDCS interelectrode distances and electric current distributions.
Appendix
Weak Formulation
We multiply (1a) by a test function and integrate over to obtain and using Green’s theorem we have Expanding the surface integral givesand substituting the boundary conditions (1c) and (1d) yields For these integrals to exist, we require and , where and we enforce the Dirichlet boundary condition (1b) on our solution space by further stipulating that for
Consequently, we have the following weak formulation.
Given , find such that
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge Frank Vogel and the entire inuTech team for their assistance with Diffpack. The authors would also like to acknowledge the financial support received from Virginia Tech’s Open Access Subvention Fund.