In silico study of the effects of cerebral circulation on source localization using a dynamical anatomical atlas of the human head

Objective. This study focuses on the effects of dynamical vascular modeling on source localization errors in electroencephalography (EEG). Our aim of this in silico study is to (a) find out the effects of cerebral circulation on the accuracy of EEG source localization estimates, and (b) evaluate its relevance with respect to measurement noise and interpatient variation. Approach. We employ a four-dimensional (3D + T) statistical atlas of the electrical properties of the human head with a cerebral circulation model to generate virtual patients with different cerebral circulatory conditions for EEG source localization analysis. As source reconstruction techniques, we use the linearly constraint minimum variance (LCMV) beamformer, standardized low-resolution brain electromagnetic tomography (sLORETA), and the dipole scan (DS). Main results. Results indicate that arterial blood flow affects source localization at different depths and with varying significance. The average flow rate plays an important role in source localization performance, while the pulsatility effects are very small. In cases where a personalized model of the head is available, blood circulation mismodeling causes localization errors, especially in the deep structures of the brain where the main cerebral arteries are located. When interpatient variations are considered, the results show differences up to 15 mm for sLORETA and LCMV beamformer and 10 mm for DS in the brainstem and entorhinal cortices regions. In regions far from the main arteries vessels, the discrepancies are smaller than 3 mm. When measurement noise is added and interpatient differences are considered in a deep dipolar source, the results indicate that the effects of conductivity mismatch are detectable even for moderate measurement noise. The signal-to-noise ratio limit for sLORETA and LCMV beamformer is 15 dB, while the limit is under 30 dB for DS. Significance. Localization of the brain activity via EEG constitutes an ill-posed inverse problem, where any modeling uncertainty, e.g. a slight amount of noise in the data or material parameter discrepancies, can lead to a significant deviation of the estimated activity, especially in the deep structures of the brain. Proper modeling of the conductivity distribution is necessary in order to obtain an appropriate source localization. In this study, we show that the conductivity of the deep brain structures is particularly impacted by blood flow-induced changes in conductivity because large arteries and veins access the brain through that region.


Introduction
Electroencephalography (EEG) source localization poses an ill-posed inverse problem, the solution to which is non-unique and sensitive to any uncertainties in the measurements and/or forward modeling [1,2]. Specifically for EEG, the measurement device records surface voltages that reflect the state of the brain and the mental condition of the patient, but also records artifacts either caused by the measurement device itself, such as electronic noise and interference (e.g. faulty electrodes and parasitic capacitances), and undesired physiological sources (e.g. muscular and cardiac activities) [3]. Finite precision measurements are another source of problems in source localization, especially in deep brain structures due to the signal-to-noise ratio necessary to record brain activity in these regions.
In particular, modeling errors in the lead field matrix, i.e. in the linear operator relating the electrical activity of the brain to potentials on EEG electrodes, can result in large reconstruction errors. Electrode position, head geometry, and the electrical properties of the tissues in the human head play an important role in the quality of the reconstruction [4][5][6][7][8]. Electrical conductivity distribution in the head is a central parameter affecting the accuracy of the lead field matrix as a forward map. The human head is a complex conductor volume and must be simplified to create a numerical model. The head is usually split into segments and an average conductivity is assigned to piecewise constant distribution-one value per segment [9,10]. In addition to the geometrical complexity, the electrical properties of the tissues change from subject to subject and even within the same subject due to physiologic changes or pathological conditions.
In this work, our aim is to investigate how conductivity changes due to the pulsatile blood flow influence the accuracy of EEG source localization. The relationship between blood flow velocity and electrical properties has previously been identified, both in vitro and in vivo [11][12][13]. To approximate the dynamical blood flow, we use an interpolated atlas [14] which follows from a cylindrical Navier-Stokes equations (NSE) model of a blood vessel combined with a statistical approach to take into account the effects of the circulatory system and inter-subject variations. Incorporating such an atlas into the forward model necessitates using a volumetric forward solver.
Our hypothesis is that EEG source localization error increases if conductivity changes due to blood flow are not taken into consideration, especially in the deepest structures of the brain. The rationale behind this hypothesis is based on the combination of two factors: (a) The deep structures of the brain are the farthest regions from the scalp electrodes, making any activity in this region the hardest to measure using surface electrodes. (b) Arterial blood enters the brain cavity via the foramen magnum (vertebral arteries) or carotid canals (internal carotid arteries). These are the largest arteries in the brain, all located near the brainstem, base of the thalamus, and hypothalamus. They form the circle of Willis at the interpeduncular fossa, and from there, the middle cerebral arteries branch to the sides along the lateral sulci, the anterior cerebral arteries branch along the longitudinal fissure, and the posterior cerebral arteries branch along the basal surface of the temporal and occipital lobes. Because of the disposition of the main arteries, any electrical activity near these vessels will be affected by the conductivity of the blood flowing nearby, influencing surface measurements. Neglecting conductivity changes in deep brain structures caused by arterial blood flow, combined with small measured signals from these regions, leads to significant source localization errors.

