# On a Fast Solution Strategy for a Surface-Wire Integral Formulation of the Anisotropic Forward Problem in Electroencephalography

Carlo Baronio, Giulio Cosentino, Paolo Ricci, Clément Henry, Maxime Y. Monin, Adrien Merlini, Francesco P. Andriulli

► To cite this version:

Carlo Baronio, Giulio Cosentino, Paolo Ricci, Clément Henry, Maxime Y. Monin, et al. On a Fast Solution Strategy for a Surface-Wire Integral Formulation of the Anisotropic Forward Problem in Electroencephalography. IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting (AP-S/URSI 2022), Jul 2022, Denver, CO, United States. 10.1109/AP-S/USNC-URSI47032.2022.9886398. hal-03762860

HAL Id: hal-03762860

https://hal.science/hal-03762860v1

Submitted on 24 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On a Fast Solution Strategy for a Surface-Wire Integral Formulation of the Anisotropic Forward Problem in Electroencephalography

Carlo Baronio<sup>(1)(2)</sup>, Giulio Cosentino<sup>(1)(2)</sup>, Paolo Ricci<sup>(2)</sup>, Clément Henry<sup>(3)</sup>, Maxime Y. Monin<sup>(2)</sup>, Adrien Merlini<sup>(3)</sup>, and Francesco P. Andriulli<sup>(2)</sup>

(1) Galileo Ferraris High School, Turin, Italy

(2) Politecnico di Torino, Turin, Italy

(3) IMT Atlantique, Brest, France

*Abstract*—This work focuses on a quasi-linear-in-complexity strategy for a hybrid surface-wire integral equation solver for the electroencephalography forward problem. The scheme exploits a block diagonally dominant structure of the wire self block—that models the neuronal fibers self interactions—and of the surface self block—modeling interface potentials. This structure leads to two Neumann iteration schemes further accelerated with adaptive integral methods. The resulting algorithm is linear up to logarithmic factors. Numerical results confirm the performance of the method in biomedically relevant scenarios.

## I. INTRODUCTION

Several neuro-pathologies require precise functional brain imaging as part of their diagnostic or therapeutic protocols (see [1] and references therein). Among non-invasive strategies, high resolution electroencephalography (HR-EEG), that images the electric activity of the brain from scalp potentials, is widely used. In HR-EEG the volume currents are retrieved from the measurements of the electric potentials on the scalp by solving the EEG inverse problem. Solving this inverse problem requires multiple solutions of the EEG forward problem (FP) in which the surface potential generated by a known current configuration is computed. Boundary element methods (BEMs) are very popular in the biomedical community to model the FP and a recent hybrid formulation [1] has introduced the possibility of modeling white matter anisotropies by coupling surface BEM with an integral equation for partially conducting wires. In this work we present a fast matrix-vector multiplication algorithm for this hybrid formulation which, by exploiting the block diagonal dominance structure (induced by the presence of neuronal fibers in the model) and coupling this matrix structure with adaptive integral methods, obtains a scheme with $O(N \log N)$ complexity in the $N$ degrees of freedom. Theoretical and algorithmic considerations will be complemented by numerical experiments showing the impact of the formulation on medical scenarios.

## II. BACKGROUND AND NOTATION

Consider a sequence of nested compartments $\Omega_i$, $i = 1, \dots, C$ modeling the different layers of the head medium, characterized by homogeneous and isotropic conductivities $\sigma_i$. The boundary of each compartment is denoted by $\Gamma_i$. Following the strategy in [1], the inhomogeneity and anisotropy

