HArtMuT—modeling eye and muscle contributors in neuroelectric imaging

Nils Harmening; Marius Klug; Klaus Gramann; Daniel Miklody

doi:10.1088/1741-2552/aca8ce

1. Introduction

Magnetoencephalography (MEG) and electroencephalography (EEG) are methods to non-invasively measure the electrical activity of populations of cortical neurons in humans and other species. EEG provides high temporal but only limited spatial resolution, because of the sensors being placed on the scalp outside the low-conductive skull. Due to volume conduction, only a spatially smeared linear mixture of all cortical sources, that are active in parallel, arrives at the sensor level (Haufe et al 2014). MEG has better spatial resolution than EEG, since the magnetic fields are less distorted by the mixture of different conductive tissues (brain, cerebrospinal fluid (CSF), skull, scalp), but suffers from low signal-to-noise ratio in particular for deeper sources. In addition, MEG is sensitive primarily to tangential-oriented sources and is more demanding in its recording procedure.

Both MEG and EEG, referred to as M/EEG in the following, record volume conducted signals composed of cortical sources that are active in parallel and mixing in a linear fashion. Dissociating the different sources can be solved by source separation techniques (Pearson 1901, Makeig et al 1996a, Blankertz et al 2011), that introduce a number of assumptions to solve the ill-posed inverse problem (Sarvas 1987). Mathematically dissociating cortical sources in M/EEG signals is further complicated by other sources contributing to the recorded signal, like physiological non-brain sources (eye movement, muscle activity) and mechanical (cable sway in EEG) or electrical sources (noise stemming from stimulation monitors and other devices).

While mechanical and electrical sources often show invariant patterns, that can be easily dissociated, eye movement and muscle activity of facial and neck muscles are often less easy to identify (Winkler et al 2011a). Importantly, these classes of non-brain sources contribute significantly to the signal as they are located outside the low-conductive skull and show higher power compared to the minuscule EEG signal. Signals from eye movements and blinks can be orders of magnitude larger than brain-generated electrical potentials (Joyce et al 2004). Therefore, activity stemming from eye movements, the closure of the lid, and its necessary muscular activity together with neck and face muscle activity are the most common types of artifacts in EEG recordings (Croft and Barry 2000).

This kind of non-brain activity is traditionally considered to be artifactual, even though eye movement and facial mimicry significantly contribute to or are a result of cognitive and affective processes (König et al 2016). As a consequence, a number of artifact rejection methods have evolved to identify and remove such kind of activity from the recorded signal. Rejecting contaminated trials, however, causes substantial data loss, and restricting eye movements and blinks might limit the experimental design and could impact the cognitive processes under investigation. Thus, newer algorithms attempt to remove artifact specific features from the signal without deleting contaminated time periods to preserve as much data as possible (Jung et al 2000, Dimigen and Ehinger 2021).

Among these algorithms are blind source separation techniques (Urigüen and Zapirain 2015) that decompose the signal mostly based on its spatio-temporal statistics into different sources after which neural sources are identified and separated. The most commonly used approach is independent component analysis (ICA) which separates the sources based on the principle of maximal statistical independence (Makeig et al 1996a). The assumption that neural and artifactual sources are independent implies the possibility of sorting them into sources of different origins. This can be done manually by an expert or by applying automatic routines (Winkler et al 2011b, Pion-Tonachini et al 2019).

Eventually, the sources that have been identified as originating in the brain, are further localized using source localization approaches (Oostendorp and Van Oosterom 1989, Mosher et al 1992, Oostenveld et al 2011). Finding a source in general requires computing forward solutions by a head model, so-called leadfields. Their difference from the original data is subject to minimization, where the exact source position is the optimization parameter. The same approach to localize brain sources could thus be used to also localize other physiological sources. However, to our knowledge, no head model actually provides leadfields that incorporate the eyes and face and neck muscles in the forward solution.

Furthermore, the simulation of different signal sources including those of eyes and muscles is useful for educational and validation purposes, as it has been done for neural origins (Krol et al 2018). This also involves a biophysical forward model, that mimics the field propagation within the head, and has so far been done almost exclusively for neural sources in existing models.

In order to improve the source localization accuracy, to offer a better guide for users to distinguish brain from non-brain sources and to provide a model, which is suitable for eye and muscle simulations, we developed the Head Artifact Model using Tripoles (HArtMuT), incorporating eyes, facial muscles and muscles in the dorsal neck region into EEG leadfields.

In this study, exact source positions and orientations were extracted from the open-source, high-resolution multimodal imaging-based detailed anatomical (MIDA) model (Iacono et al 2015). Dipolar, tripolar and symmetric dipolar structures for realistic artifact modeling were physiologically motivated and placed onto these positions. A four-shell boundary element method (BEM) head model of the anatomy of Colin27 (Holmes et al 1998) with neck extended scalp mesh was used to evaluate the artifact model on both simulated highly detailed finite element method (HD-FEM) data and ICA patterns from EEG recordings.

The final HArtMuT is publicly available as BEM and FEM model for different commonly used standard head anatomies and is ready to be implemented into standard approaches of all kinds.

2. Source modeling

In the lower frequencies ( ${\lt}1\,$ kHz), the basis of common electrical models of volume conduction in the human head with conductivity σ is the quasi-electrostatic assumption (Häamäläinen et al 1993), which leads to a Poisson equation that has to be solved in order to model the field distributions (see e.g. Sarvas 1987):

$\begin{equation} \nabla \cdot \left(\sigma \nabla \phi\right)) = \nabla \cdot \vec{J}_\mathrm{P}. \end{equation} \tag{ 1 }$

$\vec{J}_\mathrm{P}$ are the sources of primary currents that lead to the volume conduction $\sigma \nabla \phi$ , hence the electrical potential φ. $\vec{J}_\mathrm{P}$ are usually modeled to be evoked by a superposition of point sources like single dipoles called the equivalent current dipole (ECD). The sources produce a local current flow which leads to a return current in all connected conductive tissues. The infinite homogeneous solutions of equation (1) for a given geometry with the ECD as $J_\mathrm{P}$ are commonly used to simulate the electrical potential of a firing population of neurons. These solutions are unique and can be numerically determined. The most common numerical methods nowadays are BEM, FEM, and finite differences method (FDM). In general, FEM/FDMs can incorporate more detail in anatomy and conductivity, such as anisotropy or certain local inhomogeneities, but are computationally much more expensive.

In the following, the physiology of the different source types, their electrical fields, and corresponding source models are examined in detail, demonstrating that the current densities of dipolar ( $\vec{J_\mathrm{d}}$ ) and tripolar ( $\vec{J_\mathrm{t}}$ ) source models are reasonable approximations for $\vec{J_\mathrm{P}}$ , when modeling EEG contributors.