Method
In this study, the EEG source localization error resulting from the interpolated conductivity atlas is evaluated using numerical simulations. For this objective, we use a forward model to simulate surface electrode measurements and inverse solvers to reconstruct the sources. Figure 1 presents an overview of the method we employed to assess the reconstructions.
We consider volumetric FEM based forward modeling to solve Maxwell's electromagnetic field equations for the non-invasive EEG [9,15,16], where 128 contact electrodes are attached on the scalp (figure 2) with electrode impedance of 2 kΩ uniformly. The multi-compartment head model contains superficial and deep brain structures using realistic geometries provided by an openly available anatomical T1-weighted Magnetic Resonance Imaging (MRI) data (https://brain-development.org/ixidataset/).
For the calculation of the lead fields, we have used the Zeffiro Interface toolbox for electromagnetic brain imaging [17]. Zeffiro utilizes highly adaptive and automated tetrahedral mesh generation for multi-compartment head models [18]. For source modeling, we use FEM and the current preservationbased H(div) approach [19][20][21].
To investigate the role of blood in EEG source localization, we created conductivity distributions σ f and σ p for the forward and inverse solvers, respectively, based on an anatomical atlas. The process to construct the conductivity distributions is detailed in the following subsections.

Conductivity distributions
Three conductivity distributions were considered: (a) σ 0 the conductivity of the main tissues in the human head without blood vessels and flow effects; (b) σ d  containing the conductivity of the main tissues with arterial flow effects in the blood vessels and vascular territories at the end of diastole; and (c) σ s similar to the second case but at peak systole. With this choice, σ 0 represents a static model that does not take into account the effects of blood and its flow, whereas σ d and σ s represent a dynamic model in two instants with minimum and maximum blood flows, respectively. It is worth emphasizing that choosing end diastole and peak systole instants can provide bounds to pulsatility effects. Conductivity of biological tissues depends also on the frequency. For this investigation, we will use a frequency of 3 kHz.

Anatomical atlas
In Moura et al [14], the authors develop a statistical atlas of the electrical properties of the human head based on 3D-MRI of 107 healthy human subjects. The atlas takes into consideration the natural variability of the internal structures and the electrical properties of the main tissues. The proposed atlas also comprises the dynamic effects of pulsatile arterial blood flow in the main vascular territories of the brain and their influence on the electrical properties of the brain. The relation between blood flow velocity and its electrical properties have been identified before, both in vitro and in vivo [11][12][13] and modeled as by Visser and collaborators [25,26]. This expression relates the longitudinal conductivity change ∆σ ℓ with the cross-sectional average velocityv in the blood vessel. In this expression, σ ref is the reference conductivity (still blood), H is the hematocrit (volumetric percentage of red blood cells),v is the average cross-sectional velocity, and R is the radius of the artery. The velocity is computed via the solution of the Navier-Stokes equations for the pulsatile blood flow in the main arteries. The blood flow model and simulation are described in the next subsection.
The atlas is presented as a Gaussian random vector S(t) given by where A is the static component of the atlas and A ′ (t) is the dynamic component, assumed to be independent of each other. We selected the averages in different situations to create σ 0 , σ d , and σ s to reflect the overall distribution of the population. From figure 5 in the results section, it is evident that the end of the diastole occurs at t = 0 s while peak of systole happens at t = 0.25 s and used to define σ d and σ s

