Cerebrospinal fluid dynamics coupled to the global circulation in holistic setting: Mathematical models, numerical methods and applications

Abstract This paper presents a mathematical model of the global, arterio‐venous circulation in the entire human body, coupled to a refined description of the cerebrospinal fluid (CSF) dynamics in the craniospinal cavity. The present model represents a substantially revised version of the original Müller‐Toro mathematical model. It includes one‐dimensional (1D), non‐linear systems of partial differential equations for 323 major blood vessels and 85 zero‐dimensional, differential‐algebraic systems for the remaining components. Highlights include the myogenic mechanism of cerebral blood regulation; refined vasculature for the inner ear, the brainstem and the cerebellum; and viscoelastic, rather than purely elastic, models for all blood vessels, arterial and venous. The derived 1D parabolic systems of partial differential equations for all major vessels are approximated by hyperbolic systems with stiff source terms following a relaxation approach. A major novelty of this paper is the coupling of the circulation, as described, to a refined description of the CSF dynamics in the craniospinal cavity, following Linninger et al. The numerical solution methodology employed to approximate the hyperbolic non‐linear systems of partial differential equations with stiff source terms is based on the Arbitrary DERivative Riemann problem finite volume framework, supplemented with a well‐balanced formulation, and a local time stepping procedure. The full model is validated through comparison of computational results against published data and bespoke MRI measurements. Then we present two medical applications: (i) transverse sinus stenoses and their relation to Idiopathic Intracranial Hypertension; and (ii) extra‐cranial venous strictures and their impact in the inner ear circulation, and its implications for Ménière's disease.


| INTRODUCTION 2 | MATHEMATICAL MODELS
In this section, we present the mathematical models for describing the human circulatory system and the CSF dynamics in the craniospinal space. As already pointed out, the present paper is based on two independently developed pieces of work. The circulatory system part emerges from References 59 and 69, with some significant improvements. For the 1D representation of the blood vessels we include viscoelasticity of the vessel wall. 39,70 The resulting partial differential equation system is solved numerically using the high-order Arbitrary DERivative Riemann problem (ADER) framework 71 with a solver that allows for local time stepping (LTS). 72 The microcirculation and heart models are also partially modified with respect to Reference 59. Representation of pulmonary circulation follows the same models as in References 59 and 69 wherein a simple three-compartment (arteries, capillaries, veins) description for systemic microcirculation based on the work by Sun et al 73 is adopted. Venous valves and SR are modeled following. 41 The second piece of work underpinning this paper is the coupling of the circulation to a refined mathematical description of the CSF dynamics in the craniospinal cavity, for which we follow the model proposed by Linninger. 62 The Linninger model 63 was chosen because it integrated results of extensive MR imaging studies [47][48][49] with detailed two and 3D mathematical models into a comprehensive mathematical description of the major intracranial dynamics with fluid structure interactions of blood, CSF in the ventricular system, as well as cranial and the SSASs and the deformable brain parenchyma.
Since in this article we are particularly interested in the cerebral dynamics, we incorporate into our model one of the most important physiological control systems, the cerebral autoregulation. This process aims at maintaining adequate and stable CBF during changes in blood pressure working on dilatation or contraction of arterioles and capillaries. 77 The model used to account for this phenomenon is based on References 75 and 76.
F I G U R E 1 Schematic representation of the global model for the full human circulatory system coupled to the craniospinal fluids and brain parenchyma. The arterial 1D network is represented in the right dotted box with red vessels while the venous 1D network is displayed in the left dotted box with blue vessels. The terminal arteries of the arterial network are connected to draining veins of the venous circulation through 0D models representing arterioles, capillaries, small arteries (red boxes) and venules, small veins (blue boxes). The dotted arrows indicate the connection between 1D network and terminal vessels, depicted for simplicity as dots in the yellow bar. Left and right cardiac chambers are displayed by red and blue boxes, respectively, connected to green atrioventricular valves; the left ventricle is connected to aortic root and venae cavae are linked to right atrium through black arrows. Cardiac chambers are connected to the pulmonary circulation, comprising arteries, capillaries and veins, represented by the pink box, through the aortic and pulmonary valves (green arrows). The CSF compartments are represented by cyan boxes. The arrows between cardiovascular system and CSF circulation represent the fluid exchange between these systems through production and absorption while the green dashed boxes indicate the mechanical interaction between between all components in the cranial cavity through the Monro-Kellie hypothesis

| Equations for blood flow in major vessels
In this section, we review the partial differential equations representing blood flow in major vessels, along with closure laws and reformulations.

| Conservation laws and closure conditions
The flow of blood in major arteries and veins is represented through 1D cross-sectional averaged models resulting in time-dependent systems of partial differential equations. We start from the classical laws of conservation of mass and of balance of momentum For details on the derivation of (1) see Reference 29, for example. The 2 Â 2 system (1) contains three unknowns, namely A x, t ð Þ, the cross-sectional area of the vessel lumen; q x, t ð Þ, the blood flow rate and p x,t ð Þ, the crosssectionally averaged internal pressure. Parameters in the equations include b α, the Coriolis coefficient, ρ the blood density, assumed constant, and f the friction force per unit length of the tube. The Coriolis coefficient depends on the assumed velocity profile; here we take b α ¼ 1, which corresponds to an assumed parabolic velocity profile.
The system of differential Equations (1) has more unknowns than equations; hence, one extra closure condition is required. Such extra condition, usually called tube law, attempts to couple the internal blood flow distribution with the mechanical properties of the solid moving vessel wall. A comparative analysis of various mathematical descriptions of elastic properties of vessel walls in modern 1D models of hemodynamics can be found in Reference 77. In the existing versions of our model, 59,69 we used elastic tube laws for both arteries and veins. In the present paper, we improved upon this by adopting viscoelastic tube laws for both arteries and veins in the entire circulation. To this end, we follow recent works concerned with approximating time-dependent parabolic systems with hyperbolic balance laws with stiff source terms, 78,79 The approach was applied in Reference 70 to a simplified arterial network. In the present paper, we deploy the framework to the full human circulatory system, major arteries and veins. Following Reference 62, the following pressure-area relation (tube law) is adopted Here, the first term ψ A x, t ð Þ;A 0 x ð Þ, K x ð Þ, P 0 ð Þ is the elastic part of the tube law, which in turn depends on the reference pressure P 0 and the parameters A 0 x ð Þ and K x ð Þ and can be written as The term φ A x, t ð Þ; A 0 x ð Þ ð Þ ∂ t A represents the viscoelastic part of the tube law, while p ext x, t ð Þ denotes the external pressure. As usual, the transmural pressure is defined as We now define geometric and mechanical parameters in the tube law. A 0 x ð Þ defines the vessel cross-sectional area at equilibrium; K x ð Þ represents the vessel wall stiffness, wile m and n are two real numbers, to be specified. Note that A 0 x ð Þ and K x ð Þ are variable parameters, they depend on distance x along the vessel.Throughout this work, we assume m ¼ 0:5 and n ¼ 0 for arteries, while m ¼ 10 and n ¼ À1:5 for veins. Moreover, K x ð Þ is a positive function that was obtained from the reference wave speed c 0 assumed for each vessel, distinguishing arteries, veins and dural sinuses.
The viscoelastic term of the tube law depends on the time partial derivative of the cross-sectional area of the vessel and is defined as Γ is related to the viscoelastic properties of the vessel wall, which following 39 is chosen as where γ is the wall viscosity. The wall viscosity is evaluated as the product of the viscoelastic parameter K M and the volume fraction of smooth muscle. K M is chosen such that hysteresis behavior of pressure-area plots in peripheral arteries and veins lies within the physiological range. For arteries we take K M ¼ 3 Â 10 5 dyn s/cm 2 and a percentage of smooth muscle of 10%. 80 For veins we take K M ¼ 5 Â 10 4 dyn s/cm 2 and a smooth muscle fraction of 8%. Concerning the wall thickness h 0 x ð Þ we follow 80 and express it in relation to the vessel radius at equilibrium. For arteries h 0 ¼ 10%r 0 , while for veins h 0 ¼ 5%r 0 .
The momentum equation in (1) contains the friction term f x, t ð Þ, which depends on the local velocity profile (assumed a priori). Here, we take with μ being the blood dynamic viscosity.