2.1. Cortical sources

The main neuronal activity is believed to be based on membrane potentials and electrical currents through ions. Whether being action potentials (APs) along the axons or postsynaptic potentials (PSPs), they produce a current flow through the membrane, that leads to a secondary return current flow through the neuron's surrounding, which can be measured outside the head (Buzsáki et al 2012). In particular scalp EEG is assumed to be primarily elicited by the mean field produced by PSPs of larger assemblies of neocortex layer pyramidal cells. The ECD is a widely accepted approximation for the electric field of such a population of parallelly aligned and synchronously firing neurons. For a detailed derivation of the ECD from the contribution of principal neocortical neurons to the EEG, please refer to Murakami and Okada (2006).

The far-field of an ECD is described by an analytic expression, which is then used in combination with the head model to describe the resulting volume conduction within the whole head. In Cartesian coordinates, this electrical far-field potential of a dipolar current source in an infinite homogeneous conductive medium is

$\begin{equation} \phi_\mathrm{d}(\vec{r},\vec{r}_q)=\frac{1}{4\pi\sigma}\frac{\vec{Q}\cdot\left(\vec{r}-\vec{r}_q\right)}{\left|\vec{r}-\vec{r}_q\right|^3}\hspace{5mm}. \end{equation} \tag{ 2 }$

And the current density is calculated as

$\begin{equation} \vec{J}_\mathrm{d}(\vec{r},\vec{r}_q)=\sigma \nabla \phi_\mathrm{d}(\vec{r},\vec{r}_q)=\frac{1}{4 \pi } \left(3 \frac{\left(\vec{r}-\vec{r}_q\right) \cdot \vec{Q}}{\left|\vec{r}-\vec{r}_q\right|^5} \cdot \left(\vec{r}-\vec{r}_q\right)-\frac{\vec{Q}}{\left|\vec{r}-\vec{r}_q\right|^3}\right)\hspace{5mm}, \end{equation} \tag{ 3 }$

where $\vec{Q}$ is the local current density produced by the source, $\vec{r}$ the observed point, and $\vec{r}_q$ the position of the source.

2.2. Ocular sources

Ocular artifacts in the EEG are caused by eye movements or blinks and are associated with higher amplitudes and lower frequencies than those of cortical EEG signals. The underlying sources are mainly permanent charge distributions in the eyeballs, but additionally include the eye muscles, that produce an electrical field, when active, and the closure of the eyelids, which modifies the volume conduction. The eyes are frequently the components of strongest amplitudes and can be easily identified by their specific frontal patterns. Eye movements of both eyes are generally linked in healthy humans which leads to a joint field produced by both eyes equally. Since the potentials are correlated among both eyes, single eye sources cannot be separated by ICA/principal component analysis (PCA). Hence, this symmetry has to be incorporated into the model, which is taken care of by investigating both single and symmetric eye models.

2.2.1. Physiology

The eyes are electrically charged at different structures as depicted in figure 1. Firstly, the photoreceptors in the retina maintain a standing negative resting potential in order to turn incoming light into an excitation. In the pigment epithelium, this results in an electrical signal, which then propagates through the optic nerve to the visual cortex. Contrary to this strong negative charge distribution, fewer positive potassium and sodium ions are located in the outer segment of the photoreceptor cells, forming the overall membrane potential of around $-40\,$ mV. When stimulated, the photoreceptors become even more hyperpolarized. Opposing the retina's negativity is the overall positivity of the cornea (Iwasaki et al 2005). In its function to protect the frontal eye from physical and chemical agents, it transports positively charged sodium ions inward and negatively charged chloride ions outward. To this end, a permanent gradient is maintained, such that the surface, the corneal epithelium, is constantly positively charged, opposed by less positivity inside. During eye movements or blinks, the described charge distributions move and their electrical field changes. Since it is these changes, that are recorded by EEG (Berg and Scherg 1991), an equivalent ocular model is expected to describe the charge distribution difference during movement rather than the static charge distributions. Additionally, the eyeball rotation is driven by small muscles around the eyes, that produce an electrical field when active. The closure of the eyelid is thought to additionally change the local field by altering the geometry of the volume conductor (Plöchl et al 2012).

2.2.2. Previous approaches to model eye activity

Some studies already investigated modeling the contribution of eye activity to EEG signals, where the strong corneo-retinal potential was mostly approximated by one equivalent ocular dipole with its negative pole posteriorly (directed towards the retina) and its positive pole anteriorly (directed towards the cornea). Using electrooculography (EOG) recorded signals during vertical and horizontal eye movements and eye blinks, Berg and Scherg (1991) observed, that 'a reasonable fit was only obtained if the equivalent dipoles were allowed to take up different locations and orientations depending on the type of eye activity'. A single dipole was determined for each movement across subjects. However, the difference dipole of a single eye movement was assumed to be approximately stationary due to asking the subjects to move their eyes always by the same fixed angle. Therefore, one might conjecture, that even more equivalent dipoles for the same movement direction but different angles are needed. One specific dipole per eye was also determined to approximate the EOG signal during eye blinks. This dipole bears analogy to that of vertical eye movements, which reflects Bell's phenomenon: 'Spontaneous blinks produced small eye movements directed down and inward, whereas slow or forced blinks were associated with delayed upward eye rotations', as stated by Iwasaki et al (2005). They undertook a detailed EOG analysis of eye and lid movements, reasoning that 'EOG signals during vertical eye- and lid-movements are greatly influenced by the eyelids', since a closed eyelid conducts the ocular dipole field to the frontopolar scalp region within the scalp tissue.

2.2.3. Conclusions regarding modeling eye activity

In order to mirror the actual biological charge distribution differences correctly, dipolar sources are placed into retina and cornea, respectively. Accompanied by eyeball rotations, activity of the extraocular muscles is expected, as well as Muller's muscle activity due to eyelid closures during blinks. This muscular activity is included in the muscular model.

Additionally, we hypothesize that both eyes move simultaneously and synchronously during regular usage in healthy subjects, which leads to one electric field and therefore patterns stemming from the superposition of two sources—one in each eye. To this end, we included symmetric dipolar fields into HArtMuT, that consist of a summation of leadfields of left and right eye positions with identical orientation vectors. The performance of the different ocular models (single and eye-symmetric dipoles together with tripolar sources in the eye muscles) was investigated. The same dipolar model as used for the brain (equations (2) and (3)) and the same tripolar model as used for the muscles (equations (6) and (7)) are obeyed, respectively. The modeling of the eyelid closure is incorporated by the inclusion of the eyeball meshes into the scalp mesh, which in turn, however, excludes the open eyelid state. Moreover, this approach lacks modeling eyeballs (and muscles) as their own tissue(s) with different conductivity than the remaining scalp.

2.3. Muscular sources