Blood flow model
Visser's model (1) requires the cross-sectional average velocity in each vessel of the arterial tree. The velocities are obtained from the solution of the Navier-Stokes equations in the superior aortic arterial tree model, presented in figure 3, using the openBF solver [27,28]. Blood is assumed Newtonian with density ρ = 1050 kg m −3 , dynamic viscosity µ = 4.5 × 10 −3 Pa s, and hematocrit H = 0.5 moving in a network of narrow long circular vessels (arteries) with fully developed flow. The vessels are straight, have linearly elastic compliant walls, and are connected to other vessels at the ends. The following system describes the 1D blood flow, together with the constitutive equation of the compliant walls [27,28] ∂A ∂t where z is the longitudinal coordinate, A(z, t) is the cross-sectional area of the vessel, Q(z, t) is the volumetric flow rate, α is the Corioli's coefficient, P(z, t) is the static blood pressure, γ is a velocity profile shape parameter, P ext is the external pressure, E(z) is the vessel wall Young's modulus, and ν is the Poisson's ratio of the wall. A 0 (z) is the reference cross-sectional area, and h 0 the wall thickness at rest. Geometrical and mechanical properties of the arteries under consideration (figure 3) are presented in [14]. The boundary conditions at the terminal vessels are three-element windkessel models, and the input is a half-sinusoidal wave model of the cardiac output flow in one cycle where K = 485 ml s −1 , τ = 0.3 s, and the cardiac cycle period T = 1 s (60 bpm), following [29].

EEG localization error assessment
Source localization is conducted using three different methods: LCMV beamformer, sLORETA, and DS. All of them were implemented as plugins in the Zeffiro Interface. The methods were selected based on their widespread usage and methodological differences. In this case, we have the following mathematical linear presentation to relate the measured EEG data to neural activity sources modeled as 3D dipoles where y ∈ R m denotes the measurements obtained and L ∈ R m×n is the lead field matrix. Each triplet of the unknown vector x ∈ R n represents the coefficients of the basis current vector functions that form the estimated current distribution in 3D. The vector n represents the measurement noise that is assumed to be zero-mean Gaussian. The following subsections describe the compared inverse methods in detail.

LCMV beamformer
The beamforming is a widely used spatial filtering approach imported from signal processing to brain imaging [22,30]. In this study, we are using a vectortype implementation of LCMV beamformer, where the weight matrix is defined as The weight matrix is defined as the minimizer of the following optimization problem where C is the measurement noise covariance matrix and the vector p represents a single source position to which the constraints' submatrices correspond. The optimized weights are then used to estimate the source activity as a weighted sum of the measured data.

sLORETA
Standardized LORETA, also called sLORETA, is based on the idea of standardizing current density, which means that both noise variances in measurements and biological variances are assumed to be taken into account by scaling the minimum norm estimated reconstruction [31] by the resolution matrix [23]. In the case of noiseless measurements, this results in linear source localization with respect to data with zero localization error [32][33][34].

DS
DS is based on the best-fit solution that minimizes the residual variance between the measured data and a dipole forward mapped to measurement space [24]. The goodness of the fit for jth source location is defined as where L j ∈ R m×3 is a sub-matrix of the lead field and (·) + is the pseudoinverse operator.

