Pressure-Poisson Equation in Numerical Simulation of Cerebral Arterial Circulation and Its Effect on the Electrical Conductivity of the Brain

This study considers dynamic modelling of the cerebral arterial circulation and reconstructing an atlas for the electrical conductivity of the brain. While high-resolution 7-Tesla (T) Magnetic Resonance Imaging (MRI) data allow for reconstructing the cerebral arteries with a cross-sectional diameter larger than the voxel size, electrical conductivity cannot be directly inferred from MRI data. Brain models of electrophysiology typically associate each brain tissue compartment with a constant electrical conductivity, omitting any dynamic effects of cerebral blood circulation. Incorporating those effects poses the challenge of solving a system of incompressible Navier-Stokes equations in a realistic multi-compartment head model. We postulate that circulation in the distinguishable arteries can be estimated via the pressure-Poisson equation, which is coupled with Fick's law of diffusion for microcirculation. To establish a fluid exchange model between arteries and microarteries, a boundary condition derived from the Hagen-Poisseuille model is applied. The relationship between the estimated volumetric blood concentration and the electrical conductivity of the brain tissue is approximated through Archie's law for fluid flow in porous media. Through the formulation of the PPE and a set of boundary conditions based on the Hagen-Poisseuille model, we obtained an equivalent formulation of the incompressible Stokes equation. Thus, allowing effective blood pressure estimation in cerebral arteries segmented from open 7T MRI data. As a result of this research, we developed and built a useful modelling framework that accounts for the effects of dynamic blood flow on a novel MRI-based electrical conductivity atlas. The electrical conductivity perturbation obtained in numerical experiments has an appropriate overall match with previous studies on this subject.


Introduction
This study focuses on the dynamic modelling of the cerebral arterial circulation [1] and reconstructing an atlas for the electrical conductivity of the brain [2].Electrical conductivity is a governing parameter in several electrophysiological modalities targeting the brain, for example, electroencephalography (EEG) [3], transcranial electrical stimulation (tES) [4], and electrical impedance tomography (EIT) [5,6,7].Brain models of electrophysiology conventionally associate each head tissue compart-ment with a constant electrical conductivity [8,9,10,11], omitting the dynamic effects of cerebral circulation.While highresolution 7-Tesla (T) MRI data allows for the reconstruction of cerebral arteries with a cross-sectional diameter greater than the voxel size [12,13], electrical conductivity cannot be directly inferred based on MRI data.Advancements in dynamic conductivity modelling have addressed this critical gap in our understanding; in recent years, there has been growing interest in utilizing dynamic conductivity modelling to enhance our understanding of cerebral blood flow (CBF) regulation and its applications.Incorporating blood flow effects is an important and timely topic that has most prominently been approached via rheoencephalography [14], EIT [15,16], and diffusion-weighted MRI [17] measurements.Alternatively, a dynamic model can be constructed in silico through blood flow simulation, as recently suggested in [18,6,7].To this end, recent studies have, for instance, involved predicting the hydraulic conductivity of the capillary network, which is a critical factor in understanding microvascular dynamics [19].Modelling CBF poses the challenge of solving a system of Navier-Stokes equations (NSEs) numerically for an incompressible fluid flow in a complex structured blood vessel model, which has been demonstrated in several different contexts, e.g., in [20,21,22,23].Finding a numerical solution for incompressible NSEs includes, among other things: (1) the ability to incorporate incompressibility conditions into the system; (2) calculating pressure equivalently from flow velocity by assuming pressure boundary conditions (BCs) for the system; and (3) guaranteeing the numerical stability of the system [24].In addition, since a full set of arteries cannot be perfectly reconstructed from the MRI data, only approximate solutions can be obtained.Hence, a simplified model is well-motivated.In particular, the model of [6] has been built upon a one-dimensional (1D) approximation of NSEs in the main arteries [20].In this study, we propose using the pressure-Poisson equation (PPE) [25] in conjunction with Fick's law of diffusion for the microcirculation [26,27] to estimate blood circulation in an arbitrary set of arteries distinguishable in MRI data.In incompressible flows, the pressure is in equilibrium with a time-varying divergence-free velocity field, and its gradient is a significant physical quantity since it represents a force per unit of volume.We derive PPE for a three-dimensional (3D) steady state blood flow and domain.We do not directly discretize the NSEs [28]; instead, we (i) model the blood flow in the arteries as the simplest possible system that can be derived from Stokes equation (SE), (ii) derive BCs assuming that the pressure gradient can be obtained from a Hagen-Poiseuille flow in arterioles, which is in accordance with [29], and then (iii) discretize and solve the resulting PPE system.We utilize PPE to estimate pressure fields within the arterial domain by discretizing this system using a numerical method such as FE analysis while considering Neumann BCs.This process involves breaking down the arterial system into discrete elements, which enables us to solve the equation numerically.Subsequently, we utilize Fick's law to estimate the excess volumetric blood concentration within the microcirculation domain [26,27].We discretize Fick's law using the FE method, to numerically solve it within the domain.Our hypothesis posits that there is a gradual reduction in this concentration with increasing distance from the arteries.This presumption is consistent with physiological principles, as it aligns with the process of diffusion, wherein blood diffuses from areas characterized by higher concentrations to those with lower concentrations.The parameters determining the BCs are based on the knowledge of average blood viscosity, cerebral blood flow rate, pressure, and diameter of the arterioles constituting an interface between arteries and capillaries [1].The diffusion and decay parameter of the microcirculation model are, additionally, affected by experimental data of the microvessel density in different brain tissues [30].
The excess volumetric blood concentration plays a central role in determining the electrical conductivity σ within the tissue.We incorporate this critical information into the electrical conductivity atlas, which is a key component of the proposed model.This atlas quantifies how the electrical conductivity of the tissue varies in relation to the volumetric blood concentration.The relationship between microcirculation estimates and electrical conductivity is derived from Archie's law, a well-known and widely used two-phase mixture model.Archie's law allows us to determine the effective electrical conductivity of a mixture comprising fluid and porous media.Noteworthy studies [31,32,33,34] have contributed to our understanding of this relationship via electrical measurements.As outlined in this study, the impact of the volumetric blood concentration on electrical conductivity remains significant within a specific distance, typically ranging from 10 to 20 mm, from the arteries.To place our work in a broader scientific context, we explore the proposed model's connection with recent advancements in modelling dynamic blood flow along the capillary bed and their effect on the electrical conductivity distribution of the brain tissue.Our results obtained with a finite element (FE) discretization of a multi-compartment head segmentation suggest that, given a high-resolution 7T MRI dataset, PPE, together with Fick's and Archie's laws, allows us to approximate blood pressure effects on the electrical conductivity in the brain.We compare the results to a tissue-wise constant distribution [8] which provides the background model for Archie's law.Furthermore, we investigate potential future directions and applications of the proposed model, ultimately concluding this paper with an extensive discussion that emphasizes its importance and potential influence on future research.Our numerical implementation is openly available [35].

