A comparison of two mathematical models of the cerebrospinal fluid dynamics

In this paper we provide the numerical simulations of two cerebrospinal fluid dynamics models by comparing our results with the real data available in literature (see Section 4). The models describe different processes in the cerebrospinal fluid dynamics: the cerebrospinal flow in the ventricles of the brain and the reabsorption of the fluid. In the appendix we show in detail the mathematical analysis of both models and we identify the set of initial conditions for which the solutions of the systems of equations do not exhibit blow up. We investigate step by step the accuracy of these theoretical outcomes with respect to the real cerebrospinal physiology and dynamics. The plan of the paper is provided in Section 1.5.


INTRODUCTION
The present paper introduces a mathematical analysis of two simplified models describing the cerebrospinal fluid dynamics. In order to clarify the framework to the reader, we shall provide the preliminary notions of anatomy and physiology on the processes underlying the craniospinal dynamics modelized by the systems of equations we are going to study. Indeed, the cerebrospinal fluid (CSF) circulation network is an intricate system which surrounds the central nervous system and is incorporated into it. It has been the subject of debate since its first description in the 18 th century, emphasized by the complex vascular network of the choroid plexus that has been conventionally considered the most important structure in the production of CSF through a variety of transporters and active channels. Recent outcomes in scientific methodology and outstanding improvements in the experimental tests represent fundamental tools in order to clarify the mechanisms of the CSF circulation mechanisms.

Cerebrospinal fluid
Cerebrospinal fluid is a clear fluid mainly composed of water (99%), with the remaining 1% accounted for by proteins, ions, neurotransmitters, and glucose ( [5]). It fills the surrounding spaces of the central nervous systems (CNS) of mammals and is a versatile wonder which sustains continuously the entire nervous system through the life of the organism. In the average adult man, the amount of CSF circulating at any given time is about 150 ml: in the ventricular compartments it is possible to detect the 17% of the total fluid volume, the rest perfuses the cisterns and the subarachnoid space (SAS). CSF is produced at a rate of about 0, 3 − 0, 4 ml/ min, i. e. 430 − 530 ml per day ( [3]). Classically the CSF flow has been thought as driven by the forces generated by heartbeats and pulmonary respiration. The original theory about the CSF production establishes that the 75% of all the fluid is produced by the choroid plexus epithelium, while the remaining 25% arises from other CNS structures such as the ependymal wall, the cerebral parenchyma and interstitial fluid ( [20]). Recently, however, there have been criticisms around the design of experiments on the choroid plexuses, questioning the thruthfulness of what we know about CSF.

Cerebrospinal fluid production
The secretion of CSF from the choroid plexus takes place as a process characterized by two different steps ( [2]). In the first stage, the plasma is passively filtered through the capillary endothelium into the choroidal interstitial space by means of the osmotic pressure gradient between the two surfaces. The ultrafiltrate fluid undergoes progressively active transport via the choroidal epithelium in the ventricular compartments. An alternative hypothesis on the production of CSF supported by Oresković and Klarica ([35]) clarifies new data about the choroid plexus as the main compartment in which occurs the CSF secretion. The authors state that there are not any experiments that prove the ability of the choroid plexus of producing the expected volume of CSF. The new theory proposed describes the CSF formation as an active process that is not influenced by intracranial pressure. In balanced physiological conditions, the rate of CSF production must be equal to the rate of absorption. The authors postulate that this could extend to the flow rate, since production and absorption take place in different compartments of the circulation system. Hence, it is immediate to assert that CSF secretion is the driving force of flow and circulation if there is a constant volume of CSF. The new theory moves from a more systematic approach, focusing on the perivascular spaces, which lie between the point where the cerebral vasculature descends from the SAS in the CNS through the pia mater by perforations ( [6]). It is at this junction level that production and reabsorption of both interstitial and CSFs occur, due to the differences in hydrostatic and osmotic pressure between the CSF circulation system and the surrounding tissue. This means that CSF is continuously produced by means of the bloodstream and not in isolated organs involved in secretion, and any variation in CSF volume is influenced by CSF osmolarity ( [35]). Although there exists evidence which supports these statements on the mixing and production of CSF, there is also extensive literature describing CSF ebbs and flows, and net flow, as well ( [41]). Indeed, according to Spector ([41]), the proposed active process in CSF formation and absorption by the entire CSF circulation system, ignores the mixing of CSF which is corroborated by the mobile cilia on the ependymal wall and by the transport of growth factors towards certain regions of the brain.

