A multiple-network poroelastic model for biological systems and application to subject-speciﬁc modelling of cerebral ﬂuid transport

Biological tissue can be viewed as porous, permeable and deformable media inﬁltrated by ﬂuids, such as blood and interstitial ﬂuid. A ﬁnite element model has been developed based on the multiple-network poroelastic theory to investigate transport phenomenon in such biological systems. The governing equations and boundary conditions are adapted for the cerebral environment as an example. The numerical model is veriﬁed against analytical solutions of classical consolidation problems and validated using experimental data of infusion tests. It is then applied to three-dimensional subject-speciﬁc modelling of brain, in- cluding anatomically realistic geometry, personalised permeability map and arterial blood supply to the brain. Numerical results of smoking and non-smoking subjects show hy- poperfusion in the brains of smoking subjects, which also demonstrate that the numerical model is capable of capturing spatio-temporal ﬂuid transport in biological systems across different scales. the


Introduction
Poroelastic theory is widely used in civil, petroleum and biomedical engineering. Two directions are commonly acknowledged for the development of the poroelastic theory ( Coussy, Dormieux & Detournay, 1998 ;de Boer, 1992 ): one is based on mixture theory while the other one is founded on macroscale theories, mainly represented by the work of Biot. Mixture theory owes much of its current structure to the early works of Truesdell (1957aTruesdell ( , 1957bTruesdell ( , 1962 . Extensive reviews of the literature on mixture theory can be found in the papers by Atkin and Craine (1976a) , Bowen (1976) and Bedford and Drumheller (1983) , and the books by Samohyl (1987) and Rajagopal and Tao (1995) , Truesdell (1984) . Mixture theory assumes that the domain of a mixture can be viewed as a superposition of several single interpenetrating continua, each representing a different constituent and following its own motion; at any time, each position in the mixture domain is occupied simultaneously by one particle from each constituent, in a homogenised sense ( Atkin & Craine, 1976a ). Based on mixture theory, the governing equations for poroelasticity are derived through averaging procedures, generally using a Eulerian description. In contrast to mixture theory, the macroscale theories (e.g. Biot's consolidation theory) assume that the stress and other related variables hold at the macroscale, as well as most of the classical field equations of continuum mechanics, generally using a Lagrangian description ( Coussy et al., 1998 ). It has been demonstrated that the Biot's theory is a special case of the more general formulation derived within the content of mixture theory, when appropriate assumptions are made ( Bowen, 1982 ;Rajagopal, 2007 ;Rajagopal & Tao, 1992. The coefficients used in these two groups of theories also have different physical interpretations ( Schanz & Diebels, 2003 ).
Based on mixture theory, the governing equations for poroelasticity are derived on a strict mathematical ground. For example, the Darcy's equations in a generalised form can be deduced from mixture theory ( Atkin & Craine, 1976b ;Kannan & Rajagopal, 2008 ;Munaf, Lee, Wineman & Rajagopal, 1993 ), whereas they were assumed in Biot's work as an approximation of the balance of linear momentum for the fluid that is flowing through a porous solid ( Rajagopal, Saccomandi & Vergori, 2010 ;Subramanian & Rajagopal, 2007 ). However, the Biot theory -a heuristic model, provides a straightforward continuum theory for materials with microstructures and leads to satisfactory results in describing porous media models ( de Boer et al., 1988 ;Rohan & Lukeš, 2017 ). One of the main advantages is that the constitutive relationships involve well defined and measurable quantities at the macroscale, such as the increment of fluid content, which makes it more tractable for engineering applications ( Coussy et al., 1998 ).
The consolidation theory describing a porous medium at the macroscale was proposed by Biot (1941) and then further developed by Biot (1955Biot ( , 1956) and other researchers ( Prevost, 1980 ;Thigpen & Berryman, 1985 ;Zienkiewicz, Chan, Pastor, Paul & Shiomi, 1990 ;Zienkiewicz & Shiomi, 1984 ). A traditional engineering field to apply the poroelastic theory is soil and rock mechanical problems in geological systems, where it is mainly used to understand subsidence, estimate hydrocarbon volumes, and predict stresses around boreholes ( Wang,20 0 0 ). For complicated engineering applications, the finite element method is often used to solve poroelastic problems. Its application to poroelasticity began in 1969 ( Sandhu & Wilson, 1969 ) on a special case of incompressible fluid and solid constituents in soil. Since then numerous research activities have been devoted in applying various spatial and temporal discretisation schemes and in the investigations of numerical characteristics, e.g. stability, convergence and error estimation, of existing algorithms. Some examples are the work of Vermeer and Verruijt (1981) to investigate the numerical instability at early stages of consolidation, the work of Wheeler and Gai (2007) , which used mixed finite element spaces for pore pressure and continuous Galerkin methods for displacements, and that of Ferronato, Castelletto and Gambolati (2010) , which used the lowest order Raviart-Thomas mixed space for the fluid pore pressure and flux, and linear tetrahedral elements for the displacement.
A more complicated poroelastic theory considers the multi-porous nature of geological systems, i.e. characterised by several distinct families of diffusion or flow paths ( Wilson & Aifantis, 1982 ). To capture diffusion or infiltration processes in such media, the multiple-network poroelastic theory (MPET) was developed from the classical form of Biot's consolidation theory ( Biot, 1941 ;Detournay & Cheng, 1993 ), which describes an isotropic and nearly incompressible solid matrix and homogeneous porous medium. This concept of multiple-porosity was first introduced by Aifantis (1977) and then developed and used for geomaterials ( Aifantis, 1979 ;Aifantis & Hill, 1980 ;Bai, 1999Bai, , 1993Berryman & Wang, 1995 ;Choo & Borja, 2015 ;Kim, Sonnenthal & Rutqvist, 2012 ;Mehrabian & Abousleiman, 2014 ), where fissured rocks or aggregated soils exhibit more than one scale in terms of fluid flow and solid deformation. One example is shown in Fig. 1 , which compares two types of porous media representing naturally fractured reservoirs -the one on the left-hand side has only matrix pores and the other one has both matrix pores and fissures. In the latter case, the connected fissure network provides conduits for fluid flow, which can be more dominant, compared with the generally low permeability of the rock matrix. Therefore, for such a medium with discrete fractions of varying solid compressibilities and permeabilities, a multiple-porosity/multiplepermeability model appears more appropriate to represent the realistic system ( Bai et al., 1993 ).
Porous structures also exist in many biological systems. One example is biological tissue, such as brain, kidney and bone. Regardless of the specific and distinctive nature of the microstructures within the tissue, from a mechanical point of view, matrix (biological tissue) is permeated by i = 1, 2, 3, . . . , n different fluid networks -each having a unique porosity ( φ i ), density ( ρ i ), permeability ( k i ), pressure ( p i ), filtration velocity (Darcy velocity) ( v i ) and increment of fluid content ( ζ i ).
In the cerebral environment, which is the application this paper aims at, the multiple-network poroelastic model can by illustrated by Fig. 2 . The 4-MPET (four-network poroelastic theory) model was proposed to conduct mechanistic modelling of fluid transport through the brain parenchyma ( Tully & Ventikos, 2011 ). Biologically, in a porous medium representing the cerebral environment, the solid matrix represents the brain parenchyma, and the communicating fluid phases considered are: an arterial network ( a ), an arteriole/capillary network ( c ), a cerebrospinal fluid/interstitial fluid (CSF/ISF) network ( e ) and a venous network ( v ) ( Fig. 2 ). The separation of arterial and arteriole networks is based on the consideration of wildly different resistances between large and small arteries. Similar consideration was adopted in the modelling of coronary blood flow in the heart ( Lee & Smith, 2012 ;Smith, Pullan & Hunter, 2002 ), where the arterial tree consists of several compartments. In general, arterioles are defined as the primary resistance vessels that enter an organ to distribute arterial blood into capillary beds, which provides more than 80% of the resistance to blood flow in the body ( Christensen & Mulvany, 2001 ;Martinez-Lemus, 2011 ;Mulvany & Aalkjaer, 1990 ). Therefore the arterial blood is further segmented into a high pressure arterial network and a lower pressure arteriole/capillary network ( Tully & Ventikos, 2011 ).
This model allows for the simultaneous solutions of continuity and momentum conservation equations, in four interconnected fluid compartments, within a deformable solid matrix (parenchymal tissue). In this paper, this 4-MPET model is chosen as the theoretical template for the development of the finite element model.

Governing equations
The formulation of the MPET theory comprises three components: a mechanical equilibrium equation governing elastic deformation; Darcy's equations for modelling fluid flow; and a mass conservation expression. The governing equations use the parenchymal tissue displacement ( u ) and the pore pressures of the fluid networks p i as the primitive variables. The general formulations are listed below.
Eq. (1) is the equilibrium equation, which describes the momentum balance in the porous medium. Here, u is the displacement vector; p i the scalar pore pressure in each fluid compartment ( i = a, c, e, v ); G the shear modulus; λ the Lamé's constant; ε the dilatational strain; α i the Biot-Willis coefficient for each fluid compartment which satisfies φ ≤ α a + α c + α e + α v ≤ 1 ( Berryman, 1992 ; Wang, 20 0 0 ), where φ is the total porosity. Body forces on the bulk material, such as gravity, are neglected in Eq.
(2) is the continuity equation, which describes the mass balance. Here, S i is the specific storage (a measure of the released fluid volume per unit pressure in the control volume at constant strain for each fluid compartment); S ij ( i = j ) is the coupling term describing the effect of excess pore pressure of fluid compartment j on the fluid content increment of another fluid compartment i ; k i is the permeability tensor for each of the four fluid compartments, which reduces to k i = k i I , with k i a constant and I the unit tensor, for an isotropic medium; and μ i is the viscosity of each fluid. The s ji terms on the right-hand side of Eq. (2) define spatially varying source s ji > 0) or sink ( s ji < 0) densities (rate of fluid transfer between networks) ( Tully & Ventikos, 2011 ;Vardakis, Tully & Ventikos, 2013 ).
It should be noted that the cross-porosity storage effect ( S ij ), which means the microscopic coupling effect between the volumetric deformation of different pore systems ( Khalili, 2003 ), is not considered in this paper. In engineering applications, this effect can be considered ( Mehrabian & Abousleiman, 2014 ) or neglected ( Bai, Meng, Elsworth, Abousleiman & Roegiers, 1999 ) based on the specific assumptions adopted. In this paper, the cross-porosity storage effect is not considered mainly due to the lack of experimental data to quantify the coupling terms in a physiological sense ( Vardakis, Chou, Guo, Tully & Ventikos, 2017 ). Therefore, the specific formulations of the governing equations for a 4-MPET model presented in this paper are: The transfer between fluid networks is considered to be driven by a hydrostatic pressure gradient of the form, s ji = ω ji ( p j -p i ), where ω ji is the transfer coefficient scaling the flow from network j to network i . The fluid flow between networks is derived from physiological considerations ( Tully & Ventikos, 2011 ) and required to obey the law of continuity for the entire domain; hence, directionality between fluid compartments must be accurately specified. These are listed as follows: Compared with the general framework of mixture theory, the main assumptions and approximations used in the 4-MPET model are summarised as follows: 1. Body forces on the bulk material are neglected. 2. The inertia terms are neglected in the momentum balance equation. 3. Darcy's equations are used for fluid flow through a porous medium. 4. The pressure gradients are low in biological systems. 5. The fluid viscosity is constant. 6. The perfused tissue undergoes small deformations. 7. The transfer between fluid networks is driven by a hydrostatic pressure gradient.

Boundary conditions
A serious difficulty associated with solving boundary value problems within the context of mixture theory is the specification of appropriate boundary conditions ( Tao & Rajagopal, 1993 ) as the assignment of boundary conditions to overlapping continua is not usually straightforward ( Rajagopal, Yalamanchili & Wineman, 1994 ). For example, if traction boundary conditions are involved, then the difficulty is how to prescribe individual partial tractions when only the total traction is known. A method for overcoming this difficulty for some special mixtures has been provided by Rajagopal, Wineman and Gandhi (1986) . In this paper, the boundary conditions are mainly defined from physiological considerations. The cerebral environment itself is complicated; therefore, assumptions and simplifications are introduced to implement boundary conditions for numerical modelling, which are discussed in detail in the rest of this section. Interpreting the boundary conditions in a physiological sense also makes the partition of boundary conditions between different fluid networks less complicated. Another factor that needs to be taken into account is the availability of data -either in vivo measurements or in vitro experiments, which are necessary to make the modelling results clinically relevant. Based on the considerations mentioned above, a set of boundary conditions applicable for the human brain has been derived for the 4-MPET model.
There are two boundary surfaces to define the volumetric domain that represents the brain parenchyma; the outer boundary is the cortical surface and the inner boundary is the ventricular wall ( Fig. 3 ). Both surfaces are closed, and the cavity enclosed by the ventricular wall lies entirely within the cortical surface, i.e. no overlapping or penetration. The volumetric domain used in modelling is the space between these two surfaces. The boundary conditions for both the solid and fluid phases are listed in Table 1 .  3. The two boundary surfaces defining a typical simulated volumetric domain (brain parenchyma). The outer surface (yellow) represents the cortical surface and the inner surface (turquoise) represents the ventricular wall. It should be noted that it is a cavity inside the ventricular wall, i.e. not included in the modelling. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1
Boundary conditions used for the 4-MPET model.

No displacement constraints
Arterial blood ∇ p a n = Q a ∇ p a n = 0 (14)

Venous blood
The skull is assumed rigid to represent adult cases ( Smillie et al., 2005 ). Previous research showed a compliant behaviour of the subarachnoid space, although the estimated deformations were small due to the large parenchyma surface area ( Gupta, Soellinger, Boesiger, Poulikakos & Kurtcuoglu, 2009, 2010. In this paper, the CSF flow and mechanical response in the subarachnoid space is not considered; for simplicity, the rigid boundary condition is applied directly to the cortical surface, which assumes no displacement is allowed at this boundary ( Eq. (12) ). There are no displacement constraints at the ventricular wall, so it can expand or contract freely. For continuity of stresses, the pressure exerted by the CSF within the ventricles on the inner ependymal surface must balance the poroelastic stress in the parenchymal tissue ( Tully & Ventikos, 2011 ;Vardakis et al., 2013Vardakis et al., , 2016, which is applied as a pressure boundary condition in the numerical modelling. The arterial blood supply to the brain is mainly provided by two pairs of arteries -internal carotid arteries and vertebral arteries ( Tortora & Derrickson, 2009 ). Due to the lack of explicit characterisation of vasculature in the 4-MPET model, the arterial blood supply to the brain is simplified into a flux boundary condition (Neumann boundary condition) Q a at the cortical surface ( Eq. (13) ). More details about assigning the boundary condition for the arterial blood compartment is given in Section 6 .
For the arteriole/capillary blood compartment, the production of CSF from the blood results in a pressure drop in the arteriole/capillary blood ( Eq. (16) ), where κ c → vent is the resistance of the flow from the capillary network to the ventricles (through the choroid plexus), and Q p is the rate of CSF production. Two assumptions are adopted in this boundary condition. First, there is no separation of the two extracellular fluid compartments in the brain -the cerebrospinal fluid (CSF) and the interstitial fluid (ISF) in the 4-MPET model, which assumes that all of the CSF/ISF is produced within the ventricles from blood at a production rate Q p . However, it has been reported that approximately 20% of CSF in the human brain originates from brain ISF ( Edsbagge, Tisell, Jacobsson & Wikkelso, 2004 ;Lei, Han, Yuan, Javeed & Zhao, 2017 ). In the current 4-MPET model, this part of CSF production is implicitly embedded in the combined CSF/ISF compartment. Second, the main site of CSF production in the ventricles is the choroid plexus, which is a highly vascularised tissue located within each ventricle of the brain and develops from several locations along the dorsal axis of the neural tube ( Lun, Monuki & Lehtinen, 2015 ). The classical hypothesis involves the production of CSF at the choroid plexus of the lateral, third and fourth ventricles. However, it is still speculative as to the exact proportions of CSF production in the various choroid plexus sites ( Gupta et al., 2009 ;Vardakis et al., 2013 ). The 4-MPET model simplifies the production of CSF as a uniform distribution on the entire ventricular wall, instead of at specific locations. This simplification is also consistent with the homogenisation assumption adopted for the 4-MPET model. The CSF/ISF compartment has a Dirichlet boundary condition at the cortical surface and a mixed boundary condition at the ventricular wall. At the cortical surface, the boundary condition ( Eq. (17) ) represents the pressure rise resulted from the absorption of CSF into the venous network, where p bp is the venous blood pressure at the cortical surface; μ e is the viscosity of CSF; R is the resistance to outflow through the arachnoid granulations; and Q 0 is the out-flux of CSF at the skull (the rate of absorption Q 0 is assumed to be equal to the production rate Q p in the quasi-steady approach). At the ventricular wall, the boundary condition ( Eq. (18) ) represents the conservation of the mass of fluid in the ventricles. Within the ventricles, it is assumed that any CSF that is produced ( Q p ) and does not flow through the cerebral aqueduct (Poiseuille's law) or the parenchyma must accumulate within the ventricles, where d and L are the diameter and length of the cerebral aqueduct, respectively; r 1 is the distance from the centre to the ventricular wall; and u n 1 is the displacement at the ventricular wall.

Numerical implementation
The highly coupled governing equations of the MPET theory have been discretised using the finite element method and implemented into an in-house numerical code. First, for the equilibrium equation ( Eq. (3) ), the displacement field u is approximated in the continuous piecewise linear polynomial space, where u, v and w are the displacement components in x, y and z -directions, respectively; N i is the continuous piecewise linear polynomial ( i = 1, 2, …, n n ); n n is the total number of nodes in the finite element mesh; and u i , v i and w i are the nodal displacements in x, y and z -directions. Based on the principle of minimum potential energy, the algebraic form of the equilibrium equation is K is the stiffness matrix with B the deformation matrix and D the elasticity matrix. Q i is the load on the solid phase contributed from the i th fluid network and h is a mapping vector. N is the matrix of continuous piecewise linear polynomial functions (shape functions). F is the load vector with b the vector of body force in the three-dimensional domain and t N is the vector of external force at boundary N . The Dirichlet boundary conditions are imposed in a strong way.
For the continuity equations of the fluid networks ( Eqs. (4) - (7) ), the method of weighted residuals and the continuous Galerkin formulation are applied to derive the integral form (weak form) of these mass conservation equations. For example, the continuity equation of the i th fluid network can be written as, where is the three-dimensional domain; 2 is the boundary where Neumann boundary condition is applied, and q i is the flux prescribed in the Neumann boundary condition. The pore pressure of the i th fluid network in the MPET model p i is approximated in the continuous piecewise linear polynomial space, where N j is the continuous piecewise linear polynomial ( j = 1, 2, …, n n ); n n is the total number of nodes in the mesh. Then the discretised continuity equation for one of the four fluid networks is The elements in matrices A, C and vector P are N i and N j are continuous piecewise linear polynomial functions at node i and j in the mesh, respectively. It can be seen from the right-hand side of Eq. (31) that the strain rate calculated from the displacement field acts as an additional source or sink term. The natural boundary condition, i.e. Neumann boundary condition, is applied at boundary 2 . It is known that finite element solutions to consolidation problems may exhibit oscillation in the pressure field and locking when the permeability is very low, or at the very early stage in the consolidation process when small time-steps are used ( Murad & Loula, 1994 ;Phillips & Wheeler, 2009 ;Vermeer & Verruijt, 1981 ;Yi, 2017 ;Zienkiewicz et al., 1990 ). This problem cannot be avoided without smoothing algorithms especially if the pressure and displacement fields are discretised by equal order elements ( Choo & Borja, 2015 ;Reed, 1984 ). In this paper, the finite element model is developed for application in biological systems with effects taking place over long-term periods (such as chronic cases of hydrocephalus), where the focus is on steady-state responses, i.e. the early stage of consolidation is not important for using this model. For example, in the numerical application of Section 6 , the numerical results at the very early stage of the simulation are disregarded and only the final states are considered in the analysis.
The temporal discretisation is conducted by using the method of weighted residuals. If the time-dependent pore pressure p is assumed to vary linearly between time t and t + t , the intermediate value can be estimated by where p n and p n + 1 are the values of p at time t and t + t , respectively. θ is a scalar value between 0 and 1. The continuity equation can therefore be expressed as where the right-hand side can be approximated using the same scheme, In this paper, the time discretisation uses an implicit backward Euler scheme, i.e. θ = 1. The discretised governing equations are solved by the standard KSP linear equation solver in the PETSc library ( Balay et al., 2018a( Balay et al., , 2018b. In the present work, the discretised governing equations are solved sequentially in a tightly coupled manner, i.e. the pressure and displacement solutions are solved sequentially during a time-step until a convergence tolerance of 10 −3 is reached. Previous comparisons showed that the iteratively coupled method and the fully coupled monolithic method gave the same accuracy if sufficiently tight convergence tolerance was used for iterations ( Dean, Gai, Stone & Minkoff, 2006 ;Wheeler & Gai, 2007 ). The inner iterations within a time-step can be illustrated in Fig. 4 . p 1 i is the initialised pore pressure of the i th fluid network at the first time-step and p 0 i is the initial condition of pore pressure. The current time-step is denoted by n + 1 and the previous time-step by n. ξ is an interpolation coefficient used on pressure solutions between inner iterations, which is chosen to be 0.9 in the current numerical code (0 < ξ < 1).
In the finite element formulation of the MPET model, the primitive variables are displacement and pore pressures of the fluid networks, i.e. a p / u formulation. Other variables of practical interests, such as Darcy velocity and increment of fluid content, can be derived from these primitive variables in the post-processing. For example, Darcy's equations are used to calculate Darcy velocities of the four fluid compartments. where where B i is Skempton's coefficient of the i th fluid network and K is the bulk modulus. A positive ζ i corresponds to a 'gain' of fluid by the porous medium. It should be noted that the accuracy of the Darcy velocity and the increment of fluid content depends on the accuracy of the pressure and displacement fields because they are not primitive variables in the formulation. An alternative way is to use a three-field formulation, where the Darcy velocity vector is a primitive variable, e.g. the work of Ferronato et al. (2010) . However, it will increase the computational cost in terms of more degrees of freedom and longer computational time. Therefore, from a practical point of view, it is not necessary to include velocity in the discretised governing equations. Its accuracy can be accepted if the solutions of pressure and displacement are sufficiently accurate.

Model verification
To ensure the formulation is consistent regardless of the number of fluid networks, an n -MPET model should give accurate results if reduced to a single-network poroelastic model ( Berryman, 2002 ). In this section, the 4-MPET finite element model is reduced to one fluid compartment (i.e. n -MPET, n = 1), and verified against analytical solutions of two classical consolidation problems -the Terzaghi's problem ( Terzaghi, 1925 ) and the Mandel's problem ( Mandel, 1953 ).

Terzaghi's problem
The Terzaghi's problem ( Terzaghi, 1925 ) describes the consolidation process of a finite layer. A constant stress is applied suddenly on the upper surface of a fluid-saturated layer ( Fig. 5 a). The only permeable boundary is where the load is applied, i.e. fluid can only drain from the upper surface. As the constant load continues, the sample keeps draining and gradually consolidates.
In numerical modelling, a column of length L = 15 m is used to represent the fluid-saturated finite layer; the crosssection is a 1 m × 1 m square. The column is discretised by 4-node tetrahedral elements with uniform size distribution; the average length of the element edge is h = 0.5 m. The positive z -direction is pointing downwards. Detailed mesh sensitivity tests have revealed that this mesh results in fully mesh independent solutions. For the solid phase, the only constraint is applied to the lower boundary ( z = 15 m), where no displacement is allowed. For the fluid phase, non-permeable boundary conditions are applied to all the boundaries, except the upper surface ( z = 0), where free drainage happens to maintain a constant pore pressure of zero.
The material used for this numerical test is assumed to be sand ( Ferronato et al., 2010 ) and the material properties are given in Table 2 . A time-step of t = 0.1 s is used in the modelling. The analytical solutions of Terzaghi's problem are given by Eqs. (37) and ( 38 ) ( Ferronato et al., 2010 ;Wang & Hsu, 2009 ). where p is the pore pressure and u is the vertical displacement (positive downward). The consolidation coefficient c is defined as where c M is the vertical uniaxial compressibility, and M is the Biot modulus, where φ is the porosity, β the fluid compressibility, and c br the solid grain compressibility. It should be noted that for sand it is reasonable to assume that the sand grain is incompressible, which means c br = 0. Therefore, the Biot modulus M can be expressed as where S is the specific storage coefficient defined in the governing equations.
In the analytical solutions, p 0 is the initial overpressure generated by the load applied instantaneously at time t = 0, given as and an initial vertical displacement where K u is the undrained bulk modulus, given by  The constant load applied at the upper surface of the column ( z = 0) is P L = 10 4 Pa ( Ferronato et al., 2010 ;Wang & Hsu, 2009 ).
The numerical results are compared and plotted with the analytical solutions in Figs. 6 and 7 . It can be seen that both match the analytical solutions very well. It should be noted that the vertical displacement at the long time limit used for normalisation is ( Wang,20 which is the fully drained response at the upper surface z = 0.

Mandel's problem
The second test for verification is a simulation of the Mandel's problem, which describes a poroelastic medium compressed in one direction by two identical plates. The plates are assumed to be rigid, frictionless and impermeable. It can be seen from Fig. 8 a that two constant loads are applied vertically on the plates. The fluid can only be drained from the two lateral surfaces, where there are no external forces.
In the numerical simulation, due to symmetry a cube with edge of 1 m is used to represent the poroelastic layer sandwiched between the plates (the shaded area in Fig. 8 a). This cubic domain is discretised by 4-node tetrahedral elements, with an average edge of h = 0.0625 m. A time-step of t = 0.01 s is used in the modelling. The sample is assumed to be made of sand as well, the same as the one used in the Terzaghi's problem ( Table 2 ).
The poroelastic layer in Mandel's problem is assumed to be infinitely long in the y -direction (the direction perpendicular to the paper), so it can be simplified into a two-dimensional plane strain problem ( Fig. 8 a), in which case the analytical where α n are the positive roots of the following equation In this analytical solution, c is the consolidation coefficient as defined in Eq. (39) ; a is the edge length of the cubic sample, a = 1 m; and p 0 is the initial overpressure generated by the loading suddenly applied to the plates at time t = 0, and ν u is the undrained Poisson's ratio. It can be derived from the drained Poisson's ration ν by The load F in Eq. (49) is set to be F = 10 4 N/m. Once again it can be seen that the numerical results are in very good agreement with the analytical solutions ( Figs. 9 and 10 ). Importantly, the non-monotonic response of the pore pressure, i.e. the Mandel-Cryer effect ( Cryer, 1963 ;Gibson, Gobert & Schiffman, 1990 ;Mandel, 1953 ;Schiffman, Chen & Jordan, 1969 ), has been accurately captured in the simulation ( Fig. 10 ). At the very early stage, a uniform distribution of pore pressure is generated by the Skempton effect ( Abousleiman, Cheng, Cui, Detournay & Roegiers, 1996 ;Skempton, 1954 ). As the consolidation progresses, the pore pressure near the lateral sides ( x close to a ) starts to decrease due to the drainage boundary conditions. In a poroelastic medium, both pore pressure and solid skeleton can contribute to the overall compressive stiffness. Therefore, the sample behaves in a more compliant manner at the lateral boundaries due to the reduction of pore pressure, while the central part of the domain ( x close to 0) carries more compressive stress by the compatibility requirement, i.e. the compressive stress is transferred from the compliant region to the stiff region. The result is that the pore pressure continues to rise under this excessive compressive stress. In the long time limit, pore pressure will vanish everywhere in the domain.

Convergence analysis
The convergence of the 4-MPET finite element model was tested on a three-dimensional anatomically accurate geometry of human brain to investigate mesh independence and provide guidance on mesh resolution for subject-specific modelling. Twelve meshes were created from the same brain geometry, with total element numbers ranging from 104 075 to 8 688 950. All the simulations used the same set of parameters and boundary conditions. As it can be evidenced from Fig. 11 , the  method exhibits good convergence characteristics and meshes of around two million tetrahedral elements are sufficient to capture salient features with good accuracy. Additional mesh independence testing results can be found in Guo et al. (2018) , where plots of mean values in the entire parenchymal domain, as well as at selected representative positions, e.g. the white matter, the cortical surface and the ventricular wall, were given.

Model validation
In this section, the numerical model is validated against experimental data from infusion tests on six C57BL/6 J mouse brains ( Vindedal et al., 2016 ), where the focus is on the circulation of the cerebrospinal fluid (CSF). The classical physiological theory believes that the CSF is produced in the choroid plexus and then flows through the cerebral ventricles and the subarachnoid space, until being reabsorbed into the bloodstream via arachnoid villi of the dural sinuses, along cranial nerve sheaths or through the nasal lymphatics ( Abbott, 2004 ;Koh, Zakharov, & Johnston, 20 05 ;Praetorius, 20 07 ). Any breakdown in this circulation can cause serious consequences -one of them is communicating hydrocephalus, where obstacles exist in the cerebrospinal fluid (CSF) flow between the ventricular system and the subarachnoid space. In this case, infusion test is often used for the diagnosis because imaging alone is not sufficient to identify the obstacles ( Boon et al., 20 0 0 ;Juniewicz et al., 2005 ). In clinical practice, the most frequently used procedure is to inject CSF-like fluid (e.g. saline) into the subarachnoid space at a constant rate, usually via a lumbar puncture. The infusion test has been simulated before either by mathematical modelling ( Cieslicki, 2007 ) or poroelastic modelling ( Sobey et al., 2012 ), but never using a highly physiologically specific multi-compartmental poroelastic model. During an infusion test, the injected CSF causes the intracranial pressure (ICP) to rise in order to accommodate the additional fluid, which also results in increased drainage. The ICP value recorded before infusion is called the baseline ICP. In the test, the ICP will reach a peak and then decline after the infusion is completed, but not return to the baseline value; this value is called post-infusion resting ICP. Two measurements are important here: The first one is the rate of rise to the peak providing an estimate of cerebral compliance and the second one is the difference between the baseline ICP and the peak value giving information about the resistance to CSF outflow ( Andersson, Malm & Eklund, 2008 ). A mathematical model to describe such phenomenon was first proposed by Marmarou, Shulman and Rosende (1978) and later further developed by Avezaat and Eijndhoven (1984) . The Marmarou model says that during an infusion test, the flow of CSF satisfies the following relationship according to the continuity requirement, where Q P is the production of cerebrospinal fluid (CSF), Q i the external infusion, Q R the reabsorption, and Q S the storage. In the numerical modelling of this section, the external infusion Q i = 2 μl/min (over a period of five minutes), as previously described ( Iliff et al., 2012 ) and it is the same value as the one used in the experiments ( Vindedal et al., 2016 ). The value of CSF production is obtained from literature ( Simon & Iliff, 2016 ), Q p = 0.35 μl/min. The reabsorption of CSF happens at the cortical surface at a rate of where p ss is the sagittal sinus pressure, p ss = 2 mmHg; μ is viscosity and R is resistance to CSF outflow, μR = 80 mmHg • s/ μl ( Jones, Deane & Bucknall, 1987 ;Moazen et al., 2016 ). Three flow rates are used in the numerical modelling as boundary conditions ( Eqs. (54) and ( 55 )) and the storage effect Q p is modelled by the poroelastic model itself because the storage capacity is an inherent characteristic of porous media. Details about the experiments can be found in the paper by Vindedal et al. (2016) . The geometry used for the validation is a spherical shell, bounded by two concentric spheres. This is a common practice if the focus is on the mechanical characteristics of the system ( Tully & Ventikos, 2011 ). The ICP is almost constant spatially so a relative coarse mesh is used for the model -the total number of 4-node tetrahedral elements is 1105. The radius of the outer sphere ( r = c ) is 4.95 mm and the radius of the inner sphere ( r = a ) is 1.05 mm ( Fig. 12 ). The volume that lies between these two spheres is the domain for numerical modelling, which represents the brain parenchyma of a mouse. The dimensions of the model are derived to make the volume of the parenchyma equivalent to the volume of a mouse (C57BL/6 J) brain ( Badea, Ali-Sharief and Johnson, 2007 ;Vindedal et al., 2016) .
The four-compartmental poroelastic model described in Section 3 is used for the modelling. The boundary conditions for the arterial compartment, the arteriole/capillary compartment and the venous compartment are the same as listed in Table 1 , where the outer spherical surface represents the cortical surface and the inner surface represents the ventricular wall. For the CSF/ISF compartment, modifications of boundary conditions are made to mimic the setup of an infusion test. More specifically, Robin boundary conditions are applied at both inner and outer spherical surfaces ( Sobey et al., 2012 ). At and it is assumed that at this surface the displacement is zero. At the inner surface boundary ( r = a ), which represents the ventricular wall, the boundary condition also defines the flux conservation of CSF, 4 π a 2 ∂u ∂t = − π d 4 128 μL p − p | r= c + 4 π a 2 k μ ∂ p ∂r and there are no displacement constraints at this boundary. From the geometry of the model, it can be derived that the length of the aqueduct (a single thin canal connecting the ventricles to the subarachnoid space) is L = c -a = 3.9 mm, and the diameter of the aqueduct d = 0.04 mm ( Feng et al., 2009 ). The flow in the aqueduct is assumed to be Poiseuille flow ( Sobey et al., 2012 ), based on which the first terms on the right-hand sides of Eqs. (54) and ( 55 ) can be decided. The permeability for the CSF/ISF compartment k e is 1.5 × 10 −8 m 2 and the viscosity μ e is 8.0 × 10 −4 N • s/m 2 ( Binder, Papadopoulos, Haggie & Verkman, 2004 ;Oshio, Watanabe, Song, Verkman & Manley, 2005 ); other parameters used in the validation test are listed in Table 3 ( Combe et al., 2016 ;Wilde et al., 2017 ). The time-step used in the modelling is t = 1 s.
In the experiments ( Vindedal et al., 2016 ), six wild types of C57BL/6 J mouse were tested. The intracranial pressure (ICP) recordings in response to injections of artificial cerebrospinal fluid are plotted in Fig. 13 . The results from numerical modelling are plotted in the same figure by the black line.
It can be seen that the numerical results fall within the range of experimental data. This validation test demonstrates that the change of ICP during an infusion test can be successfully modelled by the numerical simulation. Both the infusion and decay stages, which are separated by the peak of ICP, are distinctively captured. It should be noted that the highfrequency fluctuations of ICP from experiments are mainly caused by arterial pressure pulsations, which are not considered in the numerical modelling due to the lack of experimental data to quantify the input parameters. Therefore, the plot from numerical results is smooth without pulsatile fluctuations.

Model application to three-dimensional subject-specific human brain
The numerical solver of the three-dimensional MPET model can be applied to biological systems to investigate fluid flow, transport phenomena and deformation of tissue, etc. A successful simulation generally requires three categories of data to be fed into the MPET solver: finite element mesh, material parameters and boundary conditions. As an example, a schematic illustration of using the n -MPET model for human brain is shown in Fig. 14 , where n represents the number of fluid networks. The finite element mesh reflects the geometric characteristics and also controls the spatial resolution of the numerical results; the material parameters are required to completely define the governing equations Eqs. (3) -(7) of the  For biological systems, the individual variability and inherent complexity make it necessary to conduct numerical modelling in a subject-specific manner so quantitative understanding of their structures and functions in health and disease can be obtained ( Taylor et al., 1999 ;Winslow, Trayanova, Geman & Miller, 2012 ). To achieve such a goal by using the MPET model, ideally, every element in the ellipse in Fig. 14 should use subject-specific expressions or values. However, the reality is that many of them are not available on a subject-specific basis; this limitation is further discussed in Section 7 . In the examples presented in this section, the subject-specific requirements are considered in the modelling from three aspects 2.9 × 10 −4 m 2 N −1 Q p 5.8 × 10 −9 m 3 s −1 : First, the finite element meshes are generated from T1-weighted magnetic resonance (T1w-MR) images, which can represent anatomically accurate geometry of the individual brain. Second, the permeability of the cerebrospinal fluid/interstitial fluid (CSF/ISF) compartment is extracted from diffusion-weighted imaging (DWI), which can provide a heterogeneous and anisotropic permeability distribution for the brain parenchyma. Third, the boundary condition for the arterial compartment at the cortical surface is obtained from clinical measurements. The details about these three aspects of the subject-specific brain modelling pipeline can be found in the paper by Guo et al. (2018) . All the other material parameters and boundary conditions ( Tables 1 and 4 ) are either from literature or previously tested in numerical modelling ( Smillie et al., 2005 ;Tully & Ventikos, 2011 ;Vardakis et al., 2013 ).
The arterial blood supply to the entire brain, including carotid arteries (left and right) and vertebral arteries (left and right), was obtained through a combination of ambulatory blood pressure measurements, clinical ultrasound flow measurements and mathematical modelling . A lumped parameter circulation model (LPCM) ( Ursino, 1998 ) was used to translate spot measurements collected at 15-min intervals to continuous waveforms of arterial blood flow, which are used as boundary conditions for the arterial compartment at the cortical surface . Each waveform has duration of about one second and their amplitudes vary during the day. Here two specific moments are chosen as indicative of the subject's activity in a period of 24 h -the first one is high activity (e.g. exercise) identified by the highest peak value of the arterial blood supply; and the second one is low activity (e.g. rest or sleep) identified by the lowest peak value. To ensure that the numerical results represent steady states, in every simulation the same waveform is repeatedly used 50 times and only the last output is regarded as the final result. This criterion of running 50 cycles is determined by comparing the final results after different numbers of cycles and it is found more than 50 cycles only generate negligible differences. All the simulations in this section use a fixed time-step of 0.1 s.
Three variables are selected from the output files to link mechanistic modelling with clinical practice: Blood perfusion is defined as the filtration velocity (Darcy velocity) of the arteriole/capillary compartment, i.e. v c ; CSF/ISF clearance is defined as the Darcy velocity of the CSF/ISF compartment, i.e. v e ; and the deformation of tissue is measured by the displacements of finite element nodes, i.e. u . In this section, the results from subject-specific modelling are compared between age and gender matched subjects in the context of investigating the effect of smoking on cerebral fluid transport. It should be noted that the comparisons do not indicate any statistical significance but focus on individual differences, which demonstrate the ability of the MPET model to generate clinically relevant output from subject-specific input data. All of the numerical results presented in this section were obtained from the subject-specific brain modelling pipeline of the VPH-DARE@IT project , which was implemented on the MULTIX platform ( Kasztelnik et al., 2017 ). The study 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.
Chronic cigarette smoking is generally recognised as a lifestyle factor that can cause neurobiological deficits and disturbances in brain perfusion ( Durazzo, Mattsson & Weiner, 2014 ), which are related to stroke and vascular dementia, etc ( Rogers et al., 1983 ). To test the ability of the MPET model to capture the differences between smokers and non-smokers, two subjects are simulated using the 4-MPET model; both are healthy, female and 68 years old. The geometry of the smoking subject's brain is discretised by 2 294 536 tetrahedral elements, and 2 533 529 elements are used for the non-smoking subject's brain. Both high and low activities are modelled for each subject so in total there are four scenarios.
First, the total arterial blood supply to the brain in these four scenarios is plotted in Fig. 15 . It shows that at both high activity (e.g. exercise) and low activity (e.g. rest or sleep), the total arterial blood flow to the brain is lower in the case of the smoking subject than the non-smoking subject. This is in agreement with the physiological observations that show cigarette smoking impairs cerebral blood flow ( Rogers, Meyer, Judd & Mortel, 1985 ;Yamashita, Kobayashi, Yamaguchi, Kitani & Tsunematsu, 1988 ;Toda & Okamura, 2016 ), which is due to cerebral vascular dysfunction and has been identified as a major clinical feature of Alzheimer's disease ( de la Torre, 2012 ).
Next, the differences between high and low activities are compared between smoking and non-smoking subjects. The contours in Figs. 16-18 are obtained by conducting the following calculation at every node in the entire domain, where v diff is the difference between high activity ( v high ) and low activity ( v low ); here v can be any vector (or scalar) from the output variables. In this part, the filtration velocity (Darcy velocity) of the arteriole/capillary compartment v c is used to Fig. 15. The total arterial blood supply to the brain at high and low activities between a smoking subject and a non-smoking subject.   show brain perfusion; the magnitude of the displacement vector u is used to show the tissue deformation; and the Darcy velocity of the CSF/ISF compartment v e is used to show CSF/ISF clearance. One important difference between smokers and non-smokers is brain perfusion -smokers show lower perfusion than non-smokers (hypoperfusion), both globally and in multiple regions ( Durazzo, Meyerhoff & Murray, 2015 ;Siennicki-Lantz, Reinprecht, Wollmer & Elmståhl, 2008 ), which can be seen from the numerical results. In high and low activities, the perfusion of the smoking subject has a maximum value of 3.05 × 10 −4 m/s and 2.89 × 10 −4 m/s, respectively ( Fig. 16 a,  c); while the non-smoking subject has a higher maximum value of 3.06 × 10 −4 m/s and 2.98 × 10 −4 m/s, respectively ( Fig. 16 b, d). Moreover, it can be seen that the smoking subject shows greater difference between high and low activities, i.e. the brain perfusion of the smoking subject is not as stable as the non-smoking subject in the 24-hour period ( Fig. 16 e, f). The contours exhibit a similar pattern, e.g. the closer to the posterior horns of lateral ventricles, the greater the difference. However, the magnitude of the difference of the smoking subject is more than twice the value of the non-smoking subject. The differences of brain parenchyma deformation between high and low activities show similar magnitudes and patterns between the smoking and non-smoking subjects ( Fig. 17 ), where the values are at the same order of magnitude and largest differences are around the lateral ventricles.
CSF of the central nervous system (CNS) is critical for the clearance of solute from the brain. The variables related to the CSF/ISF compartment show marked differences between smoking and non-smoking subjects, which is mainly due to the subject-specific permeability maps used for the CSF/ISF compartment. In the modelling, every fluid network is assigned a base permeability value k i ( Table 4 ). For the arterial compartment ( k a ), arteriole/capillary compartment ( k c ) and the venous compartment ( k v ), the permeability is assumed to be an isotropic and homogeneous property due to the lack of data from clinical measurements; however, for the CSF/ISF compartment ( k e ), a three-dimensional permeability map was generated from processing diffusion-weighted imaging (DWI) for every subject Ravikumar, Gooya, Çimen, Frangi & Taylor, 2018 ), which assigns an anisotropic permeability tensor (dimensionless) for each tetrahedral element in the finite element mesh. In the simulations, the value of permeability for the CSF/ISF compartment is the base permeability multiplied by the dimensionless permeability tensor.
For the differences of CSF/ISF clearance between high and low activities, it can be seen from Fig. 18 that both the patterns and magnitudes are different between the smoking and non-smoking subjects. The non-smoking subject shows a more uniform field, while the smoking subject has a greater gradient. This can also be seen from the magnitude perspective, where it shows the maximum difference of CSF/ISF clearance is about 3.5 times higher in the smoking subject. This indicates that the smoking subject has greater fluctuations of CSF/ISF clearance in the 24-hour period. The dimensionless fractional anisotropy (FA) derived from diffusion-weighted imaging of the two subjects is shown in the right column of Fig. 18 . The fractional anisotropy indicates the deviation from pure isotropic diffusion and is used to quantify the microstructural integrity of the white matter . Previous studies showed that smoking was associated with reduced microstructural integrity of the white matter, such as white matter hyperintensities (WMH) ( Gons et al., 2011 ;Habes et al., 2016 ). From an imaging point of view, fractional anisotropy is believed to decrease with impaired structural integrity of the white matter ( Le Bihan et al., 2001 ;Pierpaoli, Jezzard, Basser, Barnett & Di Chiro, 1996 ;Wardlaw, Hernández & Muñoz-Maniega, 2015 ); it can be seen from Fig. 18 that the value of fractional anisotropy is indeed smaller in the smoking subject, especially in the parietal lobes.

Discussion
The computational platform of the three-dimensional multiple-network poroelastic model ( n -MPET) can be used to conduct steady-state static and transient quasi-static analyses of porous biological systems, such as brain parenchyma, kidney tissue and bone tissue. The in-house numerical code is written in Fortran 95/2003, and will later be released as open source. It has a flexible structure in the following aspects. First, the number of fluid networks n is not a fixed number but can be customised to specific applications. On the one hand, latest development in the understanding of physiology, such as the glymphatic system, which postulates that waste products are cleared from the brain via a glia-dependent mechanism, analogous to the lymphatic system in other organs ( Iliff et al., 2012 ;Nedergaard, 2013 ), can be added to the numerical model by increasing the value of n . On the other hand, the number n can also be reduced for validation purposes or applications where additional fluid networks have negligible influences on the problem of interest. Second, the material properties and boundary conditions can use pre-defined default values in the code or subject-specific data if available. The difference in execution is to import extra data files into the solver. Third, the pre-processing and post-processing support various formats and open-source software, such as Netgen ( Schöberl, 1997 ), ISO2Mesh ( Fang & Boas, 2009 ), ParaView ( Ayachit, 2015 ) and Mayavi ( Ramachandran & Varoquaux, 2011 ), which gives users the freedom to choose the one that they are most familiar with.
Several limitations still exist in the current model. First, the values of material properties are not the most accurate from a physiological perspective. For example, there was research showing the properties (e.g. Young's modulus) of the white and grey matter are different ( Budday et al., 2015 ;Christ et al., 2010 ); however, many numerical models ( Levine, 20 0 0 ;Tully & Ventikos, 2009 ;Wirth & Sobey, 2006 ), including the one presented in this paper, still treat the brain parenchyma as homogeneous material due to the lack of quantitative experimental data of human brain. In this aspect, magnetic resonance elastography (MRE), which is a non-invasive imaging technique to quantify mechanical properties of tissue ( Feng et al., 2018 ;Green, Bilston & Sinkus, 2008 ), can offer more realistic values for numerical modelling. Recent research also suggested that the typical permeability currently used in modelling was too high by one or two orders of magnitude ( Holter et al., 2017 ), which needs further investigation as well. Moreover, the parameters used in numerical modelling need to be validated. This has two-fold meaning: First, the values themselves must by physiologically reasonable; and second, the numerical results obtained from these values must also be validated against experimental or clinical data.
Second, more effort needs to be put into refining the underlying biological and mechanical mechanisms of the multiplenetwork poroelastic model used for biological systems. One example of such physiological mechanisms is the transfer of fluids between different fluid compartments. In the current model, this transfer purely depends on hydrostatic pressure gradients, while in the biological context it is more likely to be driven by more complicated biological processes, such as via osmotic pressure -a fact for biological fluids that have a lot of minerals including ions dissolved in them. This limitation can be overcome by implementing empirical relationships once they are validated ( Tully & Ventikos, 2011 ). Another example in the aspect of mechanical mechanisms is the cross-porosity storage effect ( Mehrabian & Abousleiman, 2017 ;Vardakis et al., 2017 ). In the future, it is worth exploring the differences in numerical results between models neglecting and employing such effects.
Third, as already pointed out in Section 6 , the subject-specific modelling at the current stage has not reached an ideal level of personalisation in order to differentiate between subjects. The higher the complexity of the model, more factors need to be specified. The main reason for this limitation is that the data are not available from experiments or in vivo measurements. It can be seen from Fig. 14 that for a complete subject-specific simulation, every element in the ellipse should be personalised, which itself is a non-trivial task. A solution that can improve this situation is to translate more physiological understanding into mechanical relationships and mathematical expressions, and implement them into the numerical model; then subject-specific data can be fed into the modelling. Similar to other interdisciplinary research, the difficulty here is to find connections across the boundaries of disciplines. Ideally, the entire workflow to develop the model should form a virtuous circle -better understanding of physiology helps improve numerical models and more advanced numerical models provide testing beds for testing hypotheses and developing theories.

Conclusions
A finite element model for biological systems has been developed based on the multiple-network poroelastic theory ( n -MPET). The governing equations were defined for the mathematical model, and boundary conditions applicable to the cerebral environment were given to provide supplementary physical meaning from a physiological point of view. The model was implemented into an in-house finite element code. The numerical model was verified against the analytical solutions of the Terzaghi's problem and the Mandel's problem, and validated against experimental data from infusion tests on mouse brain. The model was successfully applied to three-dimensional subject-specific modelling of water transport in the cerebral environment, which included anatomically realistic geometry, personalised permeability map for the diffusion of cerebrospinal fluid, and subject-specific arterial blood supply to the brain. The numerical results were compared between smoking and non-smoking subjects, which showed hypoperfusion in the brain of the smoking subject.
The three-dimensional n -MPET model has been demonstrated to be capable of capturing the spatio-temporal fluid transport in biological systems across different scales. Clinically relevant mechanisms, such as blood perfusion, cerebrospinal fluid clearance and tissue deformation, can be visualised in fully three-dimensional anatomically realistic models. The numerical code will later be released as open source, which will create an open and flexible platform for researchers to tailor the numerical tool for their specific biomechanical applications, which can help translate computational modelling paradigms into clinical practice.

Declaration of Competing Interest
None.