Methods
In this study, a domain for a human brain vasculature model is defined as the union Ω ∪ Ω of a microcirculation domain Ω and a domain Ω composed of distinguishable arteries.In Ω, the total pressure p in the arteries is assumed to be of the form where p (D) is a time-average of a dynamic arterial pressure and p (H) is a hydrostatic venous pressure distribution following from a (constant) gravitational force field f, blood density ρ and position x.Given a Riemannian metric g in Ω, p (H) (x) can be expressed as where C(x, x 0 ) is a geodesic, i.e., the shortest path on the surface, from a reference point x 0 to x and dr its differential.The Riemannian metric utilized in formulations preserves local differences in shape and size across various head regions while maintaining the domain's shape across coordinate transformations.We include the Riemannian metric in our mathematical model for generality but do not explicitly discuss the curvature of the domain in this study.We solve PPE to obtain an approximation for p, after which volumetric blood concentration in Ω is approximated using Fick's law of diffusion.Finally, electrical conductivity atlases are obtained based on the concentration via Archie's law.This section briefly reviews the theoretical grounds of PPE, Fick's law of diffusion, and Archie's law of two-phase electrical conductivity mixtures.

Circulation in arteries
In this study, blood is modelled as a non-homogeneous, incompressible viscous fluid moving through blood vessels as a Newtonian flow with constant absolute dynamic viscosity of µ = 0.004 Pa s, which can be considered a typical value in vessels with diameter from one to few millimeters with hematocrit between 45 and 60 % [36,37].The corresponding 3D time-dependent Stokes equation is: where Ω is the physical domain of the problem and [0, T ] is the time domain.The blood velocity and pressure are defined as u = u(x; t) and p = p(x; t), respectively, in which x ∈ Ω and t ∈ R + .The specification of the diffusion force, denoted as Lu, can be found in Appendix A. We identify ρ(x; t) = ρ as a constant blood (mass) density, and the term f = f(x; t) on the right-hand side accounts for the possible action of external forces.It is assumed that the initial velocity field u 0 is divergence-free.The vector ∇p is a function of given velocity data and depends on the constitutive properties of blood, ∇p = ρ f − ρu ,t + µLu , whose divergence leads to PPE.Derived from the momentum equation applying the incompressibility (3b), PPE of laminar flow is of the form Because any harmonic function with a vanishing mean can be added to the above equation, it is clear that this does not define a unique p.As a result, we must consider the specific BC for the system (4).Now, we face two critical questions: 1) Can equation (4) be used to calculate pressure p, and 2) does it imply incompressibility (3b)?It seems that the answer to both questions is yes, if the divergence of Ricci curvature term ∇•(Ri(u)) is zero or small enough, i.e., if the geometry is locally flat.Incompressibility for a static pressure field is implied by (4) under ∇ • (Ri(u)) = 0, as justified in Appendix B.

Pressure-Poisson equation
We solve the pressure in system (3) applying PPE and determining the proper BC based on that.Previously, FEs have been used to solve the pressure fields in the arteries driven by the socalled Neumann BCs, which are often sensitive and challenging to determine [38,39].It is theoretically sufficient to provide a BC for the system (3) by projecting (3a) onto the boundary in either a normal or tangential direction.Thus, when Lemma A1 and the incompressibility requirement (3b) are applied to the formula (3a) together with the assumption that f is a constant gravity field with ∇•(ρf) = 0. Utilizing formulas (1) and ( 2), the 3D-PPE model representing the dynamical aspect of the pressure field p (D) takes the following form: where parameter ζ is assumed to be constant, affecting the level of total blood flowing into the arterioles, and p (B) is a distribution enforced on the boundary, determining the contribution of the incoming flow.Note that the domain Ω is flat with zero or close-to-zero curvature, the boundary condition is set on the common boundary ∂Ω ∩ ∂ Ω of Ω and Ω, and ⃗ n is the normal unit vector that is defined on the artery wall.The pressure drop on the boundary is denoted by −g(∇p, ⃗ n) = −∂p/∂⃗ n which is the inward normal derivative of the blood pressure and characterizes the behavior of the fluid near the boundary.The BC follows from the assumption that the flow is laminar in arterioles, which leads the blood through the boundary wall to the microcirculation domain.The total level of this flow is scaled by the pressure p and λ = ξ/ξ, which is defined as the ratio between the length density ξ of microvessels [30] per unit volume (m −2 ), i.e., the total number of cross-sections per unit area, and the integral mean where |∂Ω| = ∂Ω dω ∂Ω .We discuss the Hagen-Poiseuille model [1] motivating the BC in Section 2.3 as a means of determining the BC.Assuming, for simplicity, a steady state flow, the blood velocity can be approximated based on the solution of PPE, implying the following formula: with u = 0 on ∂Ω (Appendix B).