Numerical simulation protocol and localization error mapping on the brain surfaces
Our aim with the analysis is to assess (a) whether we benefit from taking into consideration the dynamic effects of the arterial blood flow in source localization when we have a personalized conductivity model of the patient and (b) whether we benefit from the more accurate model when we employ an average conductivity model for various patients and with measurement noise of various signal-to-noise (SNR) ratios.

Personalized conductivity model
In this experiment, we have sampled 1500 uniformly randomized source locations within the brain model with diastolic conductivity σ d , representing data collection of a patient at the end of the diastole. We picked the corresponding location from the opposite hemisphere to ensure symmetry from a synthetic data perspective. Following that, a dipole was placed at each location one at a time to generate synthetic EEG data caused by the dipole adding 30 dB SNR Gaussian noise. Thus, we used a total of 3000 different data realizations in our experiment. We used the same random number generator seed to generate noise for each EEG measurement data set. In this way, the random effect created by the randomized noise has been removed between the models for comparison. The dipoles were then reconstructed using LCMV beamformer, sLORETA, and DS for three conductivity distributions used as prior information: (a) σ d indicating the case where the prior information about the conductivity is perfectly known, (b) σ s indicating a dynamic model but at the wrong instant, and (c) σ 0 neglecting the effects of blood.
From the reconstructions, we evaluated the location of a dipole using the maximum principle. The maximum principle is known to be a reliable estimation technique with dipole scanning [35][36][37] and different beamforming methods, e.g. oscillation [38,39] and standard [40,41] techniques. The method is also used to determine the location estimation via sLORETA [23].
The dipole localization error (DLE) was determined as the Euclidean distance from the true and estimated dipole locations. The DLEs obtained were placed in the head model for the location of each true dipole and extended to cover the entire model by using smooth nearest-neighbor interpolation, giving a visual representation for each case. Moreover, the localization discrepancies, i.e. distance between estimated locations, between models based on different prior information were calculated.

Interpatient variation and measurement noise level
The effect of conductivity is analyzed further by simulating 225 virtual patients σ i (.) . The patients are divided into 3 sets of 75 patients. In the first set, denoted as set B (Brain), the statistical atlas (2)-(4) is utilized. More precisely, we sampled 25 conductivity distributions for σ d , σ s and σ 0 each while fixing the conductivity of the skull, scalp and the electrode positions. This set represents a case where the geometry of the electrodes and skull conductivities are known, but the conductivity of the brain tissues are unknown. In the second set, denoted as BS (Brain+Skull), the state of the brain and the conductivity of the skull are unknown. In the third set, denoted as BSE (Brain + Skull + Electrodes), electrode positions are also not accurately known. Skull conductivity and electrode positions are Gaussian distributed. Standard deviation of the skull conductivity is set to 0.0075 S m −1 while the mean is 0.02 S m −1 , following [10]. The skull conductivity is not allowed to go under 0.0042 S m −1 , and scalp conductivity was set to 0.33 S m −1 . Standard deviation of electrode positions is 2.5 mm. Those conductivity values are then used to create simulated measurements of a deep activity from the left entorhinal cortex, shown in figure 4, where the conductivity difference between the average σ 0 and both average σ d and σ s is the largest. In this experiment setup, we used ten noise realizations with 30 dB SNR for each simulated patient. Dipole strength is selected to be 10 nA m, a typical dipolar primary current amplitude in the brain [42,43]. The localization is done in each of the virtual patients, assuming the standard model without vessels σ 0 and the arterial blood flow based model with σ d as prior information to invert (equations (5) and (6)). These priors emulate the situation of using an average head model on various patients.
The effect of measurement noise is studied for the case of a single patient from the Brain set, σ d using either the correct prior σ d or σ 0 and measurements with SNR ranging from 5 dB to 30 dB. We simulated 50 realizations for each SNR value to evaluate its effects over the localization of the source presented in figure 4.
To further analyze the effect of the vessel model in different locations of the brain and the presence of uncertainties in the blood, skull conductivity, and electrode positions, we performed a numerical experiment with 25 patients σ i d of set BSE by sampling 200 random source locations and adding a counter part of each of them from another hemisphere for each patient, resulting in 10 000 randomized source locations. The experiment and its visualization is done likewise in section 2.5.1. Briefly, we solved the inverse problem with both σ d and σ 0 priors, computed the DLEs, and determined the discrepancies between the solutions and the true locations.

