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.

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 INSIST 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 arteries5,54 and a three-dimensional (3D) porous perfusion model39 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–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.

Computational porous and poroelastic models have been applied over a number of decades and have been proposed to describe aspects of cardiac36,35,50 and cerebral biomechanics.24–26,41,71,72 Such organ-scale computational models have been applied to capture both healthy and pathophysiological states, such as brain injury,21,28,31 oedema,73 conditions associated with dementia,25,26,72 and ischaemic stroke.39 The majority of these models utilise the Finite Element (FE) method21,24–26,33,39,71–73,81 with relatively few exceptions.17,31,58 The present study focuses on a porous perfusion FE model developed specifically to capture blood flow changes in ischaemic stroke39 as part of the in silico trial pipeline of INSIST.40 Even though organ-scale porous and poroelastic brain models are becoming increasingly popular,12,13,21,24–26,31,33,39,41,81,71–73 efforts regarding their verification, validation,24 comprehensive sensitivity analysis,31,64 and uncertainty quantification12,13 remain relatively rare and there is thus a need to fill this gap before the models are applied for clinical trial design and decision making.

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 model39 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 heart35,36,50 and brain24–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 (Figs. 1b-1d). 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.

Figure 1

Schematic of the microcirculation along a slice in a cortical column: descending arterioles (red), capillaries (black) and ascending veins (blue) are shown in both grey (GM) and white matter (WM) (a). The porous cerebral blood flow model relies on the idea of separating these vessels into coexisting arteriole (b), capillary (c), and venule (d) 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 system24–26,35,36,50,71,72 is written in a compact form1 asHere, denotes the Darcy pressure in the compartment and is the corresponding permability symbolising flow conductance through pores (lumen of vessels) with well-separated length scales. Considering n compartments, and is the intercompartment coupling coefficient matrix with elements representing conductance between compartments i and j. It follows that is arbitrary if . A volumetric source term in the compartment is defined by .

Considering a domain of interest denoted by with a boundary , a set of Dirichlet type boundary conditions (DBC) and Neumann type boundary conditions (NBC) are formulated asHere, 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 , integrating over and applying the divergence theorem:Equation (4) is discretised by the continuous Bubnov-Galerkin method.34 The Sobolev space contains both the test and the trial functions so thatNumerical solutions of are obtained with the open source FEniCS library2,44,45 using first or second order Lagrange elements denoted by and , respectively. The permeability tensors and the coupling coefficients are captured by piecewise constant, zeroth order discontinuous elements () unless stated otherwise. For simplicity, we restrict ourselves to a model with arteriole (), capillary (), and venule () 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 () and capillary-venule () coupling coefficients, respectively.

Numerical solutions of linear equation systems originating from the FE discretisation are obtained iteratively based on the BiCojungate Gradient STABilised 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.

Three-Dimensional Synthetic Data

To derive manufactured solutions,61 we assume solution functions of Eq. (1) of the form The domain of interest chosen here is a unit cube with periodic boundary conditions applied on four sides:Equations (8) and (9) define four DBCs in every compartment and another one is added in the arteriole () and the venule () compartments based on the manufactured solutions such thatThereafter, the equation set is closed by prescribing homogeneous NBCs on the remaining surfaces (in Eq. (2) ).

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 , and 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 “Three-dimensional 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 () 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 () and white matter () subdomains, so that as depicted in Fig. 2c. The bounding surface regions () include a transverse cut-plane of the brainstem , the ventricles and the pial surface so that 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 . 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.

Figure 2

Coronal views of boundary regions and subdomains of the human brain model19: (a) pial surface; (b) ventricles and cut-plane at brainstem; (c) GM and WM along a slice.

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: . Setting the zero level of the pressure at the outlet of the venous compartment () leads toPressure values in the case of 3D brain simulations are presented relative to the venous outlet pressure . 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 () so that the resulting mixed boundary conditions become

The model is parametrised based on several simplifications as discussed in Ref. 39. Considering that the arteriole-capillary () and capillary-venule () coupling coefficients are known in the grey matter, it is assumed that the ratio of grey and white matter coupling coefficients () 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 .16 It is hypothesised that a reference coordinate system , , can be found at every point so that is parallel to the local axis of the descending arterioles and ascending venules. In this coordinate system, the anisotropic arteriole () and venule () permeability tensors are modelled as tensors with a single non-zero diagonal element ( and ) because they represent vessel bundles:The , , 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 ofHere, the permeability of each compartment () is a scalar. Introducing the dimensionless variables and based on a length scale () and a pressure scale () leads to a dimensionless form of Eq. (16) asThe dimensionless matrix group 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 is introduced. Thereafter, the unknown vectorcontains the pressure and pressure gradient functions for each compartment. Finally, a matrix differential equation based on Eq. (17) can be formulated asHere, the prime () denotes spatial differentiation and the matrix iswhere every submatrix has a dimension of , and is the identity matrix. The connection between the coupling coefficient matrix and the submatrix is given byHere, the n dimensional vector 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) isHere, and are the jth eigenvalue and eigenvector of , respectively. is used to take the real part of the expression and to handle complex eigenvalue and eigenvector pairs.