CSF circulation and absorption
After secretion, CSF generally (see Fig. 1) flows through the ventricular system, and the circulation is partially combined with the ciliated ependyma that beat in synchrony ( [39]). In general it is assumed that the CSF net flow perfuses the ventricular system, by starting from the lateral ventricles ( [20]). From the lateral ventricles, the CSF flows through the left and right foramen of Monro to the third ventricle. Then, it flows downwards along the aqueduct of Sylvius in the fourth ventricle. From the fourth ventricle, the CSF is allowed to exit laterally through the foramen of Lushka, or medially through the foramen of Magendie into the SAS. The displacement through the foramen of Magendie determines the filling of the SAS. The CSF that flows along the foramen of Lushka, circulates in the subarachnoid space of the cisterns and the subarachnoid space above the cerebral cortex. Finally, the cerebrospinal fluid from the subarachnoid space is reabsorbed by arachnoid granulations which are a sort of outpouchings in the superior sagittal sinus (SSS). Arachnoid granulations act as a pathway for the reabsorption of CSF into the bloodstream through a pressure-dependent gradient ( [2], [10]). Arachnoid granulations appear as outpouchings in the SSS due to pressure in the SAS which is greater than the pressure detected in the venous sinus (nevertheless, in vivo examination of arachnoid granulations would reveal the inverse). As well as new theories on CSF formation, recently new hypothesis about CSF absorption theories have been moved ( [24] and references therein). Studies on rabbit and ovine models have revealed that CSF can also be significantly absorbed by means of cervical lymphatics ( [2]). CSF not reabsorbed through arachnoid granulations can reach cervical lymphatics throughout two potential pathways. The first is the SAS corresponding to the outlet of the cranial nerves ( [2]). This provides a direct path in which CSF can be shifted from cisterns to extracranial lymphatics. The second pathway across which CSF is allowed to reach the lymphatic system is along the perivascular space of the arteries and veins that penetrate the brain parenchyma ( [7]). The perivascular space (Virchow-Robin space) is the potential space which surrounds the arteries and veins of the cerebral parenchyma, that can assume different sizes depending on the pathology. When the CSF is not absorbed through the classical route, it can flow into the perivascular space or may be moved into the interstitial fluid (ISF) which is a compartment with bidirectional flow to the perivascular space and SAS. If CSF fills the ISF, finally it will be reabsorbed into the bloodstream, flows into the Virchow-Robin space or passes back to the subarachnoid space. From the perivascular space, the CSF can flow back into the subarachnoid space or be absorbed by the cervical lymphatic vessels dependent on the forces exerted by the heartbeats and pulmonary respiration. Moreover, there have been conducted studies concerning the reabsorption of CSF into the dural venous plexus. At birth, the arachnoid granulations are not fully developed, and the absorption of CSF is due to the venous plexus of the inner surface of the dura which is more robust in children ( [29]). Although not very wide in adult subjects, the dural venous plexus is still believed to have a role in reabsorption. Nevertheless, the exact mechanism of CSF absorption has not been completely clarified yet ( [29]).

CSF and cerebral vasculature
Cerebral vasculature has the important role of providing blood supply to the brain. As the brain is one of the most important organs of the human body, an extensive network of arteries and veins assolves the task of providing oxygenated blood to the organ and to drain deoxygenated blood from it. The arterial system supplying the brain presents many collateral vessel that  insert redundancies into the system. The cerebral venous compartment is believed to be the main responsible of the intracranial compliance (defined as the ratio of volume and pressure change). In fact, the possibility of emptying part of the venous compartment to make up for an increased volume of the other compartment is a key phenomenon in the dynamic equilibrium of intracranial pressure (ICP). Moreover, an increased ICP can cause an unstable collapse of veins and so induce a change of resistance that may alter the cerebral blood flow (CBF). Such modification of the resistance may trigger autoregulation mechanisms, which act to ensure a constant blood flow, and it may be responsible for the amplification of the CSF pulsewave observed during infusion studies. However, no definite conclusions have been drawn on the latter point.
One of the main feature of the cerebral vasculature involves the coupling with the cerebrospinal fluid circulation. The CSF bulk flow associated to the production and reabsorption is small compared to its pulsatile component. The oscillations of pulsatile CSF are assumed to be managed by systolic vasculature dilatation followed by diastolic contraction. During the normal cardiac cycle, approximately 750 ml of blood are pushed into the head every minute. The increase in systolic blood pressure inflates arterial blood vessels, therefore an increase in the volume of brain blood during systole is detected. Since the cranial vault is characterized by a rigid skull in adults, the vascular expansion of the major arteries which pass through the spaces perfused by the cerebrospinal fluid activates the CSF motion. MRI techniques show evidences that the total volume of cerebral blood expands and comprises itself in every cardiac cycle of about 1 − 2 ml, the same volumetric quantity corresponding to the exchange of CSF between the cranial and spinal SAS ( [24] and references therein). The expansion of the vascular volume could be transmitted from the cortical surface through brain tissue, whose compression causes the contraction of the ventricular space. Otherwise, the movement of the ventricular wall could derive from inside the ventricles by systolic expansion of the choroidal arteries, such that the ventricular walls pulsate against the periventricular ependymal layer. Ventricular dilatation due to the expansion of the choroid was proposed in the theoretical model of Linninger et al. (see [25]).
Alternatively, Buishas et al. ([4]) proposed a model of water transport through the parenchyma from the microcirculation as driven by Starling forces. This model investigates the effect of osmotic pressure on water transport between the cerebral vasculature, the extracellular space, the perivascular space and the CSF and predicts the effects the osmolarity of ECS, blood, and CSF on water flux in the brain, establishing a link between osmotic imbalances and pathological conditions such as hydrocephalus and edema.