Results
The wave forms of flow rates, cross-sectional average velocities, and longitudinal conductivity changes ∆σ ℓ /σ ref from the Navier-Stokes system (8)-(10) along the cardiac cycle are shown in figure 5. Simulations show that the conductivities of blood are 12%-16% higher than the reference (still blood). The variability of conductivity in different vessels is the result of different flow rates and vessel diameters. Figure 6 present slices of the conductivity distributions σ 0 , σ d , and σ s obtained from the anatomical atlas. In order to emphasize the differences, σ d and σ s are presented as differences between these distributions and σ 0 . The regions with the largest differences are located where the main arteries of the brain are located.
Relative differences between the lead fields of conductivity distributions σ s and σ 0 with respect to the lead field of σ d are presented in figure 7, together with histograms. The relative difference is defined (element-wise) as where L j ∈ R m represents a column of the lead field matrix L in (12). The figure shows that the relative difference ϵ(σ s ) is much smaller (Q 1 = 1.3 × 10 −3 , Q 2 = 1.8 × 10 −3 , and Q 3 = 2.6 × 10 −3 ), than the relative difference ϵ(σ 0 ) (Q 1 = 6.3 × 10 −2 , Q 2 = 8.4 × 10 −2 , and Q 3 = 1.2 × 10 −1 ). This shows that modeling the vessels brings considerable changes to the lead field, but the effects of pulsatility are much smaller. Figures 8-10 present exploded-view plots of the brain colored accordingly with the DLEs for each source localization method. The first three rows in each figure show the DLEs using σ d , σ d and σ 0 priors, respectively. The fourth and fifth rows present localization discrepancies between the reconstructions using σ d versus σ s and σ d versus σ 0 priors, respectively.   In the case of the LCMV beamformer and DS, we see an improvement in source localization using either σ d or σ s , i.e. prior models with vessels, in comparison to the solution with σ 0 (no vessels). DLEs are  larger with σ 0 as the prior, corroborating the fact that σ 0 is the worst prior among the options considered.

Personalized conductivity model
sLORETA errors are visibly poor in the deep structures of the brain and locations where conductivity values are highest. On the surface level of the brain and brainstem, however, we can obtain less than 10 mm localization errors. Additionally, we can obtain small differences between the static model and the ones derived from the arterial blood flow model. The obtained differences, however, are concentrated in areas close to the main blood vessels.