| Variable material properties and augmented equations
As already pointed out, the material and geometric parameters K x ð Þ, A 0 x ð Þ and p ext x,t ð Þ are in general functions of distance x along the vessel length. Computationally, in order to deal with this situation we adopt the variable-parameter formulation of Toro and Siviglia, 81 admitting now, viscoelastic tube laws for arteries and veins. 70 System 1, along with the tube law (2), is then written as the following extended 5 Â 5 system In succinct form system (9) reads where Q is the vector of unknowns A Q ð Þ is the coefficient matrix , ð12Þ S Q ð Þ is the source term vector and is a higher-order differential term emanating from the viscoelastic part of the tube law. This last differential term in the advection-diffusion-reaction system defines a parabolic problem, no longer hyperbolic, as in the case of a purely elastic tube law. 81 2.1.3 | Hyperbolic approximation of a parabolic system Toro and Montecinos,78,79 proposed a method to approximate time-dependent parabolic problems by hyperbolic systems with stiffsource terms, by extending the Cattaneo relaxation approach. 82 Montecinos and coworkers 70 applied the procedure to a network of arterial blood vessels, thereby setting the bases for its extension to the global, closed-loop circulation model of this paper, including arteries and veins.
To start with, a new variable ζ and a relaxation parameter ϵ > 0 are introduced such that ζ ! ∂ x q as ϵ ! 0: One then replaces the spatial gradient of the flow variable in (9) by the new variable ζ in (15), which is constrained to satisfy an additional evolutionary PDE, namely the constitutive Cattaneo's law, which reads We now have an augmented 6 Â 6 system. Keeping the same notation for the vector of unknowns, the coefficient matrix and the source term we may write with 0  0  0  0   2   6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  5 , ð19Þ Here, c is the wave velocity (analogous to the sound speed) associated to theelastic tube law. The 6 Â 6 system (17) with the state vector (18) and coefficient matrix (19) is hyperbolic, 70,78 provided the relaxationparameter is chosen so as to satisfy All eigenvalues of the coefficient matrix (19) are real and given as where now e c denotes thewave speed associated to the complete tube law.
The corresponding right eigenvectors are These eigenvectors associated to the real eigenvalues (23) can be shown to be linearly independent, and hence system (17) is hyperbolic.
We now consider two fundamental properties of system (17), namely the nature of the characteristic fields and the generalized Riemann invariants. It can be shown that the characteristic fields associated to eigenvectors R 1 and R 6 are genuinely non-linear, while the remaining characteristic fields are linearly degenerate. The generalized Riemann invariants associated to R 1 and R 6 are respectively. We also note that for constant p ext , A 0 and K, the generalized Riemann invariants for the linearly degenerate fields associated with R 2 , R 3 , R 4 and R 5 are and where More details about the hyperbolic reformulation of the problem and its eigenstructure are found in References 79, 70. Next, we present the 0D mathematical models, which consist of systems of ODEs.

| Equations for lumped-parameter models
In the previous section, we described mathematical models for major arterial and venous blood vessels, consisting of systems of partial differential equations. Here, we present compartmental, or 0D, models consisting of systems of ODEs, for other districts of the circulation. These include the microvasculature (arterioles, capillaries and venules/veins), the heart, the pulmonary circulation, venous valves and SR. We also present a mathematical model for cerebral autoregulation, which works on the terminal portion of the cerebral arteries and on the cerebral vascular beds. Lumped-parameter models for CSF compartments 62 will be introduced in section 2.3.

| The microvasculature
Physiologically, the arterial system is connected to the venous system through arterioles, capillaries and venules. To describe this connection, the microvasculature is represented by lumped-parameter, or 0D, models. This connection can be simple, such as between one artery and one vein, or entail numerous compartments. The generic vascular bed model for all microvasculature beds, depicted in Figure 2, is inspired in the three-element Windkessel model. The model is characterized by • Characteristic impedances. These couple any number of connecting 1D arteries/veins to lumped-parameter models for the microvasculature (R da or R vn ) and regulate the pressure drop between 1D domains and vascular beds; • Peripheral resistances and compliances divided into arterioles R al ,C al ð Þand capillaries R cp ,C cp À Á ; • Venous compartment with related compliance C vn .
As illustrated in Figure 2, each connecting artery can be linked to one or both venous capacitors, while each venous capacitor can be connected to any number of terminal veins. Note that the second artery splits into both venous capacitors, while the other arteries supply only one of them. Moreover, any number of veins can be connected to each venous capacitor.
For each vascular bed, the variables to be computed are pressure and flow, denoted as follows: for arterioles P al , Q al ; for capillaries P cp , Q cp and for venules P vn and Q vn . Thus, for each element of the vascular bed one has where C is compliance and P ext is the external pressure, generally taken as zero (relative to atmospheric pressure) or equal to the intracranial pressure in the case of intracranial peripheral beds. Q in and P out are flow and pressure in neighboring compartments or obtained from the 1D models through boundary conditions. Well-matched coupling to the connecting 1D arterial and venous segments is achieved via the characteristic impedances, as suggested in Reference 83.

| Valves, Starling resistors and stenosis
Here, we describe a valve model based on Reference 41 that predicts valve motion on the basis of the instantaneous difference between upstream and downstream pressures. In the present paper, the model is applied to describe cardiac valves, venous valves, SR in the cerebral circulation and stenosis. Here, we illustrate the general model, while for each specific application, more details are given in sections 2.2.3, 2.2.4 and 6. The valve dynamics is described by means of the function A e t ð Þ defining the effective cross-sectional area, expressed by F I G U R E 2 Schematic representation of the generic vascular bed model. The red boxes include the terminal 1D arteries and the corresponding arterioles and capillaries compartments; the green boxes represent the venous capacitors while the light blue boxes refer to terminal venules/1D veins where A a is the annulus area and the time-dependent function ξ t ð Þ is the valve state, a dimensionless number in the range 0, 1 ½ and representing the rate of opening and closing of the valve. ξ t ð Þ is given by the solution of the following variable-coefficient ODEs Here, K o and K c are rate coefficients for opening and closing respectively; Δp o and Δp c are opening and closure threshold pressures. In (32), M r and M s are parameters representing regurgitating and stenotic valves respectively. Specifically, a healthy valve corresponds to M r ¼ 0 and M s ¼ 1, an incompetent valve is described by M r > 0, while a stenotic valve is represented by M s < 1. Flow variation in time across the valve is given by a first-order, variable coefficients, ODE Δp t ð Þ is the pressure difference across the valve length, defined as where p up t ð Þ and p down t ð Þ are the upstream and downstream pressures at time t, with respect to valve direction. p up t ð Þ and p down t ð Þ are evaluated from other compartments of the global model, specified later. L t ð Þ and B t ð Þ are timedependent coefficients; L t ð Þ is blood inertance, which accounts for the component of the pressure difference related to blood acceleration; B t ð Þ is the Bernoulli's resistance, which governs pressure differences related to convective acceleration and dynamic pressure losses due to diverging flow. They are expressed by Here, ρ is the constant blood density and l e is the effective length. More details on the valve model are found in References 41,59, 75 and in the appendix of Reference 84.

| Heart and pulmonary circulation
In the present paper, we consider a heart model for the dynamics of the four chambers and the cardiac valves. The chambers are denoted as ch ¼ RA, RV , LA,LV f g , where RA and RV are the right atrium and ventricle, while LA and LV represent the left atrium and ventricle. For the chambers we follow thetime-varying elastance model in References 58,73, while cardiac valves are modeled as presented in section 2.2.2 following. 41 Briefly, blood pressure in each cardiac chamber is calculated as where V ch is the current cardiac volume; V ch,0 is the dead volume (assumed to be 0); γP ch t ð Þ is the viscoelasticity coefficient of the cardiac wall and P peri t ð Þ is the external pericardial pressure defined by where V H t ð Þ is the sum of the volume of each heart chamber and V PC , Φ PC are constant parameters. P ext t ð Þ is the external pericardial pressure and E t ð Þ in (38) is the time-varying elastance, defined as where the constants E A and E B are the active and passive elastances, respectively, while e t ð Þ is the normalized timevarying elastance, taken as in Reference 73 as eðtÞ e a ðtÞ ¼ for atria and as e t ð Þ e v t ð Þ ¼ for ventricles. The modeling of cardiac valves is described in what follows. As each chamber of the heart contracts, blood is pushed through a valve either into another chamber or out of the heart into an artery (aorta or pulmonary). The four cardiac valves (tricuspid, pulmonary, mitral and aortic) ensure one-way blood flow by (a) opening to let blood through and (b) closing to prevent backflow. The mechanism that leads to opening or closure of a valve is driven by blood pressure changes as the heart contracts and relaxes. Such mechanism is modeled here following section 2.2.2. The pressure drop across each cardiac-valve length is defined from pressure data from the neighboring cardiac chamber, from the aortic root and from the pulmonary arterial compartment. In particular, the tricuspid and mitral valves are atrioventricular valves that prevent backflow of blood from the ventricles into the atria; in these cases, the upstream pressure is the pressure of the atrium while the downstream pressure is that of the ventricle. The pulmonary valve is located at the opening between the right ventricle and the pulmonary trunk, therefore its upstream and downstream pressures are the right ventricle and the arterial pulmonary pressure, respectively. Finally, the pressure drop across the aortic valve is determined by the difference between the left ventricle and the aortic root pressure. In Equation (33), Δp o and Δp c are set to 0 for all the cardiac valves. Other parameters are defined later in section 4.1.

| Venous valves and Starling resistors
The Müller-Toro global model 59,69 is equipped with submodels for venous valves and SR consisting of ideal diodes. In the present article, these elements are replaced with the model described in section 2.2.2.

Venous valves
Venous valves are placed at different locations of the venous network; each venous valve is located between two venous vessels and governs flow across this interface. Location of valves are reported in Table 5 and are related to the venous network depicted in Figures 5 and 6. In this case, the pressure drop which governs the flow rate across the valve is determined by the pressure difference between the upstream and downstream vessels with respect to the valve direction. Moreover, the annulus area A a is assumed to be equal to the mean area between the reference areas of the connecting vessels on the right and left side. In the same way, the effective length l e is taken as the mean diameter at equilibrium between the upstream and downstream vessels. Finally, Δp o and Δp c are set to 0.

