# Looking for the Signals of the Missing Baryons in the Extragalactic
Background Light

Wei Zhu1111Corresponding author<EMAIL_ADDRESS>and Rong Wang2
1 Department of Physics, East China Normal University, Shanghai 200241, China
2 Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000,
China

###### Abstract

The missing baryons in the universe are assumed to be hidden in the whole
space as a warm-hot intergalactic medium (WHIM). Finding them is one of the
important subjects in modern cosmology. In this paper, we point out that the
very high energy electron/positron rays may light up the WHIM due to the
anomalous bremsstrahlung according to an improved Bethe-Heitler formula. The
resulting excess of the extragalactic background light (EBL) can be observed
by the direct measurement method. A possible explanation on the difference
between the direct and indirect measurements of EBL is also proposed. Thus, we
find a new window to probe the WHIM properties via the EBL.

keywords: Missing Baryons; Improved Bethe-Heitler Formula; Extragalactic
Background Light

##  1 Introduction

The standard cosmological model (SCM) considers that most $(95\%)$ of the
universe is composed of dark energy $(70\%)$ and dark matter $(25\%)$, while
the universe’s remaining $5\%$ is the visible matter partly in the form of
stars, planets, interstellar dust… However these observed matter so far can
only account for half of the visible matter in the SCM predictions (see Refs.
1 and 2 for details). It implies that at least $1/2$ of the baryons are now
“missing”.3,4 Cosmological simulations show that these “missing” baryons in
the intergalactic medium may condense into a filamentary web and undergo
shocks that heat the baryons up to the temperatures $T\sim 10^{5}-10^{7}K$.5,6
Finding this warm-hot intergalactic medium (WHIM), which contains the missing
baryons, is crucial to validate the SCM. The WHIM is hot and highly ionized,
so that it could absorb or emit far-ultraviolet and soft X-ray photons.
However, these characteristic lines are formed in the atomic bound state. This
spectral method is invalid if the WHIM is too hot to show up in the Lyman-
absorption, or too diffused to emit the detectable X-rays.7 Recently, the
dispersion of $X$-ray afterglows from fast radio bursts (FRBs) is used to
prove the existence of the WHIM.8,9 Although the preliminary results from
these studies tend to support the WHIM model, the conclusion has a large
uncertainty since the FRB sources were produced occasionally and their
positions in the universe are undetermined. Still, all mentioned measurements
above are not enough to give a complete physical picture of the WHIM.

The purpose of this article is to try to find the possible signals of the
missing baryons in the visible light band, which relates to the extragalactic
background light (EBL). The EBL is the diffuse background radiation
accumulated over the cosmic history from ultraviolet (UV) to infrared (IR)
wavelengths. A peak in the EBL spectrum around $1~{}\mu m$ is formed by the
direct emission from stars. Measurements of both the EBL spectral intensity
and its evolution are important to study both the star formation and the
galaxy evolution processes.10 We find that the cosmic rays may light up the
WHIM and lead to the observed excess of the EBL.

Direct measurements of the optical EBL are made by observing the total sky
brightness and with the measurements or modeling of all the foreground
components.11,12 In this method, these foreground components must be very
accurately known. However, the bright zodiacal light (ZL), which is created by
the scattered sunlight off the interplanetary dust particles, may arise a
large uncertainty for the weak EBL intensity.

To get rid of this problem, the indirect techniques were developed: the
counting galaxies and the gamma-ray attenuation measurements.13 These indirect
measurements are complementary for the direct measurements and they are
important to constrain the large uncertainty of EBL. At the same time, the
direct surface photometry measurements have been improved. For example, the
detectors are set at the vantage points, where the Earth’s atmosphere and the
light from interplanetary dust are not the appreciable components of the
diffuse sky brightness. These researches show an excess brightness above the
indirectly measured EBL intensity. In particular, the NASA’s New Horizons
spacecraft is presently over 50 AU away from the sun, where the ZL almost
disappears. A recent report from New Horizons shows a new component of unknown
origin of EBL of $8.06\pm 1.92~{}nWm^{-2}sr^{-1}$.14 This discrepancy is
interesting because it might point to the presence of a truly diffuse emission
component in the universe that is not resolvable into point sources.

Let us imagine a given small range in the sky, where very high energy (VHE)
cosmic electrons and positrons interact with the Coulomb field in the WHIM and
emit a lot of soft bremsstrahlung photons. According to quantum
electrodynamics (QED), the contribution of these soft photons to the EBL
intensity is only $n_{\gamma}\sim 10^{-20}~{}nWm^{-2}sr^{-1}$ (see Sec. 3),
which is much smaller than the observed $n_{\gamma}\sim 10~{}nWm^{-2}sr^{-1}$.
Therefore, the contributions of the normal bremsstrahlung photons can be
neglected.

However, the above conclusion should be reconsidered if we use an improved
Bethe-Heitler formula.15 The reason is as follows. The electromagnetic shower
is a common phenomenon in astronomy, which comprise the bremsstrahlung and the
pair production in the nuclear Coulomb field (Figs.. 1b and 1c). The quantum
theory of the electromagnetic shower is based on the Bethe-Heitler formula,16
which has been widely applied in many branches in physics including
astrophysics. However, the recent work [15] pointed out that the Bethe-Heitler
formula should be modified to a new form in the dilute ionized gas, where the
target recoil can be neglected for the incident VHE cosmic rays. In this case
both the bremsstrahlung and the pair production cross sections have an
unexpected big increment, as $\sigma_{brem}\sim 1/m^{2}_{e}$ in the Bethe-
Heitler formula Eq. (2.1) is replaced by $\sigma_{brem}^{MD}\sim 1/\mu^{2}$ in
the modified version Eq. (2.2) and the screening length (or the Debye length)
is $\mu^{-1}\gg m_{e}^{-1}$. This leads to a big screening-length effect,
which was neglected for a long time.

The ionized gas of extreme low density is an ideal place for testing the
anomalous bremsstrahlung and pair production since the ionized atoms have a
large screening length. The nuclear Coulomb potential may spread into such a
broad space, where the bremsstrahlung or pair production induce almost no
recoil of the target since the incident high energy electrons or photons are
far away from the nucleus. This result can explain the strange difference
among the precise measurements of cosmic electron-positron fluxes at the GeV-
TeV band by Alpha Magnetic Spectrometer (AMS), Fermi Large Area Telescope
(Fermi-LAT), DArk Matter Particle Explorer (DAMPE) and Calorimetric Electron
Telescope (CALET).17 We find that the anomalous cross section has about eight-
orders of magnitude higher than the normal cross section in this example.

We will discuss the corrections from the anomalous bremsstrahlung to the
directly measured EBL intensity in Sec. 3. We find that the distribution of
the EBL intensity will be moved upward if adding the bremsstrahlung
contributions. The computed result is close to the observation data.

Then we discuss the anomalous pair production in the indirect measurement
method using gamma-ray attenuation with the same theoretical framework in Sec.
4. We noticed that the data of gamma-ray attenuation coincide with that of the
counting galaxies. However, the processes in Figs. 1b and 1c do not contribute
to the counting galaxies since the WHIM is diffusely distributed in the space.
Therefore, we consider that the screening length $\mu^{-1}$ takes a smaller
value in the gamma-ray attenuation method. We will give an explanation for
that in Sec. 5. We find a strong dependence of the directly measured EBL
intensity on the screening length in the electronic structure of the WHIM.
This result may open a new window to expose the electronic properties of the
WHIM during the universe evolution.

##  2 Improved Bethe-Heitler formula

We sketch and review some results of the improved Bethe-Heitler formula from
the work [15]. When charged particles scatter off the electric field of proton
or nucleus, they can emit real photons in such interactions. This is the
famous bremsstrahlung (braking radiation). A quantum-mechanical description of
bremsstrahlung emission at the Born approximation is the Bethe-Heitler
formula. This formula for the high-energy incident electrons and at the soft
photon limit reads

$\sigma_{brem}(E,\omega)d\omega=\frac{4\alpha_{s}^{3}}{m^{2}_{e}}\ln\frac{2E}{\mu}\frac{d\omega}{\omega},$
$None$

where $\mu^{-1}$ is the screening length. Equation (2.1) presents a strong
reduced bremsstrahlung cross section since $1/m^{2}_{e}$ is much smaller than
the geometric cross section $1/\mu^{2}$ of atom, where $\mu^{-1}\geq
10^{-8}~{}cm$.

The work [15] proved that if the target recoil can be neglected for a VHE
incident electron, the Bethe-Heitler formula should be modified as,

$\sigma^{MD}_{brem}(E,\omega)d\omega=\frac{4\alpha_{s}^{3}}{\mu^{2}}\ln\frac{2E}{\mu}\frac{d\omega}{\omega}.$
$None$

It predicts a big enhanced bremsstrahlung cross section. This is called the
anomalous bremsstrahlung effect.

The deeply ionized gas in the cosmic space is an ideal place to see the
anomalous bremsstrahlung effect since the dilute ionized atoms have a large
screening length, where the nuclear Coulomb potential may spread into such a
broader space, and the recoil of the bremsstrahlung event can be neglected
since the incident high energy electrons are far away from the nuclei.
Bremsstrahlung under this environment exhibits a big cross section. On the
other hand, the bremsstrahlung probability also decreases with the reduced
density. The above two opposite effects lead to a limited kinematic window,
where the bremsstrahlung cross section becomes anomalous. According to the
analysis of the work [17], the TeV-band of electron energy is such a possible
window.