Interpatient variation and measurement noise level
Considering now the effects of conductivity uncertainties and measurement noise on source localization, figure 11 show boxplots of the localization errors from the virtual patients σ i (.) (equations (15)-(17)), adopting either the averages σ d or σ 0 as priors (equations (5) and (6)) and 30 dB SNR values. The results of all three sets of patients are presented in the figure.
The boxplots show improvement for all methodes in the cases that there is a match between the virtual patient σ i (.) and the average used as prior for the inversion, even in the presence of that level of noise, although the improvement becomes more modest as the amount of uncertainty sources increases (B → BS → BSE). This is expressed either with smaller median values or smaller interquartile ranges, specially for groups σ i d and σ i s indicating the correct prior is more robust in the presence of perturbation of various kinds. This observation is corroborated by Wilcoxon's rank sum test to check the null hypothesis H 0 that the localization errors obtained from σ d or σ 0 as priors come from the same distribution. The results of the tests are summarized in table 1 for each group of virtual patients and methods.
The table shows that set B benefits from using an accurate prior in general, except for LCMV and σ i 0 , where p-value does not reach 0.05. The trend is the same for BS and BSE sets, but the slight increase in p-value, specially for LCMV beamformer and DS indicates that the effect of the correct blood conductivity prior diminishes when more uncertainty Values are in mm. Left column: inversion is done using a lead field with averaged σ d conductivity. Right column: inversion is done using a lead field with averaged σ0 conductivity. factors are introduced for both the inversion and the measurement models aspects. Boxplots of DLEs as functions of noise SNR for σ d and σ 0 priors are shown in figure 12. The left column shows the results with the correct prior σ d . The boxplots show that reducing the SNR causes an increase in DLE. Noise also reduces the differences between   Figure 13 presents the localization discrepancies of 25 patients σ i d of set BSE when averaged σ d and σ 0 models are used as priors. The brain plots show differences up to 15 mm for sLORETA and LCMV beamformer, and 10 mm for DS. The localization discrepancies are spread everywhere in the brain model, but the highest differences are located on the brainstem and entorhinal cortices where the differences on blood conductivities are highest in the inverse models. The highest differences can be observed near the blood vessels as well as in the cerebellum, thalami and the posterior part of the brain using sLORETA, however, as presented before, sLORETA errors are visibly poor in the deep structures of the brain. Overall, in regions far from the main blood vessels, the discrepancies are smaller than 3 mm.