Starling resistors
Venous CBF is ensured by a SR mechanism, a fluid dynamic construct which governs the flow in collapsible tubes exposed to variable external pressure. The SR act at the confluence of cortical veins in the dural sinuses; these are located in the dura mater and are more rigid than cerebral veins. During the large physiological fluctuations of the intracranial pressure, the SR mechanism prevents the vein collapse maintaining the blood pressure upstream the collapsed segment higher than the intracranial pressure.
In previous work, we described SR simply via ideal diodes. 69 Here, we adopt a model that allows for a richer description of opening/closing dynamics, as well as accounting for a better description of underlying physical processes. In our venous network, Starling resistor behavior is represented through the model proposed in section 2.2.2. The pairs of vessels between which each SR element is placed are reported in Table 5, where the left vessel index represents the number of the cortical vein while the right vessel index indicates the corresponding venous sinus. The annulus area A a and the effective length l e are assumed to be the mean area and diameter, respectively, between the reference area and diameter at equilibrium of the pair of vessels connected to the SR element. The flow across SR is limited to that given by the pressure difference between the cerebral vein (upstream vessel) and the larger of intracranial pressure and the downstream pressure. When the intracranial pressure is much higher than the downstream pressure, the flow rate through the vessel becomes independent of the downstream pressure and the driving pressure difference is given by the upstream and the external pressures.As for other intracranial compartments, the external pressure is the intracranial pressure, that in this work is taken as the pressure in the fluid part of the brain parenchyma. Therefore, the valve state in Equation (33) is determined by Δp À Δp o ¼ Δp À Δp c ¼ p down t ð Þ À p ext t ð Þ, where p ext is the intracranial pressure. If p down t ð Þ < p ext t ð Þ, in Equation (34), the driving pressure difference Δp is given by Δp ¼ p up À p ext ; on the other hand, if the downstream pressure is higher than the external pressure, the flow across the SR element is determined by the pressure difference between the upstream and downstream pressures Δp ¼ p up À p down . Other parameters used in the valve and SR models are reported in.

| Control system: Cerebral autoregulation
We consider a model of the cerebrovascular regulation mechanisms, which acts by modifying resistances and compliances of the arterial microcirculation; changes in these parameters are not independent but are related through biomechanical and geometrical laws. The model is based on References 75 and 76, with appropriate modifications. Only one control mechanism is considered in this work, the myogenic response, which is linked to changes in arterial pressure and CBF. The original autoregulation model proposed in Reference 75 reproduces also the metabolic response of cerebral autoregulation, which is linked to carbon dioxide reactivity and the amount of oxygen reaching the brain tissue. This mechanism is not included in the present paper, as the current version of our model does not yet include a submodel for the transport of CO 2 and oxygen in the brain. It is an aspect to be considered in the near future.
Cerebral myogenic autoregulation is activated by changes in CBF; its action on the arterial microvasculature (arterioles and capillaries) includes a static gain G and first-order low-pass dynamics with the time constant τ. An increase in CBF causes vasoconstriction and, consequently, a decrease in compliance and an increase in resistance. The regulatory response is modeled by a sigmoidal static relationship with upper and lower levels to account for the limits of vasodilatation and vasoconstriction capacities.
The following differential equations describe the actions of autoregulation by means of first-order low-pass dynamics, time constants and gains x i is the state variable of cerebral autoregulation of the ith cerebral terminal artery that responds to alteration of CBF. G i is the static gain for the ith cerebral terminal artery; it is evaluated from the total static gain of autoregulation G (which valued will be defined in section 4.1) according to the flow distribution inside the brain. Q i is the time averaged flow and Q T i is the reference flow at the ith terminal artery over the period t À T, t ½ , where t is the current time and T is the cardiac cycle duration.
Once the control actions x i of each ith terminal artery are found solving the ordinary differential Equation (43), arterial compliances are changed through a sigmoidal relationship.
with upper and lower saturation levels. In this section, C i , R i and V i stand for compliance (C al , C cp ), resistance (R al , R cp ) or volume of arteriolar and capillaries compartments that are in the vascular beds linked to the ith terminal artery, k i is a constant parameter, inversely proportional to the central slope of the sigmoidal curve, C i and ΔC i are the central value and the amplitude of the sigmoidal curve. ΔC i depends on whether vasodilation or vasoconstriction is considered and it is chosen for each terminal artery as follows where sat 1 and sat 2 are constant parameters that define the upper and lower saturation levels of the sigmoidal curve. According to the literature, 85 the sigmoidal curve is not symmetrical; the increase in blood volume induced by vasodilation is higher than the blood volume decrease induced by vasoconstriction; therefore, two different values must be chosen for the parameter ΔC i . From (44) and (45), it follows that the upper and lower saturation levels of the sigmoidal curve are C i þ sat 2 2 and C i À sat 1 2 , respectively. The cerebrovascular control mechanisms affects not only the compliances, but also the arterial resistances. The variation of compliance in time changes the arterial volume. Since blood volume in a vessel, as a first approximation, is proportional to the radius squared (V / r 2 ), for a given vessel length, and vessel resistance is proportional to the inverse of radius to the power four (R / 1=r 4 ), volume varies according to the inverse square root of the resistance (V / 1= ffiffiffi R p ). Therefore, the following relationship is used to update the resistances of regulating arteries where V i is the mean volume of the ith arterial compartment (arterioles and capillaries) in the interval t À T, t ½ , while V T i is the mean baseline condition volume, R i is the current resistance of the arteriolar-capillaries compartment and R T i is the resistance under baseline conditions.

| Equations for cerebrospinal fluid and brain dynamics
A major aspect of the present paper is the coupling of the circulatory system to a more refined description of the CSF dynamics than in the first version of our model. 59,69 We depart from the model proposed by Linninger and collaborators. 62 The version of the Linninger model we present here differs somehow from the original version and includes nine CSF compartments: the lateral (LV and RV), the third (3 V) and the fourth (4 V) ventricles; the cerebral aqueduct (AoS), the CSAS, the SSAS and the bi-phasic brain parenchyma, comprising the left and right hemispheres. Each compartment is spatially idealized as a cylinder of length l and variable cross-sectional area A t ð Þ. Each brain parenchyma hemisphere is treated as an incompressible, deformable medium composed of two phases, the solid cell matrix, representing neurons, glial cells and axon fibers (70% of its total volume), and the extracellular fluid (remaining 30%); the model assumes that the volume of the solid matrix does not change and therefore brain parenchyma size changes depend only on the extracellular fluid content variations, that is, changes in porosity. All CSF compartments are interconnected and contain CSF assumed to be a Newtonian and incompressible liquid, with a constant viscosity of 0:001kg= ms ð Þ and a constant density of 998:2kg=m 3 .

62
In our model, we assume that CSF is secreted by the choroid plexuses, a highly vascularized region from the microcirculation of the anterior and posterior cerebral arteries, into the lateral ventricles. Also included is a constant CSF production from arterioles to the ventricles and the diffuse capillary production throughout the brain parenchyma to the ventricles. Then, CSF flows from the lateral ventricles to the third ventricle and, through the cerebral aqueduct, to the fourth ventricle. Then CSF is assumed to enter the subarachnoid space. Here, CSF is absorbed into the venous system through arachnoid villi into the superior sagittal sinus. Moreover, from the CSAS, CSF is displaced into another CSF compartment, namely the SSAS. Figure 3 shows the craniospinal compartments involved in the CSF system and their connectivity. The arrows indicate fluid exchange between compartments driven by pressure differences, while dashed arrows denote constant production of CSF, from the cerebral arterioles into the lateral ventricles q Al!LVs,const and from the brain capillaries into the extracellular space of the parenchyma. The type of arrows identifies whether the exchange of CSF between different compartments is unidirectional or bidirectional.
Flow of CSF in CSF compartments is governed by mass conservation and momentum balance. Such equations are accompanied by a tube law, relating deformation state and pressure, as for 1D vessels. Other equations are included in superior saggital sinus Brain: fluid part of brain parenchyma; Al: cerebral arterioles; Cp: cerebral capillaries. Solid double arrows denote fluid exchange between different compartments driven by pressure differences, while dashed arrows describe constant CSF production. The combination of a single dashed arrow and a solid double arrow between the brain parenchyma and the capillaries indicates that there are both fluid exchange driven by pressure differences (q in br driven by p Cp À p br ) and constant CSF production (q Cp!br,const ). Single solid arrow denotes CSF reabsorption into the venous system (SSS) the model to account for fluid exchange between different compartments of the CSF system or with the vasculature. The full CSF model is composed of 36 equations and 36 unknowns. The continuity equations read l br,L dA br,L dt ¼ q in br,L þ q Cp!br,const À q out br,L À q br!LV,const : Equations (47) to (55) are continuity equations that ensure that CSF is neither gained nor lost. Equations (47) to (53) refer to continuity equations for ventricles, Aqueduct of Sylvious, cranial and SSAS. Each equation guarantees that the volume change is given by the difference between the volumetric flow rates in and out of that compartment. Equations (54) and (55) are the continuity equations for the right and left fluid part of the brain parenchyma, respectively. In this case, the right-hand-side of the equations considers both the volumetric flow rate in and out of the compartment that is driven by pressure differences and the constant mass transfer. Flow into the the brain parenchyma is the sum of a constant CSF production from the brain capillaries into the extracellular space of the parenchyma, q Cp!br,const and the pressure driven seepage from the capillaries to the brain parenchyma, q in br . Flow exiting the brain parenchyma is the sum of a constant seepage from the extracellular space of the parenchyma into the ventricles, q br!LV,const , and a pressure driven exchange between brain parenchyma and lateral ventricles q out br,L . The momentum equations are effectively Darcy's law of flow and relate the pressure difference between two compartments to the volumetric flow q exchanged between them and a resistance to flow R. For the brain parenchyma compartments, there are two momentum equations, one refers to CSF exchange between the lateral ventricles and the extracellular fluid matrix of the brain, while the other one relates to the secretion of CSF from cerebral capillaries. As shown in Figure 3, these exchange pathways are bi-directional, depending on the hydrostatic pressure differences. When intracranial pressure exceeds the capillary pressure, reverse flow occurs, that is, in the present model capillaries are a pathway for CSF drainage. The equations for CSF flow are The notation for pressures is obvious for most compartments; for example p CSAS denotes pressure in the cerebral subarachnoid compartment. Just for clarity, in the last four equations p Cp,R is pressure in the capillary compartments of the right side of the brain, while p Cp,L is pressure in the capillary compartments of the left part of the brain; p br,R is pressure in the extracellular fluid part of the right brain parenchyma and p br,L is pressure in extracellular fluid part of the left brain parenchyma. We note that Equation (62) attempts to account for the interacting dynamics of two major CSF compartments. We are currently investigating these aspects as it has clearly some limitations, particularly regarding the omission of inertial terms that are known to influence the timing of flow exchange between CSAS and SSAS. 87 As already pointed out, the distensibility equations play the role of the tube law; they relate the internal pressure with the cross-sectional area of the compartment in a linear manner. They describe the dilation and compression of a compartment; if the pressure of the compartment exceeds the external pressure, the compartment is dilated with respect to the reference state ϵ0ϵ; in the opposite case, the compartment is compressed. For each compartment inside the cranial cavity, the external pressure is that of the brain parenchyma; for the SSAS, the external pressure is taken equal to zero. For a generic compartment z, the distensibility equation expresses pressure p z as a function of crosssectional area A z in the compartment and three additional parameters, namely an external pressure p ext,z , baseline cross-sectional area A 0 z and a coefficient E z denoting elastance, that is Therefore, for each specific compartment the equations are Additional equations connecting different compartments are required to complete the description of CSF flow. Specifically, for the right lateral ventricle RV, the amount of CSF that enters the right lateral ventricle RV is equal to the amount of CSF exiting the fluid part of the brain parenchyma plus the constant production rate from arterioles q Al!RV ,const plus the constant production rate from capillaries q Cp!br,const , namely Similarly for the left lateral ventricle LV, Then, CSF flows from lateral ventricles to the third ventricle from the third ventricle to the AoS from the AoS to the fourth ventricle and from the fourth ventricle to the CSAS From the cerebral subarachnoid space CSAS, CSF is temporarily displaced into the spinal cavity and reabsorbed into the superior saggital sinus (SSS) through the arachnoid granulations. 62,87 Reabsorption is represented by a mass transfer flux, which is a function of the pressure difference between the CSAS and CSAS and a reabsorption constant coefficient k We take the maximum value between zero and the mass transfer flux to enforce the unidirectional flow from the cranial SAS to the venous sinus. Previous MRI measurements 88 have shown that in a normal subject CSF reabsorption in the spinal cavity is negligible, since almost the total flow into the spinal SAS goes back into the cerebral SAS; for this reason, we set the CSF outflow from the SSAS equal to zero, Finally, the Monro-Kellie hypothesis is enforced: all compartments, except the SSAS, are enclosed inside the cranium and the volume of each cerebral hemisphere remains constant over time The volume of each compartment is evaluated as product of its length and its cross-sectional area. So far, we have presented the complete model for the circulatory system and the craniospinal dynamics. In the next section, we briefly present the numerical methods to solve all the equations.

