Fluid–structure interaction for highly complex, statistically defined, biological media: Homogenisation and a 3D multi-compartmental poroelastic model for brain biomechanics

human brain (specifically, we target the hippocampus), to exemplify the qualities and efficacy of this methodology during the course of Alzheimer’s Disease. The methodology we present has been implemented through the Finite Element Method, in a general manner, allowing for the co-existence of an arbitrary number of compartments. For the applications used in this paper to exemplify the method, a four-compartment implementation is used. A unified pipeline is used on a cohort of 35 subjects to provide statistically meaningful insight into the underlying mechanisms of the neurovascular unit (NVU) in the hippocampus, and to ascertain whether physical activity would have an influence in both swelling and drainage by taking into account both the scaled strain field and the proportion of perfused blood injected into the brain tissue. A key result garnered from his study is the statistically significant differences in right hemisphere hippocampal NVU swelling between males in the control group and females with mild cognitive impairment during high and low activity states.


a b s t r a c t
Numerous problems of relevance in physiology and biomechanics, have at their core, the presence of a deformable solid matrix which experiences flow-induced strain. Often, this fluid-structure interaction (FSI) is directed the opposite way, i.e. it is solid deformation that creates flow, with the heart being the most prominent example. In many cases, this interaction of fluid and solid is genuinely bidirectional and strongly coupled, with solid deformation inducing flow and fluid pressure deforming the solid. Although an FSI problem, numerous cases in biomechanics are not tractable via the traditional FSI methodologies: in the internal flows that are of interest to use, the number and range of fluid passages is so vast that the direct approach of a deterministically defined boundary between fluid and solid is impossible to apply. In these cases, homogenisation and statistical treatment of the material-fluid system is possibly the only way forward. Such homogenisation, quite common to flow-only systems through porous media considerations, is also possible for FSI systems, where the loading is effectively internal to the material. A prominent technique of this type is that of poroelasticity. In this paper, we discuss a class of poroelastic theory techniques that allow for the co-existence of a multitude of -always statistically treated -channels and passages of widely different properties: termed multiple-network poroelasticity (or multicompartmental poroelasticity). This paradigm is particularly suitable for the study of living tissue, that is invariably permeated -perfused -by fluids, often different in nature and across a wide range of scales. Multicompartmental poroelasticity is capable of accounting for bidirectional coupling between the fluids and the solid matrix and allows us to track transport of a multitude of substances together with the deformation of the solid material that this transport gives rise to or is caused by, or both. For the purposes of demonstration, we utilise a complex and physiologically very important system, the human brain (specifically, we target the hippocampus), to exemplify the qualities and efficacy of this methodology during the course of Alzheimer's Disease. The methodology we present has been implemented through the Finite Element Method, in a general manner, allowing for the co-existence of an arbitrary number of compartments. For

Introduction
Fluid-Structure Interaction (FSI) methods allow for the coupling of solid and fluid mechanics phenomena, and enable the study of processes that involve the exchange of loads between the fluid and solid (Dowell and Hall, 2001). An example of an internal flow where such interactions are at play is shown in Fig. 1: In the first subfigure, Fig. 1a, a deformable solid is shown, which is permeated by a series of channels, of arbitrary shape and with non-trivial interconnections. Traditional FSI methods can capture the flow field in the channels/passages and the deformation field of the solid matrix, Fig. 1b, since it is straightforward and computationally feasible to mesh the two domains (in a wide variety of ways) and solve the coupled problem.
The challenge necessitating a different viewpoint is illustrated in Fig. 1c which of course can be further complicated almost ad infinitum, to the limit of real biological materials, as depicted in Fig. 1d. For the latter figure, it becomes extremely difficult, and in most cases computationally intractable, to discretise the multitude of fluid channels, as well as the fine solid ligaments and strands that define them; it is clear that an alternative is needed. An approach that is often used, originally in geotechnical engineering and groundwater studies and as of recently in increasing frequency within the field of biomedical engineering, is that of poroelasticity.