The similar argument also satisfies for the electron-positron pair creation
from the real high energy photon. We consider a high energy photon traversing
the atomic Coulomb field. This photon has a probability of transforming itself
into a pair of electron-positron. The corresponding cross section in the
leading approximation is given by,

$d\sigma_{pair}=\frac{\alpha^{3}}{m_{e}^{2}}\ln\frac{4\omega^{2}}{\mu^{2}}(1-z)[(1-z)^{2}+z^{2}]dz,$
$None$

where $z=E_{e}/\omega$. In a situation where the target recoil can be
neglected, we have the anomalous pair production cross section,

$d\sigma_{pair}^{MD}=\frac{\alpha^{3}}{\mu^{2}}\ln\frac{4\omega^{2}}{\mu^{2}}(1-z)[(1-z)^{2}+z^{2}]dz$
$None$

The above improved Bethe-Heitler formula has been used to explain the
difference of VHE electron-positron fluxes in the ionosphere about $400\sim
500~{}km$ height above the Earth ground, where the oxygen atoms are strongly
ionized and have extremely dilute density.

We take AMS-Calet-DAMPE-FermiLAT data to test the improved Bethe-Heitler
formula since their measurements are on orbit. There are other data of cosmic
electron-positron flux, for example, HESS and ATIC are ground-based
measurements and balloon observations. They relate to the propagating models
of electrons in atmosphere and both them have a large uncertainty. Therefore,
we did not use them.

The anomalous bremsstrahlung effect relates to the degree of ionization and
density of atoms. We only give an order of magnitude estimation of the
screening length, which is much larger than the atomic scale. This is
consistent with the characteristics of the physical environment in plasma. On
the other hand, at altitudes 400-500 km the plasma density and ionization
degree are dependent from solar activity considerably. It would be result to
lepton flux variations with time. These need to be further tested.

##  3 Anomalous bremsstrahlung in the WHIM

The direct measurement of the EBL intensity was carried out by the special
camera in the space. After the subtraction of the contributions from the ZL
and other foreground components, the data show an obvious excess of EBL
photons compared to the indirectly measured EBL intensity. This indicates the
presence of an unknown diffuse emission component. However, the conclusion
remains uncertain since the estimations of the ZL contributions have large
uncertainties.

The NASA’s New Horizons spacecraft is an excellent platform for the EBL
observation since it is presently over 50 AU away from the sun, where the
zodiacal light is almost disappeared. A recent report from the New Horizons
confirms an EBL component of unknown origin of $n_{\gamma}=8.06\pm
1.92~{}nWm^{-2}sr^{-1}$.14 Although some of the authors have argued that the
EBL includes a substantial component of light from stars tidally removed from
galaxies, or a population of faint sources in extended halos. None of these
hypotheses may be correct.14 We try to explain this excess of the EBL photons
and to point out that it could be arisen from the contributions of the
anomalous bremsstrahlung when the VHE cosmic electron/positron fluxes pass
through the WHIM.

We consider a small range in the deep space, where the incident VHE cosmic
electron-positron rays interact with the Coulomb field of the WHIM and emit
soft bremsstrahlung photons. The screening length in Eq. (2.2) depends on the
ionization state and spatial distribution of the missing baryons in the WHIM.
If the probing range is small enough, the resolution of camera can not
distinguish the WHIM structure. Thus, we take an average screening length
$\mu^{-1}$ in Eq. (2.2) to describe the probed WHIM. The contributions of the
anomalous bremsstrahlung to the EBL intensity is

$n_{\gamma}(\omega)d\omega=\int_{E_{min}}^{E_{max}}n_{\gamma}(E,\omega)dEd\omega=\int_{E_{min}}^{E_{max}}J_{e}(E)P(E,\omega)dEd\omega.$
$None$

where $J_{e}(E)$ is the electron-positron flux at energy $E$. We take
$E_{min}=0.1~{}TeV$ and $E_{max}=10~{}TeV$ according to the Ref. [17].
$P(E,\omega)d\omega/\omega$ is the photon number emitted by the electrons in
the range $\omega\rightarrow\omega+d\omega$. The relation between the emission
power spectrum $P(E,\omega)d\omega$ and the bremsstrahlung cross section is
written as [18],

$\sigma_{brem}^{MD}(E,\omega)d\omega=\frac{P(E,\omega)d\omega}{n_{p}v_{e}\omega},$
$None$

where $n_{p}\sim 1/m^{3}$ is the average density of the missing baryons in the
WHIM and $v_{e}\simeq c$.

We use the parameterized distribution function

$J_{e}(E)=10^{-2}\left(\frac{E}{100GeV}\right)^{-2}[m^{-2}s^{-1}sr^{-1}],$
$None$

to model the electron-positron flux in the energy range
$100~{}GeV<E<1~{}TeV$17 and assume that this estimation is available in the
intergalactic space, although there is a large uncertainty of the spectra
index in the TeV range.

Note that the azimuth angles in both $n_{\gamma}(\omega)$ and $J_{e}(E)$ are
almost coincide, since the soft photons are concentrated in the moving
direction of the high energy electrons. Thus, we have,

$n_{\gamma}(\omega)=\int_{E_{min}}^{E_{max}}J_{e}(E)\frac{4\alpha_{s}^{3}}{\mu^{2}}\ln\frac{2E}{\mu}n_{p}v_{e}dE.$
$None$

This result shows a flat distribution and its height sensitively depends on
the value of $\mu$. The intergalactic electron intensity may be different from
Eq. (3.3), but it will not cover up a large anomalous bremsstrahlung effect.

In Fig. 2 we show our resulting $n_{\gamma}(\omega)\sim\lambda$ (dashed curve)
with $\mu=10^{-15}~{}GeV$ or $\mu^{-1}=10~{}cm$. Dotted curve is a global
fitting of the EBL intensity data from the indirect method.19 The sum of the
above two components (solid curve) is consistent with the data from New
Horizons spacecraft (solid point).14 We also present some other directly
measured EBL intensities (all other points),20-22 where the ZL contributions
are deducted by using the special models rather than the measurements in the
observations. Although our estimation is still preliminary, one can find that
the contribution from the anomalous bremsstrahlung to the EBL density seems
possible. On the other hand, the large excess around $\lambda\sim 1~{}\mu m$
was suspected to be from the highly redshifted ultraviolet radiation of the
first generation of stars. However, the evolution models based on the
observations of high-redshift galaxies indicate that the brightness caused by
first stars are much lower than the observed excesses.23 The increment of the
EBL due to the anomalous bremsstrahlung may compensate this deficiency.

Note that the solution of Eq. (3.4) can not be applied without restrictions.
The applications of the anomalous bremsstrahlung requests very low photon
energy so that neglecting the recoil effect is safe. On the other hand, the
radiation corrections are needed when the photon energy $\rightarrow 0$ in Eq.
(3.4).24 The infrared cutoff should be used to stop the infinite increasing of
$\sigma^{MD}\propto 1/\omega$ due to the negative corrections of the virtual
processes. The radiation corrections depend on the experimental conditions
(energy and angular resolutions). Thus, we assume that the intensity
$n_{\gamma}\propto\omega\sigma^{MD}$ is suppressed when the photon energy
$\omega$ is below the resolution of the detector. Therefore, we restrict the
solution of Eq. (3.4) to be in the range of $\lambda\approx 0.1-3~{}\mu m$,
which is shown in Fig. 2.

According to the Bethe-Heitler formula (2.1), Eq. (3.4) should be replaced by

$n_{\gamma}(\omega)=\int_{E_{min}}^{E_{max}}J_{e}(E)\frac{4\alpha_{s}^{3}}{m_{e}^{2}}\ln\frac{2E}{\mu}n_{p}v_{e}dE,$
$None$

and we have $n_{\gamma}\sim 10^{-20}~{}nWm^{-2}sr^{-1}\ll
10~{}nWm^{-2}sr^{-1}$. This is why the contributions of the normal
bremsstrahlung photons are neglected.

##  4 Anomalous pair production in the WHIM

An indirect measurement of the EBL intensity is based on the absorption of
gamma rays via the pair production off real EBL photons. This effect leads to
the distance-dependent suppression of the gamma-ray flux from active galactic
nuclei (AGNs).10,13 AGNs are the extragalactic objects characterized with
extremely luminous electromagnetic radiations. The missing baryons in the WHIM
may deform the gamma-ray distributions due to the anomalous pair production
process described by Eq. (2.4). We begin from a case without the WHIM. The EBL
may modify the propagation of VHE gamma-rays traveling through intergalactic
space due to pair production in the channel
$\gamma_{VHE}+\gamma_{EBL}\rightarrow e^{-}+e^{+}$ (Fig. 1a). The observed VHE
gamma-ray flux on the Earth is related to the initial flux at the source as,

$\Phi_{\gamma}^{ob}(E_{\gamma})=\Phi_{\gamma}^{in}(E_{\gamma})e^{-\tau_{\gamma\gamma}},$
$None$

where the optical depth $\tau_{\gamma\gamma}$ relates to the pair production
cross section of two photons $\sigma_{\gamma\gamma}$ and the density
$n_{\gamma}$ of the soft photon in the EBL as,

$\tau_{\gamma\gamma}(E_{\gamma})=\int_{0}^{R}\int_{\lambda_{min}}^{\lambda_{max}}\sigma_{\gamma\gamma}(E_{\gamma},\lambda)n_{\gamma}(\lambda,r)d\lambda
dr,$ $~{}~{}~{}~{}around~{}E_{\gamma}\sim 1~{}TeV.$ $None$