| NUMERICAL METHODS
Much of the numerical methodology utilized here has already been applied in the original Müller-Toro model for the human circulatory system. 59 Therefore, in the present paper we limit ourselves to a succinct presentation of the numerical methods, with emphasis on the new aspects, along with relevant references.

| Overview of methods for PDE system
The hyperbolic system of blood flow Equations (17) is solved numerically using the ADER high-order numerical framework, originally reported in Reference 64; some of the numerous, subsequent ADER developments are found in References 89-96. An up to date review of ADER is found in 97 and references therein. An introductory presentation of the ADER methods is given in Chaps. 19 and 20 of Reference 98. A key ingredient of the ADER method is the solution of the generalized Riemann problem (GRP). In the present paper, we use the Dumbser-Enaux-Toro (DET) method, 95 extended to nonconservative systems in References 99 and 100. This is a locally implicit method that is able todeal with stiff source terms, in particular, the stiff source terms resulting from the hyperbolic approximation, 70,78,79 of the parabolic system incorporating the viscoelastic nature of the vessel wall mechanics. The initial conditions for the GRP are piece-wise smooth and result from a high-order nonlinear spatial reconstruction; here we use the WENO method, as presented in Reference 101. The DET method makes succesive use of a solver for the conventional piece-wise constant data Riemann problem. As already pointed out, the governing Equations (17) are written in non-conservative form, for which we deploy the path-conservative framework 102 to compute the numerical fluctuations, the analogues of the numerical fluxes for conservative methods. As is well known, the presence of source terms in the equations requires a suitable modification of the scheme so as to make it well-balanced. The present numerical scheme is indeed well-balanced, as proposed in Reference 60; this is a modification of the original Dumbser-Osher-Toro (DOT) 103,104 for constructing well-balanced fluctuations for a first-order nonoscillatory scheme in the framework of path-conservative schemes. Full details of the schemes for the present application are found in References 105, 60.
An enhancement to the numerical methodology for the present model of the human circulatory system is the LTS technique, 72 adapted from Reference 106. The technique is akin to the adaptive mesh refinement methodology, 107 in which different patches of the full domain are identified by a common spatial discretization length; then each patch is updated with a constant time step for that patch. Explicit methods, such as the present ADER method, are constrained to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition. This is normally implemented on the basis of the maximal wave speed present in the full domain and a minimal length emerging from the spatial discretization. Then, the time step for the updating of each cell is constant, though the local Courant numbers are variable. By allowing a variable time step for each cell, or patches of cells, according to its length and maximal wave speed one may reduce local numerical diffusion and increase efficiency. It is here where the LTS technique can be very useful. The procedure requires a careful global time matching at selected time levels. In our previous works on the human circulatory system 59,70 we have employed the common technique of a fixed time step per time level. In the present paper, we adopt the LTS technique, 72,106 in a vessel-wise manner. It seems as if the first example of a LTS solver applied to blood flow is in Reference 108. In the present paper, we follow, 72 allowing an explicit local time-stepping temporal discretization of the underlying finite-volume type ADER scheme. The net benefits of the LTS technique are two fold: reduced numerical dissipation and enhanced computational efficiency by orders of magnitude.
Considering the introduction of viscoelasticity and the LTS technique, it is essential to adopt a numerical methodology that considers the viscoelastic effects and solves the GRP at the coupling of several 1D vessels (junctions) accurately enough in order to preserve the formal accuracy of the scheme used for the resolution of conservation laws within the 1D domain. For the numerical treatment of the junctions, we follow. 59,72,109,110 Briefly, the adopted coupling strategy enforces mass conservation and total pressure continuity among vessels sharing a node, while generalized Riemann invariants are used to ensure that coupling conditions and states within 1D domains belong to a smooth solution of the original PDE system.
A succinct description of the numerical methodology for solving the system of PDEs (17) for the human circulatory system has been presented. Next, we describe the solution methods for the equations for the CSF system in the craniospinal space, along with the coupling procedure between the circulation and the CSF and brain dynamics.