Muscular activity is usually measured using bipolar electromyography (EMG) recordings. Several muscles are located in the face, on the human skull surface, and in the posterior and anterior neck regions. The relative strength of their activity patterns on the scalp surface is related to their location outside of the poorly conductive skull. This is why EMG activity is often observed in EEG recordings, although participants are restrained in their movements in traditional M/EEG protocols. Supracranial and facial muscles contribute to the EEG signal mainly in the higher frequency range starting at around 20 Hz. However, newer mobile EEG and mobile brain/body imaging (MoBI) studies (Jungnickel and Gramann 2016, Gramann et al 2021) demonstrated a strong contribution of neck muscles in moving participants with a contribution of lower frequencies in addition to the known higher frequency range from other muscles. Thus, it is even more important to dissociate this non-brain activity from the activity of interest, also in lower frequencies. Most of these muscle artifacts are larger in amplitude than the neuro-electric activity measured at the scalp. Their spatial scalp patterns differ from those of neural origin, since their sources are located directly in the scalp in close proximity to the EEG cap, and their field is thus more local. Additionally, due to their physiological properties, they have a mix of dipolar and tripolar field distributions, which has an impact on the patterns as well.

2.3.1. Physiology

Skeletal muscles are composed of bundles of muscle fibers that are each connected with the corresponding motor neuron of the nervous system at neuromuscular junctions (NMJs), so-called motor end plates. Many nerve impulses for muscle contraction in form of motor neuron's APs reach the end plates and at some point lead to a depolarization of the resting membrane potential of approximately $-70\,$ mV. If the threshold potential (often between around −55 and $-50\,$ mV) is exceeded, the contraction of all sarcomeres of the muscle fiber is triggered (Merletti and Farina 2016). This depolarization and the subsequent re-polarization phase, called hyperpolarization, is followed by the afterhyperpolarization (AHP), in which the membrane potential falls below the resting potential until it returns to its normal resting voltage (Merletti and Farina 2016).

The muscle fibers innervated by a single motor neuron are known as the muscle unit. Those muscle units innervated by the same motor neuron form together with its motor neuron a single motor unit, such that a motor unit action potential (MUAP) is produced by the summation of the accompanying muscle unit's end plate potentials. The generated MUAPs propagate in opposite directions from the muscle end plate to the tendons, that lead to a recordable EMG signal (and as unwanted EEG artifacts) at the scalp surface, which is originating from multiple distributed source locations, and thus rather a mean field signal. However, at the end of the single muscle fibers the propagating MUAPs progressively overlap and cancel out (end-of-fiber effect), that in sum yield non-propagating signals, that mostly dominate over the propagating signals (Merletti and Muceli 2019).

2.3.2. Previous approaches to model muscle activity

While modeling ocular contributors to the EEG have already been discussed among researchers, muscular sources, to our knowledge, have not been added to any head model so far.

At the tendons the gradual extinction of the MUAPs result in a positive potential field at the side of the muscle fiber and at the side of the tendon a negative potential field, slightly attenuated by the positive AHP bump (Roeleveld et al 1997). Hence, this non-propagating part of the EMG signal is commonly modeled as dipole.

To approximate the propagating MUAP's curve some researchers used a simple dipole (Boyd et al 1978, Winter et al 1994), while others used (two) balanced tripoles pointing opposingly into the direction of the muscle fiber ends (Griep et al 1982, Roeleveld et al 1997, Merletti et al 1999, Farina and Rainoldi 1999). This is because the potential curve appears more tripolar with a skewness related to the tripolar asymmetry in propagation direction, when additionally taking the AHP into account, i.e. the negative sink after the AP peak (figure 2). As shown by Merletti and Farina (2016), the field of an equivalent tripolar source model has a similar curvature to the more realistic analytic waveform, that was originally developed by Rosenfalck (1969) and later modified by Nandedkar and Stålberg (1983).

**Figure 2.** Simulated transmembrane current per unit length (based on Nandedkar and Stålberg (1983)) and the equivalent tripole source model as approximated by Merletti and Farina (2016). The electric field of the first positive (300 nA) and the negative (420 nA) sources corresponds to the depolarization and the subsequent re-polarization phase. This is covered by an equivalent dipolar current source model. Adding the second positive charge (120 nA) a bit apart additionally incorporates the following small AHP bump. This constitutes the equivalent tripolar current source model *Tripole A*. Adapted from Merletti and Farina (2016) John Wiley & Sons. Copyright © 2016 by The Institute of Electrical and Electronics Engineers, Inc.
Download figure:
Standard image High-resolution image

Used MUAP models vary in charge distribution and locations of the three monopoles, that form the tripole. Furthermore, the line, on which the three monopoles are placed, is not necessarily straight but can also be curved. All this has effects on the appearance of their surface potentials.

2.3.2.1. Thoughts on MUAPs in the human head

In the human head, these muscular sources are found within what is usually modeled as scalp: a thinner sheet of conductive tissue bordering the non-conductive air and the much less conductive skull (the estimated conductivity ratio ranges from 1:40 to 1:80 (Goncalves et al 2003, Clerc et al 2005)). These rather sudden changes of conductivity also have an effect on the inner distributions of electrical potential and current flow. In figure 3 one can observe that depending on the orientation of the tripole relative to the measured surface potential, it can have one, two, or three focal patterns on the surface. Moreover, the central (in our case positive) peaks of tripoles decay faster with distance than dipolar peaks in amplitude. Compared to brain sources, however, the missing spatial smearing effect of the CSF and skull additionally leads to more focal patterns of these scalp sources.

**Figure 3.** Comparison of scalp potential behavior for different source model types (dipole, sym. tripole, *Tripole B*) with three different orientations (tangential, radial, 45^∘) at similar source positions. Different orientations lead to different surface potential patterns: a tripole can look similar to a dipole or a monopole depending on its orientation, however, certain characteristics like the fall-off with distance differ. The dipole (upper row) consists of a double negative sink and a double positive source. A symmetric tripole (middle row) consists of two negative sinks at the same distance in opposite directions apart from a double positive source. *Tripole B* additionally has different distances leading to a rather dipolar potential distribution. The surface at x = 0 is to a non-conductive domain for x > 0, which simulates the scalp surface touching air. The lines depict different voltage levels (equipotential lines), while the flow lines represent the current flow (the density of the lines is arbitrary).
Download figure:
Standard image High-resolution image

If tripoles are modeled similar to in Merletti and Farina (2016) (see figure 2), the resulting field is in transition between dipole and tripole as the closer and stronger pair of sources (300 nA/420 nA) dominates the field.

2.3.3. Conclusions regarding modeling muscle activity

Due to existing evidence for both dipolar and tripolar approaches, we decided to test both.

