Anisotropic conductivity tensor imaging for transcranial direct current stimulation (tDCS) using magnetic resonance diffusion tensor imaging (MR-DTI)

Transcranial direct current stimulation (tDCS) is a widely used non-invasive brain stimulation technique by applying low-frequency weak direct current via electrodes attached on the head. The tDCS using a fixed current between 1 and 2 mA has relied on computational modelings to achieve optimal stimulation effects. Recently, by measuring the tDCS current induced magnetic field using an MRI scanner, the internal current pathway has been successfully recovered. However, up to now, there is no technique to visualize electrical properties including the electrical anisotropic conductivity, effective extracellular ion-concentration, and electric field using only the tDCS current in-vivo. By measuring the apparent diffusion coefficient (ADC) and the magnetic flux density induced by the tDCS, we propose a method to visualize the electrical properties. We reconstruct the scale parameter, which connects the anisotropic conductivity tensor to the diffusion tensor of water molecules, by introducing a repetitive scheme called the diffusion tensor J-substitution algorithm using the recovered current density and the measured ADCs. We investigate the proposed method to explain why the iterative scheme converges to the internal conductivity. We verified the proposed method with an anesthetized canine brain to visualize electrical properties including the electrical properties by tDCS current.


Introduction
Transcranial direct current stimulation (tDCS) was originally developed for brain injuries such as strokes or major depressive disorder. tDCS has reported to be effective in treating a number of neurological and psychiatric disorders. Recent findings show novel possibilities by confirming significant changes in the pain-related neural networks among patients with chronic pain [1,2]. tDCS can improve performance on a variety of cognitive tasks in humans, in both healthy and in patients, despite its exact mode of action remains unclear [3][4][5] memory, a necessary condition before it can be used as a therapeutic tool. For the long-term heroin-addicted subjects, the tDCS over bilateral frontal-parietal-temporal (FPT) region has been shown to significantly reduce subjective craving score induced by heroin queues of heroin addicts [6]. tDCS is a therapeutic modality as an add-on treatment that showed promising results for major depressive disorder [7][8][9]. tDCS typically uses the intensity of injection current between 1-2 mA for a duration of up to 20 minutes. However, it is insufficient to directly estimate the internal electrical properties (current density and/or conductivity) when tDCS operates by sending low direct current through the attached electrodes on the head. Up to now, computational modeling and simulation techniques have been used to characterize and understand the effects of tDCS, including electrode size, the number of anode and cathode locations, current density, and electrical field.
Using an MRI scanner, combining with the external current injections, electrical properties (current density and/or conductivity distribution) have been widely studied [10,11]. Recently, direct visualization of magnetic field changes by tDCS current was studied [12]. Magnetic resonance electrical impedance tomography (MREIT) technique aims to visualize the internal current density and the isotropic conductivity by measuring only one component of magnetic flux density using an MRI scanner [11,13]. Using MREIT modality, for the human brain cases, the internal current density distribution in the human brain from the measured magnetic flux density by tDCS current was studied [14]. However, no direct techniques exist to visualize the electrical properties including anisotropic conductivity tensor, electric field, and apparent concentration of ions by tDCS current. To visualize the isotropic conductivity, MREIT requires at least two independent injected currents because an infinite number of isotropic conductivity distributions can generate the same internal conductivity [11,15]. However, biological tissues show anisotropy due to asymmetric cellular structures. A few studies for anisotropic conductivity image reconstructions based on the projected current density have been studied. Even though they reported the experimental data using multiple independent injection currents, but still have difficulty determining the other three components of the conductivity tensor [16,17].
Diffusion-weighted imaging (DWI) is an imaging modality which measures the random Brownian motion of water molecules within a voxel of tissues such as brain white matter and skeletal muscle by applying the diffusion-sensitizing magnetic field gradient [18][19][20][21]. The statistical relationship between the mobility and the diffusivity was established by Einstein: where hΔr 2 i is the average squared displacement of a particle over the time interval, Δt, and D is the diffusion coefficient which incorporates temperature and media viscosity properties.
In this paper, we aim to investigate the electrical properties including the current density and the anisotropic conductivity tensor by the external current injection through attached electrodes during tDCS. Adopting the effective macroscopic anisotropic tensor modeled by a two-phase anisotropic medium, the electrical conductivity tensor, the mobility of water molecules, and the water diffusion tensor have been analyzed in terms of the intra-and extra-cellular transport coefficients [22]. Based on the model, the anisotropic conductivity tensor map has been investigated using an approximate linear relationship between the conductivity tensor and those of the water diffusion tensor [23][24][25].
The electrical conductivity in the extracellular space (ECS) can be decomposed into the concentration of ions and the mobility of charge carriers such as ions and charged molecules. Using the linear relationship between the conductivity tensor and those of the water diffusion tensor, the conductivity tensor is a product of scalar scale parameter and the diffusion tensor, where the scalar scale parameter mainly reflects the concentration of ECS ions. We propose an update scheme to reconstruct the scale parameter. By analyzing the update scheme and investigating how the updated conductivity converge to the true conductivity tensor, we provide more accurate insight into the current flow patterns and the conductivity characteristics.
Despite various applications in tDCS clinical and cognitive researches, the electrodes position on the head is sensitive to the internal current flow delivered to the brain [26]. As the potential applications of the proposed method for electrical properties (current density, conductivity, and electrical field), we can directly simulate the internal current flow for arbitrary electrodes configuration by solving the forward elliptic partial differential equation with the reconstructed internal conductivity distribution. The proposed method can design the electrodes configuration that is able to target specific areas of the brain with the reconstructed electrical properties.
Through an animal experiment, we verify that the proposed method can provide the electrical properties including the current density, ion-concentration information, and the conductivity tensor map due to tDCS.