| Numerical treatment of CSF equations and its coupling to the circulation
The description of CSF and brain dynamics leads to a system of 36 equations with 36 unknowns. The unknowns for each CSF compartment are pressure, cross-sectional area (which defines a volume since each compartment has an assigned length), inflow and efflux. The coupling between blood flow and CSF dynamics is explicit. The two systems are solved in a sequential manner. As we are using a LTS technique in a vessel-wise fashion, each 1D vessel is allowed to evolve in time according to a local time step given by its local stability criterion. All vessels have a common synchronization time defined by the prescribed maximum time step Δt max allowed by the LTS procedure. Therefore, the coupling between blood flow and the cranio-spinal systems is performed every synchronization time t n ¼ t 0 þ nΔt max , with t 0 the initial time. Figure 4 describes the coupling procedure from time t n to time t nþ1 . At the beginning of each time step, the vectors S n and B n are known. The vector S n represents the unknowns for the blood circulation system that includes area, flow and pressure in 1D vessels, as well as other 0D blood compartments. The vector B n represents the unknowns for the CSF and brain dynamics models. In the first step, the equations for the blood circulation models are solved. First, the system of ODEs for the cerebral autoregulation model are solved by the explicit Euler scheme, in order to find the new cerebral resistances and compliances, as described in section 2.2.5. Then we solve the system of partial differential equations for blood flow in 1D vessels and the 0D blood compartments for the heart and pulmonary circulation, microvasculature, SR and venous valves. Each jth vessel is evolved using the ADER scheme according to its local time step until it reaches the next time step t nþ1 ¼ t n þ Δt max . The 0D blood compartments are solved by an explicit Euler scheme and coupled to the 1D vessels. The external pressure for the intracranial 1D vessels and vascular beds, as well as for the SR models, is given by the mean pressure between the left and right sides of the brain parenchyma (p br,R and p br,L ) at time t n . Once the blood circulation equations have been solved, the cerebral capillary pressures (p Cp,R and p Cp,L Þ, the superior sagittal sinus pressure (p sinus ) and the intracranial blood volume (V blood,R , V blood,L ) are provided to the CSF models. This determines the CSF production, reabsorption rates and the blood volume inside the skull for the Monro-Kellie hypothesis. At this point, the system of differential and algebraic equations for the CSF and brain dynamics are solved by an implicit Euler scheme.
At the beginning of the simulation, given the initial conditions for the blood circulation models, the ODEs and systems of equations for the CSF and brain dynamics models are solved. In this way, the initial intracranial pressures are found and used as external pressure in the first time step update of the blood circulation.

| PARAMETRIZATION OF THE MODEL
In this section, we present all the parameters needed for the implementation of the global closed-loop model. The parameters assignment correspond to a healthy young male subject. We underline that the parameters are, unless otherwise specified, the ones proposed in References 59,69.

| Arteries and veins
The 1D vascular network contains 323 1D segments, of which 98 are arteries and 209 are veins, all linked by 143 junctions. Figure 5 shows a schematic representation of such networks, while Figure 6a,b provide a detailed description of head and neck arteries and veins. Compared to the network used in References 59,69, in the present article we propose a detailed description of the highly vascularized regions of the posterior part of the brain. Blood supply to the brainstem is crucial for the function of sensory and motor pathways, as the nerve connections of these systems from the main part of the brain to the rest of the body pass through it. More importantly, the brainstem plays a key role in maintaining cardiac and respiratory functions, such as heart rate and breathing. Three main arteries supply blood to the cerebellum: the superior cerebellar artery (SCA), the anterior inferior cerebellar artery (AICA) and the posterior inferior cerebellar artery (PICA). AICA (No. 304, 305, 306, 307) was previously included in our model for the ear circulation network 68 ; the other arteries are added here using data from the ADAN network 111 and from the literature. 112   I G U R E 4 Schematic representation of the coupling between blood circulation and CSF and brain dynamics models. S n and B n are the vectors of unknowns in the blood circulation models (1D vessels and other 0D blood compartments) and in the CSF and brain dynamics model at time t n . From B n , the pressures of the left and right fluid part of the brain parenchyma p n br,R and p n br,L at time t n are used to find the solution S nþ1 for the hemodynamics equations and the cerebral regulation. From S nþ1 , the superior sagittal sinus pressure p nþ1 sinus , the pressure of capillaries p nþ1 br,R , p nþ1 br,L and the total cerebral blood volume V nþ1 blood,R , V nþ1 blood,L are used to find the solution B nþ1 for the CSF and brain equations terminal part of the vertebral artery and the pontine arteries (No. 316, 317, 318, 319), lateral branches from the basilar artery that supply the pons. All these arteries end in the vascular beds of the cerebellum and brainstem.
According to experimental observations, [113][114][115] the mean value of blood flow to the cerebellum and the brainstem is about 10% of the total CBF. Following Reference 113, we estimate that flow to the cerebellum and to the brainstem are 1.01 and 0.13 ml/s, respectively. The posterior part of the brain is drained by the group of cerebellar veins, such as the superior cerebellar veins and the inferior cerebellar veins.  101,102). The addition of the posterior brain vasculature is essential Arterial and venous network composed of 114 arteries and 209 veins; numbers refer to those used in Tables A1 and A2 to explore some medical conditions that affect the brainstem and the cerebellum, such as the effect of vertebral artery hypoplasia in the ipsilateral PICA. 116 Moreover, in order to analyze better the implications of venous strictures in the pathophysiology of Ménière's disease, we redefine the ear vasculature previously included in References 68. The ear is mainly supplied by the labyrinthine artery, 117-119 which arises from the AICA, passes through the internal acustic meatus and then perfuses the inner ear. More details about the complete vessel network can be found in the Appendix. The coefficient K present in tube law (5) is obtained from the reference wave speed c 0 assumed for each vessel; in this work, we estimate its value, distinguishing arteries, veins and dural sinuses. For arteries, this wave speed is computed as proposed by Olufsen et al, 120 namely, where r 0 is the artery's radius at the reference configuration, k 1 , k 2 and k 3 are empirical constants and are taken to achieve normal wave speeds in the large vessels for a young adult human and a reasonable increase in smaller vessels. We set k 1 ¼ 3 Â 10 6 g=s 2 =cm, k 2 ¼ À7 cm -1 and k 3 ¼ 40 Â 10 4 g=s 2 =cm. Following Reference 59, the venous reference wave speed is estimated as follows where c 0,max ¼ 400 cm=s, c 0,min ¼ 150 cm=s and r max ¼ 0:8 cmr min ¼ 0:08 cm are the maximum and minimum vein radii in the network. Due to the physiological rigid nature of the dural sinuses, for them we set a constant reference wave speed equal to 1500 cm=s.  Tables A1 and A2 4.1.2 | Vascular beds As a consequence of the addition of vessels to the original network presented in References 59,69, the present work includes four new terminal models representing the microvasculature of inner ear and brainstem-cerebellum. Table 1 summarizes the simple connections between one artery and one vein while Figure 7 shows the complex vascular beds. Parameters corresponding to the microcirculation are not always retrievable from the literature; therefore, we use the strategy proposed by References 58 and 61. Total arterial resistance, arterial compliance and venous compliance are fixed up to a constant, according to literature data. Total arterial resistance is fixed to 0:85 mmHg=ml while arterial and venous compliances are 1:7 ml=mmHg and 146 ml=mmHg, respectively. For each terminal artery, we set a total arterial resistance R T (taken as the equivalent resistance of the circuit formed by distal arteries, arterioles and capillaries, using data from References 59,69, and then modified to match the fixed arterial resistance). Then, each resistance is distributed between vascular subsystems according to the general pressure distributions among vascular segments. Arterialcharacteristic impedance R da is set to be 15% of R T while the remaining part is partitioned between arterioles and capillaries as 70% for R al and 30% for R cp ; if an artery splits into two venous capacitors (as for the second artery in the Figure 2), we divide the resistance of the capillaries part according to the flow distribution into venous capacitors. In order to approximate the flow distribution of each venous capacitor, we use the Murray's law, that is, we assume that the flow rate of a vessel is proportional to the cube of its radius; therefore, the flow rate of each venous capacitor is proportional to the sum of the cube of the radii of 1D terminal veins draining from this capacitor. As for the total arterial resistance, the total arterial compliance is fixed to a value taken from Reference 61 and then it is distributed among the respective vascular bed compartments (C art ) according to Reference 58. Finally, C al is set equal to C art and C cp is set to be 0:15C art . In Table A1 we report the total arterial resistance R T and compliance C art for each terminal artery.
Concerning the venous compliances, we redistribute the venous compliance of each territory following Liang et al. 58 If a vascular bed is composed of two venous capacitors, we divide the venous compliance of the entire vascular beds according to flow divisions among them. The venous impedance R vn is taken as in References 59,69.  Figure 2. Arrows show if an artery is linked to one or both capacitors and the veins connected to each capacitor. Vessel numbers refer to those used in Tables A1 and A2  Tables 2 and 3. Parameters for the pulmonary circulation are the same as in the Müller-Toro model, previously taken from References 73, and reported in Table 4. Finally, concerning the pericardium parameters, we set V PC equal to 400 ml and Φ PC equal to 100 ml.

| Venous valves and Starling resistors
According to the vessel network extension, this work presents additional valves and SR, for a total of 17 valves and 21 resistors. Table 5

| Autoregulation
The cerebral autoregulation model works on 12 terminal cerebral arteries; the baseline haemodynamic parameters of these arteries are set from the periodic solution obtained for a baseline simulation and reported in Table 6. Other parameters of the model are taken from Ursino and Giannessi work 75 and adjusted to match our cerebral vessel network values (Table 7). For each terminal artery, the gain of autoregulation G aut,i is computed from the total G, according to the flow distribution inside the brain.
T A B L E 5 Location of venous valves (on the left) and Starling resistors (on the right). Vessels numbers refer to those used in

| CSF model parametrization
Parameters for the CSF model are based on Linninger et al. 62 Table 8 shows length and area at rest of the cylindrical volume representing each cerebral compartment, and different values of elastance taken from the literature. Table 9 reports flow resistance values; they account for the pressure drop in the fluid along the length of a compartment due to viscous forces and they are obtained from the dynamic fluid viscosity μ, the length of the compartment l and the square of the compartments' cross-sectional area. As already written in section 2.3, the CSF model adopted here accounts for constant production of CSF, from capillaries q Cp!br,const and from arterioles to lateral ventricles q Al!LV ,const . Almost two-thirds of the total CSF production takes place in the choroid plexus of the lateral ventricles; it was found clinically that this process is almost invariant to pressure changes suggesting an active transport process. 121 As in Reference 62, we fix a constant mass transfer independent from pressure equal to q Al!LVs,const ¼ 0:00583 ml/s. Moreover, there is CSF mass transfer from capillary beds into the brain parenchyma; the constant diffuse CSF production is set equal to q Cp!br,const ¼ 0:0005 ml/s. The active exchange between capillaries and brain parenchyma is governed by Equations (62) and (63), where CSF seepage is governed by pressure differences. CSF reabsorption is described in (80) by a mass transfer flux that is a function of the pressure difference between the subarachnoid space and the superior sagittal sinus and a reabsorption constant k. In this work, we use k ¼ 0:0027 mmHg/ml/s. We underline that variation of reabsorption constant could simulate pathological situations; for example, an increase of the reabsorption resistance may be due to inflammation of meninges while acute communicating hydrocephalus could be simulated by reducing k .