Noteworthy, little is known about EMG signals fall-off in current flow power within the muscle fibers (Kuiken et al 2001), that, strictly speaking, one would additionally need to take into account. Moreover, tripole models might play a more important role in the near-field, whereas dipoles dominate the far-field potential. Also orientation and distance of muscles to the electrodes strongly influence the appearance of the signal at the electrodes.

Since the potential form of EMG signals at different fiber positions is not simply a delayed versions of each other (Merletti and Farina 2016), we place several myogenic sources along the signal's two propagation directions from the NMJ to the tendons in the model.

Because the distance between the sources is in a similar range as the distance to the surface, a far field approximation, such as the ECDs used for neurons inside the brain, does not hold for sources outside the skull. Instead, the optimal model would be to directly model the current flow out of the muscle along its fibers using analytic solutions comparable to Nandedkar and Stålberg (1983). However, such an implementation in existing head model software would not allow the usage of most existing inverse fitting techniques, that optimize location and orientation of point sources. We, therefore, model the tripolar structure as three single monopolar sources, as described by Merletti and Farina (2016). In Cartesian coordinates, the electrical field potential of a single (monopolar) current source $\phi_\mathrm{m}$ in an infinite homogeneous conductive medium is calculated as

$\begin{equation} \phi_\mathrm{m}(\vec{r},\vec{r}_q)=\frac{I}{4\pi\sigma}\frac{1}{\left|\vec{r}-\vec{r}_q\right|}. \end{equation} \tag{ 4 }$

The corresponding current density $\vec{J}_\mathrm{m}$ is given by

$\begin{equation} \vec{J}_\mathrm{m}(\vec{r},\vec{r}_q)=\sigma \nabla \phi_\mathrm{m}(\vec{r},\vec{r}_q)=\frac{I}{4 \pi } \frac{\left(\vec{r}-\vec{r}_q\right)}{\left|\vec{r}-\vec{r}_q\right|^3}, \end{equation} \tag{ 5 }$

where I is the total current of the source. The field is radial symmetric.

In order to model the field of a tripole, three monopolar sources with different parameters (location, amplitude) are used and their fields are summed to receive the field of the according tripole:

$\begin{align} \phi_\mathrm{t}(\vec{r},\vec{r}_{q_1},\vec{r}_{q_2},\vec{r}_{q_3},\vec{a})=&\ a_1 \phi_\mathrm{m}(\vec{r},\vec{r}_{q_1}) + a_2 \phi_\mathrm{m}(\vec{r},\vec{r}_{q_2})\nonumber\\ &\ + a_3 \phi_\mathrm{m}(\vec{r},\vec{r}_{q_3}). \end{align} \tag{ 6 }$

Accordingly

$\begin{align} \vec{J_\mathrm{t}}(\vec{r},\vec{r}_{q_1},\vec{r}_{q_2},\vec{r}_{q_3},\vec{a})=&\ a_1 \vec{J}_\mathrm{m}(\vec{r},\vec{r}_{q_1}) + a_2 \vec{J}_\mathrm{m}(\vec{r},\vec{r}_{q_2})\nonumber\\ &\ + a_3 \vec{J}_\mathrm{m}(\vec{r},\vec{r}_{q_3}), \end{align} \tag{ 7 }$

where $r_{q_i}$ is the location of source/sink i and a_i the amplitude. In order to fulfill the necessity of closed current loops, the amplitudes a_i have to sum to one, i.e. are in total normalized during modeling.

The tripolar model approaches tested within this paper are the following:

Tripole A: One negative sink with amplitude $-300\,$ nA was placed at $-1\,$ mm and one with $-120\,$ nA at $2.8\,$ mm away from a central $420\,$ nA source (Merletti and Farina 2016).
Tripole B: Two negative sinks ( $-210\,$ nA) were set $-1.6\,$ mm and $3.2\,$ mm away from a central double positive source ( $420\,$ nA) (Roeleveld et al 1997).
Tripole C: Two negative sources ( $-210\,$ nA) were set $-2.1\,$ mm and $4.8\,$ mm away from a central double positive source ( $420\,$ nA) (Kuiken et al 2001).

As described above and shown in figure 3, a tripoles scalp potential is not necessarily consisting of three peaks, even in proximity to the source. When using a tripole model with a non-equal distance between the monopoles such as the used Tripoles A–C, in particular one of the sources' surface pattern peaks is weakened (the closer one). The resulting muscular Tripoles A–C are actually very similar to dipoles in the sense that the potential and current distributions are dominated by the stronger source-sink pair. This stronger source-sink pair was reached by either a larger distance between the sources or by larger source amplitudes.

3. Deriving artifact model positions and orientations

For a realistic artifact source model, detailed muscle and eye positions within the head are essential. However, common MRI scans lack the necessary resolution and segmentation techniques to deduce the relevant anatomical details of muscles and eyes. Therefore, the open-source MIDA model of the human head and neck (Iacono et al 2015), an open-source, high-resolution MRI scan, segmented manually by experts, has been utilized.

3.1. MIDA open-source Atlas

MIDA has been obtained by segmenting high-resolution T1- and T2-weighted MRI scans, magnetic resonance angiography, and diffusion tensor imaging data of one healthy 29-year-old female volunteer. It was acquired at the Institute for Biomedical Engineering of the ETH Zurich (Switzerland) and contains 153 structures, segmented at $500\,\mu$ m isotropic resolution.

3.2. Surface mesh corrections

The segmented face and neck muscle meshes as depicted in the center of figure 4 were corrected for self-intersections, unconnected nodes, and duplicated elements. Unconnected mesh parts were detected and removed. The triangulated surfaces were simplified, coarsened to an appropriate number of mesh vertices, and qualitatively enhanced (uniformity) while preserving the manifold.

3.3. Artifact source model creation

Evenly spaced grids of $1\,$ mm size were constructed within every muscle surface, building the artifact model's source positions, of that a random sample of size 3.9k was drawn³ . According to the propagation direction of the MUAPs, the source's pole orientations should follow the muscle's fiber directions for simulation. The fiber directions were approximated by using a PCA on close neighboring grid points of a muscles grid point cloud, where the closeness criteria have been varied depending on the underlying muscle shape.

For eye activity, source grids were built for six tissues in total. Firstly, the retina of the left and right eye was modeled with orientations pointing to the positive y-axis of the anterior and posterior commisura (ACPC) coordinate system (i.e. forward). Secondly, two source grids each were placed into the cornea of the left and right eye, where one of them is oriented horizontally and the other vertically. To every position, a mirrored (on the x-axis) counterpart for the calculation of the symmetric dipoles was determined.

3.4. Warping to MNI coordinates

A warping procedure for transferring the artifact source model from MIDA into the Montreal Neurological Institute (MNI) coordinate system was developed, which interpolates every point within the scalp mesh of one person into the scalp mesh of another person.

