Decomposition of high-frequency electrical conductivity into extracellular and intracellular compartments based on two-compartment model using low-to-high multi-b diffusion MRI

Background As an object’s electrical passive property, the electrical conductivity is proportional to the mobility and concentration of charged carriers that reflect the brain micro-structures. The measured multi-b diffusion-weighted imaging (Mb-DWI) data by controlling the degree of applied diffusion weights can quantify the apparent mobility of water molecules within biological tissues. Without any external electrical stimulation, magnetic resonance electrical properties tomography (MREPT) techniques have successfully recovered the conductivity distribution at a Larmor-frequency. Methods This work provides a non-invasive method to decompose the high-frequency conductivity into the extracellular medium conductivity based on a two-compartment model using Mb-DWI. To separate the intra- and extracellular micro-structures from the recovered high-frequency conductivity, we include higher b-values DWI and apply the random decision forests to stably determine the micro-structural diffusion parameters. Results To demonstrate the proposed method, we conducted phantom and human experiments by comparing the results of reconstructed conductivity of extracellular medium and the conductivity in the intra-neurite and intra-cell body. The phantom and human experiments verify that the proposed method can recover the extracellular electrical properties from the high-frequency conductivity using a routine protocol sequence of MRI scan. Conclusion We have proposed a method to decompose the electrical properties in the extracellular, intra-neurite, and soma compartments from the high-frequency conductivity map, reconstructed by solving the electro-magnetic equation with measured B1 phase signals.

characteristics of the soma size and density belonging to the intracellular compartment. The SANDI model needs to include the direction-averaged DWI signal at high b-values ( ≥ 3000 s/mm 2 ) to detect the apparent soma size and density. Combining with the estimated DWI data for the higher b-values, we apply the SANDI model to estimate the micro-structures of extracellular compartment and use it to identify the characteristics of low-frequency conductivity.
To determine the micro-structural parameters of biological tissues using the measured Mb-DWI data, we use the random forest regression, an ensemble machine learning algorithm by constructing a multitude of decision trees at training time [20,21]. The applied machine learning method builds a forest of uncorrelated trees, combined with randomized node optimization and bootstrap aggregating [22].
To demonstrate the electrical property decomposition from the recovered high-frequency conductivity, we generate the high-frequency conductivity using the convection-reaction partial differential equation with a small regularization parameter [1]. The recovered high-frequency conductivity reflects the combined electrical properties by the intra-and extra-cellular compartments using a routine protocol sequence of MRI scan.
To validate the proposed method, a conductivity phantom including a giant vesicle suspension was conducted as a model for the membranes of biological cells. The giant vesicles were cell-like materials with thin insulating membranes to validate the two compartment model including both extracellular and intracellular spaces. We conducted human experiments by comparing the results of reconstructed conductivity of extracellular medium and the conductivity in the intra-neurite and intra-cell body. We extracted the apparent total ion concentration from the high-frequency conductivity map and the estimated micro-structural parameters. The phantom and human experiments indicate that the total ion concentration and the extracellular diffusion tensor can predict the extracellular electrical properties without externally injected currents. The accuracy and precision of the reconstructed low-frequency conductivity distribution in the extracellular compartment were evaluated.