Discussion
It is known experimentally that blood flow affects the electrical properties of blood. In this in silico study, we assessed the effects of cerebral circulation on EEG source localization inverse problem. For that, a 3D time-varying anatomical atlas of the electrical properties of the human head is used to create virtual patient heads for EEG. The atlas contains a model of the cerebral circulation, connecting the Navier-Stokes equations with the electrical properties of blood via Visser's model. We considered virtual patients in three situations: (a) at the end of diastole and (b) at peak systole to evaluate the extremes of the pulsatile flow, along with (c) a static model without blood flow that is usually used for EEG source localization.
The results from the Navier-Stokes solver show that the blood in the main arteries of the brain is between 12% (diastole) and 16% (systole) more conductive than the reference (still blood). These changes are significant for source localization in deep structures near the brainstem, base of the thalamus, and hypothalamus because of the presence of the circle of Willis and main arteries of the brain.
EEG source localization was done using the LCMV beamformer, sLORETA, and DS. Dipole localization reconstructions using σ d (end of diastole) and σ s (peak sistole) as priors are very similar (figures 8-10), which might be due to the fact that the relative forward model difference is small ( figure 7). DLEs using σ 0 are larger in most cases. This difference indicates that the average flow rate plays an important role in source localization, while the pulsatility effect is very small and can be neglected. The results suggest that one could adopt average flow rates to set the conductivity of blood, substantially simplifying the analysis. It is important to stress that the Navier-Stokes equations are still needed to find average flow rates in each blood vessel.
Based on the visual representation of the DLEs, we can see the greatest differences at locations with the highest discrepancies in conductivity values, i.e. the base of the brain, brainstem, thalamus, hypothalamus, lateral sulci, longitudinal fissure, and basal surface of the temporal and occipital lobes. These are all regions near the main arteries of the brain. This result can be clearly seen in figures 8-10.
When noise is introduced together with varying conductivities (virtual patients), the results highlight the importance of an accurate forward modeling scheme in figure 11 and table 1. The aim was to simulate patient-wise conductivity deviations that can be witnessed in a real scenario. We studied the effects of the blood vessels with increasing number or sources of uncertainty, namely, brain tissue electrical properties, skull electrical properties, and electrode positions, in addition to measurement noise. The results show improvement in the cases where there is a match between the virtual patient and the prior average for moderate noise levels, specially for sLORETA and DS in all scenarios under consideration. However, when the amount of measurement noise increase above a certain level, these dominate the DLEs, masking the effects of conductivity priors.
The results indicate also that the effects of conductivity mismatch are detectable even for moderate measurement noise for DS, sLORETA, and LCMV beamformer ( figure 12 and table 2). Noting that the EEG measurement model assumes dipolar activity, it is somewhat expected to obtain good results with dipole fitting methods like DS and beamformer. Even with SNR noise of 30 dB, measurement noise seems to dominates conductivity prior mismatch in DS. The linear growth of the median DLEs of DS can be explained by the fact that it does not consider noise as its own modeling parameter, unlike the compared methods. Noise weakens the goodness of the fit and thus the confidence in the accuracy of the estimation.
Considering the inversion via the lead field with the correct conductivity distribution σ d and the static conductivity distribution σ 0 , we can state from the results that the conductivity model σ d is valid to 15 dB and after with sLORETA and LCMV beamformer and to 30 dB and after with DS to reconstruct deep activity. The table of the Wilcoxon rank sum test (table 2) further demonstrates this for sLORETA and LCMV. The DLE difference of DS shown by the box plots cannot be neglected because it shows numerically the difference between the compared models. On the other hand, even if the relative difference is high, DS is so accurate in both cases that the differences are small in clinical settings and in terms of EEG source localization accuracy and its limits [44].
Previously, Fiederer et al [45] studied the effects of blood vessels in a high-resolution volume conductor model for EEG. The authors observed large errors near the main cerebral arteries and piercing vessels, e.g. the carotid foramina and intraosseous veins. Their results show that the piercing vessels act as shunting paths for electrical signals, causing large errors if the model does not take them into account. Our work is consistent with their results and complements part of the limitations of the referred work as mentioned by the authors, namely taking in consideration interpatient variations, inhomogeneous segments, and the relation between the electrical properties of blood and its speed inside the blood vessels. The main novelty of our approach is that we employ a numerical solver of the Navier-Stokes equations for pulsatile flows in the superior aortic system, associated with Visser's model to connect these equations with the electrical properties of the blood. We also consider the main foramina to build the anatomical atlas. Although the averaging process to determine the statistics of the atlas makes these openings not very well resolvable in the images, their effects are statistically present.
Effects of conductivity uncertainties on localization of deep activities have not been studied much, but comparing to other studies about the influence of tissue conductivity to localization accuracy of dipole fitting methods on EEG [46][47][48] and in deep or otherwise inferior sources in a case of a static conductivity distribution [5,49], the obtained DLEs of DS are consistent with the reported ones. The accuracy of the LCMV beamformer was slightly worse than it has been reported with a cortical source in low and high noise cases [41] which could be explained by the depth of the true source. However, the discrepancy is so small that it can be explained by any modeling difference or by a combination of those.
The obtained measurement noise robustness of sLORETA has previously been demonstrated [50,51]. As a comparison to localization errors obtained with real measurements and similar SNR levels, DS in the case of σ d -based inverse model and sLORETA have lower median errors than those obtained at 20-30 dB (over 11 mm for an inferior source) [52,53]. However, it is worth noting that the SNR levels are estimated and not claimed to be defined, and the head models in both cases are simplified in the studies of Cuffin et al. On the other hand, our settings are synthetic, and the Simplex method used is not compared in this study. Blood flow prior model has an effect on source localization, even when skull conductivity and electrode locations uncertainty are added to the problem ( figure 13). The results show differences up to 15 mm for sLORETA and LCMV beamformer and 10 mm for DS in the brainstem and entorhinal cortices regions. In regions far from the main arteries vessels, the discrepancies are smaller than 3 mm.