| SAMPLE NUMERICAL RESULTS AND VALIDATION
In this section, we present computational results obtained with the presented model in order to perform a comprehensive validation of the model's outputs. 1D domains are divided into computational cells with a reference length of Δx ¼ 1 cm, imposing a minimum of one computational cell in each vessel. Once that the mesh spacing of a vessel is fixed, the respective relaxation time ϵ for each vessel is computed in order to ensure that the accuracy criterion for the hyperbolic reformulation proposed in Reference 70 is satisfied. The CFL coefficient is set to CFL ¼ 0:9 according to the linear stability limit of ADER finite volume schemes for 1D problems. A maximum time step of Δt max ¼ 1 Â 10 À3 s is allowed. All computations are run using a second-order accurate version of the numerical scheme previously described. Other parameters linked to the blood characteristics are the blood viscosity taken as μ ¼ 0:045 P and the blood density  3 . The reference pressures taken as initial conditions are reported in Table 10. Given the closed-loop nature of our model, such pressures are important since they determine the periodic solution that the system will reach by defining the stretched blood volume. All the computational results shown in this section are obtained with simulations of 2000 cardiac cycles. A periodic state is reached after approximately 1600 cycles; compared to References 59, the time used to reach the periodicity of the simulation is higher due to coupling between two systems (blood and CSF) that have different time scales. Therefore, the verification of convergence of the solution is mainly based on the equality of the CSF production and CSF reabsorption, since CSF production involves the arterial pressure, and the CSF reabsorption rate is related to intracranial venous pressure. While future work will regard a more efficient treatment of coupling for the two systems under investigation, it is interesting to note that this difference in time scales poses severe constrains as to the mass conservation properties of the numerical schemes used to solve this problem. In fact, a discretization that is not able to enforce mass conservation at a discrete level would result in inability to reach a periodic solution due to mass conservation error accumulation. In particular, Figure 8 shows the computed waveforms along the aorta and major arteries of the lower limb. We can notice that, as the wave travels away from the heart toward the periphery, the systolic peak pressure increases according to physiological patterns, with a pulse pressure from 25 mmHg in ascending aorta to 45  literature data and results obtained with the previous version of the model 59 ; the corresponding bar plots are reported in Figure 9 (left). We note that when using the terminology literature data we always make reference to experimental data gathered in vivo and published by other research groups. Main cardiovascular indexes, such as systolic, diastolic and mean blood pressure and pulse pressure are computed and compared to literature data in Table 11. We conclude that waveform patterns in the arterial system are in accordance with general physiological data and that blood flow distribution along the aorta is reasonable. Concerning the venous circulation, Figure 10 shows the pressure and flow rate along the main systemic veins while Figure 9 (right) depicts a comparison of the predicted flow at different locations of the systemic venous circulation with literature data and results obtained in Reference 59. It is well known that blood flow in large to medium vessels is a convection-dominated process; therefore neglecting viscoelasticity of vessel walls in 1D models is often chosen as compromise. However, the viscoelastic behavior of arterial and venous walls is well-known. It has an impact on fundamental hemodynamic characteristics of the cardiovascular system and plays a determinant role in setting the functional level of the cardiovascular system under physiological and under pathological conditions. Previous works 39,134,135 have shown the benefits of considering viscoelastic properties of vessel walls in arterial circulation comparing model predictions and in vivo measurements of pressure and flow at different location. The effects of viscoelasticity become more significant in the periphery, especially on the flow wave. The same happens in the venous circulation. Zocalo et al 136 showed the importance of the dynamic process of veins walls to understand venous functioning under normal and pathological conditions; pressures and diameters of anterior cava, jugular and femoral veins from sheep were registered during cyclical volumepressure pulses. The vein viscosity was higher in the peripheral segments and this could be important in the response to acute overloads and in haemodynamic control. In this work, we introduce a viscoelastic tube law not only for the arterial tree but also in the venous circulation. Both arteries and veins are represented as a Voigt-type visco-elastic material (Equation 7). Figures 11 and 12 compare computed pressure and flow rate when vessels are represented with elastic and viscoelastic behavior of their walls. Figure 11 refers to thoracic aorta, femoral and carotid artery while Figure 12 depicts superior vena cava, femoral and jugular vein; the effect of the viscoelastic tube law is more evident in peripheral vessels (femoral artery and vein, carotid artery and jugular vein) with respect to central vessels. In particular, it can be seen that the solution obtained for viscoleastic vessels presents significantly less highfrequency components with respect to the solution obtained for elastic vessels. This fact is consistent with the dissipative capacity of real vessels, which, at least for physiological states, do not display pressure and flow waveforms with very high-frequency components. Figure 13 shows the computed pressures during a cardiac cycle for three different compartments. The pressure values display a physiologically behavior in all compartments; from arterioles to venules, the pressure slowly decreases. It ranges between 40 and 80 mmHg for arterioles, between 20 and 25 mmHg for capillaries and between 5 and 15 mmHg for venules. In particular, 13a refers to a simple connection in the kidney, 13b is a vascular bed in the abdominal region with four supplying arteries and one draining vein while 13c represents the microcirculation pressures in the left part of the posterior brain. Figure 14 shows the computed pressures and volumes for the four cardiac chamberswhile Figure 15 displays the pressure-volume relationship for the left and right ventricles. The heart model well represents the physiological   Figure 17 (left). The functioning of cerebral autoregulation is verified by changing arterial resistances of all but the cerebral arteries in order to cause a mean arterial pressure change, which would cause an increment in cerebral flow if peripheral cerebral resistance would not adapt. Figure 18 shows the computed autoregulation curve compared with literature data from Reference 75. The autoregulation curve relates mean arterial pressure (pressure of vessel No. 1) and the percentage change in CBF with respect to the baseline situation (evaluated as sum of mean flow over a cardiac cycle of internal carotid arteries and vertebral arteries).

| Validation of cerebral haemodynamics
Particular attention is given to the head and neck veins; in this case, PC-MRI flow measurements were available from Reference 59; these data were collected by the MR Research Facility at Wayne State University,  Figure 19). Figure 20 depicts computed intracranial pressure and pressures of a dural sinus and a cerebral vein; the effect of the SR is evident: pressure in cerebral veins is always higher than intracranial pressure while dural sinus pressure is governed by downstream conditions. Figure 21 shows the changes in volume of the main blood cerebral compartments: arteries, arterioles, capillaries, venules and veins.

| Validation for CSF and brain dynamics
Here, we analyze pressure, flow and volume of different CSF compartments and their interaction with arterial and venous blood within the brain; see  Physiological intracranial pressure (ICP) values have been investigated extensively, ranging from 7 to 15 mmHg for an adult in supine position. 142 Generally, ICP refers to the CSF pressure, regardless of where it is measured. In our simulation, slight differences are found between the mean pressure of the extracellular fluid part of the brain parenchyma and other cranial and spinal CSF compartments. In this context, we call ICP the pressure of CSF located in the brain parenchyma compartments; specifically, ICP is defined as 1 2 p br,R þ p br,L À Á , that is the average between the CSF pressure of the right and left extracellular fluid part of the brain parenchyma compartments. Figure 22 reports the pressures of the CSF compartments. The ventricular pressure ranges from 7.7 to 9.7 mmHg with a pulse pressure of about 2 mmHg. The cranial and SSASs pressures values lie between 8.1/8.2 and 9.3/9.2 mmHg and the pulse pressure is around 1.0 mmHg. Comparing the left and right hemisphere, there are no distinct differences concerning pressures. The brain parenchymal pressure varies from 7.8 to 9.8 mmHg, with a mean value of 8.42 mmHg on the right and left sides. Figure 23 shows an analysis of the intracranial pressure waveform according. [143][144][145][146] We can notice three physiological peaks: the first one, P1 or percussion wave, is the highest, followed by P2 or tidal wave and finally there is P3 or dicrotic wave, which appears after the dicrotic notch. The peaks come from the arterial pulse wave from the heartbeat on the brain which essentially floats in CSF; the ICP waveform can usually be seen in time-synchronized fashion relation to the arterial waveform. The three typical peaks of the intracranial pressure can be observed in all CSF compartments but in the spinal canal they are less pronounced. P1 is caused by the arrival of the arterial flow pulse. P2 is caused by the second arterial flow peak, which happens before the dicrotic notch, while the third peak, P3, is related to the increased flow occurring right after the dicrotic notch. Figure 23 compares the intracranial pressure waveform of the current CSF model based on Reference 62 and that evaluated considering the simple CSF model proposed by Reference 141 and adopted in Reference 69. While mean ICP depends on initial conditions for both models, in this second case, the waveform peaks are less well defined. Figure 24 reports changes in volume of the CSF model compartments. The major part of CSF is contained in the cranial and SSASs, about 30 and 90 ml respectively, where the changes over the cardiac cycle tale place. Figure 25 shows the time variation of the volumes occupied by different compartments of the cranio-spinal system, stressing the effect of the Monro-Kellie hypothesis. During systole, intracranial arterial blood increases and arterial pulsations are transmitted directly into the incompressible CSF filled SAS. This evokes a chain of events in the following temporal order: CSF shifts out of the skull into the spinal canal; venous blood from the sinuses flows out of the brain  Figure 26 underlines the relation between blood and CSF compartments' flow. In Figure 26 (left), arterial inflow, venous outflow through internal jugular veins, flow in AoS and inflow in SSAS are depicted during a cardiac cycle. In Figure 26 (right), the CSF and blood normalized flow of the same compartments is shown over a cardiac cycle. The lag in time between the systolic peaks is reported in Table  In order to show the applicability of the presented model to situations of clinical relevance, we carry out a preliminary analysis of the effects of transverse sinus stenosis and of strictures in the main extra-cranial venous outflow routes on the cerebral circulation, CSF and brain dynamics.