Variational form
The equation determining blood flow pressure in the cerebral arteries is approximated numerically through FE discretization.We need to identify an appropriate variational formulation of PPE in order to continue with the FE discretization.Both the mathematical analysis and the numerical solution are based on weak formulations.Integration by parts results in a weak general form of the above equations for blood pressure (5) and velocity field (6) .We assume that p ∈ W with W = H 1 (Ω) and u ∈ V 0 and f ∈ V, where V 0 and V denote spaces of vector-valued functions in a physical 3D space, with V = [H 1 (Ω)] 3 and In other words, each of the three Cartesian components in V 0 and V is in the Sobolev space H 1 (Ω) of squareintegrable ( Ω |u| 2 dω Ω < ∞) functions with square integrable partial derivatives and A variational form of ( 5) can be obtained by multiplying the equation (5a) with a smooth enough test function q ∈ W and applying the divergence theorem.We arrive at the following variational problems: The continuous bilinear form b : The continuous bilinear form a :

Pressure boundary condition
The blood flows from Ω to the microcirculation domain Ω through the total cross-section area of the outlets of the arterioles on ∂Ω.Arterioles are microvessels that connect to arteries on one end and are continued by capillaries and thereon by venules in Ω.Most of the pressure decay in the blood flow, about 70 % of the total pressure [1], takes place in arterioles, which form a necessary transition zone for the blood pressure.In this study, we focus on modelling the blood flow within vessels and microvessels while excluding the circulation of interstitial fluid, which occurs outside the vessels.

Hagen-Poiseuille model
To obtain a value for ζ, we use the Hagen-Poiseuille equation of laminar flow [1], which, written for a single arteryarteriole interface (Figure 1), is of the form [1]: Here, it is taken into account that p stands for a pressure drop between the inlet and outlet of an artery.As a result, it is scaled by the constant ϑ, which determines the relative pressure drop in arterioles.In the equation ( 9), Q a represents the blood flow rate through a single arteriole, A a = π/4 D 2 a is its cross-sectional , i.e., a laminar balance between the inertial force and viscous drag in a cylindrical tube with a diameter D a and a length L. Right: A schematic of brain tissue cross-section with a set of blood vessels, including cross-sections of arterioles, capillaries, and venules.The length density of arterioles can be approximated by the total observed length density of microvessels [30], the cross-section areas of individual microvessels, and the total cross-section area fractions between the different microvessel types [40].It is assumed that, due to the randomness of the vessel orientations, there is no orientational dependence in the length densities.
area with a diameter of an arteriole D a , and L represents its length.The value of Q a can be approximated as where Q is the total flow through ∂Ω and ξ a is the average length density of arterioles on ∂Ω (Appendix C).It follows that where p is a given average pressure.Since the inward normal derivative of p represents force per surface area, we additionally scale the approximation following from Hagen-Poisseuille model ( 9) locally by the relative area covered by the arterioles, i.e., the product ξ a A a between the length density ξ a and per vessel cross-sectional area A a of arterioles (Figure 1), resulting in Here it has been taken into account that λ = ξ/ξ = ξ a /ξ a (Appendix C).

Microcirculation model
When blood flows from the arteries in Ω to the microcirculation domain Ω, it produces an excess blood volume concentration c = c(x; t) in comparison to the equilibrium state where the pressure in the head is constant.This can be counted as local blood supply upregulation, which varies the cross-sectional diameter [41].To approximate c, we apply advection-diffusion equation from Fick's second law of diffusion [26,27,29] and mass conservation.Fick's second law of diffusion and mass conservation have the following form: where u ⊗ c and J stand for the advective and diffusion flux, respectively.The flux J = J(x; t) is a vector pointing in the direction of movement, and the three-dimensional flux amplitude distribution |J| is proportional to the amount of blood flowing in the direction of J/|J| per unit time.The following equation is in accordance with Fick's first law under the assumption, that volumetric blood concentration and diffusive flow are proportionate to the tissue's relative microvessel density: Here ς is an effective diffusion coefficient and is referred to as diffusivity.In system (11), the term ĉ stands for a decay term.This results in the following governing differential system: where the Fickian and non-Fickian fluxes are represented by the terms u • ∇c and ςλ∆ B c, respectively.The concentration is maximized on the boundary ∂Ω∩∂ Ω due to the inward flux over ∂Ω ∩ ∂ Ω from Ω to Ω.The decay term εc accounts for the flow from the microcirculation domain to the venous circulation system, which happens at a rate proportional to the coefficient ε.
We consider a steady state solution, i.e., lim t→∞ c(x; t) = c(x), in which concentration reaches a constant state or remains stable over time, c(x) = c, and the flow's macrovelocity u vanishes.The following simplified form results from applying the conservation of mass condition (13b) to the limiting equation (13a) where s = ςλ g(∇c, ⃗ n) .