of the head medium is modeled by populating the white matter with wires of finite anisotropic conductivity contrast $\bar{\bar{\chi}}(\mathbf{r}) = (\sigma_{i_w} \mathbf{I} - \bar{\sigma}(\mathbf{r})) \bar{\bar{\sigma}}^{-1}(\mathbf{r})$ with respect to the background conductivity $\sigma_{i_w}$ of the white matter's compartment. In this setting, the EEG FP consists in finding the electric potential $\phi(\mathbf{r})$ on the scalp surface $\Gamma_C$ generated by a primary current $J_p(\mathbf{r})$. To do so, the surface $\xi$ and wire $J_{eq}$ unknowns (see [1] for their physical definition) are expanded with discrete basis functions, i.e. $\xi \approx \sum_{i=1}^{N_s} (\alpha_s)_i p_i$ and $J_{eq} \approx \sum_{i=1}^{N_w} (\alpha_w)_i h_i$ where $p_i$ and $h_i$ are the 2D and 1D linear Lagrange interpolants, respectively. Following a Galerkin approach leads to a linear system of $N = N_s + N_w$ unknowns

$$Z\alpha = \begin{bmatrix} v_s \\ v_w \end{bmatrix} \quad \text{with} \quad Z = \begin{bmatrix} -\mathbf{G}_{ss} + \mathbf{D}_{ss}^* & -\mathbf{D}_{sw}^* \\ -\mathbf{S}_{ws} & \mathbf{G}_{ww} + \mathbf{S}_{ww} \end{bmatrix}, \quad (1)$$