Giant vesicle phantom experiment
Giant membrane vesicles dispersed in aqueous solution were prepared as described in [23]. Giant vesicles are a model membrane system for the investigation of lipid membranes. In a 1L round-bottom flask containing 3 mL of chloroform and 400 μm of methanol, 2 mL of phospholipids dissolved with a chloroform solution of 30 mg/mL concentration under argon atmosphere. A 20 mL volume of distilled water or 0.75% NaCl solution was carefully added not to disturb the interface between the aqueous phase and organic solution phase. The flask was installed to a rotary evaporator to remove organic solvent at 47 °C under vacuum for 20 min at 10 rpm and then followed by another 20 min at 60 rpm. During evaporation of organic solvents, phospholipids were assembled to form giant vesicles. The resultant aqueous solution containing giant vesicles were centrifuged at 1500 rpm for 10 min. The volume fraction of the giant vesicles was about 80 to 90% by visual observations. The mean and standard deviation of the diameters of the giant vesicles were 13 ± 4.7 μm.
The conductivity phantom was constructed, which comprised two compartments of electrolyte (Background) and giant vesicle suspension (ROI) (Fig. 1a). The background was an NaCl solution of 3 g/L and ROI was filled with the giant vesicles suspension. The giant vesicle suspension was a mixture of the same NaCl solution of background and giant vesicles. In the giant vesicle suspension, the ratio of NaCl solution and giant vesicle was the same, which means that the volume fraction of the intracellular space was about 0.4-0.45.
The conductivity values of the giant vesicle suspension in ROI (Fig. 1a) were measured by using an impedance analyzer (SI1260A, AMETEK Inc., UK). The conductivity values were 0.57 S/m at 1 kHz and 1.05 S/m at 5 MHz. Figure 1b shows the estimated high-frequency conductivity, σ H , extracellular volume fraction, f ec , extracellular diffusivity, D ec , and extracellular conductivity, σ ec , a b c d Fig. 1 a Magnitude image. b-d Estimated high frequency conductivity,σ H , extracellular volume fraction, f ec , extracellular diffusivity, D ec , and extracellular conductivity, σ ec for different noise levels of ( SNR P , SNR D ). b ( ∞, ∞ ), c (100,50), d (50,10)) Page 5 of 20 Lee et al. BioMed Eng OnLine (2021) 20:29 respectively. Table 1 shows the mean and standard deviation values of the estimated high frequency conductivity, extracellular volume fraction, extracellular diffusivity, extracellular ion concentration and extracellular conductivity within ROI. The estimated high-frequency conductivity was σ H = 1.01 ± 0.02 S/m in ROI, which is close to the conductivity value 1.05 S/m measured using the impedance analyzer at 5 MHz in ROI. We used the diffusion term c = 0.05 to stabilize the convection-reaction partial differential equation in (3). Since the volume fraction of the giant vesicles after the centrifugation was about 80% to 90% and the ratio of NaCl solution and giant vesicles was same in ROI, the controlled extracellular volume fraction was about 0.55-0.60 in ROI. The value of f ec in ROI was 0.59±0.14 as expected.
To estimated the extracellular ion concentration, c ec , in (12), we assume that the intra-neurite ion concentration and soma ion concentration are same, and further assume that c in = βc ec for some constant β . Since we used the same electrolyte for both inside and outside of the giant vesicle, the ratio of ion concentrations in the intracellular and extracellular compartments, β , was 1. The extracellular ion concentration c ec can be estimated as With the estimated high-frequency conductivity and the parameters of the SANDI model for brain micro-structure, we recovered the extracellular ion concentration, c ec , extracellular conductivity, σ ec . σ ec = 0.54 ± 0.15 S/m in ROI is close to the measured conductivity value 0.57 S/m at 1 kHz.
To verify the proposed method, we conducted two cases. We artificially destroyed the measured signals by adding random noise and motion artifacts. For noisy data, we added the zero-mean Gaussian random noise to both the phase data and diffusion weighted signal. The noise standard deviation were calculated by S /( √ 2 SNR), where S is the average signal (phase signal and the signal obtained without diffusion gradient) amplitude and SNR is the signal-to-noise ratio in MR magnitude images. For the phase data, SNR P values 200, 100, and 50 were employed. For the diffusion weighted signal, we assumed that SNR D of the signal obtained without diffusion gradient, S 0 , were 100, 50, and 10. The random noise added to S 0 was added to each S b . Figure 1c, d shows images of the estimated high frequency conductivity, σ H , extracellular volume fraction, f ec , extracellular diffusivity, D ec , and extracellular conductivity, σ ec for ( SNR P , SNR D ) = (100,50) and (50,10). Table 2 shows the mean and standard deviation values of the estimated high frequency conductivity, extracellular volume fraction, extracellular diffusivity, extracellular ion concentration and (1) extracellular conductivity within ROI for the chosen noise levels. The proposed method stably reconstructed the conductivity values when SNR P and SNR D were above 100.
To investigate the impact of subject motion, we simulated a sinusoidal motion in the phase-encoding direction with amplitude 1 and a frequency of 2 Hz and obtained a motion-corrupted complex MR signals which were the convolution of the measured complex MR signals with the corresponding point spread function [24].  Figure 2b-d shows images of the estimated high frequency conductivity, σ H , and extracellular conductivity, σ ec for the different diffusion term c in (3). Since the diffusion term c acts as a low-pass filter in the reconstruction algorithm, it was possible to reconstruct the motion-corrected high-conductivity map using a high c value. For c = 1 , the values of the estimated high frequency conductivity and extracellular conductivity within ROI were 0.99±0.02 and 0.53 ± 0.14, respectively. Table 2 Mean and standard deviation values of the high frequency conductivity, σ H , extracellular volume fraction, f ec , extracellular diffusivity, D ec , extracellular ion concentration, c ec , and extracellular conductivity, σ ec within ROI for different noise levels of ( SNR P , SNR D ) (

Human experiment
MRI measurements were performed with three healthy volunteers without a documented history of any disease were recruited. The participants were located inside the bore of a 3T MRI scanner with the head coil in transmit and a 32-channel RF head coil (Achieva TX, Philips Medical Systems, the Netherlands). All experimental protocols were approved by the institutional review board of Kyung Hee University (KHSIRB-16-033). All methods were carried out in accordance with the relevant guidelines and regulations and all participants provided written informed consent. For MREPT imaging experiments, the multi-spin-echo pulse sequence with multiple refocusing pulses was adopted to minimize the measured noise. Before the data acquisition, we applied a volume shimming method with the volume defined to cover the brain region. Imaging parameters were as follows: TR/TE = 1500/15 ms, NE = 6, NEX = 1, slice thickness = 4 mm, FOV = 260×260 mm 2 , imaging matrix size = 128×128× 5, and scan time = 16 min. After MREPT scans, we performed DWI scans using the single-shot spin-echo echo planner imaging (SS-SE-EPI) pulse sequence. We applied the diffusion weighting gradients in 15 directions with 4 b-values of 1000, 2200, 3000 and 3600 s/mm 2 , respectively. Imaging parameters were as follows: TR/ TE = 2000/70 ms, δ/� = 21/33 ms, NEX = 2, slice thickness = 4 mm, and acquisition matrix size = 64×64× 5. The scan time was about 6.2 min. The matrix size of 64×64× 5 was extended to 128×128× 5 to match the spatial resolution (2.03×2.03× 4 mm 3 ) of MREPT experiment. To increase the spatial resolution with the matrix size of 64 × 64 to the spatial resolution of 128 × 128 by reducing scanning times without much loss in resolution or SNR, we used the zero-filling interpolation (ZIP) by zero filling the high spatial-frequency components of the raw k-space data, which places the acquired data in the central regions on k -space and fills the data of zero in the outer regions. An additional conventional T 1 weighted scan of 2 min was included for anatomical reference. Figure 3 shows the estimated micro-structural parameter maps of SANDI model (5) for the first human subject: the estimated extracellular volume fraction f ec , intra-neurite volume fraction f ne , soma volume fraction f so , extracellular diffusivity D ec , and intra-neurite diffusivity D in from the first subject. Figure 4a, b shows the MR magnitude and the B1 phase images at the third imaging slice of the first subject. For quantitative analyses, T 1 image was segmented into cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM) using the segmentation tool of Statistical Parametric Mapping (SPM 12) [25]. These regions are shown in Fig. 4c. Figure 5 shows the details of reconstructed conductivity images from the first subject. To estimate the high frequency conductivity, σ H , with the acquired transceiver phases of B1 maps, we solved the convection-reaction partial differential equation (PDE) in (3) with the diffusion term c = 0.02.
To estimate the extracellular ion concentration, c ec , in (12), we set that the ratio β for ion concentrations in the intracellular and extracellular compartments to be 0.41 as suggested in [6] by adopting reference values of intracellular and extracellular ion concentrations of four predominant ions ( Na + , Cl − , K + , and Ca 2+ ). Using the reference ratio value β = 0.41 , the extracellular ion concentration c ec can be estimated by (1). With the estimated high-frequency conductivity and the parameters of the SANDI model for brain micro-structure, we recovered the extracellular ion concentration, c ec , extracellular conductivity, σ ec , intra-neurite conductivity, σ ne , and soma conductivity, σ so . Figure 6 shows the 3 × 3 extracellular conductivity tensor, C ec , and intra-neurite conductivity tensor, C ne , images. To estimate conductivity tensors, we used the water molecule diffusion tensors with the b value of 1000 s/mm 2 . We fixed the principal diffusion direction (eigenvector corresponding to the maximum eigenvalue of the diffusion tensor) and solved the equations (16) and (17) to obtain the extracellular diffusion tensor, D ec , and intra-neurite diffusion tensor, D ne , images. Using these diffusion tensors, we reconstructed the extracellular conductivity tensor, C ec , and intra-neurite conductivity tensor, C ne , in (18) and (19), respectively. Table 3 summarizes the estimated values of high frequency conductivity, extracellular ion concentration, extracellular conductivity, and intra-neurite conductivity of each subject. In Table 3, longitudinal (L), transverse (T), and average (A) white matter conductivity values were found by computing the principal eigenvalue, the mean of the other two eigenvalues, and the mean of all eigenvalues of the conductivity tensor over all region, respectively. The high-frequency conductivity values in CSF regions were between 1.43 and 1.58 S/m and the extracellular conductivity values were between 1.35 and 1.49 S/m. Note SANDI model does not take into account CSF compartment. For this reason, a  The extracellular conductivity values in GM regions were higher than the average extracellular conductivity values in WM regions, whereas the average intra-neurite conductivity values in WM regions were always higher than the intra-neurite conductivity values in GM regions. Figure 7 shows the estimated extracellular conductivity tensor, C ec , and intra-neurite conductivity tensor, C ne , images represented by tri-axial ellipsoids, respectively, in the rectangular ROIs shown in extracellular conductivity, σ ec , and intra-neurite conductivity, σ ne , images. Since all of the conductivity tensors shared the same eigenvectors from the diffusion tensor, their orientations were same. The radii of each ellipsoid are proportional to the eigenvalues and their axes are oriented along the directions of eigenvectors. As expected, the volume of intra-neurite conductivity ellipsoids appeared larger in WM regions.

