Anatomical atlas of the upper part of the human head for electroencephalography and bioimpedance applications

Volume conductor problems in cerebral electrophysiology and bioimpedance do not have analytical solutions for nontrivial geometries and require a 3D model of the head and its electrical properties for solving the associated PDEs numerically. Ideally, the model should be made with patient-specific information. In clinical practice, this is not always the case and an average head model is often used. Also, the electrical properties of the tissues might not be completely known due to natural variability. The objective of this work is to develop a 4D (3D+T) statistical anatomical atlas of the electrical properties of the upper part of the human head for cerebral electrophysiology and bioimpedance applications. The atlas is an important tool for in silico studies on cerebral circulation and electrophysiology that require statistically consistent data, e.g., machine learning, sensitivity analyses, and as a benchmark to test inverse problem solvers. The atlas was constructed based on MRI images of human individuals and comprises the electrical properties of the main internal structures and can be adjusted for specific electrical frequencies. The proposed atlas also comprises a time-varying model of arterial brain circulation, based on the solution of the Navier-Stokes equation in the main arteries and their vascular territories. The atlas was successfully used to simulate electrical impedance tomography measurements indicating the necessity of signal-to-noise between 100 and 125dB to identify vascular changes due to the cardiac cycle, corroborating previous studies.


Introduction
Electrophysiology is the branch of physiology that investigates the electrical properties of biological tissues. The analysis is based on electrical measurements, voltages, or electric currents, generated by the tissue or in response to external electric stimuli.
One special group is clinical neurophysiology, where the bioelectrical activity is recorded to assess central and peripheral nervous systems. Electroencephalography (EEG) is an important monitoring and diagnostic method in this speciality to record brain electrical activity that can be used to diagnose thalamocortical rhythms, such as assessing seizure disorders, epilepsy, sleep disorders, coma, schizophrenia, Parkinson disease, and brain death (Michel andBrunet 2019, Jatoi andKamel 2017).
Electroencephalography measures voltage variations using multiple electrodes typically placed along the scalp of the patient. The measured voltages are the result of ionic currents inside the brain, therefore they are caused by spontaneous electrical activity.
In the case of epilepsy diagnostics, EEG is employed to determine the type, location, and extension of the lesion causing seizures. This is a challenging task since it depends on measurements taken on the surface of the scalp to infer the internal source of the disorder, a process called source reconstruction (Hallez et al. 2007, Grech et al. 2008.
Bioimpedance analysis is another group of methods used to assess the electrical properties of biological tissues. Measurements are made in response to external stimuli, such as measuring voltages caused by external current sources attached to the surface of the body or vice-versa.
Electrical impedance tomography (EIT) is a medical image technique in which electrical measurements on the surface of the body are used to create an image of conductivity distribution (or admittivity, the complex-valued equivalent) within the body. The image can then be associated with the physiological conditions of the patient. EIT has been used successfully in many areas, such as in lung applications to monitor acute respiratory distress syndrome, obstructive lung diseases or perioperative patients (Martins et al. 2019), monitoring mechanical ventilation, heart activity (Frerichs et al. 2016), cardiac function and detecting cancerous tissues (Adler and Boyle 2019). It is also being investigated to monitor brain activity and to distinguish between ischemic and hemorrhagic stroke (Adler and Boyle 2019).
Many electrophysiology applications require solving a nonlinear ill-posed inverse problem associated with the volume conductor. Reliable and stable solutions depend on prior information about the geometry and electrical properties of the tissues and knowledge about measurement uncertainty. In the case of source reconstruction, prior information about the electrical properties of the tissues is required for composing the forward volume conduction model needed in the process. Small errors in the electrical properties inside the head can obfuscate the effects of deep brain activity. In the case of EIT, anatomical and electrical prior information is also required to restrict the solution search space.
In addition, the brain is not a static structure. The flow of blood varies periodically during each cardiac cycle (Wagshul et al. 2011). Intracranial pulsatility has been eval-uated using magnetic resonance imaging (Vikner et al. 2019, or transcranial Doppler sonography (Kneihsl et al. 2020), and tissue pulsatility imaging (Kucewicz et al. 2008, Desmidt et al. 2018.
Blood flow centre-line velocity and artery radius influence the electrical impedance of blood in that artery (Gaw et al. 2008, Shen et al. 2016, 2018. This makes it promising to use electrical conductivity measurements, such as in EIT or impedance cardiography, in haemodynamic monitoring using surface electrodes on the skin. In fact, there have been several works recently aiming to monitor blood flow and/or pressure waveforms using them, e.g., common carotid arteries (Zhang et al. 2020), aortic artery (Badeli et al. 2020), pulmonary artery (Braun et al. 2018), radial artery , and cerebral arteries (Beraldo and Moura 2020). There are works about impedance cardiography to determine stroke volume (Bernstein 2010), electrical bioimpedance sensing to determine the central aortic pressure (CAP) curves , and pulmonary artery pressure estimation using EIT (Proença et al. 2020). There is an increasing interest in brain monitoring using electrical measurements, such as to monitor ventricular volume (Wembers et al. 2019), rheoencephalography to assess cerebral blood flow (Bodo et al. 2018, Meghdadi et al. 2019), brain perfusion of rats (Dowrick et al. 2016, Song et al. 2018, and stroke identification (Goren et al. 2018, Agnelli et al. 2020, Candiani and Santacesaria 2020, Candiani et al. 2019. Volume conductor problems in electrophysiology and bioimpedance do not have analytical solutions for nontrivial geometries and rely on numerical methods, e.g., finite element method (FEM) to discretize the head in small elements and solve the associated PDEs. Ideally, the FEM model should be built with patient-specific information, taken from MRI or CT scans to capture precisely the geometry of the head, its internal structures, and electrode positions. Unfortunately, in clinical practice, this is not the case. Often, an oversimplified geometry is employed for all patients due to the lack of computational tools and time.
The effects of mismodelling have been investigated before. EEG source localization errors increase substantially when individual-specific head models are not at disposal (Acar and Makeig 2013). The authors also show that errors in the conductivity of the skull cause large estimate errors. The latter is especially challenging because the electrical properties of the skull are highly heterogeneous and have large variability inter-individual. Cerebrospinal fluid (CSF) has a big impact on the results due to its high conductivity that forms a conductive layer surrounding the brain effectively shielding the interior (Vorwerk et al. 2014, Cho et al. 2015. Also, the authors show that distinguishing white and grey matters also impact the head volume conductor model. The objective of this work is to develop a statistical anatomical atlas of electrical properties of the upper part of the human head for electrophysiology and bioimpedance applications. The atlas is constructed based on MRI images of human individuals and comprises the electrical properties of the main structures for electrophysiology. The proposed atlas also comprises a time-varying model of the brain circulation, based on the solution of the Navier-Stokes equation for blood flow in the main arteries , Melis 2018. The atlas can be used to generate synthetic data statistically consistent with the population to compose learning sets for machine learning methods, for sensitivity analyses, and as a benchmark to test algorithms. The atlas can also be used as statistical prior information for inverse problems in electrophysiology. Anatomybased priors are found in the literature, such as for thorax applications (Martins et al. 2019, Kaipio andSomersalo 2005).

Anatomical Atlas Description and Construction
The anatomical atlas is composed of a static component A s with the electrical properties of the main tissues found in the upper part of the human head and a dynamic component A d (t) to account for blood perfusion dynamics in the human head. The two components of the atlas are considered Gaussian and independent, therefore the final statistics of the atlas A is composed by The two components of the atlas are presented in details in the following subsections.

Static component
The static component of the atlas distinguishes five main compartments of importance for electrophysiology of the human head: grey matter (GM), white matter (WM), cerebrospinal fluid (CSF), bones (BO) and other soft tissues (OT). Figure 1 depicts the general procedure to calculate the static component of the anatomical atlas. Fifty 3D Magnetic Resonance (MR) images of healthy human individuals, made available by the CASILab at the University of North Carolina at Chapel Hill were used (Bullitt et al. 2005). MR images are of type T1, obtained in a three-tesla equipment using the Fast Low Angle Shot (FLASH) sequence with a resolution of 1 × 1 × 1 mm. An equal number of male and female individuals were used, with an average age of 40 ± 15 years old.
The Symmetric image Normalization (SyN) method was applied to the images to diminish differences due to misalignment, aspect ratio, and sizes between the heads (Avants et al. 2008). For this purpose, the Advanced Normalization Tools (ANTs) was used, under the Neuroimaging in Python Pipelines and Interfaces (Nipype) framework (Gorgolewski et al. 2011). Detailed information regarding the normalization can be found in Avants et al. (2008Avants et al. ( , 2009  Each image is normalized to a reference image, segmented into the main compartments. Each segment is assigned with the corresponding electrical property and, finally, the statistics of the atlas can be estimated. with 1 × 1 × 1 mm resolution (Grabner et al. 2006, Fonov et al. 2009). The reference image is presented in Figure 2. Using an average head geometry as reference avoids having to choose one of the images in the dataset as reference, eliminating the possibility of choosing as a reference an individual with any abnormal geometric feature. Each transformation is performed in two stages, first a rigid transformation to roughly align the geometries, followed by an affine transformation. The normalization process assumes the heads have similar shapes and proportions. Therefore, the atlas represents a head with size and proportions similar to the reference head used for the normalization. If necessary, the resulting atlas can be transformed to accommodate other geometries, for example when the geometry of the head of the patient is available or if an average head model is preferred. This procedure will be described in Section 3.3 After the spatial normalization, the images were segmented into six classes: back-ground, WM, GM, CSF, BO, and OT. The Statistical Parametric Mapping (SPM) was used to segment the images (Frackowiak et al. 2004). The method is composed of three distinct components that are combined and optimized circularly: (i) modelling of the intensities of the images using a Gaussian mixture model; (ii) normalization of tissue probability maps of the five parts of the head with the images; and (iii) a bias field correction. Further details about the implementation can be found in Ashburner and Friston (2005). At the end of this phase, each voxel of the images is assigned to the label with the highest probability. Three additional steps were also performed to improve segmentation. (i) Any segmentation holes inside the head were filled with the nearest tissue in the image. This procedure was applied to all 2D slices in the three anatomical planes of each image, (ii) Small segmentation artefacts outside the human head were removed by isolating the largest connected group in the image using a six-connected neighbourhood strategy. (iii) Four iterations of morphological opening operation to the binary mask of the head to smooth the external surface of the head.

Electrical properties of the segments
Each voxel of the segmented images was assigned to the electrical property of the corresponding tissue before computing the statistics of the atlas. Tissues were modelled as isotropic, even though it is known that some tissues are anisotropic. The electrical properties depend on the type, physiological conditions and frequency in consideration Corthout 1996, Gabriel, Lau andGabriel 1996a).
Given the angular frequency of the electrical signal ω = 2πf , the complex relative permittivityˆ (ω) of a tissue can be modelled as the sum of four Cole-Cole dispersion terms (Gabriel, Lau and Gabriel 1996b) where 0 is the permittivity of free space and all the other parameters depend on the tissue Gabriel 1996b, Andreuccetti et al. 1997). The conductivity σ and permittivity of the tissue can be obtained fromˆ (ω) Biological tissues are naturally inhomogeneous, due to their complex macroscopic and microscopic structure, function and physiological condition. To account for this, the uncertainty level of the estimates from the above equations was set to ±20% following reported results in Gabriel, Lau and Gabriel (1996a).
The Electrical properties of BO were modelled as the average between cortical and cancellous bones, while OT was modelled as muscle tissue.

Atlas statistics computation
Let u ∈ R N V be a vector representing a 3D image of the human head after normalization and segmentation, where N V is the number of voxels, excluding those representing the background around the head. Let the image be segmented into N T nonintersecting regions (tissues), each one with associated characteristic funcion χ t ∈ R N V , for t = 1, 2, · · · , N T . Also, let p t be the electrical property of each tissue under consideration. The property can be real (e.g., resistivity or conductivity) or complex valued (e.g., impeditivity or admitivity). We can write the 3D image of this property x ∈ C N T (or real with the same dimension) as where p ∈ C N T is a vector composed by the electrical properties of the tissues and X ∈ R N V ×N T is a matrix where each column is a characteristic function. Assume images from N I individuals are used to build the atlas. Formally this number should be very large to represent the statistics of the population. In practice, this number is limited by the size of the dataset. To reduce this limitation each individual will be considered N S times, each time with a different value for p, following the statistics of the tissues. This implies that the N I individuals represent the general shape of the head of the population while allowing the electrical properties of the tissues to be more diverse. In addition, we assume the same number of samples N S per individual, making them equally probable.
Using these hypotheses the average and covariance of the population can be estimated efficiently. For the covariance matrix Γ in special, the formulation allows the computation in factorized form Γ = KK T , reducing storage requirements and simplifying algorithms that depend on factorizations of Γ.
Samples of the i-th individual can be composed by sampling the properties of the tissues p s and applying (7) x s,i = X i p s , for s = 1, 2, · · · , N S , where the samples p s can be generated from data fitted models or measurements. In this work, the electrical properties of the tissues are considered Gaussian with average resulting from the model (4) and standard deviation of 20% of the average, following reported results (Gabriel, Lau and Gabriel 1996a). Let x s,i represent a sample of the i-th individual. The average over all individuals can be estimated with where, again, N S is fixed and represent the number of samples with the same head geometry X i andp is the average electrical properties of the tissues.
The covariance matrix can be computed using the usual sample estimator where Γ ∈ R N V ×N V and M H denotes conjugate transpose of M . For real valued p, the conjugate transpose is the transpose M H = M T Adding (x i −x i ) to both terms between parenthesis and rearranging the terms, Proceeding with the products, Finally, taking the limit N S → ∞, Defining ∆x i =x i −x and observing the linear relation (7) we can rewrite this last expression as where W i ∈ C N V ×(N T +1) and Γ p ∈ C N T ×N T is the covariance matrix of the electrical properties of the tissues. Furthermore, the expression can be simplified to where Note that both the average (9) and the covariance (20) estimates do not require explicit sampling procedure presented in (8).
In case of complex-valued p, the pseudo-covarianceΓ ∈ C N V ×N V can be computed in a similar wayΓ The atlas in this work is assumed Gaussian, thereforex and Γ (andΓ for complexvalued p) completely specify its probability density function.

Dynamic component
Flow in the cranial cavity is pulsatile, following the cardiac cycle. Arterial blood flows in waves, forcing part of the venous blood and CSF to move, following the Monro-Kellie hypothesis (Wagshul et al. 2011). Venous blood is drained to the jugular veins via cerebral sinuses, while CSF moves in the subarachnoid space and leaves/returns to the cavity via the foramen magnum to balance intracranial pressure waves during the cardiac cycle (Greitz et al. 1992, Sakka et al. 2011. Recently, pulsatility was also observed in small cortical veins (Driver et al. 2020), and studies in rats show that the flow in microvessels is quasi-steady laminar flow, following Hagen-Poiseuille law expected in low Reynolds and Womersley numbers (Seki et al. 2006).
Arterial blood enters the cranial cavity through its base via two pairs of arteries, the (right/left) vertebral and internal carotid arteries. After entering, these arteries form the circle of Willis, a circulatory anastomosis responsible for providing backup routes for cerebral blood supply (Bradac 2017, Chandra et al. 2017. From the circle of Willis four main pairs of arteries branch out, the (right/left) anterior, middle, posterior cerebral, and superior cerebellar arteries.
The dynamic component of the atlas comprises the circulation in the main cerebral arteries. The procedure follows the same main steps presented in Section 2.1 with a few modifications. (i) Magnetic Resonance Angiography (MRA) images of 109 healthy human individuals were used (Bullitt et al. 2005). The images were obtained in a three tesla equipment with a resolution of 0.5 × 0.5 × 0.8 mm. An equal number of male and female individuals were used, with an average age of 43 ± 14 years old. (ii) Segmentation was performed first by applying a total variation filter to the images followed by a threshold segmentation. Only the lumen of the vessels with contrast agent were segmented. (iii) Electrical property assignment follows the procedure described in the following subsection.

Electrical properties of the segments
The influence of blood flow centre-line velocity and vessel radius over the electrical impedance of blood is modelled and included in the atlas (Gaw et al. 2008, Shen et al. 2016, 2018. We simulated blood flow in the main arteries of the brain using the openBF solver , Melis 2018), a 1D blood flow solver based on monotonic upstreamcentered scheme for conservation laws (MUSCL) finite-volume numerical scheme. The solver assumes the blood is an incompressible Newtonian fluid, flowing through narrow and long circular vessels with linear compliant walls. The Navier-Stokes equations are reduced to 1D by imposing axisymmetry, linearized and solved for pulsatile flows using the finite difference method. Detailed description can be found in Melis (2017).
Brain circulation simulation encompasses the superior aortic system, from the ascending aorta to the main arteries providing blood to the brain. The arteries considered in the simulation can be seen in Figure 3.  Table 1.
Blood was assumed Newtonian with density ρ = 1050 kg m −3 and dynamic viscosity µ = 4.5 × 10 −3 Pa s. The geometry and mechanical properties of the vessels are presented in Table 1, based on Alastruey et al. (2007) and complemented with data collected by Dodo et al. (2020), Fomkina et al. (2016), Schmitter et al. (2013). The terminal vessels were coupled with 3-element Windkessel models to mimic the perfusion of downstream vessels and avoid numerical oscillations. Heart flow output in one cardiac cycle was set to where Q M = 485 ml s −1 is peak flow rate, τ = 0.3 s and the cardiac cycle period is 1 s, following Alastruey et al. (2007).  Lasting one cardiac cycle, the simulated pulsatile blood flow of each vessel must be converted to the electrical property of interest. Visser's model, a nonlinear function that relates blood resistivity changes to the average blood velocity in a cylindrical vessel, can be used for this purpose (Visser 1989, 1992, Hoetink et al. 2004) where ∆ρ is the longitudinal resistivity change to the reference (still blood) resistivity ρ 0 , H is the hematocrit (volume percentage of red blood cells in the blood),v is the average cross-sectional velocity and R is the radius of the vessel. Visser's model presents a similar expression for the conductivity, however no expression was derived for other electrical properties. We will hypothesize the conductivity expression can be applied to the permittivity of blood. Visser's model applies to blood flowing in a rigid vessel and in a defined orientation. Measurements taken from impedance cardiography studies in humans showed relative variations ∆ρ /ρ 0 smaller (15% maximum) than predicted from Visser's model in the same conditions (25% maximum), 60% reduction (Raaijmakers et al. 1996). The difference can be explained by the fact that the vessels are not straight and have different orientations. To accommodate this discrepancy, changes in blood resistivity were scaled to 60% of Visser's model (23), as reported in the literature.
For each time step, the electrical property of the blood in each main artery is calculated using Visser's model and used to compute the statistics of the atlas at that time instant.
The volume occupied by the main arteries is small compared to the volume of the brain, however its area of influence is considerable. Each artery is responsible for providing blood to specific areas of the brain, known as brain arterial vascular territories (Kim et al. 2019). Six main supratentorial vascular territories were modelled, (right/left) MCA, ACA, and PCA, also the (right/left) superior cerebellar artery (SCA). The main brain territories can be seen in Figure 4. In addition to these, the (right/left) external carotid territories were also included due to the proximity to the electrodes that can impact measurements. The dynamic model does not consider collateral circulation other than the redundancy coming from the circle of Willis. Blood supply inside each territory is assumed to be proportional to the waveform of the associated main artery. To the best of our knowledge, there are not many studies on electrical property variations of brain tissues along the cardiac cycle. The majority of the studies focus on electrical property changes in response to sensorial or motor activity or epilepsy events (Newell et al. 2002, Tidswell et al. 2001, Towers et al. 2000, Holder et al. 1996. The net electrical property change is caused by a dynamic balance between the amount of blood, extracellular fluids, and cell swelling in a given location and at a given time instant. In this study, the dynamic component of this variation is set to 0.5% of the main artery of the territory (Tidswell et al. 2001).

Atlas statistics computation
Based on the segmentation of the vessels, explained in Section 2.2, and the location of the vascular territories, presented in Figure 4, it is possible to define two characteristic functions per territory, the main vessels in a given territory χ m ∈ R N V and its area of influence χ c m ∈ R N V (the complement within the vascular territory). Let N M be the number of vascular territories. We can write a 3D image of the electrical property of the dynamic component x(t) ∈ C N T as where p m,B (t) and p c m,B (t) are the electrical properties of blood in the respective segments. The effect of blood in the complements p c m,B (t) is modelled as a percentage of p m,B (t) where α m ≥ 0 adjusts the effect. The vector χ α,m represents the influence region of each vascular territory and can be used to compose images as in (7). The statistics of the dynamic component is computed following the same procedure presented in Section 2.1.2 with the modified characteristic function χ α,m . Figure 5a shows two representative individuals in the segmentation steps to compute the static component of the atlas. The first row contains the original images, the second row contains the normalized images and the third row show the segmented tissues. Figure 5b shows the average over the characteristic functions of all individuals and for each segmented tissue. A voxel with a value equal to 1.0 indicates it was classified as the same tissue in all images. The static component of the atlas at 1 kHz is presented in Figure 6. The figure shows transversal slices of the average (Figure 6a) and standard deviation (Figure 6b). It is possible to see high resistivity regions in the forehead, caused by the thick bone and the frontal sinus, in the zygomatic bones and the petrous part of the temporal bone in the base of the skull.  Figure 7 shows slices of the atlas built in terms of conductivity, resistivity, and relative permittivity and in different frequencies. The figure shows that the average resistivity and permittivity decrease with increases in frequency while conductivity increases. Although the average process tends to eliminate small features of the images, it is still possible to see small and thin structures inside the brain, like the longitudinal fissure, third and fourth ventricles and central canal.

Dynamic component
The dynamic component of the atlas was computed following the procedure described in Section 2.2. Transversal slices of the average at 1 kHz are presented in Figure 8. The main vessels that compose the circle of Willis in the base of the cranial cavity, the dense arterial vascularization in the insular cortex, and the superior sagittal sinus are visible. Figure 8: Transversal slices of the average image of the segmented vessels filled with still blood at 1 kHz. Small values were masked in grey to emphasize the structure of the main vessels. Figure 9 presents the waveforms obtained from the Navier-Stokes solver simulating one cardiac cycle (60bpm), with hematocrit H = 0.5. From left to right, the figure presents flow rate, average cross-sectional velocity, static pressure and resistivity changes to still blood following Visser's model (23). Peak velocity occurs approximately 0.25 s after the beginning of the cardiac cycle. Resistivity changes lie within -17% and -21% to still blood, indicating the resistivity in these vessels differs considerably from the electrical properties of still blood.

Effects of the cardiac cycle on surface measurements for electrical impedance tomography
As one example of application, the atlas was employed to simulate EIT surface electrode measurements at 1 kHz during one cardiac cycle (60bpm) using the finite element method to solve the complete electrode model for EIT (Cheney et al. 1999, Holder 2005. This type of in silico study is important to investigate the possibility of monitoring blood perfusion anomalies in patients. A segmented head image of an average young adult 2 was selected to create the geometry (Hammond et al. 2017, Hou et al. 2017, Song et al. 2013). This geometry is not among those used to create the atlas, avoiding statistical biases. The boundary surfaces of the segments were extracted and cleaned to remove artefacts.
A FEM mesh was created using Gmsh software (Geuzaine and Remacle 2009). The mesh, presented in Figure 10a, is composed of 2.06 million linear tetrahedral elements, split into six segments and 32 electrodes with a diameter of 15 mm were placed in two parallel planes with 20 mm of separation. Electrode numbers can be seen in the figure.
The atlas was projected into the FEM mesh using the symmetric image normalization (SyN) (Avants et al. 2008). The method requires two 3D binary images with the characteristic functions of the same volume, one in the atlas reference system and one in the FEM mesh reference system. The volumes used for the normalization comprise the soft tissues inside the cranial cavity (GM+WM+CSF). For the atlas, the characteristic function of the cranial cavity χ C is found by taking the average of the sum of the characteristic functions of these tissues over all individuals that compose the atlas, followed by thresholding at 75% where (i, j, k) is the coordinates of each voxel. For the FEM mesh, a 3D image can be created from the segmented internal structures of the mesh, seen in Figure 10a, by defining a 3D grid of points that encloses the head and checking if each voxel belongs to the segment GM+WM+CSF. The affine transformation resulting from the normalization was applied to all voxels of the atlas, projecting them into the FEM mesh reference system. Finally, the values of the atlas were interpolated into the centroids of the tetrahedron. The projection can be seen in Figure 10b.
EIT measurement simulation was performed by imposing sinusoidal bipolar current injection of 1 mA at 1 kHz and computing the electrode voltage measurements. Current pattern follows a skip-8 scheme to allow diametral current injection (two planes with 16 electrodes each) (Silva et al. 2017). This choice mitigates the electrical shunting effect of the skull that causes the majority of the current to flow along the scalp only if the pair of injecting electrodes are too close. The simulated measurements are presented in Figure 11 in four time instants along the cardiac cycle.   Figure 11: Simulated electrode measurements. (a) reference measurements v 0 ; (b) relative differences ∆v(t); (c) first 32 relative differences ∆v(t); (d) histograms of the differences ∆v(t). Figure 11a presents the measurements at t = 0 s, used as reference measurement v 0 . Figure 11b shows relative differences ∆v(t), in dB, between measurements in three other time instants v(t) and the reference v 0 , defined as where the division is computed element-wise. Figure 11c is a plot of the same differences for the first 32 measurements and show that the largest differences are measured at t = 0.25 s, time instant when we have peak velocities, as presented before in Figure 9. Figure 11d presents histograms of the differences ∆v(t). In the same histograms, the vertical lines represent the percentiles 10%, 50% (median), and 90% together with their numeric values (p10, p50, and p90).

Discussion
Conductor volume problems in electroencephalography and electrical bioimpedance cerebral monitoring require a 3D model of the head and its electrical properties for solving the associated PDEs numerically. In many situations, a 3D model of the head of the patient is not available or an average head model is preferred. Even in cases when the model is available via MRI or CT images, the electrical properties of the tissues might not be completely known due to natural variability. This work presents a novel 4D (3D+T) statistical anatomical atlas of the electrical properties of the human head for electrophysiology applications. The atlas was built for an average head shape and can be normalized to specific geometries. This process was exemplified with one EIT application. The normalization step optimized the alignment of the cranial cavity, volume comprising GM WM and CSF segments. In our experiments, this choice produces a better match between the atlas and the FEM mesh. However, this choice can cause small artefacts in the external surface of the mesh due to the small thickness of the scalp and its proximity to the skull. These artefacts can be seen in Figure 10b where the resistivity near the top of the cranial vault and the anterior part of the frontal bone affect the resistivity near the external surface.
The electrical properties of biological tissues are frequency dependent. The atlas can be built for different frequencies as exemplified in Figure 7. This flexibility expands the applicability of the atlas, such as in multifrequency EIT (Horesh 2006, Malone et al. 2014. Cerebral circulation was also modelled and added to the atlas. The atlas is capable of simulating the pulsatile blood flow in the main cerebral arteries and vascular territories and their effects on electrical measurements. As one example of application, the atlas was employed in an in silico study to investigate the possibility of monitoring blood perfusion using EIT. Among other objectives, this type of study is important to provide information on measurement sensitivity necessary to detect perfusion anomalies and serve as a guide for future EIT equipment developments for these applications. Figure 11d predicts that EIT equipments need a signal-to-noise ratio between 100 and 125 dB to identify changes due to the cardiac cycle. This is in good agreement with a previous study (Towers et al. 2000) on clamped carotid arteries.
The statistical nature of the atlas also allows quantifying its uncertainty. Figure 6 shows that the regions with the largest standard variations are located along the boundary of the bones. This can be explained by the fact that small anatomical differences between the skulls of the individuals cause large resistivity variations due to the difference between the resistivity of bones and other tissues around them.
The proposed atlas has some limitations. It is known that ageing increases the stiffness of the vessels, however the atlas does not include ageing effects on the stiffness of the walls of the arteries. Nevertheless, the atlas can be adjusted by setting Young's modulus of the vessels accordingly with the age of the population. Except for the redundancy caused by the circle of Willis, collateral brain circulation was not modelled either caused by preexisting vascular redundancy or neovascularization. The venous side of the circulation was not modelled. Although only five main tissues were segmented, the method can be readily extended to accommodate more tissues. All tissues were modelled as isotropic, even though it is known that some tissues are anisotropic. Extending the atlas to anisotropic tissue is possible but increases the complexity substantially. Also, the scalp was modelled as a uniform tissue, however it is a multi-layer tissue, composed of skin, connective tissue, epicranial aponeurosis, and muscles that have different electrical properties. Due to its proximity to the electrodes, scalp mismodellings can impact EIT recovered images. Finally, the dynamic effect of blood circulation in each territory was modelled as proportional to the velocity of the blood in the main vessel. This approximation does not take into consideration variations in the volume of blood in a given region, for example when the brain responds to external stimuli. The current limitations of the atlas act as motivation for future research topics. These challenging limitation will be the focus of future works to further improve the anatomical atlas.

Conclusion
We presented a novel anatomical atlas of the electrical properties of the human head. To the best of our knowledge, the present model is the first model capable of simulating cerebral circulation and its effects on electrical measurements. Despite the limitations, the atlas brings important implications to cerebral electrophysiology studies. This novelty has the potential to become an important tool for in silico studies on cerebral circulation and electrophysiology, such as electrical measurement sensitivity to vascular pathologic conditions like stroke classification and monitoring, arterial vasospasms, and arteriovenous malformation. The atlas can also be used as statistical prior information for inverse problems in EEG and EIT and to create training sets for machine learning algorithms.
The MR brain images from healthy volunteers used in this paper were collected and made available by the CASILab at The University of North Carolina at Chapel Hill and were distributed by the MIDAS Data Server at Kitware, Inc. (https://www. insight-journal.org/midas/community/view/21)