4. Validation

HArtMuTs application to the inverse problem was validated with simulated EEG scalp patterns and subsequently on a real experimental EEG data set.

4.1. Head model construction

A four-shell BEM head model of the anatomy of Colin27 (Holmes et al 1998) with neck extended scalp mesh was used for both validation steps. It consists of a cortex (gray + white matter), CSF, skull, and scalp mesh and was segmented using an automatic segmentation pipeline (Harmening and Miklody 2022). After segmentation, the scalp mesh was merged with the neck part of the NYheads scalp mesh, in order to extend the scalp mesh of Colin27 by the neck part. The final surface meshes consist of 1922 vertices and 3990 triangles each. This has proven to be a reasonable compromise between model size, hence accuracy, and computational speed during inverse fit methods with inverse nonlinear fit requiring frequent leadfield recalculation (Miklody et al 2016, Miklody 2020). As BEM solver OpenMEEG (Gramfort et al 2010) was used, and the conductivity was fixed at the following values: scalp $0.465\,$ S m⁻¹, skull $0.01\,$ S m⁻¹, CSF $1.65\,$ S m⁻¹, cortex $0.201\,$ S m⁻¹ (mean of gray and white matter) (Huang et al 2013).

4.2. Inverse fitting routines

In order to make existing dipolar fitting routines work with the tripolar approach, we used the following approximation: three tripoles were modeled around a central location in each direction of a Euclidean ACPC coordinate system. The linear combination of the three tripoles, that minimizes the error as measured by equation (8), was determined. This determined orientation and moment of the final single tripole and is the same procedure as for neural dipoles in particular with generic head models, as described in the following.

Location and orientation of the sources were adjusted in an iterative procedure in an attempt to maximize the amount of variance in the data explained by the model. One criterion for evaluation and comparison of models, therefore, was the relative residual variance (RV), the relative amount of variance in the data, which is unexplained by the model (Scherg and Berg 1991). The RV is a sum-of-squares type error function, that computes the relative deviations of an estimated pattern $\hat{x}$ from the target pattern x:

$\begin{equation} \mathrm{RV}=\frac{\sum{\left(x-\hat{x}\right)^2}}{\sum{x^2}}. \end{equation} \tag{ 8 }$

The RV has the advantage of being normalized to the average amplitude of the measured signal.

In order to fit the most optimal location and orientation to a measured scalp pattern, two different inverse fitting routines were used successively: dipfit_gridsearch and dipfit_nonlinear based on FieldTrips ft_dipolefitting (Oostendorp and Van Oosterom 1989, Oostenveld et al 2011).

dipfit_gridsearch compares the RV of the potentials at the electrodes (x) with the computed HArtMuT leadfields ( $\hat{x}$ ) for every source point in a gridded full source model (cortical + artifactual). For this gridded full source model, the sources were placed at a regular grid within the brain (cortical) and scalp (artifactual) meshes. Here, no restrictions were set to positions within the scalp, as e.g. from realistic positions derived from MIDA, in order to avoid systematic modeling errors due to incorrect segmentation or generic head models. The source position combined with an orientation, that produced the lowest RV with the scalp pattern, finally comprises the linear result to the inverse problem.
dipfit_nonlinear achieves a further optimized fit by allowing an error minimization algorithm to vary orientation and source position using the dipfit_gridsearch result as starting point.

For the purpose of labeling, the name of the closest artifact mesh is in the end assigned to a fitted source.

4.3. Simulated data

In the first validation step (simulated data), an FEM was used as a 'ground-truth' head model to validate the inverse fitting routines, the artifact positions with orientations and the four-shell BEM HArtMuT of Colin27 with extended neck as used for the inverse fit. HArtMuTs source localization results were compared to those of a standard three-shell BEM model from EEGLAB (Delorme and Makeig 2004a), which is also based on the anatomy of Colin27.

The simulated 'ground-truth' scalp patterns were calculated for all artifact positions with orientations from MIDA using only the common dipolar source model type, because no FEM solution for tripoles was available and the validity of the source localization is also given this way. As anatomy, the NYhead (Huang et al 2016) was chosen, which is a high-resolution FEM in the MNI coordinate system, that includes the neck. Its mesh and 75k cortical leadfields are publicly available. The assigned conductivities for each tissue type are scalp $0.465\,$ S m⁻¹, skull $0.01\,$ S m⁻¹, CSF $1.65\,$ S m⁻¹, gray matter $0.276\,$ S m⁻¹, white matter $0.126\,$ S m⁻¹ (Huang et al 2013). Since the published leadfields were calculated using the proprietary Abaqus FEA software suite, the cortical leadfields were recalculated using the open-source software package SimBio (SimBio Development Group 2022). Comparing both the original Abaqus Software and the SimBio leadfield simulations, an average correlation of 0.947 for single cortex source location leadfields was found. The discrepancy may be explained by the use of the different solvers and an additional linear transformation of source points and electrodes applied beforehand. This transformation was needed to bring the published FEM mesh and source and electrode positions into the same coordinate system.

Using the NYhead, FEM leadfields of a 3540 large artifact source model from MIDA were calculated at the same 129 aligned electrodes as in the dataset used in the experimental data validation. Each leadfield was considered as a single scalp pattern, to which the corresponding source position had to be localized.

4.3.1. Results of the validation on simulated data

The results of the source localization on simulated dipolar cortical and dipolar artifactual source points are plotted in figure 5. The standard EEGLAB model led to significantly higher RV and source localization error for both brain and non-brain sources than HArtMuT dipoles ( $p\leqslant0.005$ ). The average source localization error for brain and artifact sources using HArtMuT dipoles was estimated to be below 10 mm with an RV mainly below 0.01. The median RV was 0.0032 for the brain and 0.0042 for artifacts. The median source localization error was 8.4 mm for the brain and 9.6 mm for artifacts. In comparison, the standard EEGLAB model has a similar median source localization error of 9.1 mm for the brain, but a much higher value of 34.2 mm for non-brain sources. The median RVs were 0.010 for brain and 0.25 for non-brain sources.

4.3.2. Remarks on the validation with simulated data

The validation against simulated data from a HD-FEM model as 'ground-truth' was successful. This confirmed in particular HArtMuTs four-shell BEM model as a whole, the used inverse fitting routines, the neck extension procedure developed and applied to Colin27, and furthermore the artifact translation procedure between the different head shapes (MIDA, NYhead, Colin27).

4.4. Experimental data

The second validation step (experimental data) was performed by using data from an EEG experiment (Gramann et al 2021). From the experimental data, different components and their scalp patterns were extracted using ICA within EEGLAB (Delorme and Makeig 2004b). During source localization using HArtMuT, the performance of the different artifact modeling types (dipole, tripoles, symmetric dipole) was evaluated. The RVs were compared to see, which source model type best describes the experimental data patterns, and thus, which source model type overall allows for the best fit for the different types of measured signals in EEG. For comparison, the source localization results of the commonly used standard EEGLab model were computed as well.