Parameter estimation
To obtain an approximation for ς without taking into account the geometry, we rely on the Hagen-Poiseuille model for the arterioles (section 2.3.1) with the additional assumptions that (i) the excess concentration is zero when the pressure is zero and that (ii) the concentration gradient is constant along the length of the microvessels.Under these assumptions, Fick's first law (12) for the flux passing ∂Ω ∩ ∂ Ω is of the form where A a ξ a is the mean volume concentration of the blood and L is the distance over which the concentration decreases from ϑA a ξ a .Sustituting |J| = λQ a ξ a and the formula (10) for L, it follows that ς = A a p 8πµ .
Parameter ε can be obtained assuming that at each interior point in Ω the blood in venules flows to a venous vessel, thereby exiting Ω.In a balanced state, the concentration loss caused by this flow equals the density |J|.We establish a balance by requiring that the sink amplitude ε integrated over a radius R sphere with a volume V max = (4/3)πR 3 of the largest element in the FE discretization matches an outward flux |J| integrated over the surface of the sphere, i.e.,

Variational form
The variational form of ( 14) can be obtained in the same fashion as in the case of PPE.Multiplying (14) with a smooth enough test function h ∈ W and applying the divergence theorem, we arrive at the following form: where the linear boundary term describing the incoming flow is on the left-hand side and ⃗ n is the normal unit vector defined in the microcirculation domain.The continuous bilinear form

Discretization
We use the Ritz-Galerkin method [42] to discretize the PPE (5) and Fick's law of diffusion (14), whose solutions are assumed to be contained by the trial function spaces respectively.In each case, the discretization error is assumed to be orthogonal to the solution.Linear Lagrangian (nodal) basis functions are utilized in this context.Specifically, we have the sets {ψ i } n i=1 and {φ h } m h=1 supported in Ω, as well as {ϕ h } m h=1 supported in Ω.These sets consist of piecewise linear functions that satisfy the conditions at the nodes of Ω.Consequently, the velocity u ∈ V, pressure p (D) ∈ W, and concentration c ∈ V take the following forms: for ℓ = 1, 2, 3.The coordinate vectors are denoted by

Blood pressure and velocity in arteries
The system of ( 7) is equivalent to the following Ritz-Galerkin discretized form: II. Find u h ∈ V 0,h such that, for all Here V 0,h ⊂ V h is obtained from V h by excluding the boundary degrees of freedom, i.e., basis functions with non-zero values on the boundary ∂Ω.Equation ( 16) possesses a solution that satisfies the following equation: Here p denotes a contribution of the incoming flow, normalized to a given pulse pressure p and p is a normotensive diastolic average, p i = p for each entry i.The matrices corresponding to equation ( 18) can be expressed as follows: where K = K i j m×m and M = M i j m×m .Following from the present modelling premises, we find an estimate for the systolic pressure distribution as a steady state solution, whose boundary restriction approximately satisfies p (D) | ∂Ω = p (B) , i.e., the boundary restriction of p (D) | ∂Ω equals that of p (B) .Given an initial approximation p(0) of p, a recursive sequence of p (D) (1) , p (D) (2) , . . . is obtained by finding p (D) (k) as a solution of (18) corresponding to p(k) , and setting p(k+1) = p (D) | ∂Ω , for k = 0, 1, 2, . ... As a stopping criterion, we use the tolerance condition with ϵ = 0.01 and fix p (D) (K) as the final estimate of p (D) .The initial distribution p(0) is selected to be piecewise constant with p(0) i = p, if i corresponds to one of two inlets (Figure 2) placed in the vicinity of the anterior and posterior end-points of Circle of Willis at the base of the brain, where the blood flow enters the brain, one in basilar artery and the other one in the junction of anterior cerebral and anterior communicating arteries.Other entries of p(0) are set to zero.Two sources were applied to ensure balanced results; namely, asymmetry of the flow in the Circle of Willis has been shown to extend to global scale [23].
The solution of equation ( 17) satisfies u ℓ = 1/µ L −1 q ℓ , where the components of matrices L and q ℓ are obtained as follows: where ℓ = 1, 2, 3. Matrix L can be obtained from K by excluding the boundary degrees of freedom from row and column indices.Namely, the test function space V h is linear and nodal akin to W h , but in V h the boundary degrees of freedom are set to zero due to the zero boundary condition for the velocity.

Volumetric Blood concentration in microcirculation
The discretized diffusion problem related to the system of ( 15) can be formulated as follows: A numerical solution c of ( 19) can be obtained via where and w = w 1 w 2 . . .w m T .

Electrical conductivity of brain tissues
To approximate how excess volumetric blood concentration perturbs the electrical conductivity distribution of the head, we apply Archie's law, a two-term linear combination of power functions that approximates the effective electrical conductivity σ for a two-phase mixture of fluid and inhomogeneous medium [31,32,33,34].For a two-phase mixture, Archie's law is of the form where σ f and σ m denote conductivities of fluid and medium, respectively, and β is so-called cementation factor [31,33], which for the cerebral cortex is between 3/2 and 5/3 [31].The lower and upper limits for β follow from spherical and cylindrical inhomogeneities, which in the cortex are represented by the somas and dendrites of the pyramidal cells, respectively.When substituted in the formula of Archie's law, β = 3/2 and β = 5/3 yield a lower and upper bound for the effective electrical conductivity, respectively.
Alternatively, the effective electrical conductivity can be estimated from above and below via Hashin-Shtrikman upper and lower bound, defined as respectively.Hashin-Shtrikman bounds have shown to be valid for coated spherical inhomogeneities of all different sizes, filling the space [43].