and where $\alpha = [\alpha_s \alpha_w]^\mathrm{T}$, $(\mathbf{G}_{ss})_{ij} = \frac{\sigma_{i_n} + \sigma_{i_n+1}}{2(\sigma_{i_n+1} - \sigma_{i_n})} \langle p_i, p_j \rangle_\Gamma$, $(\mathbf{D}_{ss}^*)_{ij} = \langle p_i, \mathbf{D}^* p_j \rangle_\Gamma$, $(\mathbf{D}_{sw}^*)_{ij} = \frac{1}{\sigma_{i_n}} \langle p_i, \mathbf{D}_v^* \bar{\bar{\chi}} h_j \rangle_\Gamma$, $(\mathbf{S}_{ws})_{ij} = \langle h_i, \nabla S p_j \rangle_\Omega$, $(\mathbf{G}_{ww})_{ij} = \langle h_i, \bar{\bar{\chi}}^{-1} h_j \rangle_\Omega$, $(\mathbf{S}_{ww})_{ij} = \frac{1}{\sigma_{i_n}} \langle h_i, \nabla S_v \bar{\bar{\chi}} h_j \rangle_\Omega$, $(v_s)_i = -\frac{1}{\sigma_s} \langle p_i, \mathbf{D}_v^* J_p \rangle_\Gamma$, and $(v_w)_i = -\frac{1}{\sigma_p} \langle h_i, \nabla S_v J_p \rangle_\Omega$, with $(Sf)(\mathbf{r}) = \int_S G(\mathbf{r}, \mathbf{r}') f(\mathbf{r}') dS', (\mathbf{D}^*f)(\mathbf{r}) = \int_S \hat{n} \cdot \nabla G(\mathbf{r}, \mathbf{r}') f(\mathbf{r}') dS'$, $(S_v f)(\mathbf{r}) = \int_V G(\mathbf{r}, \mathbf{r}') \nabla' \cdot f(\mathbf{r}') dV'$, and $(\mathbf{D}_v^* f)(\mathbf{r}) = \int_V \hat{n} \cdot \nabla G(\mathbf{r}, \mathbf{r}') \nabla' \cdot f(\mathbf{r}') dV'$. Above $\hat{n}$ denotes the unit normal vector pointing outwards $\Gamma_i$ and $G(\mathbf{r}, \mathbf{r}') = \frac{1}{4\pi \|\mathbf{r} - \mathbf{r}'\|}$ is the static Green function. Once (1) is solved, $S$ and $S_v$ can be applied to $\alpha$ to get the potential $\phi(\mathbf{r})$ on $\Gamma_C$.

## III. A FAST SOLUTION STRATEGY

With respect to a standard integral formulation for isotropic media, corresponding to the left diagonal block in (1), the inclusion of the white matter anisotropy adds a new wire-diagonal block and two coupling blocks in the system. First, the new scheme aims at decoupling the surface and wire solution via block diagonal inversion and Neumann series solution of the remainder: after separating diagonal and off-diagonal blocks: $Z = Z_{\text{self}} + Z_{\text{coupl}}$ with $Z_{\text{self}} = [\mathbf{Z}_{ss}, \mathbf{0}, \mathbf{0}, \mathbf{Z}_{ww}]$ and $Z_{\text{coupl}} = [\mathbf{0}, \mathbf{Z}_{sw}; \mathbf{Z}_{ws}, \mathbf{0}]$, we solve (1) as $(\mathbf{I} + \mathbf{Z}_{\text{self}}^{-1} \mathbf{Z}_{\text{coupl}})\alpha = \mathbf{Z}_{\text{self}}^{-1}v$ via a Neumann series approach enabled by the block diagonal dominance, in cases of practical relevance, of the original matrix (i.e. for the spectral radius $\rho_Z = \rho(\mathbf{Z}_{\text{self}}^{-1} \mathbf{Z}_{\text{coupl}}) < 1$). Thus we have
<table><caption>TABLE I
PERFORMANCE COMPARISON</caption><thead><tr><th></th><th>Setup time</th><th>Storage</th><th>Time per RHS</th></tr></thead><tbody><tr><td>Standard Iterative Solution</td><td>24 657 s</td><td>31.2 GB</td><td>89.15 s</td></tr><tr><td>This work</td><td>2498 s</td><td>0.4 GB</td><td>26.42 s</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Fig. 1. Solution time vs number of unknowns.</figcaption></figure>

$\alpha_s^{(k+1)} = \mathbf{Z}_{ss}^{-1}(v_s - \mathbf{Z}_{sw}\alpha_w^{(k)})$ and $\alpha_w^{(k+1)} = \mathbf{Z}_{ww}^{-1}(v_w - \mathbf{Z}_{ws}\alpha_s^{(k)})$ whose complexity reduces to the one of the two inversions and of the multiplication of the coupling terms. The multiplication of the coupling terms can be done efficiently if a fast matrix vector product algorithm is available. We have opted for an adaptive integral method (AIM) [2]. In other words, all kernel interactions in $\mathbf{D}_{ss}^*$, $\mathbf{S}_{ww}$, $\mathbf{S}_{ws}$ and $\mathbf{D}_{sw}^*$ between all Gaussian quadrature points are interpolated on the same Cartesian grid with a number of nodes proportional to the number of unknowns $N$ and handled via FFT in $\mathcal{O}(N \log N)$ complexity. As is standard in AIM [2], a near field precorrection is required for all kernels: a generic $\mathbf{D}^*$ and $\mathbf{S}$ (for surface, wire, and off diagonal couplings) is written as $\mathbf{D}^* = \mathbf{D}_{near}^* - \tilde{\mathbf{D}}_{near}^* + \mathbf{\Phi}_p \mathbf{\Lambda}^\top \mathbf{g}_D \mathbf{\Lambda} \mathbf{\Phi}_f$ and $\mathbf{S} = \mathbf{S}_{near} - \tilde{\mathbf{S}}_{near} + \mathbf{\Phi}_p \mathbf{\Lambda}^\top \mathbf{g}_S \mathbf{\Lambda} \mathbf{\Phi}_f$, where $\mathbf{D}_{near}^*$ and $\mathbf{S}_{near}$ are the uncompressed near fields, $\tilde{\mathbf{D}}_{near}^*$ and $\tilde{\mathbf{S}}_{near}$ are the FFT precorrections, $\mathbf{\Lambda}$ is unique for every product and interpolates the quadrature points, and $\mathbf{\Phi}_p$ and $\mathbf{\Phi}_f$ map quadrature points to basis functions; all these matrices are sparse. The FFT is applied to the Toeplitz matrices $\mathbf{g}_S$ and $\mathbf{g}_D$ that, because of the translation invariance of all Green functions involved, require $\mathcal{O}(N)$ memory storage. Since the double layer kernel is $\hat{\mathbf{n}} \cdot \nabla G(\mathbf{r}, \mathbf{r}') = \hat{\mathbf{n}} \cdot \frac{\mathbf{r}' - \mathbf{r}}{4\pi \|\mathbf{r}' - \mathbf{r}\|^3}$, the product of $\mathbf{g}_D$ with a vector is split into three scalar components.

Since the $\mathbf{Z}_{ss}$ block corresponds to the classical homogeneous multilayer BEM formulation, once a fast matrix vector product algorithm is available, it can be inverted iteratively with standard techniques (see [1] and references therein). Regarding $\mathbf{Z}_{ww}$, the near field kernel interactions are extracted with an octree and the resulting sparse matrix $\mathbf{N}$ is used as a preconditioner of the linear system $\mathbf{Z}_{ww}\mathbf{x} = \mathbf{b}$. The near field dominance of $\mathbf{Z}_{ww}$—due to the electric current flowing along the fibers, i.e. $\rho_w = \rho(\mathbf{N}^{-1}(\mathbf{Z}_{ww} - \mathbf{N})) < 1$— enables a second usage of a Neumann series from which $\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \mathbf{N}^{-1}(\mathbf{b} - \mathbf{Z}_{ww}\mathbf{x}^{(k)})$. A sparse solver is used to invert $\mathbf{N}$ in $\mathcal{O}(N)$ time complexity and the multiplication of $\mathbf{Z}_{ww}$ is done in $\mathcal{O}(N \log N)$ with the AIM.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Brain fibers current density on the MRI head model.</figcaption></figure>

## IV. NUMERICAL RESULTS

The favorable complexity scaling of the proposed scheme has been verified on a set of canonical geometries composed of spherical surfaces and orthogonal brain fibers. The total timings are reported in Fig. 1 and clearly confirm that the scheme we propose is, up to logarithmic factors, linear in complexity. The relevance of our fast solution strategy for real case scenarios has been studied on a realistic head model obtained with magnetic resonance imaging (MRI) data that includes white matter neuronal fibers with a tangential anisotropic conductivity of $1.3 \text{ S m}^{-1}$ and four layers (gray matter, cerebrospinal fluid, skull, scalp) with conductivities $0.13 \text{ S m}^{-1}$, $1.79 \text{ S m}^{-1}$, $0.01 \text{ S m}^{-1}$, and $0.43 \text{ S m}^{-1}$ respectively. The obtained current on the neuronal fibers is shown in Fig. 2. For this problem the radius of the fibers is chosen to match a total volume of $450 \text{ mm}^3$. The total number of unknowns is 63 922. The two spectral radii are $\rho_Z = 0.439$ and $\rho_w = 0.799$, both less than one, thus allowing the Neumann strategy. For this experiment we have compared in Table I the method proposed in this work with the uncompressed solution. In both cases the tolerance iterative schemes has been set to $10^{-3}$ and the results show the advantage of the new scheme.

## ACKNOWLEDGMENT AND CONTRIBUTIONS

This work and the associated Early-Research Program for talented high school students was supported in part by the European Research Council (ERC) through the European Union's Horizon 2020 Research and Innovation Programme under Grant 724846 (Project 321) and in part by the H2020 Marie Skłodowska-Curie Actions ITN-EID project COMPETE GA No 955476. The first two authors are listed in alphabetic order and have contributed equally to this work.

## REFERENCES

[1] M. Y. Monin, L. Rahmouni, A. Merlini and F. P. Andriulli, "A hybrid volume-surface-wire integral equation for the anisotropic forward problem in electroencephalography," IEEE J. Electromagn. RF Microw. Med. Biol., vol. 4, no. 4, pp. 286-293, Dec. 2020.

[2] E. Bleszynski, M. Bleszynski and T. Jaroszewicz, "AIM: Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems," in Radio Science, vol. 31, no. 5, pp. 1225-1251, Sept.-Oct. 1996.