The coefficients are determined based on combinations of DBCs and NBCs at . A DBC in compartment i reads asBased on the unknown vector defined in Eq. (18) and the solution expression Eq. (22), a DBC can be satisfied by imposingHere, denotes the ith element of and is the ith element of jth eigenvector. A NBC in compartment i is defined asThe second half of the unknown vector defined in Eq. (18) includes the pressure gradients. Therefore, a NBC can be satisfied by imposingEquations (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 (). The resulting linear equation system determines uniquely the coefficients of the solution expressions (). 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.

Results and Discussion

Reflection on the ASME V&V40 Framework

The INSIST software suite40 aims to quantify the efficacy of thrombolysis43 and thrombectomy4,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 arteries54 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.

Figure 3

The ASME V&V40 software credibility assessment framework detailed in Ref. 1.

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 Scale60,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.

Spatial Resolution

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 () and second order () 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 -norm defined over the entire domain based on the difference between the analytical and the numerical solutions. With increasing spatial resolution, the -norm 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

Figure 4

Manufactured solution in a unit cube visualised in the arteriole (a), capillary (b), and venule (c) compartments with a sampling line (0.5,  0.5 , z).

Table 1

Results corresponding to the grid convergence study using manufactured solutions.

NGS	n_edge	n_ele	Element order	n_DoF	L^2-norm	∂p_1/∂zat (0.5, 0.5, 1)	t_inv(s)
16	2	48	1	81	0.15174	1.273	0.0042
8	4	384	1	375	0.08714	1.800	0.0287
4	8	3072	1	2187	0.02857	1.867	0.2278
2	16	24,576	1	14,739	0.00771	1.927	1.8428
1	32	196,608	1	107,811	0.00197	1.965	15.0789
16	2	48	2	375	0.06215	1.824	0.0161
8	4	384	2	2187	0.01399	2.006	0.1403
4	8	3072	2	14,739	0.00266	2.008	1.4050
2	16	24,576	2	107,811	0.00060	1.999	16.3241
1	32	196,608	2	828,375	0.00014	1.998	162.1361

Normalised Grid Spacing (NGS), number of elements along each edge of the unit cube (); number of elements (); number of degrees of freedom (); wall time required for the iterative solution of the linear system ()

Figure 5

Exact manufactured solutions and FE approximations with second order (ord.) elements and normalised grid spacing (NGS) 2 along a line (a). Numerical error of the FE approximation (b), FE approximation and exact value of the arteriole gradient at a single point (c), and wall time of the linear system matrix inversion () (d) as functions of the NGS and element order.

According to Fig. 5b, solutions achieved using a second order scheme on a given mesh have about the same -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 settings39 summarised in Table 2: baseline (healthy) and an RMCA occlusion. The algorithm 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 asFrom here, the overall volumetric flow rate entering the arteriole compartment through the pial surface is simply defined asThe product corresponds to an outward-pointing vector with magnitude equal to the area of the boundary surface. (The Darcy velocity has non-zero components perpendicularly to the boundary surface only at and therefore in Eq. (28) the integration domain is replaced with for a more concise notation.)

Table 2

List of model parameters and some reference values ( deviation when available).

Parameter	Value	Reference	Unit	Optimised
k_1	1.234	431; 0.021 ± 0.01915	mm3 s kg−1	Yes
βi,jG/βi,jW	2.538	1.631; 1.024–26	–	Yes
k_2	4.28 × 10−4	4.28 × 10−416;	mm3 s kg−1	No
k_3/k_1	2	231	mm3 s kg−1	No
β1,2G	1.326 × 10−6	1.5 × 10−1924,26; 1.5 × 10−725;	Pa−1 s−1	No
		5 × 1031; (4.13 ± 5.49) × 10−515
β2,3G	4.641 × 10−6	1.5 × 10−1924,26; 1.5 × 10−725; 5 × 10−631	Pa−1 s−1	No
CPP	75	78 ± 1047	mmHg	–
β1,2G+β2,3G	5.967 × 10−6	3 × 10−1924,26; 3 × 10−725; 10−531	Pa−1 s−1	–

The distribution of the parameters reported in Ref. 16 based on microscale vessel network simulations is not normal. The last column indicates which parameters are optimised as detailed in Appendix 2

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 phase-contrast MRI85 for validation purposes.

The multi-compartment formulation offers an alternative method to compute volumetric flow rate through the brain. Introducing perfusion , and integrating it over the brain volume leads toAccording to mass conservation so the imbalance between the two formulations is present because the utilised FE method is not conservative. The present (pressure) and (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 FE formulation48,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 and pressure elements, the natural velocity pair is and 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 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.

Table 3

Results corresponding to the grid sensitivity analysis of the brain simulations.

				Healthy		RMCA occl.
Mesh	n_ele	p_i FE	v_i FE	Q^Ω	Q^Γ	Q^Ω	Q^Γ
Baseline	1,042,301	P_1	dP_0	604	464	493	380
		P_1	P_1	604	497	493	406
		P2⋆	P1⋆	600^⋆	559^⋆	478^⋆	447^⋆
		P_2	P_2	600	578	478	460
Loc. ref.	3,400,570	P_1	dP_0	603	539	489	433
		P_1	P_1	603	546	489	439
		P_2	P_1	600	582	476	461
		P_2	P_2	600	597	476	471
Unif. ref.	8,338,408	P_1	dP_0	602	511	486	411
		P_1	P_1	602	527	486	425

Due to the computational cost associated with the locally refined (loc. ref.) and the uniformly refined (unif. ref.) meshes, these simulations are run using 12 threads of an Intel Xeon E5-2640 processor and 128 GB RAM. Computations on the unif. ref. mesh with elements require more than 128 GB RAM and therefore they are omitted

highlights the case used to evaluate the impact of pressure and velocity inlet BCs, whereas the bold text indicates the case with the smallest difference between and

Figures 6a and 6b show the imbalance between and highlighting that without improved approximation the surface flux is always underestimated by . Both refinement and increased FE order can mitigate the error of . However, higher order pressure and velocity approximations are required to keep %. 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 elements) can provide a reasonable compromise between accuracy and performance if computing is important. Otherwise, perfusion and pressure distributions can be estimated reasonably well on coarse meshes with about 1 million elements.

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

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.

Table 4

Comparison of some model parameters and results with literature data: volumetric flow rates of the anterior (ACA), middle (MCA) and posterior (PCA) cerebral arteries obtained by surface integration of the velocity vector over perfusion territories; cortical surface area () and brain volume ().

	Simulation	Reference
Q ACA/MCA/PCA (mL/min)	66/126/67	75 ± 15/131 ± 23/51 ± 1085
Q total (mL/min)	600	657 ± 9485; 484 ± 8883
A_cortical(cm^2)	1005	177082
V_b(mL)	1390	869 ± 9937

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 BC24–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.

(ml/min) is preserved, indicating that the simulations with pressure and velocity inlets are equivalent regarding average brain perfusion. 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 () and perfusion () 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 % because and ) are close to 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 and in Fig. 7).

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 ; average pressure in the arteriole (), capillary (), and venule () compartments over the entire brain; minimum (min) and maximum (max) arteriole pressure () and perfusion (F); average perfusion in grey () and white () 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.

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.