Measurement of magnetic flux density
The current injected into the brain through tDCS produces the internal current density distribution. The induced internal current density induces the magnetic flux density, B = (B x , B y , B z ), by the Biot-Savart law. An MRI scanner allows to measure the z-component magnetic flux density by the injected current. A multi-echo in injected current non-linear encoding (ICNE) MREIT by currents from the end of RF pulse additionally accumulates the magnetic flux density B z as a total phase offset. By comparing the l-th complex data S l after injection current and S 0 without injection current, the z-component of magnetic flux density distribution generated by the injected current, B l z , can be measured: where γ = 26.75 × 10 7 rad/Ts is the gyromagnetic ratio of the hydrogen, N E is the number of echoes, T E l is the l-th echo time, α l and β l are the imaginary and real parts of S l and S 0 , respectively [10]. Due to a small amount of injection current magnitude by tDCS (1-2 mA), the intensity of measured B z data is within several nT ranges and the measured B z data include non-negligible noise artifact. The weighted combination B z ¼ P N E l¼1 w l B l z using the solved weighting factor w l reduces noise artifacts of magnetic flux density [27]. The noise standard deviation of B z , sd(B z ), depends on the intensity of MR signal and the time width of injection current: where S c is the measured complex signal and T c is the time width of injection current per repetition time (T R ).

Initial current density
Using the background conductivity σ 0 , we can generate an initial current density J 0 = −σ 0 ru 0 by solving the following equation: where g is the injected current density through the attached electrodes, ν is the outward normal vector, and O is an imaging object. The initial current density J 0 satisfies the divergence-free condition r Á J 0 = 0. According to the Biot-Savart law the z-component of magnetic flux density can be easily updated from the current density.