Key points of the paper
In the light of the above, the main issues in analyzing the intracranial dynamics is represented by the particular structure of the human brain that includes very complex elements and the occurrence of many different phenomena: • The flow of the CSF fluid throughout the CSF compartments.
• The mechanical interaction between the fluid and the brain. The brain parenchyma and the CSF compartments interact by exerting reciprocal stress on each other.
• Coupling with the circulatory-system. Models of the global circulatory tree may be necessary to provide pressure values in the brain with enough accuracy and to impose the Monro-Kellie doctrine on arterial and venous system, CSF and brain parenchyma. The Monro-Kellie doctrine, may be expressed mathematically through the Linninger reformulation (see [26]): The subscript, exf, refers to the extracellular fluid flow inside brain. Moreover, the volume of the brain parenchyma in (1.1) 2 is the sum of the volume of extracellular fluid, V exf , and the volume of the solid part of the brain, V solidBrain . The first quantity in (1.1) 1 is characterized by the sum of all the volumes of the vascular system components (i. e. arteries, arterioles, capillaries, veinules, veins and venous sinus), while the second quantity is the sum of the volumes related to the compartments involved in the ventricular system (i. e. lateral ventricles, third and fourth ventricles, cranial subarachnoid space).
• Production and reabsorption laws. In order to describe the CSF system, suitable laws are needed to model the production of fluid inside the ventricles and the absorption of the fluid in the sagittal sinus.
Lumped or compartmental models are particularly appreciated in the description of complex multicompartmental systems, in fact they are usually simple to develop as prototype and fast to solve. Marmarou pioneered the field of mathematical modeling of intracranial dynamics by lumped models in 1973. In his thesis [31], he introduced the concept of PVI (pressure-volume index), studied the intracranial pressurevolume relation, defined how to assess PVI and CSF reabsorption resistance by infusion studies and developed the first lumped model for intracranial dynamics focused on CSF [33]. His aim was to predict changes in ICP when the content of CSF inside the CSF system changes.
Brain parenchyma has been modelled through three dimensional poroelastic models ( [9], [18], [8]). Brain fluid systems have been modelled through lumped parameter models ( [43], [44], [15], [16]) and through multi-scale models ( [34]). Ursino ([43], [44]) proposed a mathematical model of the human intracranial hydrodynamics. The group of Linninger ([25]) created a dynamic model of the ventricular system to test the hypothesis that choroid plexus expansion drives CSF flow in this system. Later, they proposed a mathematical model of blood, cerebrospinal fluid and brain dynamics, including the Monro-Kellie doctrine ( [26]). The same group presented a mathematical model of the intracranial fluid dynamics based on the Bulat-Klarica-Orešković hypothesis ( [36], [28]). Gehlen et al. ( [16]) studied the effect of postural changes in the CSF dynamics through a lumped-parameter model of the CSF system and major compartments of the cardiovascular system. Our aim is to compare from a mathematical point of view two models which describe the simplified dynamics of the cerebrospinal fluid in the parenchyma. Both models are obtained by considering as a starting point the CSF model studied in [14], where the authors provide a comprehensive mathematical analysis of the system of equations derived from the CSF compartmental model introduced by Linninger et al. in [25]. The decision to adopt non recent CSF models ( [25], [33]) as starting point for our research, was dictated by the fact that the analysis that we propose below represents a first attempt to work from a purely mathematical point of view with models that describe cerebrospinal fluid dynamics. Because of the complexity of the mechanisms involved in CSF production, circulation and absorption, we decide to neglect in this first step some fundamental contributions of the brain vasculature or the parenchyma. Therefore, the aim of this paper is to lay the foundations for a much more detailed future study, based on models that treat in detail the interactions that emerge in the intracranial pattern among the vascular system, the CSF compartments and the brain parenchyma, and are able to predict blood alterations or potential cerebral water content variation. In this paper we want to improve the simplified model analyzed in [14] in the following way: • the cross sectional area of the CSF compartments, which is affected by the CSF flow and by the intracranial pressure in the real physiology, is not assumed anymore as a constant but it is free to change during the cerebral process modelized; • we add to the Linninger model (see [25]) an equation, developed by Marmarou et al. in [33], which describes the evolution of the intracranial pressure; • we introduce a CSF absorption rate in order to study also the equations that rule the final part of the CSF process which involves the venous superior sagittal sinus. We will define this new quantity in terms of intracranial pressure.
A model which reflects a perfect mixture of these features is the one proposed by Linninger et al. in 2009 (see [26]) that provides a complete representation of the CSF pathways interacting with the vascular system and the porous parenchyma.

Definition of equations symbols
The quantities involved in the models are the following (see Fig. 3).