4.4.1. Data and processing

The data of 19 healthy subjects (11 female, 8 male, aged 20–46 years, 30.25 years on average) were recorded in a MoBI study (Gramann et al 2021), focusing on heading computation in a stationary and a mobile experimental condition. The virtual environment consisted of a simple floor rendering without external landmarks. Each trial started with a pole (indicating zero heading) which was then replaced by a sphere moving around the participant at a constant distance but with different velocity profiles. Participants were instructed to follow the movement of the sphere and to rotate back after the sphere stopped moving unpredictably at eccentricities between 30^∘ and 150^∘ to the left or right. In the stationary condition, participants followed the sphere using a joystick while standing in front of a 2D monitor. In the mobile condition, they physically rotated wearing a virtual reality head-mounted display. Rotation eccentricity, speed, and direction were varied across 140 trials per movement condition (stationary vs. mobile). For the present study and modeling approaches, only the full-body rotation condition was taken into account, as we were interested in the contributions of neck muscles and eye movements to the sensor signal.

EEG data were recorded with a $1\,$ kHz sampling rate from 157 active electrodes (Brain Products, Gilching, Germany) with 129 channels in a cap with an equidistant layout (with two vertical EOG channels, both eyes) and 28 channels in a custom neckband. Individual electrode locations were recorded (Polaris Vicra, NDI, Waterloo, ON, Canada). Additional motion data was recorded but is not considered for the modeling approach in this study.

The data was processed in MATLAB (R2016b version 9.1; The MathWorks Inc. Natick, Massachusetts, USA) using custom scripts based on the EEGLAB toolbox (Delorme and Makeig 2004a, version 14.1.0). The 28 channels located in the neckband did not increase the decomposition quality of the ICA (Klug and Gramann 2021) and were thus removed. Subsequently, the data were re-sampled to 250 Hz, line noise was removed with Zapline (de Cheveigné 2020) using default settings, and breaks and pre-/post-experiment segments were removed from the data. We then used the clean_rawdata function of EEGLAB to detect and interpolate bad channels (disregarding EOG channels) using an 0.8 correlation threshold and no other measures. Subsequently, the data were re-referenced to the average of all channels except the EOG channels. We then applied a 1 Hz high-pass filter as suggested by Klug and Gramann (2021) and computed an ICA using the AMICA algorithm (Palmer et al 2011) with 2000 iterations and automatic rejection of samples (three times with three SD threshold).

In summarizing, the 19 subjects had on average 112 independent components (ICs) (with standard deviation of 10.35) depending on the number of valid channels after the preceding artifact channel rejection and rank of the data. In total 2132 ICs were extracted and to each IC a source position had to be localized using the above described dipfit routines from EEGLAB. Finally, the components were labeled using the ICLabel 'lite' classifier (Pion-Tonachini et al 2019).

4.4.2. Results

Of the total 2132 extracted components of all subjects, 69.4% were detected as muscle sources based on HArtMuTs source location, while 28.8% were brain and around 1.7% eye sources. The amplitudes of eye sources were overall strongest (median $8.8\,\mu$ V, maximum $419.0\,\mu$ V) followed by muscle sources (median $0.4\,\mu$ V, maximum $181.1\,\mu$ V) and brain sources (median $0.2\,\mu$ V, maximum $36.0\,\mu$ V).

Comparing different approaches for modeling the muscles in figure 6, Tripole B led to the lowest median and lower quartile errors, which is why it was also used as the main tripolar model of HArtMuT tripoles, although the differences are not significant.

Most of the muscular IC patterns were localized close to the combination of Muscles temporalis and temporoparietalis ( $21.6$ %), followed by MIDAs unspecific Muscle (M.) generalis class ( $18.1$ %), which is a collection of muscular tissue in the lower neck. Then it follows the M. splenius capitis ( $4.2$ %) and M. sternocleidomastoid ( $4.1$ %). All other muscles (26 in total) were below $4$ % in their occurrences.

For the eyes, using dipolar models that are symmetrically aligned (mirrored at the plane spanned by the x- and z-axis in the ACPC coordinate system) with equal amplitudes led to the lowest error using dipfit_gridsearch compared to any other combination (figure 7). Since fitting symmetric dipoles with dipfit_nonlinear was not implemented, the final results of the symmetric dipole fit remain those from dipfit_gridsearch.

**Figure 7.** The relative residual variance (RV) of all (a), cortical (b), muscular (c) and ocular (d) source reconstruction results for different source models using *dipfit_nonlinear* over 19 subjects: HArtMuT mix (dipoles for the brain, tripoles for muscle, and symmetric dipoles for eye sources) led to the lowest errors overall. For cortical sources, dipoles (blue) led to significantly ( $p\leqslant0.001$ ) lower errors than tripoles (red) (1/19) while the standard EEGLAB model (gray) produced the lowest errors (19/19). All types of HArtMuT (dipoles, tripoles, symmetric dipoles (grid symmD, dark blue)) lead to significantly lower errors for all subjects ( $p\leqslant0.05$ ) for muscles than the standard EEGLAB model. Tripoles lead to significantly lower errors for muscles than dipoles in 2/19 subjects ( $p\leqslant0.05$ ) and for two additional participants very close to significance ( $p\leqslant0.05$ ). Symmetric dipoles led to the lowest errors for eye sources while significant only for 2/19 subjects most probably due to the low number of sources per subject (0–5, 2 on average). The standard EEGLAB model significantly ( $p\leqslant0.05$ ) produced the highest RVs in all subjects for muscle and eye sources. Note, that the symmetric dipoles were only fitted using *dipfit*_gridsearch and could potentially be further improved by *dipfit*_nonlinear.
Download figure:
Standard image High-resolution image

**Figure 7.** The relative residual variance (RV) of all (a), cortical (b), muscular (c) and ocular (d) source reconstruction results for different source models using *dipfit_nonlinear* over 19 subjects: HArtMuT mix (dipoles for the brain, tripoles for muscle, and symmetric dipoles for eye sources) led to the lowest errors overall. For cortical sources, dipoles (blue) led to significantly ( $p\leqslant0.001$ ) lower errors than tripoles (red) (1/19) while the standard EEGLAB model (gray) produced the lowest errors (19/19). All types of HArtMuT (dipoles, tripoles, symmetric dipoles (grid symmD, dark blue)) lead to significantly lower errors for all subjects ( $p\leqslant0.05$ ) for muscles than the standard EEGLAB model. Tripoles lead to significantly lower errors for muscles than dipoles in 2/19 subjects ( $p\leqslant0.05$ ) and for two additional participants very close to significance ( $p\leqslant0.05$ ). Symmetric dipoles led to the lowest errors for eye sources while significant only for 2/19 subjects most probably due to the low number of sources per subject (0–5, 2 on average). The standard EEGLAB model significantly ( $p\leqslant0.05$ ) produced the highest RVs in all subjects for muscle and eye sources. Note, that the symmetric dipoles were only fitted using *dipfit*_gridsearch and could potentially be further improved by *dipfit*_nonlinear.
Download figure:
Standard image High-resolution image