Relationship between conductivity tensor and diffusion tensor
Let D denote a 3 × 3 positive definite symmetric diffusion tensor matrix from DTI, which measures the intrinsic diffusion of water molecules in tissue: where the column vectors of S D = {s 1 , s 2 , s 3 } are the orthonormal eigenvectors of D, the superscript T denotes the transpose and d i for i = 1, 2, 3 are the corresponding eigenvalues. The signal intensity ρ D of diffusion MRI is given by where ρ 0 is the signal obtained without diffusion-sensitizing gradient, g j is the j-th normalized diffusion-sensitizing gradient vector and b denotes the diffusion-weighting factor depending on the gradient pulse used in the DT-MRI sequence. The b-value can be expressed as where δ, Δ, and G are the duration of applied gradient, the duration between the paired gradients, and amplitude, respectively. By accepting the effective macroscopic anisotropic tensor model by Sen and Torquato [22], the effective electrical conductivity tensor in the case of macroscopically anisotropic media shares the eigenvectors with the water diffusion tensor allowing the well-known Hashin-Shtrikman bounds [28]. The eigenvalues l s i of the conductivity tensor C σ satisfy the following relation for the relatively small intra-cellular diffusion coefficient d int : where σ ext is the extra-cellular conductivity, d ext is the extra-cellular diffusion coefficient and Oðd 2 int Þ is bounded as d 2 int tends to infinity. Based on the effective anisotropic conductivity tensor model, the eigenvalues of the electrical conductivity and the water diffusion tensors approximately satisfy a linear relationship: where the coefficient η denotes a scale parameter between the eigenvalues l s i and l d i , i = 1, 2, 3. The scale parameter η is a scalar function influenced by the coefficients of d int , d ext , σ ext , and d ext . For the anisotropic conductivity tensor C σ , the divergence-free condition of current density and applied external current density satisfy The relation (1) implies that

Computation of projected current density J p
Due to the limitation of the measured B z signal which is the z-th component of magnetic flux density induced by tDCS, the vector field J p as the projected current density of J is directly calculated from the measured B z data as follows where O t is an imaging slice of O, J 0 is the background current density, and the potential ψ solvesr Here, we use the following two-dimensional terminologies: rf ≔ @f @x ; @f @y ; 0 ;r ? f ≔ @f @y ; À @f @x ; 0 : The projected current density J p is the quantity that is optimally observable from the measured B z data [13]: where k Á k O t denotes the L 2 -norm and the constant C does not depends on J.

Diffusion tensor J-substitution algorithm
From the measured diffusion tensor and the recovered current density, the water diffusion coefficient can distinguish the discontinuous regions classified according to the diffusion difference of water molecules. The proposed DT-J-substitution algorithm is based on the distinguishable diffusivity signals of diffusion tensor across the discontinuous regions. Since a single-shot echo planar imaging (SS-EPI) pulse sequence is typically used to measure the diffusion tensor map, the measured each component of D is noisy due to the fast scan MR pulse sequence. The reconstructed current density J p from the measured B z data also contains nonnegligible noise due to a small amount of injection current. Utilizing the informations of J p (% J) and D, the aim is to determine the scale parameter that satisfies the following condition: Since the noise level of B z is inversely proportional to the magnitude intensity, we introduce a neighborhood N r ¼ fr i : i ¼ 1; Á Á Á ; Ng of a voxel r and set a weighting factor as where S is the magnitude of the MR image and h is a function of the noise level of B z .
We propose an iteratively formula for the anisotropic conductivity tensor C σ : The n-th updated voltage potential u n satisfies Here, the voltage potential u n is uniquely determined up to a constant. The proposed algorithm for the current density J p , and the conductivity tensor C σ , using only tDCS can distinguish the discontinuous regions.

Characteristics of DT-J-substitution algorithm
For a discontinuous anomaly D σ and a point ξ 2 @D σ , the relation between the conductivity and the current flow satisfies where u + and u − denote The normal component of the current density J is continuous across the interface of the discontinuous anomaly. The relation (6) implies that only one current flow for visualizing discontinuous conductivity tensor C σ is insufficient to distinguish the edge of anomaly.
Since the updated scale parameter η n+1 is determined by satisfying the condition, −η n+1 Dru n % J = −ηDru, we have the following identity: Using the identity (7), the integration by part leads to The right hand side of (8) can be rewritten as Z O ZððDruÞ Á ru À ðDru n Þ Á ru n Þdr The applied external injection current intensity g in (5) and the integration by part yield Z O ðZDru À Z n Dru n Þ Á ðru n þ ruÞdr By the relations (8), (9), and (10), we have the following identity The identity (11) shows that the updated scheme is deeply related to the directions of Dru n and ru. From the measured diffusion tensor and the current density, for an anomaly region D σ with different electrical concentrations of ions, we have the following relationship at a point ξ 2 @D σ : The distinguishable signals of D on the edges of anomalies can updated the conductivity tensor even with the measured one component of magnetic flux density.