Discussion
Literature conductivity values including those obtained from direct impedance measurements (IM), MREPT and diffusion tensor MR electrical impedance tomography (DT-MREIT), are summarized in Table 4. The conductivity values of brain tissues heavily depend on the biological tissue structures, participant's age and pathology, frequency, and the measurement conditions (in vivo, ex vivo, and in vitro) [26][27][28]. Without external injection currents and using conventional MR pulse sequences minimizing magnetic field inhomogeneity, MREPT is a promising research area for practical clinical medical devices. Comparing to MREPT, by injecting a dc current into the imaging subject, MREIT can reconstruct images of the internal low-frequency conductivity distribution.   20:29 In this paper, we proposed a new way to decompose electrical properties in each compartment (extracellular, intra-neurite, and soma compartments) from the reconstructed high-frequency conductivity using MREPT technique and the micro-structural parameters using SANDI model. From the intrinsic noise in the MR measurements, MREPT reconstruction techniques have been proposed to improve the quality of the high-frequency conductivity map [33][34][35][36]. The high conductivity values in Table 3 using the reconstruction algorithm in [1] also depend on the first term in (3), which is the diffusion term to stabilize the solution. The diffusion term acts as a low-pass filter, leading to some blurring of the final highconductivity maps.
Because SANDI model does not take into account CSF compartment, a slight difference between the high-frequency conductivity values and the extracellular conductivity values was found in CSF region as was described in Result section. SANDI model estimates the orientation-independent features of microstructure using the directionaveraged DWI signals. In future studies, it is worth exploring a new model that take into account CSF compartment. In the SANDI model, the direction averaged signals using many uniformly distributed directions are advantageous to estimate the microstructural parameters since the DWI data is noisy at a high-b value. Although the accuracy of estimated parameters can be improved as the diffusion directions are increase, the numbers of diffusion directions were 30 and 15 for the phantom and human experiment, respectively, to recover the low-frequency electrical properties within a feasible MR scanning time.
At the low frequency, the internal electrical current flow caused by external current stimulation occurs only in the extracellular space and CSF, excluding the intracellular space, due to the insulation properties of cell membrane [27,37,38]. Despite extensive researches for the electrical properties in ECS, the reported low frequency conductivity values do not match and show considerable standard deviation. Recently, a meta-analysis of reported human head electrical conductivity values at low frequency (< 1 kHz) provides a recommended value estimated under suitable and realistic conditions, in which data acquisition techniques were categorized into five groups including directly applied current and MREIT. In Table 4, the reference values for low-frequency conductivity were 1.71 ± 0.3 (CSF), 0.47 ± 0.24 (GM), and 0.22 ± 0.14 S/m (WM) with broadly similar results to ours [30].
At fairly low frequencies, the conductivity of brain tissue, particularly white matter, is known to be anisotropic [32]. Comparing to the anisotropicity of intra-neurite compartment, the average ratio between longitudinal WM extracellular conductivity values and transverse WM conductivity values of 1.87-2.00 found here were lower than 3 and 5.9 reported in [31,32], respectively. The extracellular conductivity values in GM and WM regions exhibited much stronger frequency dependencies compared to CSF region, because of their complicated tissue structures. More rigorous analysis including geometric and viscous components of the tortuosity of ECS will be needed in the future work.
In this paper, the acquired DWIs are combined because SANDI model assume an isotropic diffusion to distinguish the intra-and extra-cellular diffusion signals. We focus on separating the microscopic parameters for extracellular space from the multiple DWI data. The recovering procedure for the low-frequency electrical properties is very Page 13 of 20 Lee et al. BioMed Eng OnLine (2021) 20:29 complex and requires further developments, including measurement techniques, noise artifact reduction, main magnetic field inhomogeneity artifact correction, and more reliable models to separate the extracellular and intracellular compartments.
To compensate the difficulties of measuring DWI data for higher b-values, we adopt the hypothesis that the DWI data corresponding to the b-value range (1000-3600 s/mm 2 ) reflects distinguishable diffusion signals between the soma and the extracellular space. To avoid the ill-posedness to determine the six unknowns, f in , f ec , D in , D is , D ec , and r s , in the two compartment model (5) from the smoothly decayed exponential curves, the parameter of D is was fixed as 2 ×10 −3 mm 2 /sec. However, the determination of parameters by matching the observed DWI data for Mb-DWI data and SANDI model was still very sensitive to measured noise artifacts.
As described in the above, there are several problems to overcome, especially including the follows • SANDI model is not yet suitable for clinical MRI scanner because it requires high b-values (higher than 3000 s/mm 2 ) DWI and relatively many diffusion directions to apply direction averaged DWI signals. • To estimate the extracellular ion concentration, c ec , it is assumed that the intra-neurite and soma ion concentrations are the same. Moreover, the ratio of ion concentrations in the intra-and extra-compartments is assumed as a fixed constant, β = 0.41 for human experiments. Since the ratio β depends on the apparent ion concentration, β cannot be fixed as a constant in the whole brain region. However, no method to experimentally estimate of the human brain is available.
Nevertheless, the method of extracting low-frequency electrical properties in a noninvasive way from the high-frequency conductivity is critical for clinical usefulness. Electrical brain stimulation (EBS) techniques, such as transcranial direct current stimulation (tDCS) and deep brain stimulation (DBS), are promising treatments for human disorders [39][40][41][42][43]. Since there is no clear explanation for the mechanism, EBS studies have relied on computational modeling using reference conductivity values in the whole brain region. The proposed electrical property decomposition from the high-frequency conductivity can be a promising work for the EBS techniques.

