On the Sensitivity Analysis of Porous Finite Element Models for Cerebral Perfusion Estimation

Computational physiological models are promising tools to enhance the design of clinical trials and to assist in decision making. Organ-scale haemodynamic models are gaining popularity to evaluate perfusion in a virtual environment both in healthy and diseased patients. Recently, the principles of verification, validation, and uncertainty quantification of such physiological models have been laid down to ensure safe applications of engineering software in the medical device industry. The present study sets out to establish guidelines for the usage of a three-dimensional steady state porous cerebral perfusion model of the human brain following principles detailed in the verification and validation (V&V 40) standard of the American Society of Mechanical Engineers. The model relies on the finite element method and has been developed specifically to estimate how brain perfusion is altered in ischaemic stroke patients before, during, and after treatments. Simulations are compared with exact analytical solutions and a thorough sensitivity analysis is presented covering every numerical and physiological model parameter. The results suggest that such porous models can approximate blood pressure and perfusion distributions reliably even on a coarse grid with first order elements. On the other hand, higher order elements are essential to mitigate errors in volumetric blood flow rate estimation through cortical surface regions. Matching the volumetric flow rate corresponding to major cerebral arteries is identified as a validation milestone. It is found that inlet velocity boundary conditions are hard to obtain and that constant pressure inlet boundary conditions are feasible alternatives. A one-dimensional model is presented which can serve as a computationally inexpensive replacement of the three-dimensional brain model to ease parameter optimisation, sensitivity analyses and uncertainty quantification. The findings of the present study can be generalised to organ-scale porous perfusion models. The results increase the applicability of computational tools regarding treatment development for stroke and other cerebrovascular conditions.

Abstract-Computational physiological models are promising tools to enhance the design of clinical trials and to assist in decision making. Organ-scale haemodynamic models are gaining popularity to evaluate perfusion in a virtual environment both in healthy and diseased patients. Recently, the principles of verification, validation, and uncertainty quantification of such physiological models have been laid down to ensure safe applications of engineering software in the medical device industry. The present study sets out to establish guidelines for the usage of a three-dimensional steady state porous cerebral perfusion model of the human brain following principles detailed in the verification and validation (V&V 40) standard of the American Society of Mechanical Engineers. The model relies on the finite element method and has been developed specifically to estimate how brain perfusion is altered in ischaemic stroke patients before, during, and after treatments. Simulations are compared with exact analytical solutions and a thorough sensitivity analysis is presented covering every numerical and physiological model parameter. The results suggest that such porous models can approximate blood pressure and perfusion distributions reliably even on a coarse grid with first order elements. On the other hand, higher order elements are essential to mitigate errors in volumetric blood flow rate estimation through cortical surface regions. Matching the volumetric flow rate corresponding to major cerebral arteries is identified as a validation milestone. It is found that inlet velocity boundary conditions are hard to obtain and that constant pressure inlet boundary conditions are feasible alternatives. A one-dimensional model is presented which can serve as a computationally inexpensive replacement of the three-dimensional brain model to ease parameter optimisation, sensitivity analyses and uncertainty quantification.
The findings of the present study can be generalised to organscale porous perfusion models. The results increase the applicability of computational tools regarding treatment development for stroke and other cerebrovascular conditions.

INTRODUCTION
Cerebrovascular diseases including stroke impose a heavy burden on society. 59 The majority of stroke cases are caused by a thrombus blocking a major cerebral artery leading to severe blood shortage in the brain (ischaemic stroke). In recent decades, ischaemic stroke treatment has been revolutionised by thrombolysis (dissolving of blood clot with a thrombolytic agent) 43 and thrombectomy (mechanical thrombus removal with a stent retriever). 4,18 Advancing the related clinical procedures has great potential to benefit stroke patients but is tied entirely to resource intensive animal experiments and clinical trials. 4,14,18,22,23 To ease this task, the in silico clinical trials for the treatment of acute Ischaemic STroke (INSIST) consortium aims to develop computational models of stroke and its treatments. 40 Once reliable models are available, in silico clinical trials could sharpen the focus of clinical trials and hence save time and money. 32,55,76 Estimating brain perfusion in both healthy and occluded scenarios is an important element of the IN-SIST pipeline and serves as a bridge between localised treatment effects and patient outcome. 40 To this end, a multi-scale haemodynamic model is under development which combines a one-dimensional (1D) network model of large arteries 5,54 and a three-dimensional (3D) porous perfusion model 39 to simulate blood flow in the entire vasculature. The models will be coupled strongly through the cortical surface. The envisaged interface will facilitate communication between the boundary conditions of the models, namely pressure and volumetric flow rate values averaged over cortical territories. 54 Previously somewhat similar weak-coupling schemes have been established between a lumped parameter model and a porous continuum model to cover arteriole boundary conditions at the brain surface. [24][25][26]71,72 This approach preserves the two models as independent functional units and avoids the positioning of volumetric sources required for volumetric coupling. 17,31,58 The anatomical connection between large arteries and the microcirculation is crucial in ischaemic stroke modelling because it determines the location and extent of ischaemic regions. Another advantage of placing the interface at the cortical surface is that anatomical connections are preserved at the cortical surface so that subcortical vessels can be homogenised to lower computational costs. In the case of volumetric coupling, subcortical arterioles must be added to the one-dimensional network model to account for such connections which can drastically increase the number of discretised branches.
Recently, the American Society of Mechanical Engineers (ASME) has introduced the Verification and Validation (V&V) 40 standard (titled ''Assessing Credibility of Computational Modeling through Verification and Validation: Application to Medical Devices''), 1 which details the requirements that a credible engineering software needs to satisfy to ensure safe applicability in the medical device industry. The present study thus sets out to establish guidelines regarding the efficient usage of a porous cerebral perfusion Finite Element model 39 both in terms of accuracy and computational cost. To this end, simulations will be compared to exact analytical solutions and a thorough sensitivity analysis will be carried out covering every numerical and physiological model parameter. The results prepare the ground for uncertainty quantification, virtual patient generation, and validation, which are essential to carry out reliable in silico trials of ischaemic stroke treatments.