For further analysis, the best models for each type of source, namely cortical dipoles, muscular tripoles (Tripole B), and ocular symmetric dipoles, were combined to one model, which will be referred to as HArtMuT mix in the following.

The performance of HArtMuT mix is demonstrated in figure 7, which compares the results of HArtMuT with those of the standard EEGLAB model for brain, muscle, and eye sources separately. For cortical sources, the standard EEGLAB model produced slightly but still significantly lower errors than HArtMuT ( $\alpha\leqslant0.01$ ). Using HArtMuT mix, the RVs were reduced from a median 0.25 to 0.11 for muscle sources and 0.18 to 0.09 for eye sources compared to using the standard EEGLAB three-shell model. With the simple dipfit_gridsearch, HArtMuT tripoles reduced the RV compared to HArtMuT dipoles for 3/19 subjects ( $\alpha\leqslant0.05$ ) only. When further optimized by dipfit_nonlinear, the differences in RV between HArtMuT dipoles and HArtMuT tripoles were not significant anymore except for one subject ( $\alpha\leqslant0.05$ ). Using dipfit_nonlinear significantly ( $\alpha\leqslant0.05$ ) reduced the errors within any source model in every subject.

4.4.3. Labeling

We also investigated the resulting classification based on source location with automatic IC classification using IClabel, an EEGLAB plugin. IClabel automatically classifies IC components into different classes based on a classifier, which was built with a crowd labeling approach (Pion-Tonachini et al 2019). It has to be noted that the classifier used in ICLabel was trained mainly on stationary data and its reliability regarding the classification of ICs from MoBI data, as used in the present study, is not well understood. The comparison of IClabel and HArtMuTs classification is summarized in table 1.

Table 1. Correspondence between IC classification based on IClabel and localization in the HArtMuT mix. Of the 15.6% brain sources labeled by IClabel, 12.9% were also found as brain sources by HArtMuT mix, which corresponds to an overlap of 83% of all brain sources. Of the 42.4% muscles even 39.6% were detected also by HArtMuT mix as muscle sources, such that the match is 93%. For the eyes, however, the 1.7% of all 3.0% that both labeled as eyes, correspond to only 40% match. IClabels line and channel noise, heart and other labels sum to in total 38.9% of all ICs, to which HArtMuT mix assigned mostly muscle (23.7%) or brain (14.8%) sources.

		IClabel
		Brain	Muscle	Eye	Heart	Line noise	Channel noise	Other	Sum
HArtMuT	Brain	12.9%	2.3%	0.1%	0.0%	0.8%	0.4%	13.6%	30.1%
mix	Muscle	2.7%	39.6%	1.7%	0.3%	0.6%	6.8%	16.0%	67.7%
	Eye	0.0%	0.5%	1.2%	0.0%	0.0%	0.0%	0.4%	2.1%
	Sum	15.6%	42.4%	3.0%	0.3%	1.4%	7.2%	30.0%

This simple approach, in which only the comparison of the scalp patterns of the forward model with IC patterns is incorporated, already leads to a high correspondence between the IC classification labels and the resulting source location based classification. Of the ICs classified as brain by IClabel, there were $83$ % equally classified as brain by HArtMuT mix and for IClabels muscles $93$ % as muscles by HArtMuT mix. The median RV of IC classification as 'other' was 0.22 opposing all remaining classes' RV of 0.096. The components classified as 'eye' by IClabel were $40$ % classified as eyes and $57$ % as muscles. It is noteworthy, that 0.9% of the by HArtMuT mix localized muscle sources were actually eye muscles and additionally 2.7% eyelid muscles (muscle orbicularis oculi), which one might also label as eye artifact.

This large correspondence between HArtMuT and IClabel is remarkable, although one has to take into account, that it is not assured that IClabel classified correctly.

4.4.4. Comparison with physiologically motivated artifact source model

Of the 1517 reconstructed artifact sources, only 324 (21.4%) were found to be actually inside an artifact mesh as derived from MIDAs MRI segmentation. However, the median distance of a source to the artifact mesh corresponding to its label was $4.5\,$ mm with a standard deviation of $52.7\,$ mm. The relatively high standard deviation is caused by a few outliers, such that 74% (54%) are closer than 10 mm (5 mm) to its mesh. Interestingly, these outliers have in most cases high RVs, meaning that they stem from rather unsuccessful fits. If one for example restricts this analysis to poor RVs, e.g. $\gt$ 0.3, which corresponds to 20% of the data, the mean distance of a source to the artifact mesh corresponding to its label is $35.7\,$ mm ( $84.8\,$ mm standard deviation).

4.4.5. Detailed view on muscular and ocular source model type performances

Looking at the exemplary source localization results and the corresponding patterns reveals several insights. First of all, the scalp patterns of muscular sources appeared as expected as common muscular artifact patterns and tripoles almost always led to lower RVs. However, it was surprising to us, that a dipole can lead to a tripolar-like pattern in the EEG like that measured in the upper row of figure 8. Investigations revealed that this was partly because of interpolation artifacts and partly because sources normally oriented and very close to the scalp surface can actually lead to similar patterns in particular in regions with heterogeneous scalp and skull thickness like the forehead. However, in most cases, a well-fitted tripole still remained the better model here. The lower row of figure 8 shows another example from the back of the head in the neck region with a slightly higher RV. Here, both dipole and tripole do not fully reproduce the shape of the positive peak but the approximation is still reasonable. What can be seen in both examples of figure 8 is, that some of the seemingly multi-poled solutions are rather dipolar but interpolation of the EEG scalp plotting from EEGLAB led to very strong (in this case negative) peaks in areas where there were no electrodes (no black dots). This means, that a comparison of BEM solver solutions in form of scalp plots, might not always be a good evaluation approach.

**Figure 8.** Two exemplary source models and their solution in source localization for muscle artifacts: the center shows the IC extracted from the data, in the left is the dipolar, and in the right the tripolar source on the scalp and their corresponding EEG patterns. Blue is the dipolar model and moment, red tripolar with Tripole B.
Download figure:
Standard image High-resolution image