Conclusion
We have proposed a method to decompose the electrical properties in the extracellular, intra-neurite, and soma compartments from the high-frequency conductivity map, reconstructed by solving the electro-magnetic equation with measured B1 phase signals. By decomposing the electrical conductivity into the product of mobility and charged carrier concentrations, voxel-wise micro-structures including the extracellular volume fraction and diffusivity were investigated using SANDI model by analyzing the Mb-DWI data based on a two compartment model. In the SANDI model, to distinguish the intra-soma compartment from the mixed diffusivity signals in ECS, the Gaussian phase distribution approximation of the tail was used. To determine the micro-structural parameters in each separated compartment using Mb-DWI data, a machine learning algorithm, random forests, was used by constructing a multitude of decision trees. Combining with the predicted DWI data and SANDI model, we separated the extracellular conductivity from the high-frequency conductivity using the decomposed micro-structural diffusion parameters. To verify the proposed method, we conducted human experiments to verify that the proposed method recovered the low-frequency electrical properties using a routine protocol sequence of MRI scan.

High-frequency conductivity at Larmor frequency using B1 phase map
The electrical conductivity of biological tissue as a function of frequency is complicated by the anisotropic nature of tissue, non-homogeneous natures in the extracellular and intracellular compartments, and randomly distributed cells sizes. The high-frequency conductivity is dominantly isotropic because the electrical current flow tends to pass through the cell membrane.
Since the reaction-diffusion equation is sensitive to the measured noise, to stabilize the formula (2), after adding an artificial diffusion term, the equation (2) leads to where c is a constant diffusion coefficient.