| Idiopathic intracranial hypertension patient with transverse sinus stenoses
The role of vascular abnormalities in the onset and course of neurological diseases has long been recognized. In the last decade, the influence of intra-and extra-cranial venous pathology as a trigger/cause of certain neurological disorders has gained attention. IIH is a neurological condition of unknown etiology, which requires prompt diagnosis and if left untreated can result in a rapidly progressive visual loss. As said before, how increased CSF pressure in the subarachnoid space influences intracranial arterial and venous fluid dynamics within the framework of the Monroe-Kellie hypothesis remains unclear. Bateman showed that in IIH, the total cerebral arterial inflow is increased by 21%. 69 On the venous front, there is evidence that almost 93% of patients with IIH harbor some degree of dural sinus stenosis. 139

| Problem setup
We investigate the impact of bilateral transverse sinus stenosis on cerebral venous flow and CSF dynamics, paying special attention to the role played by collateral flow pathways between deep cerebral vessels and extra-cranial regions. We   introduce stenoses to the reference venous network presented in Table 17 dividing the vessels affected by the stricture (the right and left transverse sinuses, No. 101 and 102, Figure 29) in two segments and putting between them a stenosis model. This model is based on References 152,153 and it evaluates the flow variation in time across the stenosis by means of a first-order ODE where Δp is the difference between the upstream and downstream pressures and L, R and B are defined by A 0 and r 0 are the mean reference area and radius of the vessels wherein the stenosis model is placed while l s and A s are the length and the minimum area of the stricture, taken equal to 1 cm and 10 %A 0 . Finally, k u ¼ 1:2, k t ¼ 1:52 and Þ are empirical coefficients. 152,153 6.1.2 | Comparison between healthy and IIH patient Figure 29 shows computed flow in the transverse sinuses for the healthy control (HC) and for the stenotic case (ST), along with MRI data for a healthy subject; there is a reduction of average flow rate of about 70%. Moreover, due to the stenosis, there is a re-distribution of flow in dural sinuses and an increase in dural pressure ( Figure 30). As a consequence, there is an increase in pressure in intracranial veins while the arterial flow and pressure are not modified. The venous hypertension due to transverse sinuses stenosis leads to decreased CSF reabsorption via arachnoid granulation which depends linearly on the pressure difference between intracranial pressure and superior sagittal sinus pressure. As a consequence, in order to maintain the balance between CSF generation and absorption, the intracranial pressure rises; moreover, following the Monro-Kellie hypothesis, a major amount of CSF is displaced into the spinal canal. For equal narrowing of the transverse sinuses, the severity of intracranial hypertension depends on the ability of the venous vascular network in developing collateral pathways to brain drainage. Occipital vein and sinus are the main collateral routes for flow limited by stenosis; the flow through these vessels increases significantly, in particular in the occipital vein (from 0.722 to 6.567 ml/s). We must consider that the venous network for head and neck used here represents a best-case scenario, with all possible collaterals present. If the collateral circulation is impaired, the consequences of a stenosis in the dural sinuses should be aggravated. In order to explore this hypothesis, we performed a simulation where blood is forced to flow exclusively through the dural sinuses. Results are shown in Figure 31. There we show the computed intracranial pressure for the HC and for the patient with stenotic transverse sinuses, when the collateral routes are activated and also when they are compromised. In the first case, there is an increased intracranial pressure from 8.42 to 9.83 mmHg in brain parenchyma and a comparable increase in other intracranial compartments; on the other hand, when the collateral routes are excluded, the averaged intracranial pressure over the cardiac cycle rises from 8.42 to 31.16 mmHg in the brain. Concerning the intracranial pressure waveform, in Figure 32 (left) we observe that the pulse amplitude between the healthy subject and the one with transverse sinus stenosis does not change, both in case of complete collateral circulation and without collaterals. According to Reference 143, if the ICP values were low, the pulse wave presents a descending saw-tooth appearance, with a clearly distinct P1 component; as the mean ICP rises, there is a progressive elevation in the magnitude of P2 and to a lesser extent of P3. Increase in the P2 component of the intracranial pressure wave is thought to represent decreased intracranial compliance. 140 144 From the numerical experiments shown here, we barely observe such changes in the ICP waveform. We attribute this fact to the linear character of intracranial compliant compartments. We consider it a limitation of our model in its current form, that will be addressed in future work. Table 13 reports the cardiac cycle-averaged cerebral arterial volume (1D arteries and arterioles), venous flow (capillaries, venules and 1D veins), cranial CSF and spinal CSF. Arterial volume is 28.98% of total blood volume; this data is in line with Reference 155 where the percentage of arterial blood volume with respect to total cerebral blood volume   was estimated to be approximately 20-30%. When the collateral routes are blocked, the spinal CSF volume increases from 90.42 to 100.90 ml without significant changes in cerebral venous blood and cranial CSF. This rather insensitive behavior of venous blood volume is thought to be related to the modeling approach adopted for representing SR, which resistance to flow is not influenced by transmural pressure, see section 7 for more details. Figure 32 (right) shows flow through the AoS. As for the case of ICP, we observe small changes in flow pulsatility for all cases considered. This observation is not in line with experimental observations of flow across the AoS in patients with venous obstructions to cerebral blood drainage, where visible increase in flow pulsatility is observed. 150,156,157 As already remarked for ICP, we believe that this behavior is due to the lack of non-linear compliance in intracranial compartments. Cerebral arterial inflow is similar between healthy and stenotic subjects with collaterals, 13.13 and 12.92 ml, showing that the cerebral autoregulation is maintaining the cerebral perfusion. The cerebral arterial pressure evaluated in the middle cerebral artery is 85.49 mmHg in healthy subject and 85.64 mmHg in the ST. However, the ratio between total jugular veins flow and arterial inflow is decreased from 93% in HC to 87% in subject with transverse sinuses stenosis; the outflow is deviated to external jugular veins and vertebral veins. When the collateral circulation is blocked in the ST, the arterial inflow is decreased to 10.11 ml and the ratio between total jugular veins flow and arterial inflow is 94%. The pressure in main cerebral arteries is increased by 3 mmHg. In this case, the CBF is significantly lower than that of HC, since the cerebral autoregulation is not able to fully compensate the drop in perfusion pressure (cerebral arterial pressure minus intracranial pressure) by reducing peripheral arterial resistances. According to Reference 63, standard MRI, MR venography and MR flow quantification studies revealed that the mean arterial inflow in stenotic patient with IIH is 21% above normal, but the SSS outflow was within the normal range; this means that the mean outflow as a percentage of the total inflow was reduced: this is evidence of collateral flow. In our mathematical model, we do not observe an increase in arterial inflow, both with and without collaterals circulation. However, when the collaterals are present, the percentage of venous outflow to arterial inflow is decreased, as in Reference 63. When the collateral circulation is reduced, the arterial inflow predicted by our model decreases. While this is the expected behavior if one analyses how our model is constructed, it contradicts observations in Reference 63. This disagreement could be due to the fact that we consider a sudden change from a baseline situation to a pathological one. There certainly are shortand/or long-term mechanisms not considered in our model that produce the above mentioned experimental observation of increase cerebral flow in IIH patients. This aspect will be investigated in future work. Table 14 shows the flow coming in and out of the left and right brain parenchyma (considering the pressure driven and the constant flow) and the brain porosity, evaluated as the ratio between the volume of the fluid part of the brain parenchyma and its total volume. The flow through the brain parenchyma decreases when there are stenotic transverse sinuses with respect to the HC; when the collateral circulation is blocked, this decrease reaches 5%, that is almost 16% if we consider the pressure driven seepage of extracellular fluid flow from capillaries into the brain. The brain porosity remains almost constant in all cases.

| Extracranial venous outflow strictures and their implication for Ménière's disease
Chronic cerebrospinal venous insufficiency (CCSVI) has been described as a chronic syndrome, characterized by extracranial venous malformations involving internal jugular veins, vertebral veins and the azygos vein. 12 The narrowing of these veins hampers the normal outflow from the brain, causing an impact on intracranial haemodynamics, as well as on CSF and brain dynamics. Some works on modeling this condition are available in the literature. In References 69 and 158, stenoses of the internal jugular veins were studied, while Reference 84 concerns stenotic venous valves. In Reference 68, the neck venous strictures are associated to Ménière's disease, a pathology of the inner ear. These works reveal that CCSVI leads to a significant increase in intracranial pressure; however, since in previous versions of our work, a single CSF compartmental model was used, no refined information about the CSF dynamics could be obtained.
Here, we investigate the impact on the CSF and brain dynamics resulting from the CCSVI condition, using a more sophisticated multicompartment model for CSF.

| Problem setup
We consider two different malformations of the extracranial venous vessels. The first one (case A) includes left and right stenotic internal jugular veins, symmetrically above the insertion of the middle thyroid vein, and also a stenosis in the azygos vein. In order to account for these strictures in the model, we introduce stenoses as represented in section 6.1.1 for vessels No. 224, 225, 244. The minimum area of the strictures is taken equal to 10% of the reference crosssectional area of the vessels where the stenosis model is placed ( Table 2). As for transverse sinuses stenosis, this area, as T A B L E 1 5 Cerebrospinal fluid exchange q in br,L , q in br,R , q out br,L , q out br,R and brain porosity Φ br (ratio between fluid part and total volume of brain parenchyma, considering a solid part of 980 ml): Comparison of results for a healthy control and subject with IJV stenoses (case 1) with and without collaterals, and subject with IJV stenotic valves (case 2)

HC
ST -case 1 -with Cols.  well as the reference cross-sectional area of the vessels, defines the parameters of the stenosis model in Equation (87). The second configuration (case B) takes into account stenotic valves, symmetric in both left and right internal jugular veins. The parameter M s of Equation (33) is taken equal to 0.25, causing an obstruction of 75% of the reference crosssectional area.