METHODOLOGY
Inspired by recent advances in the fields of heart 35,36,50 and brain 24-26,71,72 haemodynamic modelling, a multi-compartment porous continuum framework is utilised to describe time-averaged cerebral blood flow in the microcirculation embedded within brain tissue. In this approach, microscale individual blood vessels are not captured. Instead, porous media are introduced to describe characteristic units of the vessel network, such as bundles of penetrating arterioles and the capillary mesh as shown in Fig. 1. During the macroscale continuum description, the corresponding vessels are grouped into multiple compartments . In practice, the corresponding vessels can be isolated, for instance, on a geometrical basis using diameter thresholding or branching order, or on a structural basis considering layers of the vessel walls (smooth muscle cells and endothelial cells). The remaining vessels, which are not part of any of these compartments and are not shown in Fig. 1, serve as a link between the compartments.
In the porous formulation, the primary flow variable is the volume-averaged (Darcy) pressure which can be deduced from the pressure distribution of the individual blood vessels by spatial averaging. The general partial differential equation set describing a steady-state multi-compartment porous Darcy system [24][25][26]35,36,50,71,72 is written in a compact form 1 as Einstein summation convention is not utilised in this study.

BIOMEDICAL ENGINEERING SOCIETY
Here, p i denotes the Darcy pressure in the ith compartment and K i is the corresponding permability symbolising flow conductance through pores (lumen of vessels) with well-separated length scales. Considering n compartments, i ¼ 1; 2; _ s n and b i;j is the intercompartment coupling coefficient matrix with n Â n elements representing conductance between compartments i and j. It follows that b i;j is arbitrary if Considering a domain of interest denoted by X with a boundary @X ¼ C ¼ C D;i [ C N;i , a set of Dirichlet type boundary conditions (DBC) and Neumann type boundary conditions (NBC) are formulated as Here, n is the outward-pointing unit vector perpendicular to the boundary surface. DBCs can be used to describe, for example arterial blood pressure along the cortical surface as shown in Figs. 1a and 1b. By comparison, NBCs are typically used to enforce a specific volumetric blood flow rate across boundary regions, such as zero blood flux through the ventricular surface in Figs. 1b-1d. (NBC could be also used to account for cerebrospinal fluid filtration at the choroid plexus.) Thereafter, the weak form can be derived by multiplying every component of Eq.
(1) with a non-zero test function t i , integrating over X and applying the divergence theorem: Z Equation (4) is discretised by the continuous Bubnov-Galerkin method. 34 The H 1 Sobolev space contains both the test and the trial functions so that t i 2 H 1 ðXÞ : t i ¼ 0 on C D;i , and ð5Þ Numerical solutions of p i are obtained with the open source FEniCS library 2,44,45 using first or second order Lagrange elements denoted by P 1 and P 2 , respectively. The permeability tensors K i and the coupling coefficients b i;j are captured by piecewise constant, zeroth order discontinuous elements (dP 0 ) unless stated otherwise. For simplicity, we restrict ourselves to a model with arteriole (i ¼ 1 ¼ a), capillary (i ¼ 2 ¼ c), and venule (i ¼ 3 ¼ v) compartments, and then generalise our findings whenever it is possible (Fig. 1). In this case, the remaining pre-capillary and post-capillary vessels are lumped into the arteriole-capillary (b 1;2 ) and capillary-venule (b 2;3 ) coupling coefficients, respectively. Numerical solutions of linear equation systems originating from the FE discretisation are obtained iteratively based on the BiCojungate Gradient STA-Bilised method (BiCGSTAB). 77 Computations are accelerated with an Algebraic MultiGrid (AMG) preconditioner. 66 Simulations were run on a desktop with an Intel Xeon E-2146G processor and 32 GB RAM unless stated otherwise. In order to verify the resulting model and carry out a comprehensive sensitivity analysis, three test cases are considered as detailed in the following subsections.
The domain of interest chosen here is a unit cube X ¼ ½0; 1 Â ½0; 1 Â ½0; 1 with periodic boundary conditions applied on four sides: p i ð0; y; zÞ ¼ p i ð1; y; zÞ, and ð8Þ p i ðx; 0; zÞ ¼ p i ðx; 1; zÞ: Equations (8) and (9) define four DBCs in every compartment and another one is added in the arteriole (i ¼ 1) and the venule (i ¼ 3) compartments based on the manufactured solutions such that Thereafter, the equation set is closed by prescribing homogeneous NBCs on the remaining surfaces (in Eq. (2) b i ¼ 0). The permeability tensors and coupling coefficient matrices corresponding to the manufactured solutions are detailed in Appendix 1. The chosen polynomial functions (7a)-(7c) are straightforward to differentiate and are the simplest functions which satisfy the imposed boundary conditions. The source terms S 1 , S 2 and S 3 satisfying Eq. (1) can be obtained by substituting Eqs. (7) and (31)-(32) into Eq. (1). The source terms are set to zero in every other test case (details in ''One-dimensional brain tissue column'' and ''Threedimensional human brain'' Section).

Three-Dimensional Human Brain
Details of the baseline human brain simulation are provided in Ref. 39 and therefore only the most important settings are covered here. The computational domain (X) is obtained by postprocessing a tetrahedral patient-specific head model utilised in multiple recent studies. 19,20,39,79 The geometry depicted in Fig. 2a is remeshed with tetrahedral elements using Tetgen. 68 The volumetric region of interest includes both grey matter (X G ) and white matter (X W ) subdomains, so that X ¼ X G [ X W as depicted in Fig. 2c. The bounding surface regions (@X) include a transverse cut-plane of the brainstem C BS , the ventricles C V and the pial surface C P so that @X ¼ C BS [ C V [ C P as depicted in Figs. 2a and 2b. The boundary region associated with the pial surface is subdivided into eight perfusion territories corresponding to major feeding arteries as shown in Fig. 2a. Perfusion territories have been identified using a voxelised atlas created based on vessel-encoded arterial spin labelling (ASL) perfusion magnetic resonance imaging (MRI). 28,29,30,51,53 Then, the surface region that is perfused, for instance, by the Right Middle Cerebral Artery (R-MCA) is denoted by C R-MCA . This approach ensures that blood arrives to the domain through cortical surface regions mirroring anatomical connections between large arteries and the microcirculation as shown in Fig. 2a.
The imposed baseline boundary conditions prescribe zero flow through the transverse cut-plane of the brainstem and the surface of the ventricles: Flow through the pial surface in the capillary compartment is zero: Because of the incompressible fluid flow model, the results depend solely on the cerebral perfusion pressure (CPP) defined as the arteriole-venule pressure difference at the pial surface: CPP ¼ p a À p v . Setting the zero level of the pressure at the outlet of the venous compartment (p v ) leads to Pressure values in the case of 3D brain simulations are presented relative to the venous outlet pressure p v . To account for totally occluded scenarios, blood flow through the perfusion territory of an occluded vessel is set to zero whereas surface pressure is assumed to remain constant in other regions. Accordingly, a R-MCA occlusion is modelled with zero flux through the corresponding cortical territory (C RÀMCA ) so that the resulting mixed boundary conditions become @p a @n ¼ 0 on C R-MCA , and ð14aÞ The model is parametrised based on several simplifications as discussed in Ref. 39. Considering that the arteriole-capillary (b G 1;2 ¼ b G 2;1 ) and capillary-venule (b G 2;3 ¼ b G 3;2 ) coupling coefficients are known in the grey matter, it is assumed that the ratio of grey and white matter coupling coefficients (b G =b W ) is constant. It has been demonstrated based on microscale one-dimensional haemodynamic network simulations that the capillary permeability tensor is isotropic and characterised by a scalar k 2 . 16 It is hypothesised that a reference coordinate system ½n, g, f can be found at every point so that f is parallel to the local axis of the descending arterioles and ascending venules. In this coordinate system, the anisotropic arteriole (K ref 1 ) and venule (K ref 3 ) permeability tensors are modelled as tensors with a single non-zero diagonal element (k 1 and k 3 ) because they represent vessel bundles: The ½n, g, f coordinate system is determined assuming that penetrating vessels grow from the pial surface to the ventricles following a vector field. This vector field is computed as a gradient of a scalar field governed by a single diffusion equation. 39 It is worth noting that no other methods have been proposed to obtain these anisotropic permeability fields even though they play a crucial role in predicting perfusion response to vessel occlusion.

One-Dimensional Brain Tissue Column
Equation (1) can be simplified so that a one-dimensional problem with homogeneous permeabilities can be posed and defined by a set of ordinary differential equations in the form of Here, the permeability of each compartment (k i ) is a scalar. Introducing the dimensionless variables x Ã ¼ x=l s and p Ã ¼ p=p s based on a length scale (l s ) and a pressure scale (p s ) leads to a dimensionless form of Eq. (16) as The dimensionless matrix group a Ã i;j ¼ l 2 s b i;j =k i provides a similarity condition for Eq. (16).
Equation (17) describes n second order ordinary differential equations. This system can be replaced by 2n first order ordinary differential equations once the pressure gradient q Ã i ¼ dp Ã i =dx Ã is introduced. Thereafter, the unknown vector contains the pressure p Ã i and pressure gradient q Ã i functions for each compartment. Finally, a matrix differential equation based on Eq. (17) can be formulated as Here, the prime ( 0 ) denotes spatial differentiation and the A matrix is where every submatrix has a dimension of n Â n, and I is the identity matrix. The connection between the coupling coefficient matrix a ¼ a i;j and the submatrix B is given by Here, the n dimensional vector e n is filled with ones, and the ''diag'' operator is used to turn a column vector into a diagonal square matrix. The solution of Eq. (19) is Here, k j and v j are the jth eigenvalue and eigenvector of A, respectively. R is used to take the real part of the expression and to handle complex eigenvalue and eigenvector pairs. The C j coefficients are determined based on combinations of DBCs and NBCs at x b . A DBC in compartment i reads as Based on the unknown vector defined in Eq. (18) and the solution expression Eq. (22), a DBC can be satisfied by imposing Here, r i ¼ p Ã i denotes the ith element of r and v i j is the ith element of jth eigenvector. A NBC in compartment i is defined as The second half of the unknown vector defined in Eq. (18) includes the pressure gradients. Therefore, a NBC can be satisfied by imposing Equations (24) and (26) provide a single equation for each boundary condition. Considering n compartments, the 2n boundary conditions lead to 2n linear equations based on the 2n coefficients of the solution expressions (C j ). The resulting linear equation system determines uniquely the coefficients of the solution expressions (C j ). Using the DBC (23) and the NBC (25), the procedure can be extended to handle multiple subdomains with different coupling coefficients. In such cases, interface conditions have to be applied to ensures that the solution is continuously differentiable. Therefore, this analytical method is suitable to obtain solutions for a column of brain tissue stretching between the cortical and ventricular surfaces. The method can handle a domain including grey and white matter with identical permeabilities but different coupling coefficient matrices.

Reflection on the ASME V&V40 Framework
The INSIST software suite 40 aims to quantify the efficacy of thrombolysis 43 and thrombectomy 4,18 using computational models. Following the framework of the ASME V&V40 standard shown in Fig. 3, the question of interest can be phrased as ''what is the best available treatment option in terms of functional outcome using stent retrievers and tissue plasminogen activator?'' The context of use of the porous finite element model is quantifying how cerebral blood flow is altered in the microcirculation during ischaemic stroke and post-treatment compared to the healthy baseline case. The output of the FE model is the spatial distribution of brain tissue perfusion, which can be measured, for example, by ASL perfusion MRI in clinical settings. 29 The porous FE model provides (i) pressure and volumetric blood flow rate inputs for haemodynamic simulations in large arteries 54 which evaluate forces on the thrombus; and (ii) perfusion input for tissue health computation used for infarct volume estimation. 84 Therefore, the output of the FE model is linked indirectly to a statistical module estimating the functional outcome of individual virtual patients based on the computed infarct volume, and other features, such as age.
Considering model risk, the final decision consequence is high because the treatments of interest have a substantial impact on patients' abilities (evaluated according the modified Rankin Scale 60,70 ). Model importance is estimated medium because the results are not the only factors in decision making, but they impact large artery simulation directly and functional outcome estimation indirectly in the envisaged INSIST pipeline. 40 The verification and validation plan relies on comparisons with analytical solutions and clinical data, respectively. Previously, preliminary validation was presented highlighting a promising agreement with clinical data regarding perfusion and infarct volume estimations during ischaemic stroke. 39 This marked the first step towards establishing model credibility. Beyond a thorough verification, the following sections underpin a detailed validation with multiple virtual patients by determining the parameters which ensure that the pressure, volumetric flow rate, and perfusion estimates provided by the porous FE model are accurate and physiologically realistic.

Manufactured Solutions
The suitability of the finite element method is evaluated based on the manufactured solution shown in Fig. 4. These synthetic solutions are reminiscent of a cortical column: the flow is driven from high pressure arterioles to low pressure venules. A uniform tetrahedral mesh is generated with two elements along each edge of the cube and uniformly refined four times to carry out a grid convergence study using first (P 1 ) and second order (P 2 ) elements as summarised in Table 1. The analytical solution and the finite element approximation along the sampling line depicted in Fig. 4 are shown in Fig. 5a, highlighting that a very good agreement is achievable. Figure 5b presents the L 2norm defined over the entire domain based on the difference between the analytical and the numerical solutions. With increasing spatial resolution, the L 2norm decreases as expected. Superconvergence of the solution with first order elements can be observed which is a well-known feature of this FE approximation. 78 According to Fig. 5b, solutions achieved using a second order scheme on a given mesh have about the same L 2 -norm as first order approximations after a single refinement, as suggested by the number of degrees of freedoms in Table 1. However, it is important to emphasise that the convergence of pointwise variables does not necessarily follow the convergence trend of integral quantities. Figure 5c suggests that the second order FE scheme outperforms the first order scheme in terms of gradient estimation even if computations with the latter are carried out with a threefold refined mesh. Beyond accuracy, the computational cost of simulations as a function of the spatial resolution and the element order can be examined in Fig. 5d. In general, computations with second order elements take one order of magnitude longer and are comparable to the computational cost of simulations with first order elements on a refined mesh (see also Table 1).

Whole Brain Model
Next, simulations are carried out with the whole brain model to investigate how uncertainty of the gradient estimation near the boundaries impacts the results. Two scenarios are considered using previously optimised settings 39 summarised in Table 2: baseline (healthy) and an RMCA occlusion. The algorithm Normalised Grid Spacing (NGS), number of elements along each edge of the unit cube (n edge ); number of elements (n ele ); number of degrees of freedom (n DoF ); wall time required for the iterative solution of the linear system (t inv ).

BIOMEDICAL ENGINEERING SOCIETY
used for parameter optimisation is described in Appendix 2. In a porous framework, the pressure gradient is directly connected to the Darcy velocity vector which can be introduced in compartment i as From here, the overall volumetric flow rate entering the arteriole compartment through the pial surface C P is simply defined as     The Cn product corresponds to an outward-pointing vector with magnitude equal to the area of the C boundary surface. (The Darcy velocity has non-zero components perpendicularly to the boundary surface only at C P and therefore in Eq. (28) the integration domain C P is replaced with C for a more concise notation.) Therefore, any uncertainty in the pressure gradient propagates directly to the volumetric flow rate. Accurate flow rate computation is important for two reasons: (i) superficial arteriole pressure and flow rates are cornerstones of a coupling scheme between the present continuum model and a network model of large arteries as detailed in ''Introduction'' section; (ii) surface fluxes associated with major cerebral arteries can be directly compared with values obtained from phasecontrast MRI 85 for validation purposes.
The multi-compartment formulation offers an alternative method to compute volumetric flow rate through the brain. Introducing perfusion F ¼ b 1;2 ðp 1 À p 2 Þ, and integrating it over the brain volume X leads to According to mass conservation Q X ¼ Q C so the imbalance between the two formulations is present because the utilised FE method is not conservative. The present H 1 (pressure) and L 2 (velocity) spaces are selected for straightforward implementation, robustness, and relatively low memory requirements. Nevertheless, it is worth noting that numerical solutions of the Poisson equation using the mixed HðdivÞ Â L 2 FE formulation 48,62 are conservative, similarly to the finite volume method. 75 To mitigate the mass conservation error of the present model, two strategies are employed: (i) grid refinement; and (ii) increasing FE orders of both the pressure and the velocity. For P 1 and P 2 pressure elements, the natural velocity pair is dP 0 and dP 1 but higher order velocity approximations can be obtained by projection. The results in Table 3 summarise the tested cases and show that the volumetric flow rate from Eq. (29) depends on the grid size and the element order only negligibly. The error in Q X should be relatively low in the investigated cases because the computed values overlap independently from the spatial resolution. For this reason, and because of the mass conservation principle, the relative difference, is used for error measurement. Figures 6a and 6b show the imbalance between Q X and Q C highlighting that without improved u i approximation the surface flux is always underestimated by Q C . Both refinement and increased FE order can mitigate the error of Q C . However, higher order pressure and velocity approximations are required to keep jD r Q C j<5%. Figure 6c displays the wall time of the simulations emphasising the increased computational cost required for more accurate flux estimations. In summary, second order velocity and pressure elements on a relatively coarse mesh (about 10 6 elements) can provide a reasonable compromise between accuracy and performance if computing Q C is important. Otherwise, perfusion and pressure distributions can be estimated reasonably well on coarse meshes with about 1 million elements.
Once it is established under what conditions the present 3D model can provide accurate estimation of volumetric flow rate through cortical surface regions, it is possible to compute volumetric flow rate of the major cerebral arteries. Simulation results are compared with values from a clinical study in Table 4. Whereas the model is optimised for overall brain perfusion, blood flow rate through these major arteries provides a chance for independent validation. It is promising that volumetric flow rate in both the ACA and the MCA is predicted within one standard deviation of the experimental values. Volumetric flow rate in the PCA is somewhat overestimated. This simple comparison is an important milestone in the quantitative validation of the 3D brain model and provides guidance about development directions for the future.
The results emphasise the importance of accurate mapping between voxelised vascular atlases and the triangulated brain surface and draw attention to geometrical flaws. The brain model does not capture the lateral sulcus and the insula. Because of these neglected regions, the brain surface is underestimated whereas the brain volume is overestimated (see Table 4). It is anticipated that the 3D porous model will provide more reliable volumetric flow rate estimations once these geometrical issues are overcome. However, obtaining accurate discretised brain geometries suitable for FE computations remains a challenging task, partially because of labelling the boundary regions, and therefore such efforts are left to a future study.

Boundary Conditions
Results from the previous section draw attention to issues related to flux computation in a porous FE framework. Instead of constant pressure, uniform velocity has been used as an inlet BC [24][25][26]71,72 where the velocity is calculated as the ratio of a given volumetric flow rate and the inlet surface area. This NBC at the inlet might distort the results because of errors related to gradient estimation so we next investigate the difference between pressure and velocity inlets using the present model. To this end, a single simulation with constant inlet velocity is run with the same settings as the case in Table 3 marked with ? so that the total inlet volumetric flow rate remains unchanged. Q X ¼ 600 (ml/min) is preserved, indicating that the simulations with pressure and velocity inlets are equivalent regarding average brain perfusion. Q C is impacted only slightly: pressure inlet leads to 559 (ml/ min) whereas velocity inlet results in 572 (ml/min) (both values should be 600 (ml/min) based on mass conservation). The effects of the inlet boundary condition on other statistics regarding the pressure and the perfusion fields are shown in Fig. 7. Most of the volume averaged pressure (hp i i) and perfusion (hp i i) values are altered only slightly by the inlet boundary conditions. However, extrema of both the pressure and the perfusion fields change substantially when the inlet velocity BC is used. Minimum values are underestimated so that the relative difference of the minimum arteriole pressure and perfusion values are approximately À100% because min p 1 and min F) are close to unif. ref. (a) unif. ref. unif. ref.
(c) FIGURE 6. Relative error of the superficial volumetric flow rate estimation as functions of the spatial resolution and the FE order in the case of the baseline (a) and RMCA occlusion (b) scenarios. Computational time averaged between the healthy and the occluded cases (c). The ''loc. ref.'' and ''unif. ref.'' abbreviations correspond to locally and uniformly refined meshes, respectively (see Table 3 for further details).   zero when uniform inlet velocity BCs are used. By comparison, the inlet velocity BCs lead to overshot maximum values nearly by a factor of 10 compared to the case with inlet pressure BCs (see max p 1 and max F in Fig. 7). In Figs. 8a-8d, and 8b-8e, the pressure and perfusion fields corresponding to inlet pressure and velocity BCs, respectively, emphasise further the difference between the two solutions. The arteriole velocity magnitude displayed over the cortical surface (Fig. 8c) shows that in the case of pressure inlet BC the velocity is strongly inhomogeneous. It is a remarkable feature of this scenario that high velocity values are predicted along the major cerebral arteries branching from the circle of Willis because the whole brain model does not include any information about the location of these vessels. The solution obtained with constant pressure BC suggests that descending arterioles branching from the major cerebral arteries play a key role in providing well-balanced grey and white matter perfusion. The velocity distribution over the cortical surface is not known a priori and therefore it is challenging to impose feasible velocity distributions with velocity inlet BC.
The volumetric flow rates corresponding to boundary regions of the major cerebral arteries are summarised in Table 4. Based on previous clinical FIGURE 8. Pressure p a (a-d) and perfusion F (b-e) distributions along a coronal plane, and velocity magnitude (c-f) along the cortical surface. Constant pressure inlet (a-c) and constant velocity inlet (d-f). In (f), the non-uniform velocity magnitude is associated with numerical errors. Results correspond to the mesh with local grid refinement at the boundaries and P 2 pressure and velocity elements (details in Table 3).
overestimated underestimated FIGURE 7. Absolute value of the relative difference between some statistics extracted from simulations carried out with inlet pressure and velocity BCs. Average arteriole pressure over the pial surface hp 1 i CP ; average pressure in the arteriole (hp 1 i), capillary (hp 2 i), and venule (hp 3 i) compartments over the entire brain; minimum (min) and maximum (max) arteriole pressure (p 1 ) and perfusion (F); average perfusion in grey (hF i G ) and white (hF i W ) matter. The colour and pattern of each bar indicates whether results with the inlet velocity BCs overshoot or underestimate the reference values corresponding to the pressure inlet BCs.
investigations, 83,85 it is clear that each major cerebral artery contributes differently to volumetric blood flow rate to the brain. Even if we assume that superficial perfusion territories of cerebral arteries are proportional to the associated blood flow rate, it is important to recognise that the volumetric flow rate should be heterogeneous as dictated by underlying grey and white matter volumes. Although the validity of homogeneous inlet velocity is physiologically questionable, prescribing constant pressure inlets in healthy scenarios is supported indirectly by clinical data. Based on ASL MRI, 29 cerebral perfusion regions have been found to be well-separated indicating that the pressure gradient between these territories is small. Furthermore, dynamic computed tomography angiograms indicate that flow through leptomeningeal collaterals is activated as a result of cerebral artery occlusion. 69 We hypothesise that in healthy cases, leptomeningeal collaterals 8,28,42,43,46,69 as well as pial arterioles 6,7,67 act to equalise blood pressure in descending arterioles near the cortical surface. This hypothesis could be tested once the present perfusion model and the 1D model incorporating the larger arteries 54 are coupled, and it is supported by pressure data obtained in rodents which suggest small pressure drop in large arteries compared to the microcirculation. 27,65,80 In summary, there is no evidence that large arteries distribute blood flow uniformly along the cortical surface. Prescribing uniform inlet velocity as a boundary condition in the arteriole compartment is questionable, and leads to extreme pressure and perfusion values out of the expected range. Based on anatomical and physiological considerations, uniform inlet pressure appears to be a more suitable boundary condition in healthy scenarios.

Sensitivity Analysis
Determining the actual values of parameters in physiological models is a challenging task because of limited experimental data, high patient-specific variability, and the typically high dimensionality of the parameter spaces. Considering porous cerebral haemodynamics, this issue is far from being settled as indicated by the broad range of reported parameters summarised in Table 2. With increasing computational costs, mapping a high dimensional parameter space quickly becomes unfeasible. Beyond parameter optimisation, sensitivity analysis and the associated uncertainty quantification are other tasks which require large batches of simulations. Here we aim to utilise the one-dimensional model described in ''Onedimensional brain tissue column'' section in order to ease these tasks.

One-Dimensional Brain Tissue Column
The one-dimensional multi-compartment porous model enables the investigation of tissue columns including both grey and white matters so that the corresponding model parameters are identical with the ones used for the three-dimensional brain model (Table 2). In the 1D case, x ¼ 0 and x ¼ l correspond to the pial and the ventricular surfaces, respectively. Therefore, the same boundary conditions can be imposed in the 1D tissue column and the 3D brain. (Similarly to the 3D brain, pressure values in the 1D case are presented relative to the venous outlet pressure p v .) Thereafter, the one-dimensional model has only two free parameters: the lengths of the grey (l G ) and the white (l W ) matter domains defined so that l ¼ l G þ l W . The equivalent l G ¼ 13:55 and l W ¼ 7:99 (mm) is uniquely determined by the human brain model used here and can be obtained by optimisation. The corresponding cost function is zero if the one-dimensional model returns hFi ¼ 43, hFi G ¼ 56, and hFi W ¼ 21 (ml/min/100 mL) which overlap exactly with volume-averaged perfusion values of the threedimensional brain model. Figure 9 displays the analytical pressure distributions in the arteriole, capillary and venule compartments. A one-dimensional FE approximation with 100 P 1 elements is also shown which is in excellent agreement with the analytical results. The average pressure values provided by the 1D model are within AE4% of those computed using the 3D brain in every case. Onedimensional simulations are approximately 3000 times faster compared to their 3D counterparts.

Mapping the Parameter Space
The 1D model is a good candidate for parameter space mapping because of its low computational cost. A disadvantage of the analytical formulation is that with the parameters presented in Table 2, the linear equation system governing the boundary conditions is very stiff, probably because the parameters have different orders of magnitude. The condition number depends nonlinearly on the reference length scale used for the nondimensionalisation. To avoid the issue of suitable reference length scale selection, we will use the 1D finite element model verified in Fig. 9 for a detailed sensitivity analysis.
Next, each parameter is perturbed independently from others. The one-at-a-time approach is applicable and the corresponding sensitivity analysis is representative of the system behaviour because the governing equation set (1) is linear. 64 The parameter range is [10%;1000%] compared to the values in Table 2 with 101 samples along each dimension. In addition, a geometrical scaling factor is introduced based on the ratio of the original and the perturbed volume of the computational domain. In order to evaluate whether the 1D model can capture the system behaviour, 11 samples along each dimension of the parameter space are tested with the 3D brain model using P 1 pressure elements.
Results are summarised in Fig. 10 showing brain perfusion as a function of the model parameters in the healthy scenario. Based on the 3D simulations, RMCA occlusion is also considered so that the perfusion in the occluded scenario and infarcted volume change are also displayed. Infarcted volume is estimated solely based on the perfusion distribution. A region with more than 70% perfusion drop compared to the baseline scenario is identified as an infarcted zone inspired by perfusion imaging-based methods. 9,38 The obtained value estimating infarcted volume is used here solely as an indicator of perfusion drop severity because this method overestimates the ischaemic region compared to diffusion-weighted magnetic resonance imaging. 11 According to Figs. 10a-10h, the 1D solutions follow the trend of the 3D simulations in a wide parameter range in every case except the geometrical scaling factor. Considering that the present 1D model does not account for curvature effects, the results could potentially be improved with a cylindrical or spherical formulation. The sensitivity of each parameter is estimated as the normalised partial derivative of brain perfusion and infarcted volume at 100% of the abscissa (around baseline values presented in Table 2). The derivatives computed based on the 3D simulations in Fig. 10a-10h are listed above the abscissa and indicate the percentage change of brain perfusion and infarcted volume as a response to unit percentage change in model parameters.
In terms of haemodynamics, it should be noted that the present model captures solely the resistive part of the microcirculation and neglects conductance and inductance associated with vessel wall stiffness and pulsatility. The linear relationship between brain perfusion and perfusion pressure is a direct manifestation of the ''hydraulic Ohm's law'' (Fig. 10a). As other parameters are increased, brain perfusion tends to exhibit a saturating behaviour even if it is not obvious within the displayed parameter range (Figs. 10b-10g). Considering the infarcted volume fraction, the insensitivity to CPP change originates from its definition restricted to relative perfusion change (as CPP is increased perfusion in the healthy and the occluded cases increases by the same amount).
Perturbations in the capillary and venous permeabilities have a relatively weak impact on both F b and IV. Fluid flow within the capillary compartment is relatively small. This result can be interpreted as flow ''avoiding'' the capillaries because of their high resistance, which overlaps with detailed network simulation of the rat microcirculation. 65 Considering an indefinite domain with solely grey or white matter and constant but different arteriole and venous pressures, blood flow is non-zero between the compartments but it is zero inside each compartment. In the microcirculation, capillaries cannot be bypassed so it is important to recognise that the coupling coefficients must account not only for pre-and post-capillaries but also for those capillary vessels which establish the shortest route between arterioles and venules. Therefore, the results show strong sensitivity to changes in coupling coefficients, which are responsible for the majority of the system resistance and consequently most of the pressure drop. The geometry scaling factor has a strong influence on the infarcted volume fraction. If one of the major feeding arteries is occluded, it seems logical that for a given purely resistive vasculature a larger brain suffers more compared to a smaller brain. Such effects might be observable in clinical data but they are unlikely to be statistically significant compared to other factors, such as the extent of collateralisation. Nevertheless, the results suggest that a representative brain geometry is essential to keep the numerical uncertainty of the simulations low. In addition, the data presented in Fig. 10, and the one-dimensional model particularly, can help to accelerate virtual patient generation, where parameter sets must be found so that the resulting perfusion distribution is representative of a patient cohort.
The sensitivity analysis shown in Fig. 10 emphasises the importance of the (i) boundary conditions and (ii) coupling coefficients. Based on the literature data in Table 2, the coupling coefficients are burdened by the highest uncertainty which can undermine the applicability of the present model in terms of perfusion and infarct volume estimation through uncertainty propagation. (a) 1.000 0.000 FIGURE 10. Sensitivity analyses of the 1D and 3D models. Brain perfusion (left axes) and infarcted volume fraction (right axes) as functions of change in the following parameters: cerebral perfusion pressure (a); arteriole (b), capillary (c), and venule (d) permeabilities; arteriole-capillary (e) and capillary-venule (f) coupling coefficients in the grey matter; ratio of grey and white matter coupling coefficients (g); geometry scaling factor (h). The top (black) and bottom (red) numbers on each subplot stand for the sensitivity of healthy brain perfusion and infarcted volume fraction, respectively. Reference values of the parameters are listed in Table 2.

Limitations
This subsection summarises some of the most important overlooked factors in the present study. Investigations are limited thus far to a single patientspecific brain geometry because of challenges associated with mesh generation. The porous model relies on multiple scale separation even though the vessel diameter in the vasculature changes continuously. 58 The present study demonstrated promising quantitative validation based on blood flow rate in major cerebral arteries but a more comprehensive validation is essential to ensure wider, multi-purpose applicability of the model.
The study is restricted to a porous model with three compartments (arteriole, capillary, venule) which are connected by homogeneous coupling coefficients both in grey and white matter. The model cannot represent arterioles which branch relatively far away from the cortical surface and hence feed deeper tissue regions. For this reason, low perfusion regions caused by occluded cerebral arteries are always connected to the brain surface so that isolated white matter infarction near the ventricles cannot be captured. On the one hand, such effects could be resolved by model expansion, for instance, by introducing two arteriole compartments which perfuse grey and white matter separately and rely on different boundary conditions. On the other hand, models with a single or two 31 compartments might be found more suitable in certain cases. Considering that the present study focuses on the sensitivity analysis of a three-compartment brain model for acute ischaemic stroke, exploring alternative porous frameworks for other specific (patho)physiological problems is left to future investigations.
From the haemodynamics point of view, the present model is purely resistive and aims to capture statistically steady state (time-averaged) brain perfusion. Therefore, the following physiological phenomena have been neglected: vessel wall stiffness (conductance), unsteady flow phenomena (inductance), 10 and cerebral autoregulation. 57 The associated processes and mechanisms are strongly time-dependent and often rely on nonlinear processes leading to significant complexity which is beyond the limits of the present study. Nevertheless, the authors hope that the present work will contribute to the creation of a comprehensive cerebral blood flow model which can incorporate these effects, in addition to pathophysiological processes, such as cerebral oedema, 73 emboli advection and blockage of the microcirculation, 3,15,49 and spreading of ischaemic tissue damage. 63,74,84

Conclusions
The present study set out to verify and explore the sensitivity of a porous brain perfusion model and to prepare the ground for validation and uncertainty quantification as specified in the ASME V&V40 standard. 1 To this end, simulations were carried out using the finite element method and three test cases: a unit cube with manufactured solutions, a human brain geometry, and a one-dimensional brain tissue column. For the unit cube and the brain tissue column, analytical solutions were presented and used for verification. This approach allowed direct comparison between analytical and numerical results and enabled the identification of certain settings which were found crucial to minimise numerical errors.
Based on the key findings of this study, the following guidelines are established regarding the efficient usage of three-dimensional porous finite element models for steady state cerebral perfusion simulation: Blood pressure and perfusion can be estimated reasonably well on a relatively coarse grid and first order pressure elements.
In order to approximate volumetric flow rate through superficial territories accurately with minimal computational effort, higher order finite elements are essential. Volumetric flow rate estimation can be further improved with local grid refinement in the vicinity of the superficial territories.
Porous models can provide a good estimation of blood flow rate through major cerebral arteries when higher order elements and constant inlet pressure boundary conditions are used.
Using uniform velocity inlet boundary conditions does not have solid physiological foundations because blood flow is likely not distributed proportionally to the cortical surface area. However, when considering healthy scenarios, constant pressure inlets can be prescribed with higher confidence.
A multi-compartmental porous system can be solved analytically in one-dimensional cases. After determining the equivalent thickness of the subdomains, the one-dimensional system behaviour is representative of the three-dimensional organ-scale model. The one-dimensional case is suitable for parameter mapping and sensitivity analysis at a computational cost three orders of magnitude smaller compared to three-dimensional models.
The results lay down the foundations for validation, uncertainty quantification, and for the development of a multiscale model which incorporates large arteries as BIOMEDICAL ENGINEERING SOCIETY well as the microcirculation. The present work contributes to the creation of a reliable computational model which can assist in the design of clinical trials and in clinical decision making related to cerebrovascular diseases.

APPENDIX 1: PERMEABILITIES AND COUPLING COEFFICIENTS FOR THE MANUFACTURED SOLUTIONS
To obtain a synthetic pressure field, the permeability tensors are set to whereas the non-zero elements of the coupling coefficient matrix are b 1;2 ¼ b 2;1 ¼ 1:5 þ 0:5 tanhð10z À 5Þ; and ð32aÞ b 2;3 ¼ b 3;2 ¼ 3 þ tanhð10z À 5Þ: The hyperbolic tangent functions in the elements of the coupling coefficient matrix are selected to capture the sudden change in the material properties of grey and white matter. K i and b i;j are approximated by fourth order polynomials for the finite element computations.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://crea tivecommons.org/licenses/by/4.0/.