Detection of the micro-structures of biological tissues based on intracellular and extracellular compartments
For the model-based micro-structure imaging based on the tissue micro architecture in the brain, the model parameters for the intracellular and extracellular compartments are associated with specific tissue micro-structure characteristics from Mb-DWI data [7][8][9]. The signal intensity S b by applying a diffusion encoding gradient is given by where S 0 is the signal obtained without diffusion gradient and b denotes the diffusionweighting factor. To distinguish the diffusion signals from ECS and the cell bodies of any brain cell type (collectively named soma) [9], SANDI model proposes the following compartment model of brain tissue micro-structure: where f ic and f ec are the intracellular and extracellular volume fractions, f ic + f ec = 1 ; f in and f is are the neurite and soma relative volume fractions in the intracellular compartment, f in + f is = 1 ; A in and A is are the normalized signals for restricted diffusion within neurites and soma, respectively, and A ec is the normalized signal of the extracellular compartment. To investigate the complicated micro-structural model (5), some assumptions are applied to each parameter. The diffusion of water molecules associated with the extracellular compartment is modeled as isotropic Gaussian diffusion: The diffusion signals A in from neurites are assumed as a collection of sticks (long thin cylinders). The direction-averaged DWI signal A in is computed as [11] where erf is the error function. The signal contribution, A is , for the inra-soma compartment is computed from the Gaussian phase distribution (GPD) approximation of the tail: where D is is the bulk diffusivity of water in soma, α m is the m-th root of the Bessel equation The total unknown parameters to be determined are f in , f ec , D in , D is , D ec , and r s . To avoid the ill-posedness to determine the six unknowns from the smoothly decayed exponential curves, typically the parameter of D is is fixed as 2 ×10 −3 mm 2 /sec [9].

Model-parameter estimation using random forest regression
Random forest is a popular machine learning algorithm using the bootstrap aggregating (bagging) with a tree model as the base model. Random forest regression is an ensemble supervised learning method that combines bagging decision trees with random subset sampling of the predictors for constructing each node split [22]. To perform the random forest regression, each individual decision tree produces a prediction individually and then predictions of all decision trees are combined to generate a prediction of the ensemble. The number of trees can be adapted to find the desired trade-off between accuracy and computational efficiency of the detection process.
The five parameters f in , f ec , D in , D ec , and r s in (5) are estimated by random forest regression using the scikit-learn python toolkit [48] as in [9]. We generated 13 5 synthetic signals using the model (5) with 13 5 combinations of the five parameters chosen uniformly distributed within the intervals: f in , f ec ∈ [0.01, 0.99] , D in , D ec ∈ [0.1, 3] × 10 −3 mm 2 /sec, and r s ∈ [3,20] µ m. We added rician-distributed noise to the synthetic signals. We split the synthetic signals into random train and test subsets with test size of 20%. Hyperparameter selection in random forest regression was performed using a grid search for the number of decision trees (150, 180, 200, 230, 250, 280, 300 and 330 trees). The remaining (6) A ec (b) = e −bD ec