This is the foundation of the indirect method of the EBL measurement.

If we considering the WHIM, the contributions of pair production in the
nuclear Coulomb field of the WHIM (figure 1c) will be added to the observed
optical-depth data in Eq. (4.2). Therefore, the measured optical depth based
on the pair production model should include the contribution
$\gamma+p\rightarrow e^{+}+e^{-}$ due to the WHIM, i.e.,
$\tau_{\gamma\gamma}\rightarrow\tau_{\gamma\gamma}+\tau_{\gamma p}$. The
optical depth in the WHIM is evaluated with the anomalous pair production as,

$\tau_{\gamma
p}(E_{\gamma})=\int_{0}^{R}\sigma_{pair}^{MD}(E_{\gamma},r)n_{p}(E_{\gamma},r)dr,$
$None$

where $n_{p}$ is the proton density in the WHIM.

We have24 $\sigma_{\gamma p}\sim 0.1\sigma_{\gamma\gamma}$ around $E\sim
1~{}TeV$ with $n_{p}/n_{\gamma}\sim 10^{-9}-10^{-10}$ and
$\sigma_{pair}^{MD}\sim(m_{e}/m_{\mu})^{2}\sigma_{Pair}$. To have the similar
magnitude of gamma-ray attenuation from the WHIM as that from the EBL,
$\tau_{\gamma p}\sim\tau_{\gamma\gamma}$, the screening length of the WHIM is
required to be,

$\left(\frac{m_{e}}{\mu}\right)^{2}=10^{10}-10^{11}.$ $None$

Note that the counting-galaxy method is irrelevant to the WHIM structure since
the virtual photons in the Coulomb field are not counted in this method. The
EBL data deduced from the counting galaxies and the gamma-ray attenuation
method coincide with each other.13 Therefore, the correction to the optical
depth from the anomalous pair production on the WHIM should be neglected. This
implies that $(m_{e}/\mu)^{2}\ll 10^{10}-10^{11}$, or $\mu^{-1}\ll
10^{-5}~{}cm$ in Eq. (4.3). We will discuss these results in the next section.

##  5 Discussions

Both the EBL and the WHIM are the products of the universe evolution, and they
convey some important information of the cosmic formation history. Many works
studied the evolution of the optical depth $\tau_{\gamma\gamma}$.12,13 In this
work, we focus on the evolution of the properties of the WHIM with the
redshift, based on the optical-depth correction from the WHIM with the
improved Bethe-Heitler formula.

The screening length $\mu^{-1}$ is sensitively dependent on the ionization
status of the missing baryons in the WHIM. We take two extreme examples to
illustrate it. (i) If the missing baryons exist in the form of neutral atoms,
we have $\mu^{-1}\sim 10^{-8}~{}cm$ or $n_{\gamma}\sim
10^{-15}~{}nWm^{-2}sr^{-1}$ calculated with Eq. (3.4); (ii) The missing
baryons and electrons are uniformly distributed in the space as the plasma
form, where the maximum of the average screening length is
$\mu^{-1}=(1/n_{p})^{1/3}\sim 10^{2}~{}cm$, and we have $n_{\gamma}\sim
10^{5}~{}nWm^{-2}sr^{-1}$. These two examples show that as we use the average
screening length to describe the ionization property of the atom gas in the
WHIM, its value has a broadly reasonable range from the beginning of re-
ionization (i.e., $\mu^{-1}>10^{-8}~{}cm$) to a maximally ionized state of the
baryon-electron plasma.

The bremsstrahlung energy loss relates to the elemental QED processes and a
concrete ionization model. The later may provides a process-dependent
parameter $\mu$. In the above discussions we give a possible range of the
parameter $\mu$ in the WHIM according to the EBL data. We noticed that there
are numerous ionization models are devoted to describe bremsstrahlung energy
loss in interstellar medium.25 It is necessary to improve our estimation using
these models but based on the improved Bethe-Heitler formula. We will study
them in future work. On the other hand, as work [15] have pointed out that the
anomalous bremsstrahlung effect origins from a fact, that the screening
parameter in a logarithmic factor of the cross section moves to a denominator
position if the recoil effect can be neglected at the high energy
approximation. We believe that the application of the ionization model will
not inhibit a big correction of the anomalous bremsstrahlung in the WHIM,
although the estimation may become more accurate.

The standard cosmological model predicts that the missing baryons in the
universe evolution are concentrated in the WHIM due to continuously shock-
heated process.1,2 Cosmological simulations find that the WHIM contains a
large fraction of the baryons at redshift $z\simeq 0$ in the form of highly-
ionized gas.26 Theoretically, an atom which is bounded in an infinite heavy
web can completely absorb the recoil effect. However, this bounded atom can
not avoid the obvious recoil corrections due to the intense collisions at a
short impact parameter. Therefore, these bounded atoms in the WHIM do not have
the anomalous bremsstrahlung or the anomalous pair production. We consider the
atomic gas inside the WHIM at high redshift. The direct EBL intensity is
measured by using a long time exposition of the camera at low redshift $z\sim
0$. We have extracted a large screening length $\mu^{-1}\simeq 10~{}cm$ at
$z\simeq 0$ in Sec. 3. It implies that a lot of atoms in the WHIM at
$z\rightarrow 0$ are almost re-ionized in the complete baryon-electron plasma
state. This result is consistent with the observations, which gives the
evidence that a significant portion of the baryons are located at $z\simeq 0$
in the WHIM.27 On the other hand, the gamma rays in the indirect measurement
of the EBL intensity come from the AGNs at large redshift. The contributions
to the optical depth from highly-ionized gas in the WHIM at low redshift only
account for a small part of the total optical depth in Eq. (4.3). Therefore,
we have the average screening length $\mu^{-1}\ll 10^{-5}~{}cm$ for the
anomalous pair production discussed in section 4, which leads to the result
$\tau_{\gamma p}\ll\tau_{\gamma\gamma}$. Thus, we explained why we use two
different values of the parameter $\mu^{-1}$ in Sec. 3 and Sec. 4.

Usually, the absorptions or emissions of far-ultraviolet and soft X-ray
photons.28 or the hyperfine transition of neutral hydrogen (called 21-cm
line)29 are used to look for the signals of the missing baryons in the WHIM.
However, the WHIM is extremely difficult to be observed, for its tremendously
low density and high temperature are believed to elude the detections of the
absorption and emission signals. Therefore, the properties of the WHIM at
$z=0$ is still elusive with the traditional detecting methods. Comparing with
the characteristic spectra of the WHIM, the sensitive dependence of the
improved Bethe-Heitler cross section on the average screening length provides
a new method to complement the study of the WHIM structure.

In summary, we try to establish a connection between the WHIM and the EBL. The
found relation answers the following cosmological questions: Where are the
missing baryons? Is the sky absolutely dark if all star lights are removed?
According to an improved Bethe-Heitler formula, we find that the WHIM may be
lighted up by the VHE cosmic electron/positron fluxes, and this effect can be
observed by the directly measured EBL intensity. Using a sensitive relation
between the EBL intensity and the WHIM electronic structure, the improved
Bethe-Heitler formula may open a new window to study the WHIM structure.

##  6  Acknowledgments

This work is partly supported by the National Natural Science Foundation of
China No.11851303 and No.12005266.

## References

  * [1] C.L. Bennett et al.,_Astrophys. J. Suppl._ $\bf{148}$ (2003) 1.
  * [2] D.N. Spergel et al., _Astrophys. J. Suppl._ $\bf{148}$ (2003) 175.
  * [3] M. Fukugita, arXiv:0312517.
  * [4] C.W. Danforth and J.M. Shull, _Astrophys. J._ $\bf{624}$ (2005) 555.
  * [5] R. Cen and J.P. Ostriker, _Astrophys. J._ $\bf{514}$ (1999) 1.
  * [6] R. Dave et al., _Astrophys. J._ $\bf{552}$ (2001) 473.
  * [7] E. Behar, S. Dado, A. Dar and A. Laor, _IL NUOVO CIMENTO._ $\bf{33}$ (2010) 55.
  * [8] S. Walsh et al., arXiv/2007.10158.
  * [9] J.P. Macquart et al., _Nature._ $\bf{591}$ (2020) 391.
  * [10] F.A. Aharonian, Very high energy cosmic gamma radiation: A crucial window on the extreme universe, Word Scientific Co. Pte. Ltd., Singapore, (2004)
  * [11] K. Mattila and P. Vaisanen, arXiv/1905.08825.
  * [12] A. Dominguez et al., _Mon. Not. R. Astron. Soc._ $\bf{410}$ (2011) 2556.
  * [13] A. Desai et al., arXiv/1903.03126.
  * [14] T.R. Lauer et al., arXiv/2202.04273.
  * [15] W. Zhu, _Nucl. Phys._ $\bf{B953}$ (2020) 114958.
  * [16] H. Bethe and W. Heitler, _Proc. Roy. Soc., London._ $\bf{A146}$ (1934) 83.
  * [17] W. Zhu, P. Liu, J.H. Ruan, L. Feng and F. Wang, _Astrophys. J._ $\bf{896}$ (2020) 106.
  * [18] Jun Han You, Radiation Mechanism in Astronomical Physics, Science Press (China) (1998)
  * [19] W. Zhu, Z.C. Zheng, P. Liu, L.H. Wan, J.H. Ruan and F. Wang, _JCAP._ $\bf{01}$ (2021) 038.
  * [20] Y. Matsuoka et al.,_Astrophys. J._ $\bf{736}$ (2011) 119.
  * [21] K. Mattila et al., _MNRAS._ $\bf{470}$ (2017) 2152, doi:10.1093/mnras/stx1296.
  * [22] S. Matsuura et al., _Astrophys. J._ $\bf{839}$ (2017) 7.
  * [23] A. Cooray, Y. Gong, J. Smidt and N.G. Santos, _Astrophys. J._ $\bf{756}$ (2012) 92.
  * [24] W. Greiner and J. Reinhardt, Quantum Electrodynamics, Springer 1994.
  * [25] G.R. Blumenthal and R.J. Gould, _Rev. Mod. Phys._ $\bf{42}$ (1970) 237.
  * [26] P. Richter, B. D. Savage, T. M. Tripp and K. R. Sembach, _PoS BDMH._ (2004) 002.
  * [27] J. N. Bregman, Ann. Rev. _Astron. Astrophys._ $\bf{45}$ (2007) 221.
  * [28] Simon D. M. White, _Rept. Prog. Phys._ $\bf{70}$ (2007) 883.
  * [29] D. Martizzi et al., arXiv/1810.01883.