Numerical experiments
Table 1: Compartments of the head model segmentation were obtained using the FreeSurfer software suite [44] together with FieldTrip's [45] segmentation interface, which has been built upon the functions of the SPM12 package [46].The vessel segmentation was performed using the Vesselness algorithm [47,48].The piecewise constant background approximation of the electrical conductivity distribution σ m was based on [8].The subcortical active nuclei were associated with the conductivity of the grey matter [49].Vessel conductivity was chosen to match the blood conductivity [50].Skull and skin conductivity values were included in the head segmentation but not in the electrical conductivity atlases, since the segmented blood vessels were fully enclosed by the skull, inside the cranial cavity.

Segmentation
We performed numerical experiments using a realistic multi-compartment head model to assess PPE in combination with Fick's and Archie's laws in reconstructing an electrical conductivity atlas of the brain.This head model was created using the open sub-millimeter precision Magneto Resonance Imaging (MRI) dataset 1 of CEREBRUM-7T [13].The dataset has been acquired using 7T magnetic flux density and, therefore, allows distinguishing the arterial vessels as a separate compartment, as shown in [12].The FreeSurfer Software Suite 1 doi:10.18112/openneuro.ds003642.v1.1.0Table 2: The physical parameters applied in numerical simulations.Gravitational acceleration has been set to its average level and oriented parallel to the z-axis.The electrical conductivity of the blood σ f and the reference pressure were chosen according to [50] and [21], respectively.Blood density ρ, µ, total CBF Q, pressure decay in arterioles ϑ, diameters D a , D c and D v of arterioles, capillaries and venules (subtracting the total wall thickness, 2.0E-5, 2.0E-06 and 2E-6, respectively), and their relative total area fractions γ a , γ c and γ v , respectively, are based on the textbooks [40,1].Microvessel density ξ in cerebral and cerebellar grey and white matter (GM and WM), subcortical WM, and the brain stem was chosen according to the median values observed in [30].The cementation factor estimates for spherical and cylindrical inhomogeneities approximating somas and dendrites of brain tissues, respectively, are based on [31].Arteriole length has been obtained by substituting other parameter values in (10) with an appropriate correspondence to the values found in literature [ [44], FieldTrip's [45] interface for the SPM12 surface extractor [46], and the Vesselness algorithm [47,48,12] were applied to segment the arteries, i.e., domain Ω.The other 16 brain compartments (Table 1), also enclosed by the skull, constituted Ω. Skin and skull were not included in the electrical conductivity atlas since the domain of arterial vessels Ω was fully contained by the skull.The microcirculation domain Ω included microvessel-containing compartments: in addition to skin and skull, the cerebrospinal fluid (CSF) compartment and CSF-filled ventricles were excluded from Ω.

Vessel extraction
The vessel extraction process was inspired by the work of Choi et al. [51] but is not entirely based on their proposed procedure.The Frangi filter was first applied to the MRI data sliceby-slice, and then the results were aggregated to produce the final arterial model.This process was performed in the following three steps: 1. Frangi's algorithm was applied to both the preprocessed INV2 and T1w slices of the dataset separately.We used the Scikit-Image [47] package's implementation of the Frangi method with different parameters for each slice.2. After applying the filter to a specific slice of the INV2 and T1w data, we created a mask by superposing these two layers in an element-wise manner.The mask was binarized using a user-defined threshold level; every element with value less than this threshold was set to zero and, otherwise, to one.3. The segmented cerebral vessels were obtained by iterating the previous steps through an axis of the MRI image and aggregating the results.In order to reduce noise, aggregation was performed separately for sagittal, axial, and coronal slices using the following scoring scheme: if a voxel was detected as a vessel in two or three of the results, it was considered a vessel in the final vessel mask; otherwise, it was neglected.

Numerical simulations
The numerical simulations were performed using the open Zeffiro Interface [35,52] (ZI) toolbox.Solvers for PPE, Fick's law, and Archie's law were implemented as Matlab codes and included in ZI 2 .Using ZI, the volume of the head segmentation was discretized by a tetrahedral FE mesh of 6.4 M nodes and 32 M elements, corresponding approximately to 1 mm overall resolution.Of these, Ω contained 0.15 M nodes and 0.54 M tetrahedra, and Ω 2.4 M nodes and 11 M tetrahedra.Other relevant parameter values can be found in Table 2.After solving the discretized pressure in Ω from ( 18), the excess blood concentration c in Ω was obtained by solving (20).Archie's law (21) was evaluated using the excess concentration c and two alternative cementation factors β = 5/3 and β = 3/2 corresponding to cylindrical and spherical tissue inhomogeneities and constituting a lower and upper bound approximation for the electrical conductivity, respectively.In addition, the Hashin-Shtrikman lower and upper bounds (23) and (22) were evaluated as an alternative approximation.As a result, altogether five different electrical conductivity atlases were obtained: one corresponding to the piecewise constant background distribution and four effective electrical conductivity atlases following from the different mixture models.To examine differences between the background σ bg and effective σ eff distributions, the following difference measures (%) were evaluated: Of these, RDM (relative difference measure) evaluates the overall relative difference between normalized distributions σ eff and σ bg , MAG (magnitude measure) shows the average amplitude of σ eff compared to σ bg , and PRD (pointwise relative difference) is the difference between σ eff and σ bg in relation to ∥σ bg ∥ ∞ for each point.

Results
This section describes the results of our four-phase numerical simulations, where (i) the brain model was segmented to obtain Ω and Ω as well as an estimate for the background tissue concentration; (ii) the blood pressure p and velocity u following from the present PPE model were found; (iii) Fick's law was applied to find an estimate c for excess blood concentration in the microcirculation domain Ω; finally, (iv) effective electrical conductivity atlases were reconstructed by estimating the effect of the excess blood concentration on the background distribution via Archie's law and Hashin-Shtrikman bounds.The results of the numerical experiments have been included in Figure 2 showing the head segmentation results; Figure 3 with sagittal, axial, and coronal illustrations of the pressure and concentration distributions obtained as numerical solutions of PPE and Fick's law; Table 3 including RDM and MAG for the estimates of background and effective electrical conductivity atlases obtained via Archie's model; histograms showing the value distributions of the blood pressure, velocity, and volumetric blood concentration (Figure 4) and the PRD of the electrical conductivity distribution (Figure 5); and Figure 6 visualizing effective conductivity atlases for sagittal, axial, and coronal projections.

Phase (i): segmentation
The segmentation obtained (Figure 2) shows that the Vesselness algorithm found a connected set of arteries inside the skull.Finding a connected set was considered important for the continuity of the PPE solution.Therefore, it was prioritized in the segmentation process over resolution, which led to an overlap of the tightly packed vessel bundles following from cortical branches.This geometrical distortion was considered unavoidable, due to the relatively small diameter of the cerebral arteries compared to the 1 mm resolution of the discretization.

Phase (ii): PPE
As shown by Figures 3 and 4, excluding the outliers, the total pressure distribution p varies between 80 mmHg and 115 mmHg in the artery domain Ω.The values are the greatest at the base of the brain, close to the basilar artery, the deepest vessel in Ω, which is oriented nearly vertically in front of the brainstem.The pressure gradually decreases when moving towards the cerebral cortex, where branched, overlapping structures are dominant.This result is expected since the total vessel area gradually increases as the branching occurs.Thus, the pressure gradually decreases as the blood flows towards smaller vessels, eventually entering the microcirculation domain Ω.The blood velocity profile varies between 0 and 1.46 m/s, excluding the outliers (Figures 3 and 4).The greatest values were Surface meshes (Table 1) were extracted using the segmentation routines of the FreeSurfer software suite [44] and FieldTrip's [45] interface for SPM12's surface extractor [46].The cerebral arterial vessels shown on the right (violet) were segmented using Frangi's Vesselness algorithm [47,48] as suggested in [12].The clipping planes correspond to the sagittal, axial, and coronal surface slices shown in this study.Bottom: Sagittal, axial, and coronal projections of the segmentation without cerebral grey and white matter surfaces, showing the arteries and how they integrate with the subcortical structures.Two spherical surfaces show the locations of two 10 mm diameter inlets placed in the vicinity of the anterior and posterior end-points of the Circle of Willis at the base of the brain.The locations correspond to the basilar artery and the junction of anterior cerebral and anterior communicating arteries.observed in anterior and middle cerebral artery which had overall greater velocities than posterior cerebral artery or the smaller arteries in cortical branches.

Phase (iii): excess blood concentration
The excess blood concentration estimate c obtained via Fick's law (Figure 3) expectedly decays when moving away from the arteries in Ω that bring blood into the microcirculation domain Ω.The amplitude of c vanishes at a distance of 10-20 mm from the arteries.Consequently, the effect of the concentration on the electrical conductivity atlases is, within the present model, limited to this distance.

Phase (iv): effective electrical conductivity atlases
As shown by Table 3, the reconstructed atlases differ overall by 1.65-1.86% and 11.19-11.85% with respect to the MAG and RDM, respectively.The major part of the differences is limited to a few percent of the volume fraction (Figure 4), which is obvious based on the excess concentration estimates obtained in Figure 3: Approximations of systolic blood pressure, velocity, and volumetric blood concentration illustrated on a logarithmic color scale Top: Sagittal, axial, and coronal views of the pressure distribution (mmHg) solving PPE in domain Ω.The large subcortical arteries can be observed to have a higher pressure up to 115 mmHg, compared to the smaller arteries of the cortical branches, for which the value range extends down to 80 mmHg.Values larger than q 75 +1.5(q75 −q 25 ), where q 25 and q 75 denote the 25th and 75th percentiles, respectively, have been excluded as outliers.Center: The velocity distribution (m/s) in Ω.The greatest values extend up to or close to 1.46 m/s.The anterior and middle cerebral arteries can be observed to have overall greater velocities than the posterior cerebral arteries or the arteries with smaller cortical branches.Outliers larger than q 75 + 1.5(q 75 − q 25 ) have been excluded.Bottom: Sagittal, axial, and coronal surface cuts of the estimated excess blood concentration (%) in the microvessel domain Ω as predicted by Fick's law.The greatest values are obtained in the vicinity of the arterial vessel boundary ∂Ω, i.e., the boundaries of the violet subdomains.The concentration decays to zero within 10-20 mm of a ∂Ω.The visible structures not included in Ω include the arterial vessels, i.e.Ω, ventricles (blue), and cerebrospinal fluid (green).
the third phase and is verified by PRD, showing that locally the largest differences are approximately 30 % with respect to the Figure 4: Histograms for the distributions of blood pressure, velocity, and volumetric concentration.The horizontal axis shows the value of the pressure (top), velocity (center), or concentration (bottom), and the vertical one shows the corresponding volume fraction of Ω (%) on a logarithmic scale.Each visualization includes 20 bars.For pressure and velocity, outliers larger than q 75 + 1.5(q 75 − q 25 ), where q 25 and q 75 have been excluded, the concentration is shown for the values above 0 and below 100 %. maximum (1.79 S/m) of the background distribution.Those can be related to regions close to the vessels where the excess blood concentration is close to one.Compared to the other atlases, the Hashin-Shtrikman upper bound yields a greater volume fraction of PRD values slightly below 30 %. Spatial differences between atlases corresponding to different mixture models are minor, which can be observed based on Figure 6.

Discussion
This study demonstrated that a simplification of the Stokes equation (SE), namely the pressure-Poisson equation (PPE) [25] allows for the estimation of the blood pressure in cerebral arteries segmented from open 7T MRI data [13].We introduced a boundary condition (BC) based on the Hagen-Poisseuille model [1] to bind PPE with the governing physical parameters of CBF, particularly the microvessel diameters [1] and densities [30].Through the formulation of the PPE and the BC, we obtained an equivalent formulation of the incompressible SE.Based on the solution of PPE, we estimated the excess volumetric blood concentration in microvessels caused by the pressure using Fick's law [26,27], the parameters of which were likewise obtained via the Hagen-Poisseuille model.Finally, the effect of the excess concentration on the brain tissues was approximated using Archie's law as well as the upper and lower bounds of Hashin and Shtrikman [31].Our four-phase modelling process (i) first generates a multi-compartment FE mesh and a piecewise conductivity atlas of the head, then (ii) finds a solution for PPE and (iii) Fick's law, and finally, (iv) reconstructs an atlas.As the strength of our approach, we suggest the direct applicability of NSEs and their approximations to individual datasets to potentially improve the quality of electrophysiological brain modelling.Thereby, the results of this study complement the recently developed statistical approaches following from 1D NSEs [6,7].Overall, this study advances the electrical conductivity approximation techniques applicable in electrophysiological modalities, where dynamic components affecting the conductivity atlases are typically absent, e.g., EEG/MEG source localization [8], tES [4], and EIT [5,53,6].In particular, we have shown how to incorporate the dynamic blood flow effects in modelling the electrical conductivity when a high-resolution and high-intensity MRI segmentation  1.Compared to the background, each effective electrical conductivity distribution the conductivity in the vicinity of the arteries is emphasized.Center row:: Sagittal, axial, and coronal illustration of the effective electrical conductivity atlas obtained with Archies's law assuming that the tissue inhomogeneities are cylindrical.Bottom row: A comparison between the different effective electrical conductivity atlases in a sagittal region of interest.The differences can be observed to be minor.Hashin-Shtrikman upper bound shows more spread effects of the blood flow than the other estimates.
with distinguishable blood vessels is available [12].We consider PPE an appropriate approximation of SE under the present modelling framework, since the MRI data does not allow a perfect segmentation of cerebral arteries.Hence, a more advanced solution based on NSEs might at least partly suffer from the limited accuracy of the segmentation.The current results suggest a spatial pressure variation between 80 and 115 mmHg, which matches with ±5 mmHg discrepancy the normotensive systolic pressures found via numerical simulation in [21] for arteries with diameter greater than 0.5 mm (e.g., 117 mmHg for 4.839 mm internal carotid artery, 113 mmHg for 3.448 mm basilar artery, 110 mmHg for 0.545 mm distal medial striate artery, and 85 mmHg for 1.039 mm posterior parietal branch of the middle cerebral artery).The velocity distribution can be considered appropriate based on experimental transcranial doppler ultrasound studies, e.g., [54,55].A peak systolic velocity of 1.4 m/s has been suggested as a threshold criterion for mild stenosis in an intracranial vessel [54].While our results reach that threshold, the vast majority of the velocity distribution stays below 1.4 m/s, i.e., in the range found for healthy subjects.While there are some obvious artifacts, the structure of the observed velocity distribution reflects the existing literature, as the greatest values were observed in the larger vessels, of which the middle and anterior cerebral arteries have overall greater velocities than the posterior cerebral artery or the vessels with a smaller diameter in the cortical branches [56,1].The present model of excess blood flow and concentration near the arteries, resulting from arterial pressure, builds upon the utilization of Fick's law within the intricate network of microvessels [26,27].This model incorporates the assumption of a linear pressure drop along the length of microvessels and considers the adaptability of microvessel diameter and length to regulate blood flow [41].Since venous vessels or any vessels outside the brain are difficult to distinguish based on MRI data [12], we introduced a uniform sink term that covers the entire brain.Although directly validating the accuracy and reliability of numerical simulation results may be impractical, an observable correlation can be found between the estimated concentration distribution and whole-brain CBF scans obtained through MRI, positron emission tomography, and single-photon emission computed tomography [57,58,59].A potential feature affecting the accuracy of our estimate obtained for the volumetric concentration is the fluid exchange between the microvessels and the tissue interstition, which was omitted in deriving the diffusion coefficient.We, however, deem that the uncertainty related to the sink term is likely to be a dominant error source as the venous flow rate is greater than the interstitial fluid exchange.Our excess concentration estimates and the general knowledge of perfusion scans both demonstrate that the excess blood in the brain is not limited to the large arteries but is somewhat spread in the neighborhood of those.Thus, the present distributions obtained using Archie's law and Hashin-Shtrikman bounds, which, in this study, were designed to take this aspect into account, might represent an improvement compared to the piecewise constant electrical conductivity estimates when a segmentation of the arteries is available.Furthermore, as suggested in [12], an additional challenge arises due to a significant portion of the blood being limited to the subcortical region of the brain.This area is acknowledged to exhibit compromised localization and stimulation accuracy when utilizing non-invasive techniques such as EEG/MEG, EIT, and tES.We consider studying this aspect from the application point of view as important future work, while the focus of this study clearly was on establishing an appropriate modelling framework in order to take dynamic blood flow effects into account in an individualized MRI-based electrical conductivity atlas of the brain.The knowledge of the mixture models suggests [31] that while Hashin-Shtrikman bounds give valuable information about the potential modelling discrepancies, they do not achieve the accuracy of Archie's law, which is better suited for brain tissues.While the discrepancies between the models can be considered significant regarding the local value of electrical conductivity, the global differences observed in this study can be considered minor.As for the limitations of this study, based on the above reasoning, we do not expect that our model in its current form would be applicable for obtaining other than coarse estimates of blood pressure, velocity, volumetric concentration, and electrical conductivity distribution; the current results are limited to showing the feasibility of evaluating the PPE approximation for pressure and velocity to obtain estimates for concentration and electrical conductivity directly based on individualized data in electrophysiological modelling.As a governing limitation, we consider the weak distingushability of blood vessels from MRI data recorded with a magnetic flux density lower than 7T, as most datasets applied in electrophysiological head model generation comprise 3T or 4T measurements.Moreover, the present mathematical model is simplified and thus limited in its capability to approximate cerebral circulation.Features omitted in this model include any time-dependencies of the blood flow, the contribution of viscoelastic arterial walls [60] and fluid exchange between microcirculation and tissue interstitium [61].A simplified model is, however, well-motivated in this study due to the incompleteness of the MRI data with respect to obtaining blood vessel segmentation.

Future prospects
While we achieved an appropriate overall match with the existing results, such as a similar range of values and distribution of the electrical conductivity perturbation as shown in [6,7], further studies are required to validate the present approach.The focus of future work will be on approximating the velocity field of cerebral circulation, which would necessitate solving a time-dependent system, e.g., the Navier-Stokes equation.This work will involve a complex and multidisciplinary study focusing on the electrical conductivity of the brain, incorporating mathematical modelling, numerical simulations, and data analysis to advance our understanding of brain function.We will use a segmentation of arterial blood vessels in the brain and a dynamic solution of NSEs coupled with Fick's law to represent microcirculation in a specific domain.For this purpose, we will solve a discretized non-Newtonian NSEs system using a two-stage process involving pressure estimation and velocity field updates.This includes regularization techniques to ensure numerical stability.The motivation behind this research lies in its significance for EEG, tES, and EIT, particularly, in scenarios characterized by dynamic modelling.Looking ahead, the future path of this research will necessitate deep exploration of theoretical foundations, such as geometry, boundary conditions and viscosity models, as well as experimental validation.An important aspect is, for example, how computing geometry influences the results obtained.To enlighten this, for example, a modelling study can be conducted to validate the performance of the current versions of PPE and Fick's law within a simplified computational geometry, such as a cylindrical domain.Furthermore, an experimental study can be conducted to compare the results of perfusion imaging with numerically simulated volumetric blood concentration.We also aim at multi-subject studies and evaluations conducted at different scales.

Figure 1 :
Figure 1: Our model of the artery-arteriole interface shows the arteries in the domain Ω and the microcirculation domain Ω, where blood flows from the arteries into a network of arterioles, capillaries, and venules.Left:The relationship between the flow rate and the normal derivative g(∇p, ⃗ n) on ∂Ω of the pressure (with respect to the Riemannian metric) is determined by the Hagen-Poiseuille equation[1], i.e., a laminar balance between the inertial force and viscous drag in a cylindrical tube with a diameter D a and a length L. Right: A schematic of brain tissue cross-section with a set of blood vessels, including cross-sections of arterioles, capillaries, and venules.The length density of arterioles can be approximated by the total observed length density of microvessels[30], the cross-section areas of individual microvessels, and the total cross-section area fractions between the different microvessel types[40].It is assumed that, due to the randomness of the vessel orientations, there is no orientational dependence in the length densities.

Figure 2 :
Figure 2: Top: A surface segmentation was obtained for the subject sub-045 of the open CEREBRUM-7T dataset 1 [13] containing 7-Tesla (T) MRI data.Surface meshes (Table1) were extracted using the segmentation routines of the FreeSurfer software suite[44] and FieldTrip's[45] interface for SPM12's surface extractor[46].The cerebral arterial vessels shown on the right (violet) were segmented using Frangi's Vesselness algorithm[47,48] as suggested in[12].The clipping planes correspond to the sagittal, axial, and coronal surface slices shown in this study.Bottom: Sagittal, axial, and coronal projections of the segmentation without cerebral grey and white matter surfaces, showing the arteries and how they integrate with the subcortical structures.Two spherical surfaces show the locations of two 10 mm diameter inlets placed in the vicinity of the anterior and posterior end-points of the Circle of Willis at the base of the brain.The locations correspond to the basilar artery and the junction of anterior cerebral and anterior communicating arteries.

Figure 5 :
Figure 5: Histograms showing the distribution of the relative difference PRD (%) between the piecewise constant background and the approximated effective electrical conductivity in those parts of Ω, where PRD exceeds 0.1 %.The horizontal axis shows the value of PRD, and the vertical one shows the volume fraction (%) of the corresponding computing domain on a logarithmic scale.Each visualization includes 20 bars.

Figure 6 :
Figure 6: Top row: Sagittal, axial, and coronal visualization of the piecewise constant background distribution σ m in which each compartment corresponds to a constant electrical conductivity given in Table1.Compared to the background, each effective electrical conductivity distribution the conductivity in the vicinity of the arteries is emphasized.Center row:: Sagittal, axial, and coronal illustration of the effective electrical conductivity atlas obtained with Archies's law assuming that the tissue inhomogeneities are cylindrical.Bottom row: A comparison between the different effective electrical conductivity atlases in a sagittal region of interest.The differences can be observed to be minor.Hashin-Shtrikman upper bound shows more spread effects of the blood flow than the other estimates.

Table 3 :
Relative difference measure (RDM) and magnitude (MAG) difference for the mixture models applied in this study.