| Comparison between healthy and pathological patients
As expected, the CCSVI condition is related to a significant pressure drop across strictures. Pressure drops between the pre-and post-stenotic locations are negligible in the case of the HC while in case A it is almost 1.37 mmHg and in case B it is 2.25 mmHg. The pressure rise observed in the extra-cranial venous strictures is transmitted to the intracranial circulation. Figure 33 (left) depicts the computed cardiac-cycle averaged pressures in the main dural sinuses of the venous network. Moreover, Figure 33 (right) shows averaged pressures for veins of the right inner ear. As reported in Reference 68, the venous pressure in the ear veins is increased due to the extra-cranial venous stenoses; because of the functioning of the SR located in the ear circulation, this rise is not caused by the backward transmitted pressure waves from the obstructed sites, as occurs in dural sinuses, but it is due to the increased intracranial pressure. As shown in Figure 34, in all intracranial compartments the pressure is raised by 1 mmHg in case A and 2 mmHg in case B. By early animal experiments in Reference 159, it was proved that the subarachnoid space is linked to the endolymphatic space and the CSF pressure increase could be transmitted via the endolymphatic duct and sac to the inner ear, leading to the formation of the endolymphatic hydrops, one of the main anomalous conditions in Ménière's disease patients.
As in the previous pathological setting, the ability of developing collateral routes for brain drainage are important in determining the severity of intracranial hypertension. When there are stenoses in the internal jugular veins, the flow is redirected to the extracranial jugular and vertebral veins. The arterial inflow and the ratio between internal jugular vein flow and arterial inflow are 13.11 ml and 76.31% for case A, while for case B they are 13.10 ml and 62.64%. When the collateral routes are blocked, the internal jugular veins are the main drainage alternatives. We simulated such a situation where collateral pathways are blocked; in this case the arterial inflow is 12.52 ml and the ratio between internal jugular veins outflow and arterial inflow is 99.4%. The intracranial pressure increases by 5.5 mmHg (see Figure 34). This raise in pressure is of the same order of magnitude of the pressure drop through the stenosis (10.13 mmHg above the stenosis and 4.48 mmHg below in the right internal jugular vein); in fact, the increase in venous pressure is transmitted up the vessels into the superior sagittal sinus (from an averaged pressure of 5.31 mmHg in the HC to 10.97 mmHg), causing a proportional reduction in CSF absorption and then an increase in intracranial pressure until equilibrium between CSF generation and absorption is reestablished. 160 Table 15 reports data about the CSF flow in the brain parenchyma; due to the presence of stenosis, the pressure driven CSF seepage from the capillaries to the brain parenchyma is decreased. Figure 35 shows flow through the AoS. Deviations in terms of pulsatility are very small, even for the configuration featuring no collaterals. However, in Reference 150 changes in pulsatility were observed for patients with CCSVI. The fact that our model correctly captures the trend, that is, increased pulsatility, but fails to capture the magnitude of such increment, is certainly due to how well our model characterizes intracranial space compliance in the pathological case. Clearly some modeling aspects have to be improved, specifically the model of the SR (see discussion in next section), as well as mechanical properties of brain dynamics model compartments for patients with CCSVI.

| CONCLUDING REMARKS
In this paper, we have presented a global, closed-loop, multiscale model of the human circulation coupled with a multicompartmental model for the cerebrospinal-fluid dynamics. The model comprises 1D descriptions for medium to large blood vessels, arteries and veins, accounting for the viscoelastic property of the blood vessel wall. Lumped-parameter descriptions are used for other components of the full model that include the heart, the pulmonary circulation, the microvasculature, venous valves, SR and the dynamics of CSF in the craniospinal cavity, along with cerebral autoregulation. The present work departs from the Müller-Toro mathematical model 59,69 for the global systemic and pulmonary circulations in the entire human body. The main improvements with respect to the original Müller-Toro model and other works fall into three categories: (a) mathematical models for a better description of the physiology of the circulatory system; (b) the computational methods used to solve the governing equations and (c) the enlarged range of potential applications of the resulting model.
On the physiological aspects of the present paper, the improvements include the adoption of a viscoelastic tube law, not just for the arterial tree as in References 39,61,134 but also for the entire venous circulation. We find that computed pressures and flow waveforms for the arterial and venous circulations are more realistic for the viscoelastic tube laws, than for the purely elastic case. The heart model in the present paper is also an improvement over that in the original global model, 59,69 the cardiac valves are represented through a model based on Reference 41. A novelty of the present paper is the coupling of the blood circulation to a refined mathematical description of the CSF dynamics in the craniospinal cavity. 62 The CSF model comprises the cerebral ventricles, the AoS, the cranial and SSASs and the brain parenchyma. Other additions include the parametrization of the vascular beds, which together with the cerebral autoregulation model, is relevant when studying anatomical malformations of the cerebral circulation. The computational results for the arterial and venous circulation are compared with literature data, while MRI measurements are used for assessing our results for the cerebral venous circulation. Our results are in good agreement with literature and MRI data.
The physiology modeling improvements have resulted in new mathematical problems to be solved, notably, the viscoelastic nature of all major blood vessels. The associated parabolic system of equations has been approximated by a hyperbolic system with stiff source terms following a relaxation approach. 70,78,79 The resulting stiff system is solved numerically with the same high-order ADER-type numerical scheme, 71,95 as in the originalmodel. 59,69 An additional numerical improvement is the adoption of the LTS technique, 106 first introduced for blood flow in Reference 110 for solving a simplified 1D vessel network. This technique results in significant computational savings, which are more evident when coupling the blood circulation to the CSF and brain dynamics, as thesetwo systems have different temporal scales and the computational time needed to reach periodicity of the solution is considerably larger than the time scale of a cardiac cycle.
The model as presented is applicable to many pathophysiological conditions associated to the circulatory system, involving both the arterial and the venous systems. In the present paper, we have placed considerable emphasis on the CNS fluids in the craniospinal cavity, in which blood (arteries, microvasculature and veins) interact with the CNS fluids. To illustrate the applicability of the present model that couples the blood circulation and the CSF dynamics in a holistic setting, we have presented two specific medical applications, namely IIH as associated to transverse sinus stenoses, and Ménière's disease as associated to extracranial venous outflow strictures. Our results reveal that obstructions in the cerebral venous network lead to intracranial hypertension and disruption of the fluid dynamics in the entire craniospinal cavity. The severity of the consequences of intracranial or extracranial venous outflow obstructions depends on the balance between CSF generation and absorption, the displacement of CSF into the spinal cord and the ability of the venous network in developing collateral routes to respond to the venous outflow obstructions. These findings are relevant to the study of a very important function of CNS fluids, namelythe clearance of brain metabolic waste and neurotoxins from the CNS. Impairment of the cerebral venous drainage will directly disrupt this clearance function. Indirectly, CSF absorption into the venous system will be hampered, due to venous hypertension, leading to decreased CSF turnover, which will also affect the clearance function.
In spite of the progress reported in this paper, there are several limitations to be addressed in future developments. One limitation is the description of SR, which are major determinants of cerebral venous dynamics; blood flow through these compliant vessels is controlled by sphincter-like structures, which regulate discharge into the main dural sinuses. Another limitation is the absence of a model for solute transport. This limitation prevent us at present from properly describing, via Starling forces, 2 the transport of fluid and selected solutes across the blood-brain barrier (BBB), for example. Overcoming this limitation will be crucial for tackling the brain waste clearance function of the CNS, alluded to earlier. Mathematical modeling steps in this direction are outlined in References 161,162. Another limitation involves the linear distenibility equations of the CSF model. Following Reference 62, these equations link the internal pressure with the cross-sectional area of the compartment in a linear manner. A simple linear pressure-volume relationship is acceptable in the physiological pressure range, but could give under-and/or over-estimation of pressure changes in case of large volume changes, especially when addressing pathological conditions. Future work will address the non-linear behavior of the pressure-volume relationship in the CSF compartments. Another potential improvement concerns the representation of CSF in the spinal canal, which at present consists of a single 0D model; a possible improvement could be the adoption of a 1D model, as proposed in Reference 148, which is based on two coaxial compliant tubes representing the spinal cord and the CSF between the cord and the dura. Potentially, such representation admitting spatial variations could provide the bases for adding new potential routes for CSF reabsorption.
Addressing the modeling of the microcirculation is a challenging task. Here, the microcirculation was simplified as a lumped resistor capacitor system. While this simplification gave acceptable system-wide predictions, it is not able to account for biphasic blood flow phenomena 163 and network effects 164 that occur in the microcirculation. Significant progress has been made to develop realistic microvascular networks models 135,165,166 which could be integrated with the proposed system models in future work.
On balance, the mathematical model presented here is a significant improvement of the original model 59,69 published 7 years ago, represents the current state of the art and could provide the bases for realistic applications that require the representation of all extracellular body fluids in a holistic setting, along with regulatory processes. APPENDIX Tables A1 and A2 report geometrical data for arteries and veins networks. The first columns refer to number, name, length and radii of each vessel while the last one show the location of the vessel in the human body (using the same location codes of Reference 59, 1 = Dural sinuses, 2 = Extracranial, 3 = Neck, 4 = Thorax, 5 = Abdomen, 6 = Upper limbs, 7 = Lower limbs, 8 = Pelvis, 9 = Intracranial). Geometrical data are the same of References 59,74,75 while, for vessels added or changed in this work, the geometry is estimated from literature data. 111,167 T A B L E A 1 Geometrical and mechanical parameters of the arterial network. R T is given in mmHg s/ml while C art in ml (s mmHg) À1

No.
Vessel name L (cm) r 0 (cm) r 1 (cm) R T C art Loc