Clinical perspectives and study limitations
The model of the superior aortic arterial tree model requires knowledge about the geometry of the vessels, as well as about the mechanical properties of their walls and physiological conditions of the patient. We can group the parameters into three main categories: (a) Geometry: vessel lengths and diameters, (b) Boundary conditions: cardiac output flow per cycle, heart rate, and terminal resistances, and (c) Mechanical: Young modulus of the vessel's walls.
It is important to stress that physiological particularities and/or pathological conditions might affect these parameters, e.g. ageing effects on the elastic properties of the vessel's walls, atherosclerosis that can greatly reduce the lumen of the affected vessels, or even the heart rate during EEG data collection. Ideally, the superior aortic arterial tree model must be adjusted to each patient to generate accurate blood flow waveforms and, as consequence, realistic electrical properties. Otherwise, a validated averaged model can be used [54], although it would be suboptimal for the specific case.
Some parameters of the model can be easily adjusted like heart rate and, to some extent, the geometry via non invasive anthropometric measurements to adjust the scale of the model. Cardiac output flow per cycle is more difficult to assess. Echocardiography is a good candidate to determine it, but would require the use of an ultrasound probe. In contrast, other parameters are today impossible or difficult to measure in a clinical setup, like the mechanical properties of the walls and terminal resistances.
The adjustment of the model to the needs or special circumstances of an individual can be specially important. A few examples are patients with stroke induced seizures, where part of the arterial tree might be permanently obstructed, in patients with atherosclerotic vessels that affects their elastic properties and reduce their lumen, and patients with impaired cerebral autoregulation that makes the blood flow more sensitive to pressure changes.
The model can be adjusted to cases like these by setting the parameters of the models, however model individualization for pathological conditions is still an open area of research, where some results can already be found in the literature [28,55,56].
Ageing effects can also be taken in consideration by setting the vessel's wall Young Modulus, however more studies are still necessary. Sensitivity analysis of the blood flow model with respect to its parameters has been studied before [55] but to what extent each of the parameters affects EEG localization error is not yet known and requires further investigation. These would be the best candidates to focus for individualization. Once the model is parameterized, the Navier-Stokes equations can be updated in reasonable time (minutes) and the atlas can be quickly updated.
The proposed model has limitations. (a) We based our analysis on in silico models. The following findings should be supported by experimental data; (b) there are no ageing effects in the blood vessel model; (c) no collaterals other than the circle of Willis were modeled; (d) the venous side of the circulation was not modeled; and (e) all tissues were modeled as isotropic.
As a further study, a more realistic measurement model could be used to make the comparison between dipolar fitting methods and distribution methods. This can be done, e.g. by constructing volumetric activity distributions within the brain model or utilizing dynamical models, like neural mass models [57], to create the measurements. The results of this study must be validated with experimental data, where cardiovascular data must be recorded in parallel to EEG data acquisition for control of the experiment.

Conclusion
In this study, we have examined the effects of modeled conductivity caused by blood flow on source localization using a dynamical anatomical atlas of the human head.
In cases where a personalized model of the head is available, blood circulation mismodeling causes localization errors, especially in the deep structures of the brain where the main cerebral arteries are located. The average flow rate plays an important role in source localization, while the pulsatility effect is very small and can be neglected.
The relevance of the blood flow model is highlighted by the localization of a deep dipolar source located at the position where conductivity differences between the models are most pronounced. The results indicate that the effects of conductivity mismatch are detectable even with moderate measurement noise and considering interpatient variations for sLORETA, and LCMV beamformer. However, smaller SNRs dominate the error, masking the effects of conductivity priors at a certain point. In the case of DS, noise dominates even with high SNRs, but the absolute error is the smallest between the methods studied. Blood flow prior model has an effect on source localization, even when skull conductivity and electrode locations uncertainty are added to the problem. The results show differences up to 15 mm for sLORETA and LCMV beamformer and 10 mm for DS in the brainstem and entorhinal cortices regions. In regions far from the main arteries vessels, the discrepancies are smaller than 3 mm While the source localization errors tend to increase towards the deeper areas for each method, the significance of the atlas model seems to be emphasized for the methods that are known to localize deeper activity in general, and in particular, the beamformer techniques.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: 10.5281/ zenodo.7277449.