In figure 9 the clear difference between symmetric and non-symmetric dipolar eye models can be investigated. As we expected, the best eye models concerning the pattern reconstruction error were linked symmetrically. Here, the simple dipfit_gridsearch outperforms all other models that use non-symmetric tripoles and dipoles. In the source locations we see that for these models the source was mostly located somewhere between the eyes, which is not realistic. Most ocular sources were found to lie in or close to the retina ( $84$ %), while the rest lie in or close to the cornea ( $16$ %). Interestingly, the commonly found pattern on the bottom in figure 9 corresponds to the straight forward looking eye position. The approach only catches static eye positions in single ICs and eye movements will be spread over linear combinations of static patterns such as the upper for upwards movements away from the lower straight forward position. A dipfit_nonlinear routine for linked symmetric positions had not been implemented within this work but is expected to further improve the results.

5. Discussion

We investigated new approaches to modeling muscle and eye contributions to EEG and used the models in source reconstruction on simulated and experimental data.

In sum, we propose our new head model HArtMuT, which includes the neck and corresponding brain, muscle, and eye tissue source positions and labels. It is available for different head anatomies and BEM/FEM models. We showed the possibility of including artifacts into common source localization routines and validated HArtMuT on simulated data. We tested the effect of modeling sources as dipoles, tripoles, and symmetric dipoles and measured the performance on experimental data in RV. While a single dipole model is most likely sufficient for reasonably good source localization of all components except the eyes, it is more precise to use a tripolar model for muscle sources. This is in accordance with the triphasic propagating MUAP and did slightly improve the estimation of source potentials for muscles for some of the subjects, but differences were not significant. If tripoles are not at hand, a dipolar model, which is commonly implemented in FEM/BEM solvers, can also be physiologically justified, as the triphasic MUAP shows a strong dipolar behavior followed by a weaker AHP bump. Moreover, the non-propagating end-of-fiber signal is supposed to lead to a rather dipolar structure and is usually excepted to dominate EMG recordings in the far field over the propagating mean field signal. However, a significant difference in RV between dipole and tripole fits depending on its distance to the electrodes could not have been observed.

For the eyes, a symmetric dipolar model led to the lowest approximation errors, which reflects the synchronous movements of the eyes. This met our expectations since ICA decompositions cannot differentiate the signals of both eyeballs.

To summarize, the best option is a mixture of the models depending on the source's tissue type. This best combination consists of dipoles for the brain, tripoles (Tripole B) for the muscles, and symmetric dipoles for the eyes. This is condensed into one model: HArtMuT mix.

Since the used tripolar models Tripoles A–C were all motivated by literature, there is the potential to further improve the tripolar sources localization results by improved architectures. Underlining this, we achieved promising preliminary results by optimizing positions and amplitudes of three single monopolar sources freely in the head by non-linear optimization. However, the success drastically depends on the initial starting point (we used the tripolar solution) and the optimization gets easily stuck at local minima. Such solutions however often lead to lower RVs and more convincing scalp patterns. We plan to further elaborate on this approach.

Not covered in this article, we also investigated extending the scalp mesh by the upper torso to for example include the possibility of localizing signals of the heart. However, neither any meaningful source had been localized in that region nor did the RV improve. This also indicates that the chosen neck extension of HArtMuT covers neck muscular activity in real EEG data sufficiently.

It remains an open question whether the case of dipfit_nonlinear optimized three-shell RV of brain ICs being significantly below the other for all subjects means that our model is not optimal here, or that the three-shell just produces lower RVs that do not directly correspond to the actual source localization error. In general, the RV tells only about fitting quality from a numeric perspective, not plausibility. The most direct interpretation of this is that the three-shell is the better-fitted model for the purpose of source localization. The four-shell is more detailed but delivers more systematic errors either concerning the geometry or in combination with a not fully appropriate source model (dipoles). In contrast, a thorough post-hoc investigation of sources from the brain revealed that many of the ICA patterns were either a mixture of multiple dipolar sources or scalp plots suggesting that a distributed source model could make more sense. In both cases, the solutions of the EEGlab standard model were mostly the better fit due to its higher spatial spread of the single sources (more spatial smearing). Additionally, the results of the validation on an FEM model proved that the four-shell model is the better model also for brain sources and we hence believe that the main reason for the three-shell performance on experimental data was rather a non-optimal IC decomposition than head model.

Another factor might be non-optimal tissue conductivities in the four-shell head model. A post-hoc analysis with preliminary results revealed that HArtMuTs RVs can be improved for the brain to a similar level by fitted conductivities. A further reason could be the missing anatomical correspondence with the subjects. One might conjecture that the three-shell is delivering lower errors because it is more generic with its constant skull thickness. However, more investigation on this is needed.

A similar argumentation can be brought up to explain the relatively large source localization error for both models HArtMuT and standard EEGLAB in the first validation step (using the FEM leadfields as ground truth). Besides not matching conductivities, in particular different solvers (FEM vs. BEM), different anatomies (NYhead vs. Colin27), and needed transformations of electrode and source positions are to name here as origins of uncertainties, that summed up in the end can lead to such an average source localization error of approx. $9\,$ mm.

All together, HArtMuT significantly improved source localization compared to the commonly used EEGLAB model in particular for muscle and eye sources. Furthermore, muscle and eye artifacts could clearly be differentiated from brain sources by their patterns via inverse fitting routines together with an appropriate model like HArtMuT.

A comparison with IClabel indicates that HArtMuT may be a promising tool not only for source localization but also for classification of ICs.

It is noteworthy, that IClabel classifications cannot be seen as the 'ground truth' in all studies, because it was mostly trained on EEG data from stationary participants (Pion-Tonachini et al 2019). As such, training a classifier mainly based on ICs from decompositions of MoBI and mobile EEG data might improve the comparison with HArtMuT. This would also allow for an integration of HArtMuTs' artifact source space as labels for training the classifier. In this sense, the comparison of HArtMuT with automated classification of ICs is likely to improve with a classifier that is better adapted to an increased number of non-brain sources with muscle and eye-related activity. In this study, IClabels classifications were introduced to serve a basic intuition for the quality of HArtMuTs results. A further study is planned to investigate the usage of HArtMuT for artifact labels in detail.

Acknowledgment

We want to thank Benjamin Blankertz for his outstanding general and logistical support of this project. N H and K G acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation): N H from GRK 2433/1, Project Number 384950143, and K G from GR 2627/8-1.

Data availability statement

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/harmening/HArtMuT.

Contributions

N H and D M implemented the segmentation and warping routines. N H deduced the artifact atlas for head models, researched mathematical approximations and physiological origins of artifacts and implemented the FEM head model. D M took care of the BEM head model including tripolar model implementation. M K implemented the first versions of the dipfit routines in EEGLAB, which were then further improved. K G had the original idea and was closely supporting the project in the beginning and in the final manuscript phase. All authors contributed to writing of the manuscript.

Conflict of interest

The authors have declared no conflicts of interest for this article.