Poroelasticity
Poroelastic systems describe fluid flow through a porous medium coupled with deformation of the solid matrix. The fundamental poroelastic system consist of conservation of mass and momentum equations, which are derived at the macroscopic scale (Biot, 1941(Biot, , 1955Terzaghi, 1925b), effectively assuming a statistical representation of the fluid passages and the fluid/solid interface; a homogenisation approach. Building on the principles of linear elasticity, the conservation of momentum in a poroelastic system also includes a fluid pressure term (as a measure of the effect of the fluid in the medium). It is also worth noting that deformation of the medium is usually much slower than the flow rate, and therefore inertial terms are ignored in the formulation, ultimately incorporating a quasistatic assumption. The momentum equation is derived by considering the total stress,σ, and a body force, f, leading to: The total stress (see Table 1) tensor accounts for the fluid pressure in addition to the material stress tensor (or effective stress) as derived from linear elasticity.
The conservation of mass is derived from the fluid content of the medium, ζ, the volumetric fluid flux, v f , and any additional volumetric fluid source/sink terms. The final form of this conservation equation is of the form: The model is closed by assuming the constitutive relationships summarised in Table 1 (Phillips and Wheeler, 2008): These relationships relate the total stress, volumetric fluid flux and fluid content to the primary variables of the twofield formulation, namely the solid matrix displacement (u) and scalar pore pressure (p). For the fluid content relationship, the c j p term represents the amount of fluid that can be injected into a fixed material volume, whilst the α∇ · u term is a measure of the amount of fluid that can be squeezed out of the same volume. (a) Porous matrix of a deformable solid, when no flow is present, at equilibrium (b) when flow is present, a pressure field develops which loads and deforms the solid matrix, (c) direct FSI through detailed deterministic meshing captures flow and deformation fields, (d) highly complex network of fluid passages render direct FSI methods through discretisation of the exact deterministic interface and passages computationally impractical or impossible.

Table 1
Summary of constitutive relations.

Governing equations and modelling
Building on the principles of poroelasticity, the standard mathematical model for diffusive flow in an elastic porous medium is the diffusion-deformation model of poroelasticity proposed by Biot (1941). This is based on the coupling between the pore-fluid potential and the solid-stress fields. An extension of Barenblatt's double-diffusion approach (Barenblatt et al., 1960) and Biot's diffusion-deformation theory leads to the Barenblatt-Biot poroelastic model representing multiple-network diffusion in elastic porous media. This model takes the following form: In the above formulation, for a given number of networks, A ∈ N, the displacement of the solid skeleton is given by u = u(x, t), whilst the fluid potentials for each respective compartment is given by , and t ∈ [0, T ]. The Biot-Willis coefficient (which traditionally couples the momentum and mass conservation equations, see Table 1), α j ∈ (0, 1], for each compartment satisfies φ ≤ α j ≤ 1, where φ is the total porosity., c j ≥ 0 is the constrained specific storage coefficient, K j is the symmetric and uniformly positive definite hydraulic permeability tensor defined by K j = κ j ( µ j ) −1 > 0 (the ratio of the compartmental permeability tensor, κ, to fluid viscosity), ξ j→i is the intercompartmental transfer coefficient, f = f (x, t) represents a body force and h = h j (x, t) (a) Determinants of optimal brain health. The blue arrow denotes factors that can promote optimal brain health, whilst the red arrows indicate factors that hinder brain health. Neurodevelopmental, genetic and environmental factors can therefore, promote or hinder brain health. (b) An example of cerebroventricular dilatation from the cohort used in this study. The horizontal section of a cognitively healthy control subject (female, 68 years old) can be seen to possess smaller cerebral ventricles than the mild cognitively impaired subject (female, 62 years old). Similarly, enlarged cerebroventricular representations exist for hydrocephalic patients, as ventriculomegaly is an overlapping characteristic of this disorder. represent any additional compartment specific source/sink terms. In this manuscript, ε (u) denotes the small-strain tensor derived from the symmetric part of the gradient of the solid matrix displacement u: The elastic stiffness tensor, C, defines a stress tensor σ using Hooke's Law: We assume an isotropic and homogeneous linear elastic medium with elasticity tensor, C, defined by the identity: where I is the identity tensor, and µ and λ are the Lamé moduli.
It is useful at this stage to describe a field of application of the above class of methodologies, namely biomechanics and in particular the biomechanics of brain diseases. Brain disorders such as developmental and neurodegenerative diseases represent an enormous healthcare burden, not only in terms of human distress, but also economic cost (in Europe, this figure approaches 1 trillion euros (Di Luca et al., 2018)). Cognition is a group of mental processes which include memory, attention, learning, decision making, problem solving, language processing and executive functions. Dementia is classified as a progressive cognitive disorder which primarily affects memory, but additional symptoms also include aphasia, apraxia, agnosia, and lifestyle impairments (Cermakova et al., 2015). The most common form of dementia is Alzheimer's disease (AD), and the major risk factor for its development is increasing age (Wallin et al., 2013), in addition to other known factors, such as hypertension and hypotension, heart failure, low levels of physical activity and of education, obesity, and genetic factors (Cermakova et al., 2015;Kivipelto et al., 2001Kivipelto et al., , 2005Huang et al., 2004;Qiu et al., 2006).
Interestingly, there is a growing body of evidence that suggests a close association between cardiovascular risk and both cognitive impairment and dementia (Gorelick et al., 2017). Furthermore, it is also postulated that cardiovascular risks are associated with modifiable risk factors (such as physical inactivity, depression, hypertension, and smoking), so it may be possible to improve brain health and to delay the onset of dementia in later life (see Fig. 2a). Delaying the onset of AD by just 1 year would lead to an estimated 9 million fewer cases by 2050 (Brookmeyer et al., 2007).

Hydrocephalus
Although a precise definition is controversial, hydrocephalus (HCP) can be succinctly described as the abnormal accumulation (imbalance between production and circulation) of cerebrospinal fluid (CSF) within the brain (Tully and Ventikos, 2011;Rekate, 2009;Thompson, 2009;Stagno et al., 2012;Vardakis et al., 2013b). This balance of CSF production and reabsorption normally allows for the maintenance of the CSF pressure to lie within a tight range. HCP is classified with regards to whether the point of CSF obstruction or discrete lesion lies within the ventricular system (obstructive) and obstructs the flow before it enters the subarachnoid space (SAS) (Corns and Martin, 2012), or not (communicating). There is currently no definitive cure for this disorder. Dilatation of the cerebroventricular system (see Fig. 2b) can lead to loss of brain cells that ultimately results in a variety of neurological symptoms (such as AD, as described in the next section), stroke, and sometimes even death due to the pressure applied on the brain parenchyma (Kartal and Algin, 2014).
1.4. Alzheimer's disease, mild cognitive impairment and the hippocampus AD can be deemed as a heterogeneous mixture of multiple age-related neurodegenerative factors and vascular related pathologies. The hallmark pathological features of the disease are the extracellular deposition of amyloid-β (Aβ) peptide into parenchymal senile plaques or within the walls of arteries and capillaries, in addition to the aggregation of hyperphosphorylated tau into intracellular neurofibrillary tangles and neuropil threads (Tarasoff-Conway et al., 2015;Ramanathan et al., 2015). Evidence also suggests that AD may be a vascular disorder (Gottesman et al., 2017;Roher et al., 2012), caused by impaired cerebral perfusion (characterised by the reduction in both total and regional CBF (Thomas et al., 2015;Roher et al., 2012;de la Torre, 2004)), which is observable at the early stages of the disease (Drachman, 2014;Austin et al., 2011;Thomas et al., 2015). This reduction in perfusion can compromise the oxygenation of neurons, negatively affect the synthesis of proteins required for memory and learning, and subsequently lead to neuronal dysfunction or death (Mazza et al., 2011;Ruitenberg et al., 2005). Ultimately, the final consideration that needs to be made is the clearance of Aβ at the level of the blood-brain barrier (BBB). The BBB forms an important part of the neurovascular unit (NVU), a functional cellular structure that allows for the highly efficient regulation of CBF. The BBB consists of endothelial cells connected by tight junctions and a thick basement membrane which is supported by astrocytic end feet (see Fig. 3). Communication between the cells of the NVU is required to enable efficient clearing of Aβ to prevent it from accumulating in the form of plaques (Iadecola, 2004). Breakdown of the BBB ultimately results in the impaired clearance of Aβ, leading to amyloid accumulation in the brain parenchyma and in and around capillaries. The latter process is known as cerebral amyloid angiopathy (CAA), and it is defined as a major pathological insult to the NVU (Thal et al., 2008). CAA is associated with cognitive impairment (Greenberg et al., 2004), is accelerated by hypoperfusion, and is present in over 80% of patients with AD (Okamoto et al., 2012).
One of the overarching foci of neuropsychology (Bondi et al., 2017) since the turn of this century has been to better understand the prodromal stages of AD. During this early stage, AD may present itself as mild cognitive impairment (MCI), an intermediate state between normal ageing and dementia. Traditionally, MCI has been defined as a condition whereby an individual experiences memory loss to a greater extent than that expected for that age, but does not meet the criteria for dementia (Bondi et al., 2017).
The hippocampus is small structure in the brain, and it plays an important role in spatial and episodic memory. It is the region of the brain that tends to show the most rapid loss of tissue earliest in the disease course. Reduced hippocampal volume results in an amnestic syndrome, a core feature of Alzheimer's disease (Halliday, 2017).

Modelling Alzheimer's disease using poroelasticity
In the literature, several works have utilised a poroelastic approach in modelling parenchymal tissue within the realm of hydrocephalus and oedema formation in the brain and small intestine (Tully and Ventikos, 2011;Vardakis et al., 2013b;Guo et al., 2018;Vardakis et al., 2017Vardakis et al., , 2016Chou et al., 2016;Chou, 2016;Vardakis et al., 2013a;Kaczmarek et al., 1997;Levine, 1999Levine, , 2000Smillie et al., 2005;Sobey and Wirth, 2006;Shahim et al., 2012;Aldea et al., 2019;Thompson et al., 2019). The poroelastic modelling of parenchymal tissue for the purpose of investigating AD yields a narrower selection of relevant work (Guo et al., 2018;Thompson et al., 2019). Recently, Aldea and colleagues (Aldea et al., 2019) utilise poroelastic theory in a multiscale model of arteries in order to test the hypothesis that cerebrovascular smooth muscle cells drive intramural periarterial drainage. Recently, Guo and colleagues (Guo et al., 2018) introduce a pipeline that intertwines a general 3D multiple-network poroelastic model of perfused parenchymal tissue, an image-based modelling pipeline and a detailed subject-specific boundary condition model that can be used to model the influence of lifestyle and environmental factors in obtaining novel biomarkers during the MCI stage of AD.

Outline of the article
This work uses a novel consolidated pipeline that integrates three important components: a three-dimensional multiple-network poroelastic theory (MPET)-based model of perfused cerebral parenchyma; an accurate, fully automated image-based model personalisation workflow (Guo et al., 2018); and a subject-specific boundary condition model that targets the driving compartment of the MPET model. Specifically, MPET allows for the detailed investigation of spatiotemporal transport of fluid between the cerebral blood (arteries, capillaries and veins), CSF and interstitial fluid (CSF/ISF) and brain parenchyma across multiple scales. This pipeline is used on a cohort of subjects (both cognitively healthy controls (CHCs) and mild cognitively impaired (MCI) subjects stratified with respect to gender) in order to assess two novel biomarkers (swelling and drainage, which is derived from the fluid content of the capillary compartment of the four-network MPET model described in Section 2) during the early stages of AD. This is conducted by providing insight into the underlying mechanisms of the neurovascular unit in the hippocampus for both controls and cognitively impaired subjects during two states of activity (high and low) within a 24-hour period. The essential breakdown of the methodology behind the full implementation scheme follows in Section 2, which highlights the consolidated pipeline embedded within the VPH-DARE@IT research platform, the prospective data collection programme used to extract the subject-specific data (including boundary conditions) for the 35 subjects used in this study, and an outline of the statistical analyses used to analyse the MPET results (swelling and drainage). The results are given in Section 3, where two MPET simulations are presented (based on one CHC and one MCI subject) in order to depict the nature of the solution fields that are obtainable at the level of the parenchyma. Subsequently, a three-way mixed ANOVA was conducted in order to understand the effects of gender, cognitive status and activity on blood flow rate in the left and right ICA and VA (as the blood flow rate was a key driver in the MPET model), followed by a Kruskal-Wallis H-test (to determine if there were differences in NVU swelling and drainage in the hippocampus between the groups considered, during two levels of activity) and Wilcoxon signed-rank test (to determine whether there was no statistically significant median decrease in NVU swelling and drainage in the hippocampus when subjects lowered their activity level). The results are discussed in Section 4, along with limitations and perspectives for future work. The conclusions to the paper are given in Section 5.

Methodology
Modelling the transport of fluid within the brain, in a personalised manner and from first principles, is essential to help decipher some of the underlying mechanisms that are currently being investigated regarding diseases of the cerebral environment, such as hydrocephalus and AD. An MPET model for perfused parenchymal tissue is coupled with an automated image-based model personalisation workflow (Guo et al., 2018), and a subject-specific blood flow variability model. The consolidated pipeline will then be used on a small cohort of 35 subjects and used to provide insight into the underlying mechanisms of the NVU in the hippocampus.
In the previous section, we discussed the basic governing equations of poroelasticity and how they can be cast for multiple compartments. In this section we shall briefly present a formulation where four compartments are used, and we describe a computational framework where these equations are solved using the Finite Element Method.

Three-dimensional MPET model for the cerebral environment
In this paper, a MPET model was used to conduct mechanistic modelling of fluid transport through the brain parenchyma. Biologically, the solid matrix represents brain parenchyma, and the communicating fluid phases considered are: an arterial network (a), an arteriole/capillary network (c), a CSF/ISF network (e) and a venous network (v). This model allows for the simultaneous solutions of continuity and momentum conservation equations, in four interconnected fluid compartments, within a deformable solid matrix (the parenchymal tissue) (see Fig. 3).
The MPET model uses the parenchymal tissue displacement (u), and the pore pressures of the four fluid compartments (p a , p c , p e , p v ) as the primitive variables in the governing equations, which are given below: The S terms in Eq. (7b-e) define spatially varying source (S ij > 0) or sink (S ij < 0) densities (rate of fluid transfer between networks). More details can be found in Guo et al. (2018). As described in Section 1, K is the hydraulic permeability tensor for each fluid compartment. It is defined as κ(µ −1 ), where κ is the permeability tensor for each of the four fluid networks; whilst µ defines the viscosity of each fluid. In this work, three of the four fluid domains are isotropic (a, c, v), which implies κ=κI, where κ is a constant and I is the unit tensor for an isotropic medium. For the CSF/ISF compartment (e), a spatially varying permeability tensor extracted from diffusion-weighted imaging (DWI) was used, as κ can be defined in a heterogeneous and anisotropic manner (Guo et al., 2018).

Numerical implementation of the MPET solver, verification and mesh independence
The highly coupled governing equations of the MPET system have been discretised using the Finite Element Method, and implemented into an in-house FORTRAN code. Both the equilibrium equation (the displacement field u) and mass conservation equations (scalar pressure p i ), is approximated in the continuous piecewise linear polynomial space. Based on the principle of minimum potential energy, the algebraic form of the equilibrium equation is: K s is the stiffness matrix, B the deformation matrix and D is the elasticity matrix. Q i is the load on the solid phase contributed from the ith fluid network and h is a mapping vector. F is the load vector, b is the vector of body forces in the three-dimensional domain , and t N the vector of external force at the boundary, N . The Dirichlet boundary conditions are imposed in a strong way.
For the continuity equations of the fluid networks, the method of weighted residuals and the continuous Galerkin formulation are applied to derive the integral form (weak form) of these mass conservation equations. The continuity equation of the i th fluid network can be written as, where N is the boundary where the Neumann boundary condition is applied, and q i is the flux prescribed in the Neumann boundary condition. The discretised continuity equation of the i th fluid network is: where, The governing equations of the MPET system have been discretised using the Finite Element Method and implemented into an in-house FORTRAN code, which has been verified (Guo et al., 2018) against Terzaghi's (1925a) and Mandel's (1953) problems. In addition to the verification, mesh independence (12 meshes with total element numbers ranging from ∼100k to ∼9 million) of the 3D MPET outputs (displacement, scalar pressures and relevant filtration velocities from the four compartments) using a subject-specific brain geometry has been confirmed (Guo et al., 2018).

Subject-specific datasets
The subject-specific datasets used in the MPET modelling of this paper were collected as part of the VPH-DARE@IT project (www.vph-dare.eu), and prospective data collection was conducted at the Istituto di Ricovero e Cura a Carattere Scientifico (IRCCS) San Camillo, Lido di Venezia, Italy. This study, including a total of 103 people (n = 50 CHC, age 71.1 ± 7.9 years, and n = 53 with diagnosed MCI, age 75.1 ± 6.7 years), was approved by the joint ethics committee of the Health Authority Venice 12 and the IRCCS San Camillo (Protocol number 2014.08), and all participants gave informed consent prior to participation in the study. A cohort of 35 subjects (n = 20 CHC, n = 15 MCI) was picked from the 103 available subjects. This smaller cohort can be stratified into 4 groups, considering males (M) and females (F): CHC M (n = 8, age 69.4 ± 8.5 years), CHC F (n = 12, age 72.5 ± 5.7 years), MCI M (n = 8, age 75.4 ± 5.0 years) and MCI F (n = 7, age 74.9 ± 8.2 years).
Each subject had several measurement modalities collected as part of the study, such as: lifestyle questionnaires and neuropsychological tests, whole brain MR imaging, clinical ultrasound flow imaging, portable Holter recordings of blood pressure, and actigraph measured activity levels. Further details can be found in Guo et al. (2018). For the Lido study cohort, Holter recordings and ultrasound flow measurements were used to generate boundary conditions of arterial blood flow using cerebral autoregulation models and lumped parameter circulation models (Guo et al., 2018;Lassila et al., 2018), T1-weighted and diffusion-weighted MR images were processed to create accurate 3D whole-brain meshes and finally permeability tensor maps of the parenchyma were extracted using the workflow described in detail in Guo et al. (2018).

Subject-specific boundary conditions and parameters
In Fig. 4, the subject-specific modelling pipeline for acquiring personalised cerebral blood flow waveforms  that are fed into the arterial compartment of the MPET model (Guo et al., 2018) is depicted. For each subject, four waveforms were calculated at every time point, which are the ICA blood to the left and right cerebrum (ICA L and ICA R ), and the vertebral artery (VA) blood to the left and right cerebellum (VA L and VA R ). In order to apply these subjectspecific waveforms as boundary conditions for the arterial compartment in Eq. (14b), the cortical surface is divided into four perfusion regions corresponding to the four waveforms (Guo et al., 2018). Furthermore, the blood flow waveform is of 1 s duration depicting a period of the maximum (high activity) and minimum (low activity) activity during a 24-hour period. Once the Neumann boundary conditions are applied at the partitioned cortical surface, the numerical simulations are executed for 50 cycles (of cerebral blood flow waveforms), for the solution fields to reach a periodic steady state. This approach is adopted in the work of Guo and colleagues (Guo et al., 2018). The MPET solutions from the final steady-state are used to conduct the statistical analysis.
The MPET system is completed with the following boundary conditions for each of the four compartments (a, e, c, v), where ∂ s and ∂ v are boundaries at the skull and cerebral ventricles respectively, and n is the outward unit normal vector. On the cortical surface: .nds on ∂ S .
(14a-f) On the ventricular wall: Since this is an adult brain that is being taken under consideration; a rigid wall approximation can be envisaged, stemming from the elimination of layers like the dura mater and scalp. A subject-specific blood flow profile (SSBFP) depicting a period of high and low activity is used as the BC for the arterial network at the cortical surface.  CSF is assumed to be produced at a constant rate, Q p within the ventricles. p bp is the blood pressure in the sagittal sinus, κ c→vent represents the capillary network resistance to the flow from the capillary network, R is the resistance due to the presence of arachnoid granulations and finally Q o is the efflux of CSF at the region of the skull. r v represents the radius of the spherical shell encapsulating the cerebroventricular system and u 1 is the maximum ventricular displacement at each time increment. It is necessary to restrict the transfer of water between fluid networks, and Fig. 2 depicts a schematic of the setting in which the quadruple MPET model functions. Table 2 gives the complete list of all parameters used to execute the MPET model.

Statistical analysis
A Shapiro-Wilk's test (p > 0.05 for all cases) is used to assess the normality of the data. Outliers in the ICA and VA based data were assessed by inspection of a boxplot. Levene's test was used to assess for equality of variances, whilst Mauchly's test of sphericity (Mauchly, 1940) was also used to test the null hypothesis that the variances of the differences of the groups (ICA R (HA), ICA L (HA), ICA R (LA), ICA L (LA), VA(HA), VA(LA)) considered are equal. A three-way mixed ANOVA was run to understand the effects of gender, cognitive status and activity on blood flow rate in the left and right ICA and VA. All pairwise comparisons were performed for statistically significant simple main effects.
A Kruskal-Wallis H-test (Kruskal and Wallis, 1952) was conducted to determine if there were differences in swelling and drainage in the hippocampus of the brain between 4 groups (CHC M (n = 8), CHC F (n = 12), MCI M (n = 8) and MCI F (n = 7)), during two levels of activity (high and low). Visual inspection of the relevant boxplots was used to determine whether the distributions of swelling and drainage were similar for all groups. This is needed, since if the distributions are not similar, one cannot make inferences about differences in medians between groups. Instead, mean ranks are then used for the analysis. Subsequently (where a statistically significant Kruskal-Wallis H-test exists), pairwise comparisons were performed using Dunn's (1964) procedure with a Bonferroni correction (Armstrong, 2014) for multiple comparisons. Adjusted p-values are presented. A Wilcoxon signed-rank test (Gibbons and Chakraborti, 2011) was used to determine whether there was no statistically significant median decrease in swelling and drainage in the hippocampus when CHC and MCI subjects lowered their activity level. It was assessed whether the differences in swelling and drainage were symmetrically distributed (via a histogram). A Wilcoxon signed-rank test requires the distribution of the differences between the two related groups to be symmetrical in shape. Where this is not the case, a sign test with continuity correction is used, as this test does not make any distributional assumptions.
The statistical analyses were performed using IBM SPSS Statistics, Version 25 (IBM Corp., Armonk, NY).

Results
Fig. 5 depicts a typical set of results (for one 69-year-old female CHC and one 75-year-old male MCI subject) obtained from executing the consolidated pipeline. In the figure, coronal sections of the brain at the level of the hippocampus are depicted. The results pertaining to the NVU fluid content, clearance and perfusion for the two subjects are shown. The intracranial pressure (pore pressure of the CSF/ISF compartment) is also shown, and this solution field was made to overlap with the filtration velocity vectors of the capillary compartment and CSF/ISF compartment respectively. There is a reduction in peak clearance observed between the two subjects, with the female CHC subject having a peak clearance of 24 µm/s compared to 31 µm/s for the male MCI subject. Similar characteristics were observed for swelling and drainage between the two subjects, with the female CHC possessing lower peak values (1.91 and −0.13) compared to the male MCI subject (1.92 and −0.49). The level of peak perfusion was reversed, with the female CHC possessing higher level of perfusion (0.26 mm/s) compared to the male MCI subject (0.14 mm/s).
Blood flow rates in the ICA L/R and VA for males and females (both CHC and MCI) during both high and low activity were normally distributed, as assessed by Shapiro-Wilk's test (p > 0.05 for all cases) (Shapiro and Wilk, 1965). There  Mauchly's test of sphericity is not violated as there are only two levels of the within-subjects factor (high and low activity) and, therefore, there is only one paired difference. For the right ICA, the three-way interaction between gender, cognitive status and activity was not statistically significant, F(1, 31) = 0.147, p = 0.704, partial η 2 = 0.005. Similarly, for the left ICA this was F(1, 31) = 0.391, p = 0.536, partial η 2 = 0.012, and for the VA this was F(1, 31) = 0.267, p = 0.609, partial η 2 = 0.009. All two-way interactions were not statistically significant (p > 0.05). For the four groups considered, the mean CBF and standard deviation of the mean are listed in the following table and approximated to 1 decimal place.
As can be seen from Table 3, there peak flow rates that are calculated from the personalised cerebral blood flow waveform pipeline are lower during low activity for the left and right ICA, and VA. It can be observed that grouped (both males and females) CHC subjects possess higher flow rates than the grouped MCI subjects. The stratified results (with respect to gender) indicate that female CHC subjects possessed on average higher flow rates for the left ICA, right ICA, and VA compared to the male CHC cohort, whilst for the MCI subjects, this trend was reversed (MCI males possessed higher flow rates). For all gender specific groups, the flow rates were higher during high activity compared to those recorded during low activity. The parenchymal brain volume (including the cerebral ventricles) for each group is also given in the last column of Table 3, for reference. A Kruskal-Wallis H-test was conducted to determine if there were differences in swelling and drainage in the hippocampus of the brain between the four (CHC M , CHC F , MCI M and MCI F ) groups, during two levels of activity (high and low). Distributions (see Fig. 6) of swelling were not similar for all groups, whilst the distributions of drainage were. Subsequently, a Wilcoxon Signed Rank Test determined whether there was any statistically significant median decrease in swelling and drainage in the hippocampus when subjects of the four groups lowered their activity level. The results for the Kruskal-Wallis H-test is given in Table 4.
In the right hippocampus, the level of NVU swelling was statistically significantly different between the categories of cognitive status and gender during the period of high and low activity, χ 2 (3) = 8.929, p = 0.030 (for both activity states). Subsequently, pairwise comparisons were performed using Dunn's (1964) procedure with a Bonferroni correction for multiple comparisons (Armstrong, 2014). The post hoc analysis revealed statistically significant differences in NVU Table 3 Mean CBF and parenchymal tissue volume (with the standard deviation of the mean, SDµ) for the four groups in addition to the overall CHC and MCI cohorts. High activity is indicated by red cells, whilst low activity by the blue cells. In the left hemisphere, there were differences in the degree of swelling and drainage between the four groups, but these were not statistically significant.
When conducting a Wilcoxon Signed Rank Test, it was observed that when comparing the four groups during high and low activity, there were differences in the degree of swelling and drainage, but that these differences were not statistically significant.

Discussion
In the work presented here, various advancements have been made in refining the integrity of the qualitative representation of fluid accumulation and drainage, and applying it to the level of the BBB. Firstly, the consolidated pipeline has been applied to a cohort of 35 subjects (20 CHC, 15 MCI), and MPET results are statistically analysed at a region-specific level, namely the hippocampus. Secondly, it was hypothesised that physical activity would have an influence in both swelling and drainage.
In Table 1 (specifically, the fluid content relationship), the c c p term represents the amount of blood (which is also exchanging fluid with the cerebrospinal/interstitial fluid -see Fig. 3) that can be injected into a fixed material volume (in this case, the parenchymal tissue), whilst the α c ∇ · u term (the Biot-Willis constant, α, is deemed to reflect mechanical compliance (Tully and Ventikos, 2011)) is a measure of the amount of blood that can be squeezed out of the same volume of parenchymal tissue. Importantly, this second term reveals the importance of the scaled (with respect to the Biot-Willis coefficient) strain field within the parenchymal tissue, and the role this plays in determining whether the overall effect of the fluid content is either positive or negative, as it is possesses a larger magnitude in the simulations conducted in this work. Of course, this is also dependent on the boundary conditions attributed to pore pressure in the capillary compartment, as in the simulations conducted here, the range in pore pressure was approximately 2.4-2.5 kPa. Enforcing a large pressure on the boundaries would translate into a larger quantity of fluid that can be injected into the fixed parenchymal volume, which would therefore dampen the influence of the scaled strain, and in the process determining the propensity of the fluid content. It is worth noting that a balance should be struck when determining the optimal course of action here, as being able to account for the influence of the surface concavity (and therefore the accumulation of strain) is beneficial in assessing the macroscopic effects of not only the geometry under consideration (as in the case of the cerebral ventricles), but also the cumulative effect of the imposed boundary conditions on the remaining compartments and the intercompartmental fluxes.
In Fig. 6, the effect of incorporating physical activity are negligible. This result at first sight seems unexpected, as the variation in CBF between activity states (see Table 3) that is applied as the arterial boundary conditions to the MPET model (which essentially drives the MPET system) shows a clear reduction during the period of low activity. The value of the constrained specific storage coefficient for the capillary compartment underestimates the contribution of the pore pressure in this compartment and allows the strain field to play a dominant role in determining the resultant fluid content. It is important to note here that the constant strain cross-storage effects relating to the other three compartments are also not incorporated in this formulation, as the usefulness of their inclusion in our model is debatable (Vardakis et al., 2017). An additional limitation is that the current quadruple-MPET model conflates all the fluid outside the vascular tree (CSF and ISF) into a single compartment. This negates the important effects of the paravascular/perivascular spaces, glial cells and the overarching influence of the glymphatic system, which are known to play an important role in BBB breakdown and dysfunction (Iliff et al., 2012;Sweeney et al., 2018). The current MPET model has been recently extended to six compartments (in a simplified, spherically symmetric 1D formulation) and has relaxed the assumption of quasi-steady behaviour  in order to account for both aqueductal stenosis during acute hydrocephalus in addition to providing insight into oedema formation (Chou, 2016). In future work, the 3D extension to this discretisation template will be implemented within the consolidated pipeline described in this work.
From Fig. 6, both swelling and drainage is asymmetric in nature, as the distribution of the results varies between the left and right portion of the hippocampus. The mean swelling in the NVU is higher in the MCI subjects (also considering the stratified results with respect to gender), which can be postulated to signify the higher probability of BBB breakdown or dysregulation within the representative elementary volume of parenchymal tissue that the MPET system covers.
The biomarker defined by the fluid content in the capillary compartment can be deemed to qualitatively describe the physiological cascade of events that evolve due to NVU malfunction in a variety of pathological processes in both acute (such as traumatic brain injury) and chronic conditions (such as AD Busch et al., 2012;Sagare et al., 2012), as these are characterised by an inflammatory response, degradation of the extracellular matrix and loss of permeability and selectivity of the BBB (Muoio et al., 2014). More specifically, within the brain tissue, BBB breakdown leads to the perivascular accumulation of blood-derived neurotic products which ultimately leads to neuronal injury, cell death and inflammatory response (de Vries et al., 2012), and albumin which contributes to oedema formation (which can be captured by a positive fluid content in the CSF/ISF compartment of the MPET model), hypoperfusion and tissue hypoxia (Yang and Rosenberg, 2011). Detachment, degeneration and loss of pericytes (see Fig. 3) also leads to BBB breakdown (Montagne et al., 2017). As opposed to BBB breakdown, a dysregulated BBB effects the key mechanisms behind Aβ homeostasis in the brain (Koizumi et al., 2015) and can lead to inefficient drainage along the perivascular route which is signified by Aβ accumulating in the space between the astrocytic end-feet and the blood vessel walls (perivascular space -see Fig. 3) (Iliff et al., 2012). It has been noted in previous work by the same authors that when considering the CSF/ISF compartment, the fluid content was here used to describe periventricular lucency, and in doing so, conceptual links could be made with the swelling activated components of the parenchymal tissue microstructure (Guo et al., 2018;Vardakis et al., 2016).
There is evidence to suggest that a disruption to the BBB compromises Aβ clearance at the level of the NVU, which is postulated to result in a cycle between Aβ accumulation and BBB dysfunction during AD progression (Yamazaki and Kanekiyo, 2017). Aβ has been deemed a likely candidate in the disruption of tight junctions in endothelial cells, which ultimately perturbs the ability of the endothelial cells to function as an effective barrier. In the results presented for NVU swelling in the hippocampus, the mean swelling in the left and right portions of the hippocampus are lower for CHC subjects when compared to MCI subjects (1.94 and 1.86 compared to 1.95 and 1.95). The level of drainage in the NVU is asymmetric between the two sets of subjects, as the left hippocampal portion portrayed more pronounced mean drainage (−0.22 for CHCs compared to −0.42 for MCIs) compared to the right portion (−0.63 for CHCs compared to −0.58 for MCIs).
When stratifying the results with respect to gender, male and female CHC subjects displayed lower degrees of swelling and drainage in both hippocampal portions compared to male and female MCI subjects. Further analysis revealed statistically significant differences in NVU swelling between CHC males (11.19) and MCI females (26.36) (p = 0.025) during both activity states for the right hemisphere.

Conclusions
This paper describes a three-dimensional multicompartmental poroelasticity model for perfused parenchymal tissue coupled with an automated image-based model personalisation workflow, and a subject-specific blood flow variability model. This unified pipeline is used on a cohort of 35 subjects (stratified with respect to gender and disease status) to provide insight (statistically analysis) into the underlying mechanisms of the neurovascular unit (NVU) in the hippocampal region of the brain, and to ascertain whether physical activity would have an influence in both swelling and drainage by taking into account the scaled strain field and the proportion of perfused blood injected into the brain tissue. A key result garnered from his study is the statistically significant differences in hippocampal (the right portion) NVU swelling between CHC males and MCI females during both activity states.
In future work, revised estimates for the constrained specific storage coefficients need to be established, as these allow for the percolation of fluid between networks to be more accurately captured. It has been postulated that systemic hypertension causes stiffening and microvascular distortion of vessels, increased tortuosity and thickened membranes in arterioles, and is associated with a reduction of capillary density (Iadecola, 2013). The CSF/ISF and capillary compartments can be used to link microvascular permeability and overall cerebral compliance, and therefore provide further insight (via rule-based remodelling processes which will also need to account for white matter changes associated with AD Amlien and Fjell, 2014) into the influence of lifestyle and environmental factors of interest, and the extent of chronic cerebral hypoperfusion (promoting the notion of age-dependent vascular compliance Tanaka et al., 2000). Finally, tortuosity may be better served by more intricate models (as opposed to Darcy's law) of flow in porous media (such as the Blake-Kozeny-Carman (BKC) model Sochi, 2010 for a packed bed).