Figure 1: (a) The interaction of high energy photon with soft photon in the
EBL. (b) and (c) Bremsstrahlung and pair creation in the nuclear Coulomb
field, which can be described with an improved Bethe-Heitler formula in the
large-scale medium of dilute ionized matter.  Figure 2: The contribution from
the anomalous bremsstrahlung (Fig. 1b) to the EBL intensity is presented with
the dashed line. Dotted curve is taken from,19 which is fitted to the
indirectly measured EBL intensity. Solid curve is the sum of the above two
components. The filled circle with error bar is a recent data from New
Horizons,14 which does not have the ZL component. All other data points are
taken from some different direct EBL measurements20-22 but with the model-
dependent corrections for ZL.


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

Fernando S Moura Engineering, modelling and Applied Social Sciences Center,
Federal University of ABC
São Bernardo do Campo, São Paulo, Brazil Department of Mathematics and
Statistics, University of Helsinki, Helsinki, Finland Corresponding author.
Email<EMAIL_ADDRESS>Roberto G Beraldo Engineering, modelling
and Applied Social Sciences Center, Federal University of ABC
São Bernardo do Campo, São Paulo, Brazil Leonardo A Ferreira Engineering,
modelling and Applied Social Sciences Center, Federal University of ABC
São Bernardo do Campo, São Paulo, Brazil Samuli Siltanen Department of
Mathematics and Statistics, University of Helsinki, Helsinki, Finland

###### Abstract

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
$125\text{\,}\mathrm{dB}$ to identify vascular changes due to the cardiac
cycle, corroborating previous studies.

Keywords: Anatomical atlas, cerebral circulation, electrical properties, human
head model, electrophysiology, electrical impedance tomography,
electroencephalography

##  1 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 and Brunet 2019; Jatoi and Kamel 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 evaluated using magnetic resonance imaging (Vikner et al.
2019; Holmgren 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 (Pesti et al. 2019), 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
(Min et al. 2019), 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 et al. 2017; 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. Anatomy-based priors are found in the
literature, such as for thorax applications (Martins et al. 2019; Kaipio and
Somersalo 2005).

##  2 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

$\displaystyle A_{s}$ $\displaystyle\sim N(\bar{\bm{x}}_{s},\bm{\Gamma}_{s})$
(1) $\displaystyle A_{d}(t)$ $\displaystyle\sim
N(\bar{\bm{x}}_{d}(t),\bm{\Gamma}_{d}(t))$ (2) $\displaystyle A(t)$
$\displaystyle\sim
N(\bar{\bm{x}}_{s}+\bar{\bm{x}}_{d}(t),\bm{\Gamma}_{a}+\bm{\Gamma}_{d}(t)).$
(3)

The two components of the atlas are presented in details in the following
subsections.

###  2.1 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.

Figure 1: Main steps necessary to compute the static component of the
anatomical atlas. 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.

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\text{\times}1\text{\times}1\text{\,}\mathrm{m}\mathrm{m}$. An
equal number of male and female individuals were used, with an average age of
$40\pm 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. (2008, 2009).

Each of the 50 images was transformed aiming to maximize its similarity with a
reference image. The reference image is the MNI ICBM 152111Available at
http://nist.mni.mcgill.ca/?p=858, a nonlinear symmetric atlas with
$1\text{\times}1\text{\times}1\text{\,}\mathrm{m}\mathrm{m}$ 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.

Figure 2: Slices of the reference image MNI ICBM 152.

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:
background, 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.

####  2.1.1 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 (Gabriel, Gabriel and Corthout 1996;
Gabriel, Lau and Gabriel 1996a).

Given the angular frequency of the electrical signal $\omega=2\pi f$ , the
complex relative permittivity $\hat{\epsilon}(\omega)$ of a tissue can be
modelled as the sum of four Cole-Cole dispersion terms (Gabriel, Lau and
Gabriel 1996b)

$\displaystyle\hat{\epsilon}(\omega)=\epsilon_{\infty}+\sum_{n=1}^{4}\frac{\Delta\epsilon_{n}}{1+(j\omega\tau_{n})^{1-\alpha_{n}}}+\frac{\sigma_{0}}{j\omega\epsilon_{0}},$
(4)

where $\epsilon_{0}$ is the permittivity of free space and all the other
parameters depend on the tissue (Gabriel, Lau and Gabriel 1996b; Andreuccetti
et al. 1997). The conductivity $\sigma$ and permittivity $\epsilon$ of the
tissue can be obtained from $\hat{\epsilon}(\omega)$

$\displaystyle\sigma(\omega)$
$\displaystyle=-\omega\epsilon_{0}\text{Im}(\hat{\epsilon}(\omega))$ (5)
$\displaystyle\epsilon(\omega)$
$\displaystyle=\epsilon_{0}\text{Re}(\hat{\epsilon}(\omega)).$ (6)

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 $\pm 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.

####  2.1.2 Atlas statistics computation