Figure 8

Pressure (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 pressure and velocity elements (details in Table 3).

The volumetric flow rates corresponding to boundary regions of the major cerebral arteries are summarised in Table 4. Based on previous clinical 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 collaterals8,28,42,43,46,69 as well as pial arterioles6,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 arteries54 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 “One-dimensional 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, and 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 .) Thereafter, the one-dimensional model has only two free parameters: the lengths of the grey () and the white () matter domains defined so that . The equivalent and (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 , , and (ml/min/100 mL) which overlap exactly with volume-averaged perfusion values of the three-dimensional brain model.

Figure 9 displays the analytical pressure distributions in the arteriole, capillary and venule compartments. A one-dimensional FE approximation with 100 elements is also shown which is in excellent agreement with the analytical results. The average pressure values provided by the 1D model are within % of those computed using the 3D brain in every case. One-dimensional simulations are approximately 3000 times faster compared to their 3D counterparts.

Figure 9

Analytical and numerical solutions of the one-dimensional problem representing a brain tissue column perpendicular to the cortical surface: arteriole (a), capillary (b) and venule (c) compartments. Results in both grey (G) and white (W) matter are displayed.

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

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.

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

Limitations

This subsection summarises some of the most important overlooked factors in the present study. Investigations are limited thus far to a single patient-specific 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 two31 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:The results lay down the foundations for validation, uncertainty quantification, and for the development of a multiscale model which incorporates large arteries as 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.

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.