A(t, z)
cross section of the CSF compartments axial CSF flow velocity [m/s] P (t, z) CSF pressure in ventricles and SAS (ICP) N/m 2 Table 1: CSF dynamics quantities involved in Models A1 and A2.
The time t is such that t ∈ [0, T 0 ] and z ∈ [0, L] is the axial coordinate. The physical constants that will be adopted in the present paper are listed in Table 2 (see also Fig. 3). The forcing function given by represents the choroid plexus periodic motion which follows the cardiac cycle.

Description of the models
This work focuses on two systems of equation given by a CSF compartmental model (see [25]) where the CSF sections are described by a cylindrical discretization with axial symmetry and radial displacements. The foramina are treated as elastic pipes. Since the CSF flow is basically laminar based on its very low Reynolds numbers (Re < 100), the flow friction was expressed as a function of the cross-sectional area. The latter, denoted by A(t, z), is considered as a function depending on time and space, and this represents the first significant improvement on the analysis developed in [14], where it is assumed to be constant. As already explained, the models we will analyze in the present paper derive from the compartmental model presented by Linninger et al. in 2005 ([25]) which the authors improved in 2009 ( [26]), by proposing a comprehensive model of human intracranial dynamics where cerebral blood, CSF and brain parenchyma as well as the spinal canal were included. In order to manage properly our mathematical analysis, in this first approach to the CSF dynamics, we take into account the first model mentioned since it represents a good starting point for our analysis at this early stage. For the same reason, to define the behaviour of the intracranial pressure, we consider the model proposed by Marmarou et. al in [33], described by the following differential equation that rules the CSF hydrodynamic without incorporating the interactions with brain vasculature and porous parenchyma, Therefore, the first simplified CSF model we study in this paper reads as follows, (1. 2) The previous model is designed as follows.
• The first equation of Model A1 is the continuity equation, consistent with the assumption of incompressible CSF. It defines the conservation of fluid mass through the motion.
• The second equation represents the distensibility balance, indeed the stresses, strains and displacements detectable in parenchyma and in the CSF compartments tissue, is described by applying the elastodynamics laws (see [25]).
• The third equation is the axial momentum balance where βu(t, z) is the Poiseuille term (see [25]).
• The last equation is the intracranial pressure equation already defined ( [33]). All the variables and constants that appear in the equations are defined in detail in Section 1.3.
As previously explained, we know that under normal physiological conditions, the rate of CSF formation (Q p ) is balanced by an equal rate of absorption (Q a ). This condition of equilibrium results in no increase or decrease in the amount of volume stored and the initial resting volume as well as the CSF pressure (P (t, z)) are maintained at a constant level. The rate of outflow (Q a ) is given by the gradient of pressure between CSF space and the venous system of the dural sinus (we assume equal to P ) divided by the resistance to absorption (R) (see [33]): This last issue allows us to introduce the second model we want to study and compare with (1.2).
Model A2 (1. 3) The previous model, along the lines of Model A1 is defined as follows.
• The first equation is the continuity equation where we are taking into account a simplified reabsorption term, Q a .
• The fourth equation is the momentum balance.
In order to simplify the models, we reduce the order of the second equation of (1.2) and (1.3). Furthermore, in the system (1.3) we plug the equation Therefore we reformulate the models in the following way In both models we impose for the unknown u(t, z) the following initial and boundary conditions respectively where s > 5 2 and and the following compatibility conditions for which we assume Similarly for the unknown η we assume and the following compatibility conditions Concerning the no-slip conditions, we can assume, without loss of generality, the following boundary conditions Moreover for the extra unknown ζ we impose as initial condition, and the following compatibility conditions By considering the boundary conditions (1.12) of the unknown η, we get and In both systems we have a homogeneous Riccati equation that describes the evolution of the intracranial pressure. For the latter we choose the following initial condition (1.18) By using classical ODE methods we can find the general solution of the Riccati equation It represents the variation in time of the pressure during the entire CSF process from the choroid plexus to the subarachnoid space.
Remark 1.1. It is important to observe that the measurement of the CSF pressure gradients can be affected by the body position. Since the common techniques are rather invasive, we have only limited data related to the upright posture. Indeed, CSF pressure is usually measured while a person is lying in a horizontal supine position. Normal CSF pressure values, in that case, are around 7 − 15 mmHg, and the pressure is the same along the SAS and inside the cranial vault (see [13]). It is known that the change in body position (from horizontal to upright, head up or sitting position) is followed by a drop in intracranial pressure to the subatmospheric value, and it can be detected a new pressure gradient along the craniospinal axis. In [1] and [21], the authors carryed out analysis using data previously collected from healthy subjects scanned in supine and sitting positions and, after modifications for the hydrostatic component, they detected ICP values at the central cranial location of all subjects were negative, with an average value of −3, 4 mmHg.
By taking into account the previous considerations, it is clear that there are two different choices of the initial condition for the intracranial pressure and they depend on which situation we want to investigate. More precisely: • if we apply Models A1 and A2 to a subject that is lying supine, we are forced to fix an initial datum b(z) > 0. In this case, the solution (1.19) blows up at the finite time 20) and, in order to study local existence in time of the solutions of the systems (1.4) and (1.5), henceforth we will consider the time t ∈ [0, We point out that this is not a big restriction on the possible values of t, since if we plug the correct values of the physical parameters in (1.20) and we take b(z) = 7 mmHg = 932, 54 Pa, which is the average intracranial pressure, we can observe that the blow up time T 0 corresponds to an age of 124 years.
• If we consider a subject in the upright position, we need to choose the initial datum b(z) < 0. This implies there is no blow up at finite time and no restriction on the time interval is required.
The last unknown to analyze is the generic axial section which is denoted by A(t, z). We impose the following initial condition and we observe that the general solution of the third equation in both systems (1.4) and (1.5) is for Model A1, and for Model A2.
Our approach in the analysis of Models A1 and A2, supplemented by the initial and boundary conditions previously defined is as follows.
• We implement a procedure of successive approximations to the systems A1 and A2 and we analyze the regularity of the sequences u n (t, z), P n (t, z), η n (t, z), ζ n (t, z) and A n (t, z) in order to prove the existence of an approximating solution. The sequence u n requires a different analysis due to the nonlinearity of the velocity flow equation.
• We show the convergence of the sequence u n for which we resort to high order energy estimates. For the other sequences, in order to obtain convergence to classic solutions, we need to apply first the Ascoli's Theorem A.1 which guarantees that the set provided the sequences are equicontinuous and equibounded.
• We pass into the limit in the approximating systems and we complete the proof of the local existence and uniqueness of classical solutions for Models A1 and A2.
• We investigate under which conditions it is possible to obtain global existence. Therefore we study the equation (1.4) 4 with the standard characteristic method and we obtain a nonhomogeneous Riccati equation (NRE) concerning which we have to find a general solution. In constructing the particular real solution to the NRE, we are forced to employ a condition on the pressure; this result allows us to prove global existence in time and uniqueness of solutions for the systems A1 and A2 under precise restrictions on the initial data.
• In order to prove that the analysis of Model A1 and A2 is reliable, we implement proper numerical simulations by fixing first initial data which fulfill the condition imposed for the global existence of the solutions and then we choose initial data which violate it.

Plan of the paper
The plan of the present manuscript is as follows. In Section 2 we collect all the definitions and the technical results we are going to use through the paper and we state our main results which are the local existence Theorems 2.3 and 2.4 and the global existence Theorem 2.5. In Section 3, we perform the numerical simulations that validate the results on the global existence of the solutions obtained for Models A1 and A2 and we verify for each simulation how the behavior of the quantities involved in both models matches the real CSF dynamics. Section 4 is entirely devoted to the comparison of our theoretical results with respect to the real data available in literature. Finally, in the Appendix we provide the detailed mathematical analysis of both models in order to prove the local existence Theorems 2.3 and 2.4 and the global existence of solutions (Theorem 2.5) for which we are able to show the set of initial data that avoids blow up formation.

PRELIMINARIES
In this section we are going to fix the notations used in the paper, we recall the main tools we need to study our problem and we state our main result.
We will use also the following Sobolev interpolation theorem  Definition 2.2. Let X be a compact metric space, Y a Banach space, and C(X, Y ) the Banach space of continuous functions from X to Y with the sup norm. A subset Ω of C(X, Y ) is equicontinuous if, for every > 0, there is a δ > 0 such that, for each f ∈ Ω, d(x, y) < δ implies where d is the metric in X. Ω is uniformly bounded or equibounded if there exists M ∈ R + such that

Statement of the main results
In this section we state the main results of this paper. The local existence and uniqueness of a solution for Models A1 and A2 is stated in the following theorems.
The proofs of Theorems 2.3, 2.4 and 2.5 are treated in detail in the Appendix 5.

NUMERICAL SIMULATIONS
In the present section, direct numerical simulations are carried out to investigate the accuracy and completeness of the theorems we proved. In this study, the CSF is treated as a Newtonian fluid with a dynamic viscosity and density equal to 1003 · 10 −3 kg/(m · s) and 998, 2 kg/m 3 , respectively [27]. The brain tissue is considered as a linear viscoelastic material with the storage and loss moduli of 2038 and 1356 Pa for the healthy subjects, and the density of 1040 kg/m 3 ([17], [42]). The CSF flow rate in the lateral ventricles is 0.35 cm 3 / min [38]. This value is used as the amplitude in the input fluid pulsatile flow rate function for numerical models. The final section of the ventricular system after the fourth ventricle is selected for the flow output location in Model A2. We set the normal baseline CSF pressure at 500 Pa. The computation is performed in 2D under the assumption of axial symmetry of the problems. Models A1 and A2 are discretized by taking into account the approximating systems defined in Section A.1 and are solved by using the Runge-Kutta MATLAB solver ODE45, with the initial mesh size determined by dz = L/nz, where L is fixed equal to 1 for the sake of simplicity, and nz is the total number of meshes that we take equal to 100. Finally we set the time discretization with ∆t = 5 × 10 −3 . A detailed analysis of the numerical methods we shall adopt for the models analyzed is beyond the scope of the present paper, and will be subject of future investigation.
The present section is organized as follows. where in Model A1 and in Model A2.
Then we run numerical simulations in order to validate the main Theorem of the paper which proves global existence of solutions provided that proper conditions on the initial data of the CSF flow velocity and of the intracranial pressure are satisfied. We perform this task by assuming u 0 (z) = 4 sin(πz) + 1, which are not so far from a good real approximation of the CSF flow velocity in a small compartment perfused by the fluid. As a byproduct of our numerical results, we obtain interesting informations about the behaviours of the quantities involved in these models.
• As a second task, we implement again the models numerically but we fix initial data u 0 (z) and P 0 (z) which violate the conditions (2.8) and (2.9) (we show here only the case of initial intracranial pressure in upright posture) of the Theorem 2.5. We shall compare these simulations to the analysis developed in this paper and show that our approach is reliable. In order to do that we assume u 0 (z) = −(exp(z) + 1), P 0 (z) = exp(z).  In Figures 6 and 7 we can observe the evolution of all the quantities involved in the CSF dynamics in a single compartment of length L and after a certain timet > 0.
• In Model A1: the intracranial pressure, P (t, z), shows the maximum value at the inlet and then it decreases through the compartment; the tissue displacement, η(t, z) undergoes a small decrease close to the midsection and an increase at the boundary, which correspond to a compression and an expansion of the compartment respectively; Figure 5: The initial configuration of P, η, u, ζ and A for Model A2.  the flow velocity, u(t, z), shows a parabolic profile that represents exactly what we expected since the CSF flow is laminar; moreover, by taking into account the initial configuration, it is possible to notice that at timet the velocity at midsection is decreased; for the cross section, A(t, z), we can observe an increase, then an enlargement of the section, in the first part of the compartment and a narrowing in the second part.
Therefore, corresponding to the maximum value of the pressure we have an increase in cross section which shrinks in relation to the decrease of the intracranial pressure. In fact, it occurs a phenomenon called Venturi effect for which in a tube with variable section the pressure decreases as section decreases (by the continuity equation). Strangely, the compression and expansion of the tissue in the compartment does not reflect the behavior of the cross section, in fact a tissue compression occurs in relation to the increase of the section and vice versa, while we would expect the opposite.
• In Model A2 the intracranial pressure shows the same behavior as in Model A1; this is due to the fact that the Riccati equation which really rules the pressure evolution in time, in general seems to evolves independently of the other quantities involved in the models; the tissue displacement increases continuously through the compartment by starting from a minimum value at the inlet; the flow velocity behaves as in Model A1 but the value at midsection in this model is greater than the corresponding value in the initial configuration; cross section decreases up to a minimum value at the midsection and in the second part of the compartment it shows an increase.
Thus, as in Model A1, in Model A2 an unusual behavior occurs when the intracranial pressure reaches its minimum value at the outlet since the cross section shows an increase rather than a compression. On the other hand, tissue displacement and cross section behaviors are properly related, in particular at the outlet of the compartment, where the maximum value of η(t, z) corresponds to a section enlargement.
The behaviour of the parameters analyzed here, up to the final timet = 1.000, is due to the boundary conditions we imposed for both models and to the fact we are neglecting fundamentals interactions in the cerebral pattern, that are assumptions which do not reflect the real CSF physiology and represent a limitation in the intracranial detailed description. But the main result obtained with these first simulations is the regularity of the solutions displayed above which is in good agreement with the Theorem 2.5. Second case: the initial data violate the conditions of the Theorem 2.5.

Comparison to real data
Once validated the main results stated in the present paper, a comparison with respect to the real dynamics detected in the intracranial pattern is required in order to understand how mathematical assumptions work and are able to alter the effective behavior of all quantities involved in both Models A1 and A2.
It is interesting to compare in particular the CSF flow velocity and the ICP values available in literature, thanks to the modern MRI techniques, with the results obtained in our analysis and simulations. In Fig. 10 we can observe that in Model A1, with a peak of 1, 4 mm/s, CSF flow velocity is   very slow with respect to the value indicated in Table 3. On the contrary, even though the boundary conditions adopted here, the values assumed by the intracranial pressure are quite consistent with the real data since after the first second it lies in the range of values specified in Table 3.
In Fig. 11 the plot of the CSF flow velocity shows a behavior similar to the one displayed in Fig. 10 but the maximum value is about 1, 23 mm/s which is smaller than the value detected in Model A1, as well as the value listed in Table 3. This is due to the fact that Model A2 describes the reabsorption process of CSF characterized by a velocity gradient decrease and by a damping effort of the resistance to the CSF absorption. The intracranial pressure is affected by a negligible alteration and the evolution plotted is very far from the data expected in the real physiology either in a supine position or in the upright posture.
Concerning the tissue displacement and the cross section, the comparison with real data is not self-evident since the compartments perfused by the CSF are very different each other either in anatomy or in physiology. Moreover, Models A1 and A2 do not take into account important interactions among the systems involved in the CSF dynamics, therefore a more detailed CSF model (see [27]) would be required in order to have a very accurate discussion regarding the differences pointed out between model predictions and published data.

Conclusions
This paper presents an extensive analysis of two CSF models, Models A1 and A2, that represent an improvement with respect to the study provided in [14]. As a first attempt in this direction, starting by the CSF model proposed by Linninger et al. in 2005 (see [25]) we have adopted a simplified model that describes the CSF flow without considering the different and complex interactions that occur inside the skull. We have also included to both models an intracranial pressure equation with the aim of ensuring a more detailed representation than the one showed in [14], although it provides a cruder representation with respect to more recent CSF models. This choice is intended to be only the first step towards the analysis of much more realistic and anatomically complete models (see [26]), able to describe every interaction that occurs between the three main subsystems in the intracranial pattern: cerebral vasculature, CSF, porous parenchyma.
In the present work we first demonstrated the local existence and uniqueness in time of solutions to both Models A1 and A2, then we proved the global existence and uniqueness of solutions under stringent conditions on the initial data. Furthermore, the numerical simulations supported the theoretical results and show evidence that the proposed analysis is consistent as a starting point also in cases in which more parameters and physiological quantities are involved in the models in order to guarantee a more comprehensive approach to the CSF dynamics. The authors will continue the study of this problem in a future work.

APPENDIX
This appendix includes the detailed proofs of the main theorems stated in the present paper.

A. PROOF OF THE THEOREMS 2.3 AND 2.4
In this Section we are going to prove the first result of our paper. In order to do that we need to proceed in two steps: we set up an iterative process by means of approximating systems associated to (1.4) and (1.5) and then we prove the convergence of the approximating solutions to classical solutions of Models A1 and A2.

A.1 Approximating scheme
For the purpose of being more precise we define one by one the approximations systems of (1.4) and (1.5) as follows.
The approximating systems take the following form As A1 with the following initial conditions Henceforth the analysis of the approximating systems is very similar, then we will specify the model only in case of significant differences. We want to show that there exist approximating solutions to the As A1 (A.6) and As A2 (A.7). First of all we observe that the solution of equation (A.6) 5 is defined by Hence, by considering that the equation (A.6) 4 with initial condition (A.8) 1 defines a first order hyperbolic system, we can apply the Theorem 2.1 with M = [0, L] and f (x, t) = ∂ z P 1 (t, z). Since every hypothesis of the Theorem is satisfied, we can conclude that there exists a unique solution u 1 such that Moreover we know that equations (A.6) 3 and (A.7) 3 with initial condition (A.8) 5 admit the following solution , (A.14) for the system (A.1) and for the system (A.2). By using the assumptions (A.8) and the fact that  (A.19) Therefore, the previous statements allow us to infer that there exist a unique solution X 1 (t, z) = (u 1 (t, z), η 1 (t, z), ζ 1 (t, z), P 1 (t, z), A 1 (t, z)), (A.20) to the system (A.6) ((A.7) respectively) such that (n + 1)-th STEP.
Now we want to prove there exists a unique solution for the systems (A.1) and (A.2). In order to do that we recall that the n-th iteration guarantees that there exists a unique solution X n (t, z) to the approximating systems such that Moreover, as well as the basic step, we know that where and A n+1 (t, z) = h(z)e −G n (t,z)t + e −G n (t,z)t t 0 H n (s, z)e G n (s,z)s ds, (A.31) in the approximating system A1 and (A.32) in the approximating system A2. Now we are ready to outline the main result of this section in the following proposition.
Proposition A.1 (Existence of the approximating solution). There exists a unique solution for every n.

A.2 Convergence of the iterative schemes
In this section we focus on the second part of the proof of the Theorems 2.3 and 2.4. From the Proposition A.1 we know that there exists a unique solution for the approximating systems associated to Model A1 and Model A2 respectively and now we want to prove the strong convergence of these solutions to classical solutions of the systems (1.4) and (1.5). Due to its nonlinear behaviour, the equation (1.4) 4 needs to be treated carefully. To this aim we recall the comprehensive analysis developed in [14] where the convergence of the approximating sequence u n (t, z) is proved by using high order energy estimates in the space H s , s > 7/2. The convergence of the other sequences of functions will be discuss in terms of equicontinuity.
A.2.1 Analysis and convergence of the sequence {u n } {u n } {u n } As we have already stated before, the sequence u n (t, z) is widely studied in [14] where the CSF velocity equation is exactly the same that appears in both Models A1 and A2, i.e.
ρ∂ t u(t, z) + ρu(t, z)∂ z u(t, z) + βu(t, z) + ∂ z P (t, z) = 0, (A. 36) with initial condition Then, based on the results obtained in [14], we define the energy in order to perform a priori estimates. We proceed by setting the following notation for any α < s. For the sake of simplicity we will not indicate the space and time dependence in this section. Now we define and we apply the operator ∂ α ∂z α to the equation (A.36) and we get By the detailed estimates in [14] we know that is a locally bounded function of its argument. Summing up (A.43) over α < s − 1 and integrating over (0, t), where t ∈ (0, T ), T ≤ T 0 , we find Taking the sup over t ∈ (0, T ) in (A.44) we get v 2 By using the notation (A.38) this means that Moreover we can show that for every n ≥ 0 u n+1 and In order to prove the convergence of the sequence {u n } ∞ n=0 , we subtract two subsequent equations (A.1) 4 and we obtain By carrying out estimates similar to the previous ones, we obtain and we find that where we assumed T so small that In this section, we shall show the convergence of {P n }, {η n }, {ζ n } and {A n }, and due to the different structure of the equations they satisfy we cannot employ the same methods as for u n . Therefore we need to resort to other methods for the convergence proof, in particular we recall the following Ascoli's Theorem ( [37]). We want to apply the Theorem A.1 to our sequences, more precisely we have to show that every sequence of functions is equicontinuous (see (2.2)).
By using the Proposition A.1 we obtain that in (A.62), we conclude that This last result suggests the sequences {P n }, {η n }, {ζ n } and {A n } are therefore equicontinuous (see Definition 2.2); hence is trivial to prove they are uniformly bounded, as well. We can apply now the Theorem A.1 which implies that the set where T < T 0 . Then, the relative compactness of Ψ allows us to conclude that

A.3 Existence and uniqueness of the solution
In the previous sections we attained the convergence of the approximating sequences of Models A1 and A2. This means that we are allowed to pass into the limit in the systems (A.1) and (A.2) and we obtain that is a classical solution of (2.5) ((2.7)). The last result to be shown is that the solution is unique. We assume there is another solution Ψ = u, P , η, ζ, A of the system (2.5) ((2.7)) with initial condition (2.6). By denoting and after analogous calculations as performed in Section A.2.1 (see in [14]), we end up with the inequalities ds,

B. PROOF OF THE THEOREM 2.5
This section is devoted to the problem of the global existence of solutions for Models A1 and A2. This issue has been already investigated in [14] and for the simplified CSF analyzed model it has been proved existence and uniqueness of a global solution under some restriction conditions on the initial data. As well as in [14] the core of the problem here is the Burgers-like behaviour of the equation ρ∂ t u(t, z) + ρu(t, z)∂ z u(t, z) + βu(t, z) + ∂ z P (t, z) = 0, (B.1) which affects the existence time of the solutions η, ζ and A for the other equations that appear in the systems A1 and A2. This peculiarity allows us to apply in the present paper the results obtained in [14] with some proper variations related to the structure of the new models.

B.1 Global existence of solutions
In order to evaluate if it is possible to prove global existence of solutions we want to draw attention to the following points. (c) According to the analysis developed in [14], we apply the method of characteristics to (B.1) and we denote by ω = u z , which satisfies along the characteristic curve Γ λ the following ordinary nonlinear differential Cauchy problem    ω + ω 2 + β ρ ω + β ρ P zz = 0, where λ = z for t = 0. Therefore, the focus of the attention shifts to the study of a Riccati nonhomogeneous equation for which we need to find a particular solution, ω(t). This problem has been already treated in [14] and we know that the existence of ω(t) is guaranteed if and only if P zz (t, z) L 2 t,z ≤ ∞. (for more details see [14], Section 7.2).
By taking into account the previous points we proceed in two steps, as follows.

STEP 1.
First of all we recall that the solution of the equation R∂ t P (t, z) − KP (t, z) 2 − KP (t, z) RQ p + P = 0, (B.9) with initial condition P (0, z) = b(z) ∈ H s , is the following (B.10) Therefore, we deduce higher regularity for the pressure than the one achieved in [14]. In particular, if b(z) > 0, this is true provided t ∈ [0, T ), where T < T 0 < T 0 and T 0 = 1 C ln RC Kb(z) + 1 is the blow up time (see Remark 1.1). By the previous assumption, in this case we are sure that We observe that the condition (B.7) is automatically satisfied by taking into account (B.10) and the local existence if b(z) > 0.
In order to complete the proof of the Theorem 2.5, we have just to recall the following theorem ( [30]).
Theorem B.1 (A Sharp Continuation Principle). Let be u 0 ∈ H s for some s > 5 2 and T > 0 some given time. Assume that for any interval of classical existence [0, T * ], T * ≤ T for u(t) solution of the equation (B.1), the following a priori estimate is satisfied: where M is a fixed constant independent of T * .
Then for T * with 0 ≤ T * ≤ T and the constants C, M depend on the space interval [0, L] and on some physics quantities.
Since the hypothesis of the Theorem B.1 are satisfied, we can conclude that there exists a global unique solution X (t, z) = (u, P, η, ζ, A)(t, z) to the systems of equations (2.5) and (2.7).