Animal preparation
In this paper, we verify the proposed method with a canine brain experiment that used previously to visualize the conductivity tensor using two independent current injections and ADCs with 32 diffusion gradient directions [15]. We select a partial data set using one injection current and ADCs to estimate the diffusion tensor. A healthy laboratory beagle (3-4 years old; weighing 5-10 kg; Harlan Interfauna, Huntingdon, Cambridgeshire, UK) was used for imaging experiments. The dog had no signs of neurologic problems on physical examination. The dog was screened for metabolic disease by means of CBC and serum chemistry analysis and acclimated to the facility for a minimum of 1 week prior to being used in this study. The dog was housed individually and fed twice daily with a commercial dry food (Science Diet, Hill's Pet Nutrition, Topeka, KS). Fresh water was supplied continuously by an automatic dispenser at our well-ventilated facility under controlled light-dark cycles (light on, 08:00 to 20:00). An indoor temperature range of 18 to 24˚C and a humidity level of 55% ± 10% with 8 air changes hourly were maintained. Before MR imaging experiments, we injected 0.1 mg/kg of atrophine sulfate to prevent dribbling. Ten minutes later, we anesthetized the dogs with an intramuscular injection of 0.2 ml/kg Zolazepam (Zoletil 50, Virbac, France). After clipping and shaving the hair around the chosen imaging area, we attached three pairs of carbon-hydrogel (HUREV Co. Ltd, Korea) surface electrodes (the size of each electrode was 30 × 30 × 5 mm 3 ) to the skin. During MR scanning, we intubated the animal using an endotracheal tube of 8.5 mm diameter and began general anesthesia using an anesthesia system (VME, MATRX, USA). We used 2% isoflurane mixed with oxygen at 800 ml/min flow rate. Ventilation was machine-controlled

Results
We injected 2 mA currents in the diagonal direction, from the left top to the right bottom in Fig 1A, using a custom-designed MREIT current source and acquired MR images using the coherent steady state multi-gradient echo (CSS-MGRE) pulse sequence. After MREIT scans, we performed DT-MRI scans using the single-shot spin-echo echo planar imaging (SS-SE-EPI) pulse sequence to measure the diffusion coefficient for the three diffusion gradient directions. Table 1 shows the parameters for the MR experiments. Fig 2A shows the cylindrical structure by accumulating the imaging slices. Fig 2B shows the three dimensional current flow in the brain region. In Fig 2C, the recovered projected current densities were displayed by the current flow streamlines overlap the T 2 weighted image in the first, second, and third slices. See S1 Fig for the three dimensional current flow in the brain region by injecting current through another pair of attached electrodes. Fig 3A shows the estimated diagonal components of the diffusion tensor of water molecule, D ii , i = 1, 2, 3 and six regions of interest (ROIs) including a local region of 3 × 3 voxels to estimate the numerical values of the scale parameter η and the conductivity tensor. The six ROIs were selected with respect to the tissue anisotropy characteristics (white matter (ROI-1, 2), gray matter (ROI-3, 4), and CSF (ROI-5, 6)). Fig 3B shows the reconstructed scale parameter and the diagonal components of the conductivity tensor, C σ , using tDCS current. We updated 2-times by following (4). To suppress the noise artifact accumulated in the diffusion tensor and the current density, we chose a search neighborhood N r ¼ fr i : i ¼ 1; Á Á Á ; 25g around an imaging voxel r. Since the noise artifact is severe in the measured B z data, the range of B z is within several nT, the weighting factor was determined in the N r as In Fig 3C, we compared the recovered scale parameterZ using tDCS current with those using two independent injection currents. By substituting the linear relationship C s ¼ZD and utilizing the curl-free condition, r × ru = 0, we have the following direct relation between the measured current density and diffusion tensor and the scale parameter: We directly recovered the scale parameterZ using two linearly independent injections currents by solving the following matrix system [25]: where J p 2 denotes the projected current density using the measured B z data by injecting assistant injection current and (D −1 J p ) x and (D −1 J p ) y denote the x-and y-component of the vector D −1 J p , respectively. In the second columns of Fig 3B and 3C, we compared the diagonal components of the reconstructed conductivity tensor using tDCS current and two independent injection currents, respectively. Table 2 shows the estimated values of the reconstructed scale parameters and diagonal components of diffusion tensors within the marked six ROIs in Fig 3A. The proposed DT-J-substitution algorithm depends on a local behavior of the current density and the diffusion tensor from ADCs of water molecules. The reconstructed scale  Fig 4A,ĉ 2 , using two injection currents. Since the current density J p 2 compensated the weak region where the intensity of J p (r) Á ru n (r) was small, the recovered mean conductivity using the two injection currents by combining the current densities J p and the assistant current density J p 2 was definitely advantageous comparing to that using tDCS current. However, the usage of two independent injection currents is not suitable for actual tDCS situations. In that sense, the  reconstructed conductivity images in Fig 4A show the feasibility of the DT-J-substitution method.
The proposed DT-J-substitution algorithm iteratively updates the scale parameter η based on the identity (11). The identity (11) implies that the convergence of the iteratively updated η n mainly relates to the angle between the velocity vector Dru n and the electric field E = −ru: For the electric field E, we solved the forward problem (2) with the scale parameterZ in (13). Fig 5A and 5B show the intensity of the projected current density and the angle map α in (14), respectively. Comparing to the recovered mean conductivity images in Fig 4A, the contrast of the reconstructed mean conductivity relates to the characteristics of the current density and the angle map between the velocity vector and the electric field distribution. The DT-J-substitution algorithm uses the mismatch condition on the boundary of discontinuous region in (12). Without the diffusion tensor information, Fig 5C shows the apparent isotropic conductivity by using J-substitution algorithm using tDCS current: s nþ1 ðrÞ ¼ À J p ðrÞ Á ru n ðrÞ ru n ðrÞ Á ru n ðrÞ where the n-th updated voltage potential u n satisfies ( r Á ðs n ru n Þ ¼ 0 in O s n ru n Á n ¼ g on @O Due to the relationship between the conductivity and the current flow in (6), the apparent  Fig 6D and 6E compared the reconstructed conductivity tensor maps using tDCS current and two injection currents, respectively, inside the rectangular ROI shown in Fig 6A. Since both of the conductivity tensors shared the same eigenvectors from the diffusion tensor, their orientations were same. The position-dependent conductivity tensor image recovered by the proposed DT-J-substitution method were quite similar to the reconstructed conductivity tensor using the two injection currents. The elliptic radii in Fig 6D and 6E were proportional to the eigenvalues of the tensor at each voxel.

Discussion
The measured magnetic flux density induced by tDCS dominantly reflects the extracellular effects because the cell membrane is regarded as an electrical insulator. The recovered current density directly from the measured B z data also reflects only the ECS effects. However, The diffusion tensor of water molecules experiences other complex different environments even within a single voxel and includes the diffusion properties of intra-and extra-cellular compartments. The proposed method actually recovered the scalar scale parameter that connects the conductivity tensor and the diffusion tensor. The reconstructed conductivity tensor using the measured magnetic flux density induced by tDCS current and the diffusion tensor is influenced by the intra-and extra-cellular compartments due to the measured ADCs. In this paper, we assume that ADC at relatively low b-value(< 1000 s/mm 2 ) dominates the fast diffusion component [29].
Based on the fact that the fast diffusion is proportional to extracellular diffusion, we used the diffusion tensor from ADCs with multiple diffusion gradient directions as an approximated extracellular diffusion tensor. In the proposed DT-J-substitution algorithm, the updated scale parameter satisfies the following structure: where α is the volume fraction of ECS and D I denotes the intracellular component of the diffusion tensor. By assuming a relatively small intracellular diffusion coefficient, the updated scale parameter contains some error by the intracellular volume fraction and D I . Using multi-b values and multi-diffusion gradient directions, the high angular resolution diffusion imaging (HARDI) technique has been introduced to overcome the limitation of DTI [23]. In recent years, non-Gaussian diffusion methods permit the analysis of the DW signal over a various range of b-values. The neurite orientation dispersion and density imaging (NODDI) requires at least two HARDI shells with different b-value for estimating the neurite morphology, the technique requires relatively long acquisition times [30]. The extracellular diffusion coefficient can be approximately estimated as a parameter using the NODDI technique. The diffusion models have been developed by considering the complexity and heterogeneity of the neuronal tissue microstructures. The proposed anisotropic conductivity model can be a different form depending on the diffusion model.
Although some specific issues (cables and electrodes, MR pulse sequence) need to be addressed when using typical tDCS device combining with the MR scanner, using the advent of MR-compatible tDCS systems, both sequential and concurrent acquisitions are possible using the conventional MR pulse sequences. In this paper, we used the coherent steady state multi-gradient echo (CSS-MGRE) MR pulse sequence to accumulate the phase signal by the injected current through tDCS and the single-shot spin-echo echo planar imaging (SS-SE-EPI) pulse sequence to measure the diffusion coefficient for the three diffusion gradient directions. For further issues for the blood-oxygen-level dependent (BOLD) contrast functional MRI and MR spectroscopy (MRS), MR field maps with and without tDCS should be carefully designed and analyzed. Combining with the electrical properties (conductivity, current density, and electric field), the interpretation of BOLD fMRI and MRS data has sufficient potential for future tDCS applications.
As a non-invasive brain imaging, fMRI is a popular method to the treatment of patients with neurological and neuropsycohological disorders, based on the increase in blood flow to the local area. The fMRI does not directly measure electrical activity since the blood flow is an indirect measure of neural activity. However, the reconstructed electrical properties using the proposed DT-J-substitution algorithm directly reflect the neural activity. Although complex neural activities may cause loss of signal due to self-cancellation of neural currents in an MR imaging voxel and the measured magnetic flux density signal is very sensitive to noise artifact, the proposed approach has potential to monitor neural electrical activity generated by the brain at a higher spatial resolution.
Transcranial alternating current stimulation (tACS) is a neuromodulatory technique similar to tDCS, which uses a sinusoidal AC current for noninvasive investigation of brain oscillations. tDCS effects depend on the applied frequency, amplitude, and phase. tACS applied currents to the scalp in the EEG frequency range (0.1-80Hz) can interfere with oscillatory brain activity. Although the tACS effects have been reported to involve manipulation and delivery of brain vibrations and corresponding biochemical changes, the exact mechanisms are still unclear. For the alternating current frequency for tACS, we can modify the MR pulse sequence instead of the CSS-MRGE which is appropriate for the injection current procedure of tDCS. To accumulate the magnetic flus density by the oscillatory current flow, we can design the MR pulse sequence by applying synchronized 180˚rephasing pulses based on the multi-echo spin echo (ME-SE) pulse sequence [27].
We believe that the DT-J-substitution algorithm is a promising approach to obtain electrical properties in the brain by tDCS current. The proposed method reconstructs the electrical properties in a voxel-by-voxel by observing local behaviors of the current density and the diffusion tensor. Our future work will include the development of an efficient and stable algorithm for electrical property imaging by accounting for local and global characteristics of the current density and the diffusion tensor.

Conclusion
We introduce a method to visualize the electrical properties including the electrical anisotropic conductivity, effective extracellular ion-concentration, and electric field using tDCS currents in-vivo by measuring ADCs and the magnetic flux density induced by tDCS. The proposed DT-J-substitution algorithm iteratively updates the scale parameter, which connects the anisotropic conductivity tensor to the diffusion tensor of water molecules. We analyze the relationship between the current flow and the scale parameter. We found that the contrast of reconstructed conductivity tensor by tDCS current mainly depended on the intensity of the current density and the angle between the velocity vector and the electric field. Through an animal experiment, we demonstrate the feasibility of the recovered electrical properties using a single tDCS current.