Let $\bm{u}\in\mathbb{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
$\bm{\chi}_{t}\in\mathbb{R}^{N_{V}}$, for $t=1,2,\cdots,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
$\bm{x}\in\mathbb{C}^{N_{T}}$ (or real with the same dimension) as

$\displaystyle\bm{x}=\sum_{t=1}^{N_{T}}p_{t}\bm{\chi}_{t}=\bm{X}\bm{p},$ (7)

where $\bm{p}\in\mathbb{C}^{N_{T}}$ is a vector composed by the electrical
properties of the tissues and $\bm{X}\in\mathbb{R}^{N_{V}\times 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 $\bm{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 $\bm{\Gamma}$ in special, the
formulation allows the computation in factorized form
$\bm{\Gamma}=\bm{K}\bm{K}^{\mathrm{T}}$, reducing storage requirements and
simplifying algorithms that depend on factorizations of $\bm{\Gamma}$.

Samples of the $i$-th individual can be composed by sampling the properties of
the tissues $\bm{p}_{s}$ and applying (7)

$\displaystyle\bm{x}_{s,i}=\bm{X}_{i}\bm{p}_{s},\quad\text{for}\quad
s=1,2,\cdots,N_{S},$ (8)

where the samples $\bm{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 $\bm{x}_{s,i}$ represent a sample of the $i$-th individual. The average
over all individuals can be estimated with

$\bar{\bm{x}}=\frac{1}{N_{S}N_{I}}\sum_{i=1}^{N_{I}}\sum_{s=1}^{N_{S}}\bm{x}_{s,i}=\frac{1}{N_{I}}\sum_{i=1}^{N_{I}}\bm{X}_{i}\frac{1}{N_{S}}\sum_{s=1}^{N_{S}}\bm{p}_{s}=\frac{1}{N_{I}}\sum_{i=1}^{N_{I}}\bm{X}_{i}\bar{\bm{p}},$
(9)

where, again, $N_{S}$ is fixed and represent the number of samples with the
same head geometry $\bm{X}_{i}$ and $\bar{\bm{p}}$ is the average electrical
properties of the tissues.

The covariance matrix can be computed using the usual sample estimator

$\bm{\Gamma}=\frac{1}{N_{I}N_{S}-1}\sum_{i=1}^{N_{I}}\sum_{s=1}^{N_{S}}(\bm{x}_{s,i}-\bar{\bm{x}})(\bm{x}_{s,i}-\bar{\bm{x}})^{\mathrm{H}},$
(10)

where $\bm{\Gamma}\in\mathbb{R}^{N_{V}\times N_{V}}$ and $\bm{M}^{\mathrm{H}}$
denotes conjugate transpose of $\bm{M}$. For real valued $\bm{p}$, the
conjugate transpose is the transpose $\bm{M}^{\mathrm{H}}=\bm{M}^{\mathrm{T}}$

Adding $(\bar{\bm{x}}_{i}-\bar{\bm{x}}_{i})$ to both terms between parenthesis
and rearranging the terms,

$\bm{\Gamma}=\frac{1}{N_{I}N_{S}-1}\sum_{i=1}^{N_{I}}\sum_{s=1}^{N_{S}}[(\bm{x}_{s,i}-\bar{\bm{x}}_{i})+(\bar{\bm{x}}_{i}-\bar{\bm{x}})][(\bm{x}_{s,i}-\bar{\bm{x}}_{i})^{\mathrm{H}}+(\bar{\bm{x}}_{i}-\bar{\bm{x}})^{\mathrm{H}}],$
(11)

where $\bar{\bm{x}}_{i}=\sum_{s=1}^{N_{S}}\bm{x}_{s,i}/N_{S}$ is the average
of the $i$-th individual. Proceeding with the products,

$\displaystyle\bm{\Gamma}$
$\displaystyle=\frac{1}{N_{I}N_{S}-1}\sum_{i=1}^{N_{I}}(\alpha_{1}+\alpha_{2}+\alpha_{3}+\alpha_{4})$
(12) $\displaystyle\alpha_{1}$
$\displaystyle=\sum_{s=1}^{N_{S}}(\bm{x}_{s,i}-\bar{\bm{x}}_{i})(\bm{x}_{s,i}-\bar{\bm{x}}_{i})^{\mathrm{H}}=(N_{S}-1)\bm{\Gamma}_{i}$
(13) $\displaystyle\alpha_{2}$
$\displaystyle=\sum_{s=1}^{N_{S}}(\bm{x}_{s,i}-\bar{\bm{x}}_{i})(\bar{\bm{x}}_{i}-\bar{\bm{x}})^{\mathrm{H}}=\left[\sum_{s=1}^{N_{S}}(\bm{x}_{s,i}-\bar{\bm{x}}_{i})\right](\bar{\bm{x}}_{i}-\bar{\bm{x}})^{\mathrm{H}}=0$
(14) $\displaystyle\alpha_{3}$
$\displaystyle=\sum_{s=1}^{N_{S}}(\bar{\bm{x}}_{i}-\bar{\bm{x}})(\bm{x}_{s,i}-\bar{\bm{x}}_{i})^{\mathrm{H}}=(\bar{\bm{x}}_{i}-\bar{\bm{x}})\left[\sum_{s=1}^{N_{S}}(\bm{x}_{s,i}-\bar{\bm{x}}_{i})^{\mathrm{H}}\right]=0$
(15) $\displaystyle\alpha_{4}$
$\displaystyle=\sum_{s=1}^{N_{S}}(\bar{\bm{x}}_{i}-\bar{\bm{x}})(\bar{\bm{x}}_{i}-\bar{\bm{x}})^{\mathrm{H}}=N_{S}(\bar{\bm{x}}_{i}-\bar{\bm{x}})(\bar{\bm{x}}_{i}-\bar{\bm{x}})^{\mathrm{H}}.$
(16)

Finally, taking the limit $N_{S}\rightarrow\infty$,

$\displaystyle\bm{\Gamma}$
$\displaystyle=\frac{1}{N_{I}}\sum_{i=1}^{N_{I}}\bm{\Gamma}_{i}+\frac{1}{N_{I}}\sum_{i=1}^{N_{I}}(\bar{\bm{x}}_{i}-\bar{\bm{x}})(\bar{\bm{x}}_{i}-\bar{\bm{x}})^{\mathrm{H}}.$
(17)

Defining $\Delta\bar{\bm{x}}_{i}=\bar{\bm{x}}_{i}-\bar{\bm{x}}$ and observing
the linear relation (7) we can rewrite this last expression as

$\displaystyle\bm{\Gamma}$
$\displaystyle=\frac{1}{N_{I}}\sum_{i=1}^{N_{I}}[\bm{X}_{i}\bm{\Gamma}_{p}\bm{X}_{i}^{\mathrm{H}}+\Delta\bar{\bm{x}}_{i}\Delta\bar{\bm{x}}_{i}^{\mathrm{H}}]=\sum_{i=1}^{N_{I}}\bm{W}_{i}\bm{W}_{i}^{\mathrm{H}}$
(18) $\displaystyle\bm{W}_{i}$
$\displaystyle=\frac{1}{\sqrt{N_{I}}}\begin{bmatrix}\bm{X}_{i}\sqrt{\bm{\Gamma}_{p}}&\Delta\bar{\bm{x}}_{i}\end{bmatrix},$
(19)

where $\bm{W}_{i}\in\mathbb{C}^{N_{V}\times(N_{T}+1)}$ and
$\bm{\Gamma}_{p}\in\mathbb{C}^{N_{T}\times N_{T}}$ is the covariance matrix of
the electrical properties of the tissues. Furthermore, the expression can be
simplified to

$\displaystyle\bm{\Gamma}$
$\displaystyle=\begin{bmatrix}\bm{W}_{1}&\cdots&\bm{W}_{I}\end{bmatrix}\begin{bmatrix}\bm{W}_{1}^{\mathrm{H}}\\\
\vdots\\\ \bm{W}_{I}^{\mathrm{H}}\end{bmatrix}=\bm{K}\bm{K}^{\mathrm{H}},$
(20)

where $\bm{K}\in\mathbb{C}^{N_{V}\times[(N_{T}+1)N_{I}]}$. 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 $\bm{p}$, the pseudo-covariance
$\bm{\tilde{\Gamma}}\in\mathbb{C}^{N_{V}\times N_{V}}$ can be computed in a
similar way

$\displaystyle\bm{\tilde{\Gamma}}$ $\displaystyle=\bm{K}\bm{K}^{\mathrm{T}}.$
(21)

The atlas in this work is assumed Gaussian, therefore $\bar{\bm{x}}$ and
$\bm{\Gamma}$ (and $\bm{\tilde{\Gamma}}$ for complex-valued $\bm{p}$)
completely specify its probability density function.

###  2.2 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\text{\times}0.5\text{\times}0.8\text{\,}\mathrm{m}\mathrm{m}$. An equal
number of male and female individuals were used, with an average age of $43\pm
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.

####  2.2.1 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 et al. 2017; Melis 2018), a 1D blood flow solver based on
monotonic upstream-centered 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.

Figure 3: Superior aortic system considered in the simulations. The names of
the vessels are presented in Table 1.

Blood was assumed Newtonian with density
$\rho=$1050\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$$ and dynamic
viscosity
$\mu=$4.5\text{\times}{10}^{-3}\text{\,}\mathrm{Pa}\text{\,}\mathrm{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

$\displaystyle Q(t)=\left\\{\begin{array}[]{ll}Q_{M}\sin(\frac{\pi
t}{\tau})&t<\tau\\\ 0&\text{otherwise}\end{array}\right.,$ (24)

where $Q_{M}=$485\text{\,}\mathrm{ml}\text{\,}{\mathrm{s}}^{-1}$$ is peak flow
rate, $\tau=$0.3\text{\,}\mathrm{s}$$ and the cardiac cycle period is
$1\text{\,}\mathrm{s}$, following Alastruey et al. (2007).

[t]

Table 1: Geometrical and mechanical properties of the arteries. Based on
Alastruey et al. (2007) and complemented with data collected by Dodo et al.
(2020); Fomkina et al. (2016); Schmitter et al. (2013). $\ell:$ length,
$r_{0}:$ Proximal lumen radius, $E:$ Young’s modulus, $R_{T}:$ Terminal
resistance, and $C_{T}:$ Terminal capacitance.

# |  Artery1 | $\ell$ | $r_{0}$ | $E$ | $R_{T}$ | $C_{T}$
---|---|---|---|---|---|---
|  | ($\mathrm{m}\mathrm{m}$) | ($\mathrm{m}\mathrm{m}$) | ($\mathrm{k}\mathrm{P}\mathrm{a}$) | ($10^{9}$$\mathrm{Pa}\text{\,}\mathrm{s}\text{\,}{\mathrm{m}}^{-3}$) | ($10^{-10}$${\mathrm{m}}^{3}\text{\,}{\mathrm{Pa}}^{-1}$)
01 | Ascending aorta |  040 | 12.00 |  0400 | – | –
02 | Aortic arch I |  020 | 11.20 |  0400 | – | –
03 | Brachiocephalic |  034 |  06.20 |  0400 | – | –
04 | Aortic arch II |  039 | 10.70 |  0400 | – | –
05 | L CCA | 208 |  02.50 |  0400 | – | –
06 | R CCA | 177 |  02.50 |  0400 | – | –
07 | R Subclavian |  034 |  04.23 |  0400 | – | –
08 | Thoracic aorta | 156 |  09.99 |  0400 |  00.18 | 38.70
09 | L Subclavian |  034 |  04.23 |  0400 | – | –
10/13 | ECA | 177 |  01.50 |  0800 |  05.43 |  01.27
11/12 | ICA I | 177 |  02.00 |  0800 | – | –
14/17 | Vertebral | 148 |  01.36 |  0800 | – | –
15/16 | Brachial | 422 |  04.03 |  0400 |  02.68 |  02.58
18/21 | ICA II |  005 |  02.00 | 1600 | – | –
19/20 | PCoA |  015 |  00.73 | 1600 | – | –
22 | Basilar I |  025 |  01.62 | 1600 | – | –
23/24 | MCA | 119 |  01.43 | 1600 |  05.97 |  01.16
25/26 | ACA I |  012 |  01.17 | 1600 | – | –
27/28 | PCA I |  005 |  01.07 | 1600 | – | –
29/30 | ACA II | 103 |  01.20 | 1600 |  08.48 |  00.82
31 | ACoA |  003 |  00.74 | 1600 | – | –
32/33 | PCA II |  086 |  01.05 | 1600 | 11.08 |  00.62
34 | Basilar II |  001 |  01.62 | 1600 | – | –
35 | Basilar III |  003 |  01.62 | 1600 | – | –
36/37 | SCA |  086 |  00.65 | 1600 | 25.00 |  00.62

  * 1

Acronyms:Anterior cerebral artery (ACA), Anterior communicating artery (ACoA),
Common carotid artery (CCA), External carotid artery (ECA), Internal carotid
artery (ICA), Middle cerebral artery (MCA), Posterior cerebral artery (PCA),
Posterior communicating artery (PCoA), Superior cerebellar artery (SCA).

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)

$\displaystyle\frac{\Delta\rho_{\ell}}{\rho_{0}}=-0.45H\left[1-\exp\left(-0.26\left|\frac{\bar{\bm{v}}}{R}\right|^{0.39}\right)\right],$
(25)

where $\Delta\rho_{\ell}$ is the longitudinal resistivity change to the
reference (still blood) resistivity $\rho_{0}$, $H$ is the hematocrit (volume
percentage of red blood cells in the blood), $\bar{\bm{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 $\Delta\rho_{\ell}/\rho_{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 (25), 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.

Figure 4: Main cerebral vascular territories. Acronyms: Anterior cerebral
artery (ACA), Middle cerebral artery (MCA), Posterior cerebral artery (PCA),
Superior cerebellar artery (SCA).

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).

####  2.2.2 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 $\bm{\chi}_{m}\in\mathbb{R}^{N_{V}}$ and its area of influence
$\bm{\chi}^{c}_{m}\in\mathbb{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
$\bm{x}(t)\in\mathbb{C}^{N_{T}}$ as

$\displaystyle\bm{x}(t)=\sum_{m=1}^{N_{M}}p_{m,B}(t)\bm{\chi}_{m}+p_{m,B}^{c}(t)\bm{\chi}^{c}_{m},$
(26)

where $p_{m,B}(t)$ and $p_{m,B}^{c}(t)$ are the electrical properties of blood
in the respective segments. The effect of blood in the complements
$p_{m,B}^{c}(t)$ is modelled as a percentage of $p_{m,B}(t)$

$\displaystyle\bm{x}(t)=\sum_{m=1}^{N_{M}}p_{m,B}(t)(\bm{\chi}_{m}+\alpha_{m}\bm{\chi}^{c}_{m})=\sum_{m=1}^{N_{M}}p_{m,B}(t)\bm{\chi}_{\alpha,m},$
(27)

where $\alpha_{m}\geq 0$ adjusts the effect. The vector $\bm{\chi}_{\alpha,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 $\bm{\chi}_{\alpha,m}$.

##  3 Results

###  3.1 Static component of the atlas

Figure 5(a) 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 5(b) 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.

(a)

(b)

Figure 5: Static atlas computations. (a) Processing steps of two
representative individuals of the dataset. The first row shows the original
images, the second presents the results of the normalization and the third row
the result of the segmentation. GM: red, WM: blue, CSF: green, BO: white, and
OT: grey. (b) Average of the characteristic functions $\bm{\chi}_{t}$ for each
tissue.

The static component of the atlas at $1\text{\,}\mathrm{kHz}$ is presented in
Figure 6. The figure shows transversal slices of the average (Figure 6(a)) and
standard deviation (Figure 6(b)). 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.

(a)

(b)

Figure 6: Transversal slices of the static component of the atlas
(resistivity) at $1\text{\,}\mathrm{kHz}$. (a) average; (b) standard
deviation.

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.

(a)

(b)

Figure 7: Statistics of the atlas in different frequencies. (a) average; (b)
standard deviation.

###  3.2 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\text{\,}\mathrm{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\text{\,}\mathrm{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 (25).
Peak velocity occurs approximately $0.25\text{\,}\mathrm{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.

Figure 9: Waveforms of the main arteries over one cardiac cycle. From left to
right: average cross-sectional velocity, static pressure and resistivity
changes.

###  3.3 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\text{\,}\mathrm{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
adult222https://www.pedeheadmod.net/pediatric-head-atlases/ 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 10(a), is composed of 2.06 million linear
tetrahedral elements, split into six segments and 32 electrodes with a
diameter of $15\text{\,}\mathrm{mm}$ were placed in two parallel planes with
$20\text{\,}\mathrm{mm}$ of separation. Electrode numbers can be seen in the
figure.

(a)

(b)

Figure 10: Finite element mesh used to simulate EIT measurements. (a) Internal
structure and electrode locations; (b) Slices of the projected atlas into the
FEM mesh.

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
$\bm{\chi}_{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%

$\displaystyle\bar{\bm{\chi}}$
$\displaystyle=E\\{\bm{\chi}_{WM}+\bm{\chi}_{GM}+\bm{\chi}_{CSF}\\}$ (28)
$\displaystyle\bm{\chi}_{C}(i,j,k)$ $\displaystyle=\begin{cases}1&\text{if
}\bar{\bm{\chi}}(i,j,k)\geq 0.75\\\ 0&\text{otherwise}\end{cases},$ (29)

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 10(a), 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 10(b).

EIT measurement simulation was performed by imposing sinusoidal bipolar
current injection of $1\text{\,}\mathrm{mA}$ at $1\text{\,}\mathrm{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
$\bm{v}_{0}$; (b) relative differences $\Delta\bm{v}(t)$; (c) first 32
relative differences $\Delta\bm{v}(t)$; (d) histograms of the differences
$\Delta\bm{v}(t)$.

Figure 11a presents the measurements at $t=$0\text{\,}\mathrm{s}$$, used as
reference measurement $\bm{v}_{0}$. Figure 11b shows relative differences
$\Delta\bm{v}(t)$, in dB, between measurements in three other time instants
$\bm{v}(t)$ and the reference $\bm{v}_{0}$, defined as

$\displaystyle\Delta\bm{v}(t)=20\log_{10}\left(\left|\frac{\bm{v}(t)-\bm{v}_{0}}{\bm{v}_{0}}\right|\right),$
(30)

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\text{\,}\mathrm{s}$$, time instant when
we have peak velocities, as presented before in Figure 9. Figure 11d presents
histograms of the differences $\Delta\bm{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).

##  4 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 10(b) 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\text{\,}\mathrm{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.

##  5 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.

## Acknowledgements

The work was funded in part by the Jane and Aatos Erkko Foundation, project
“Electrical Impedance Tomography—a novel method for improved diagnosis of
stroke”, the Academy of Finland (Centre of Excellence in Inverse Modelling and
Imaging, decision number 312339), the São Paulo Research Foundation–FAPESP
(Process numbers: 2019/09154-7 and 2017/18378-0) and Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior (CAPES) – Finance Code 001.

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)

## References

  * (1)
  * Acar and Makeig (2013) Acar, Z. A. and Makeig, S. (2013).  Effects of forward model errors on EEG source localization, Brain Topography 26(3): 378–96.
  * Adler and Boyle (2019) Adler, A. and Boyle, A. (2019).  Electrical impedance tomography, in J. G. Webster (ed.), Wiley Encyclopedia of Electrical and Electronics Engineering, 2 edn, John Wiley $\&$ Sons, Inc, New Jersey, United States, pp. 1–16.
  * Agnelli et al. (2020) Agnelli, J. P., Çöl, A., Lassas, M., Murthy, R., Santacesaria, M. and Siltanen, S. (2020).  Classification of stroke using neural networks in electrical impedance tomography, Inverse Problems 36(11): 115008.
  * Alastruey et al. (2007) Alastruey, J., Parker, K. H. and Peiró, J. (2007).  Modelling the circle of willis to assess the effects of anatomical variations and occlusions on cerebral flows., Journal of Biomechanics 40(8): 1794–805.
  * Andreuccetti et al. (1997) Andreuccetti, D., Fossi, R. and Petrucci, C. (1997).  An internet resource for the calculation of the dielectric properties of body tissues in the frequency range 10 hz - 100 ghz.  Based on data published by C.Gabriel et al. in 1996.
  * Ashburner and Friston (2005) Ashburner, J. and Friston, K. (2005).  Unified segmentation., Neuroimage 26(3): 839–51.
  * Avants et al. (2008) Avants, B., Grossman, C. E. M. and Gee, J. (2008).  Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain, Medical Image Analysis 12(1): 26–41.
  * Avants et al. (2009) Avants, B., Tustison, N., Song, G. and Gee, J. (2009).  Ants: Open-source tools for normalization and neuroanatomy., HeanetIe 10: 1–11.
  * Badeli et al. (2020) Badeli, V., Melito, G. M., Köstinger, A. R., Bíró, O. and Ellermann, K. (2020).  Electrode positioning to investigate the changes of the thoracic bioimpedance caused by aortic dissection – a simulation study, Journal of Electrical Bioimpedance 11(1): 38–48.
  * Beraldo and Moura (2020) Beraldo, R. G. and Moura, F. S. (2020).  Time-difference electrical impedance tomography with a blood flow model as prior information for stroke monitoring, Proceedings of the XXVII Brazilian Congress on Biomedical Engineering, SBEB, Vitória, Brazil, pp. 280–5.
  * Bernstein (2010) Bernstein, D. P. (2010).  Impedance cardiography: Pulsatile blood flow and the biophysical and electrodynamic basis for the stroke volume equations, Journal of Electrical Bioimpedance 1(1): 2–17.
  * Bodo et al. (2018) Bodo, M., Montgomery, L. D., Pearce, F. J. and Armonda, R. (2018).  Measurement of cerebral blood flow autoregulation with rheoencephalography: a comparative pig study, Journal of Electrical Bioimpedance 9(1): 123–32.
  * Bradac (2017) Bradac, G. B. (2017).  Vascular Territories, Springer International Publishing, Cham, pp. 105–7.
  * Braun et al. (2018) Braun, F., Proença, M., Lemay, M., Bertschi, M., Adler, A., Thiran, J.-P. and Solà, J. (2018).  Limitations and challenges of EIT-based monitoring of stroke volume and pulmonary artery pressure, Physiological Measurement 39(1): 014003.
  * Bullitt et al. (2005) Bullitt, E., Zeng, D., Gerig, G., Aylward, S., Joshi, S., Smith, J. K., Lin, W. and Ewend, M. G. (2005).  Vessel tortuosity and brain tumor malignancy: A blinded study, Academic Radiology 12: 1232–40.
  * Candiani et al. (2019) Candiani, V., Hannukainen, A. and Hyvonen, N. (2019).  Computational framework for applying electrical impedance tomography to head imaging, SIAM Journal on Scientific Computing 41(5): B1034–60.
  * Candiani and Santacesaria (2020) Candiani, V. and Santacesaria, M. (2020).  Neural networks for classification of strokes in electrical impedance tomography on a 3d head model, arXiv preprint arXiv:2011.02852 pp. 1–20.
  * Chandra et al. (2017) Chandra, A., Li, W., Stone, C., Geng, X. and Ding, Y. (2017).  The cerebral circulation and cerebrovascular disease I: Anatomy, Brain Circulation 3(2): 45–56.
  * Cheney et al. (1999) Cheney, M., Isaacson, D. and Newell, J. C. (1999).  Electrical impedance tomography, SIAM Review 41(1): 85–101.
  * Cho et al. (2015) Cho, J.-H., Vorwerk, J., Wolters, C. H. and Knösche, T. R. (2015).  Influence of the head model on EEG and MEG source connectivity analyses, NeuroImage 110: 60–77.
  * Desmidt et al. (2018) Desmidt, T., Andersson, F., Brizard, B., Dujardin, P.-A., Cottier, J.-P., Patat, F., Réméniéras, J.-P., Gissot, V., El-Hage, W. and Camus, V. (2018).  Ultrasound measures of brain pulsatility correlate with subcortical brain volumes in healthy young adults, Ultrasound in Medicine & Biology 44(11): 2307–13.
  * Dodo et al. (2020) Dodo, Y., Takahashi, T., Honjo, K., Kitamura, N. and Maruyama, H. (2020).  Measurement of the length of vertebrobasilar arteries: A three-dimensional approach, Journal of the Neurological Sciences 414: 116818.
  * Dowrick et al. (2016) Dowrick, T., Blochet, C. and Holder, D. (2016).  In vivo bioimpedance changes during haemorrhagic and ischaemic stroke in rats: towards 3d stroke imaging using electrical impedance tomography, Physiological Measurement 37(6): 765–84.
  * Driver et al. (2020) Driver, I. D., Traat, M., Fasano, F. and Wise, R. G. (2020).  Most small cerebral cortical veins demonstrate significant flow pulsatility: A human phase contrast mri study at 7t, Frontiers in Neuroscience 14: 415.
  * Fomkina et al. (2016) Fomkina, O. A., Nikolenko, V. N. and Chernyshkova, E. V. (2016).  Morphology and biomechanical properties of cerebellar arteries in adults, Russian Open Medical Journal 5(2): e0205.
  * Fonov et al. (2009) Fonov, V., Evans, A., McKinstry, R., Almli, C. and Collins, D. (2009).  Unbiased nonlinear average age-appropriate brain templates from birth to adulthood., NeuroImage 47(1): S102.
  * Frackowiak et al. (2004) Frackowiak, R. S., Friston, K. J., Frith, C. D., Dolan, R. J., Price, C. J., Zeki, S., Ashburner, J. T. and Penny, W. D. (2004).  Chapter 31 - experimental design and statistical parametric mapping, Human Brain Function (Second Edition), second edition edn, Academic Press, Burlington, pp. 599–632.
  * Frerichs et al. (2016) Frerichs, I., Amato, M. B. P., van Kaam, A. H., Tingay, D. G., Zhao, Z., Grychtol, B., Bodenstein, M., Gagnon, H., Böhm, S. H., Teschner, E., Stenqvist, O., Mauri, T., Torsani, V., Camporota, L., Schibler, A., Wolf, G. K., Gommers, D., Leonhardt, S. and Adler, A. (2016).  Chest electrical impedance tomography examination, data analysis, terminology, clinical use and recommendations: consensus statement of the TRanslational EIT developmeNt stuDy group, Thorax 72(1): 83–93.
  * Gabriel, Gabriel and Corthout (1996) Gabriel, C., Gabriel, S. and Corthout, E. (1996).  The dielectric properties of biological tissues: I. literature survey, Physics in Medicine $\&$ Biology 41(11): 2231–49.
  * Gabriel, Lau and Gabriel (1996a) Gabriel, S., Lau, R. W. and Gabriel, C. (1996a).  The dielectric properties of biological tissues: Ii. measurements in the frequency range 10 hz to 20 ghz, Physics in Medicine $\&$ Biology 41(11): 2251–69.
  * Gabriel, Lau and Gabriel (1996b) Gabriel, S., Lau, R. W. and Gabriel, C. (1996b).  The dielectric properties of biological tissues: Iii. parametric models for the dielectric spectrum of tissues, Physics in Medicine $\&$ Biology 41(11): 2271–93.
  * Gaw et al. (2008) Gaw, R. L., Cornish, B. H. and Thomas, B. J. (2008).  The electrical impedance of pulsatile blood flowing through rigid tubes: A theoretical investigation, IEEE Transactions on Biomedical Engineering 55(2): 721–7.
  * Geuzaine and Remacle (2009) Geuzaine, C. and Remacle, J. F. (2009).  Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79(11): 1309–31.
  * Goren et al. (2018) Goren, N., Avery, J., Dowrick, T., Mackle, E., Witkowska-Wrobel, A., Werring, D. and Holder, D. (2018).  Multi-frequency electrical impedance tomography and neuroimaging data in stroke patients, Scientific Data 5(1).
  * Gorgolewski et al. (2011) Gorgolewski, K., Burns, C., Madison, C., Clark, D., Halchenko, Y., Waskom, M. and Ghosh, S. (2011).  Nipype: A flexible, lightweight and extensible neuroimaging data processing framework in python, Frontiers in Neuroinformatics 5: 13.
  * Grabner et al. (2006) Grabner, G., Janke, A. L., Budge, M. M., Smith, D., Pruessner, J. and Collins, D. L. (2006).  Symmetric atlasing and model based segmentation: An application to the hippocampus in older adults, in R. Larsen, M. Nielsen and J. Sporring (eds), Medical Image Computing and Computer-Assisted Intervention – MICCAI 2006, Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 58–66.
  * Grech et al. (2008) Grech, R., Cassar, T., Muscat, J., Camilleri, K. P., Fabri, S. G., Zervakis, M., Xanthopoulos, P., Sakkalis, V. and Vanrumste, B. (2008).  Review on solving the inverse problem in eeg source analysis, Journal of NeuroEngineering and Rehabilitation 5(1): 25.
  * Greitz et al. (1992) Greitz, D., Wirestam, R., Franck, A., Nordell, B., Thomsen, C. and Ståhlberg, F. (1992).  Pulsatile brain movement and associated hydrodynamics studied by magnetic resonance phase imaging, Neuroradiology 34(5): 370–80.
  * Hallez et al. (2007) Hallez, H., Vanrumste, B., Grech, R., Muscat, J., De Clercq, W., Vergult, A., D’Asseler, Y., Camilleri, K. P., Fabri, S. G., Van Huffel, S. and Lemahieu, I. (2007).  Review on solving the forward problem in eeg source analysis, Journal of NeuroEngineering and Rehabilitation 4(1): 46.
  * Hammond et al. (2017) Hammond, D., Price, N. and Turovets, S. (2017).  Construction and segmentation of pediatric head tissue atlases for electrical head modeling.  OHBM, Vancouver, Canada.
  * Hoetink et al. (2004) Hoetink, A. E., Faes, T. J. C., Visser, K. R. and Heethaar, R. M. (2004).  On the flow dependency of the electrical conductivity of blood, IEEE Transactions on Biomedical Engineering 51(7): 1251–61.
  * Holder (2005) Holder, D. S. (2005).  Electrical impedance tomography: Methods, History and Applications, 1 edn, IOP Publishing Ltd, Cornwall, UK.
  * Holder et al. (1996) Holder, D. S., Rao, A. and Hanquan, Y. (1996).  Imaging of physiologically evoked responses by electrical impedance tomography with cortical electrodes in the anaesthetized rabbit, Physiological Measurement 17(4A): A179–86.
  * Holmgren et al. (2019) Holmgren, M., Wåhlin, A., Dunås, T., Malm, J. and Eklund, A. (2019).  Assessment of cerebral blood flow pulsatility and cerebral arterial compliance with 4d flow MRI, Journal of Magnetic Resonance Imaging 51(5): 1516–25.
  * Horesh (2006) Horesh, L. (2006).  Some Novel Approaches in Modelling and Image Reconstruction for Multi-Frequency Electrical Impedance Tomography of the Human Brain, PhD thesis, Department of Medical Physics - University College London.
  * Hou et al. (2017) Hou, J., Turovets, S., Li, K., Luu, P., Tucker, D. and Larson-Prior, L. (2017).  Spatially resolved pediatric skull conductivities for inhomogeneous electrical forward modeling.  OHBM, Vancouver, Canada.
  * Jatoi and Kamel (2017) Jatoi, M. A. and Kamel, N. (2017).  Brain source localization using EEG signal analysis, CRC Press, Boca Raton, FL.
  * Kaipio and Somersalo (2005) Kaipio, J. and Somersalo, E. (2005).  Statistical and computational inverse problems, Springer, New York.
  * Kim et al. (2019) Kim, D.-E., Jang, J., Schellingerhout, D., Ryu, W.-S., Park, J.-H., Lee, S.-K., Kim, D. and Bae, H.-J. (2019).  Supratentorial cerebral arterial territories for computed tomograms: A mapping study in 1160 large artery infarcts, Scientific Reports 9(1).
  * Kneihsl et al. (2020) Kneihsl, M., Hofer, E., Enzinger, C., Niederkorn, K., Horner, S., Pinter, D., Fandler-Höfler, S., Eppinger, S., Haidegger, M., Schmidt, R. and Gattringer, T. (2020).  Intracranial pulsatility in relation to severity and progression of cerebral white matter hyperintensities, Stroke 51(11): 3302–9.
  * Kucewicz et al. (2008) Kucewicz, J. C., Dunmire, B., Giardino, N. D., Leotta, D. F., Paun, M., Dager, S. R. and Beach, K. W. (2008).  Tissue pulsatility imaging of cerebral vasoreactivity during hyperventilation, Ultrasound in Medicine & Biology 34(8): 1200–8.
  * Malone et al. (2014) Malone, E., Jehl, M., Arridge, S., Betcke, T. and Holder, D. (2014).  Stroke type differentiation using spectrally constrained multifrequency EIT: evaluation of feasibility in a realistic head model, Physiological Measurement 35(6): 1051–66.
  * Martins et al. (2019) Martins, T. C., Sato, A. K., Moura, F. S., de Camargo, E. D. L. B., Silva, O. L., Santos, T. B. R., Zhao, Z., Möeller, K., Amato, M. B. P., Mueller, J. L., Lima, R. G. and Tsuzuki, M. S. G. (2019).  A review of electrical impedance tomography in lung applications: Theory and algorithms for absolute images, Annual Reviews in Control 48: 442 – 71.
  * Meghdadi et al. (2019) Meghdadi, A. H., Popovic, D., Rupp, G., Smith, S., Berka, C. and Verma, A. (2019).  Transcranial impedance changes during sleep: A rheoencephalography study, IEEE Journal of Translational Engineering in Health and Medicine 7: 1–7.
  * Melis (2017) Melis, A. (2017).  Gaussian process emulators for 1D vascular models, PhD thesis, University of Sheffield.
  * Melis (2018) Melis, A. (2018).  openbf: Julia software for 1d blood flow modelling.
  * Melis et al. (2017) Melis, A., Clayton, R. H. and Marzo, A. (2017).  Bayesian sensitivity analysis of a 1d vascular model with gaussian process emulators, International Journal for Numerical Methods in Biomedical Engineering 33(12): e2882.
  * Michel and Brunet (2019) Michel, C. M. and Brunet, D. (2019).  EEG source imaging: A practical review of the analysis steps, Frontiers in Neurology 10: 325.
  * Min et al. (2019) Min, M., Kõiv, H., Priidel, E., Pesti, K. and Annus, P. (2019).  Noninvasive acquisition of the aortic blood pressure waveform, Wearable Devices - the Big Wave of Innovation, IntechOpen.
  * Newell et al. (2002) Newell, J. C., Blue, R. S., Isaacson, D., Saulnier, G. J. and Ross, A. S. (2002).  Phasic three-dimensional impedance imaging of cardiac activity, Physiological Measurement 23(1): 203–9.
  * Pesti et al. (2019) Pesti, K., Kõiv, H. and Min, M. (2019).  Simulation of the sensitivity distribution of four- electrode impedance sensing on radial artery, 2019 IEEE Sensors Applications Symposium (SAS), pp. 1–6.
  * Proença et al. (2020) Proença, M., Braun, F., Lemay, M., Solà, J., Adler, A., Riedel, T., Messerli, F. H., Thiran, J.-P., Rimoldi, S. F. and Rexhaj, E. (2020).  Non-invasive pulmonary artery pressure estimation by electrical impedance tomography in a controlled hypoxemia study in healthy subjects, Scientific Reports 10(1).
  * Raaijmakers et al. (1996) Raaijmakers, E., Marcus, J. T., Goovaerts, H. G., Vries, P. M. J. M., Faes, T. J. C. and Heethaar, R. M. (1996).  The influence of pulsatile flow on blood resistivity in impedance cardiography, Proceedings of 18th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Vol. 5, pp. 1957–8.
  * Sakka et al. (2011) Sakka, L., Coll, G. and Chazal, J. (2011).  Anatomy and physiology of cerebrospinal fluid, European Annals of Otorhinolaryngology, Head and Neck Diseases 128(6): 309–16.
  * Schmitter et al. (2013) Schmitter, S., Jagadeesan, B., Grande, A., Sein, J., Ugurbil, K. and Moortele, P. V. (2013).  4d flow measurements in the superior cerebellar artery at 7 tesla: feasibility and potential for applications in patients with trigeminal neuralgia, Journal of Cardiovascular Magnetic Resonance 15(S1).
  * Seki et al. (2006) Seki, J., Satomura, Y., Ooi, Y. and Seiyama, T. Y. A. (2006).  Velocity profiles in the rat cerebral microvessels measured by optical coherence tomography, Clinical Hemorheology and Microcirculation 34(1-2): 233–9.
  * Shen et al. (2018) Shen, H., Li, S., Wang, Y. and Qin, K.-R. (2018).  Effects of the arterial radius and the center-line velocity on the conductivity and electrical impedance of pulsatile flow in the human common carotid artery, Medical & Biological Engineering & Computing 57(2): 441–51.
  * Shen et al. (2016) Shen, H., Zhu, Y. and Qin, K.-R. (2016).  A theoretical computerized study for the electrical conductivity of arterial pulsatile blood flow by an elastic tube model, Medical Engineering & Physics 38(12): 1439–48.
  * Silva et al. (2017) Silva, O. L., Lima, R. G., Martins, T. C., Moura, F. S., Tavares, R. S. and Tsuzuki, M. S. G. (2017).  Influence of current injection pattern and electric potential measurement strategies in electrical impedance tomography, Control Engineering Practice 58: 276–86.
  * Song et al. (2018) Song, J., Chen, R., Yang, L., Zhang, G., Li, W., Zhao, Z., Xu, C., Dong, X. and Fu, F. (2018).  Electrical impedance changes at different phases of cerebral edema in rats with ischemic brain injury, BioMed Research International 2018: 1–10.
  * Song et al. (2013) Song, J., Morgan, K., Turovets, S., Li, K., Davey, C., Govyadinov, P., Luu, P., Smith, K., Prior, F., Larson-Prior, L. and Tucker, D. M. (2013).  Anatomically accurate head models and their derivatives for dense array eeg source localization, Funct Neurol Rehabil Ergon 3(2-3): 275–93.
  * Tidswell et al. (2001) Tidswell, T., Gibson, A., Bayford, R. H. and Holder, D. S. (2001).  Three-dimensional electrical impedance tomography of human brain activity, NeuroImage 13(2): 283–94.
  * Towers et al. (2000) Towers, C. M., McCann, H., Wang, M., Beatty, P. C., Pomfrett, C. J. D. and Beck, M. S. (2000).  3d simulation of EIT for monitoring impedance variations within the human head, Physiological Measurement 21(1): 119–24.
  * Vikner et al. (2019) Vikner, T., Nyberg, L., Holmgren, M., Malm, J., Eklund, A. and Wåhlin, A. (2019).  Characterizing pulsatility in distal cerebral arteries using 4d flow MRI, Journal of Cerebral Blood Flow & Metabolism 40(12): 2429–40.
  * Visser (1989) Visser, K. R. (1989).  Electric properties of flowing blood and impedance cardiography., Annals of Biomedical Engineering 17(5): 463–73.
  * Visser (1992) Visser, K. R. (1992).  Electric conductivity of stationary and flowing human blood at low frequencies, Medical & Biological Engineering & Computing 30(6): 636–40.
  * Vorwerk et al. (2014) Vorwerk, J., Cho, J.-H., Rampp, S., Hamer, H., Knösche, T. R. and Wolters, C. H. (2014).  A guideline for head volume conductor modeling in EEG and MEG, NeuroImage 100: 590–607.
  * Wagshul et al. (2011) Wagshul, M. E., Eide, P. K. and Madsen, J. R. (2011).  The pulsating brain: A review of experimental and clinical studies of intracranial pulsatility, Fluids and Barriers of the CNS 8(1).
  * Wembers et al. (2019) Wembers, C. C., Flürenbrock, F., Korn, L. C., Benninghaus, A., Misgeld, B. J. E., Shchukin, S. I., Radermacher, K. and Leonhardt, S. (2019).  FEM Simulation of Bioimpedance-Based Monitoring of Ventricular Dilation and Intracranial Pulsation, International journal of bioelectromagnetism : IJBEM 21(1): 7–20.
  * Zhang et al. (2020) Zhang, H. J., Shen, H., Li, Y. J. and Qin, K. R. (2020).  An in vitro circulatory device for studying blood flow electrical impedance in human common carotid arteries, 2020 IEEE 16th International Conference on Control Automation (ICCA), pp. 1518–22.