Analysis of solutions for a cerebrospinal fluid model

TL;DR

This paper investigates a mathematical model of cerebrospinal fluid dynamics, establishing existence and uniqueness of solutions and validating findings through numerical simulations under various initial conditions.

Contribution

It provides a rigorous mathematical analysis of the cerebrospinal fluid model, including proofs of solution existence, uniqueness, and numerical validation.

Findings

Existence and uniqueness of local solutions proven.

Global solutions exist under certain initial data restrictions.

Numerical simulations confirm theoretical results.

Abstract

The main purpose of this manuscript is to analyze an intracranic fluid model from a mathematical point of view. By means of an iterative process we are able to prove the existence and uniqueness of a local solution and the existence and uniqueness of a global solution under some restriction conditions on the initial data. Moreover the last part of the paper is devoted to present numerical simulations for the analyzed cerebrospinal model. In particular, in order to assess the reliability of the stated theoretical results, we carry out the numerical simulations in two different cases: first, we fix initial data which satisfy the conditions for the global existence of solution then, we choose initial data that violate them.

Tables

Table 1. Table 1: Tissue and fluid properties.

Property	Value
Fluid density, $\rho$	$1004-1007$ kg/m3
Tissue width, $\delta$	$\sim 5\times {10}^{-4}$ m
Fluid viscosity, $\mu$	${10}^{-3}$ Pa s
Spring elasticity, $k_e$	$8$ N/m
Brain dampening, $k_d$	$0,35\times {10}^{-3}$ (N s)/m

Equations

$$\{\begin{array}{l}{\partial }_t\eta (t,z)+{\partial }_{t}a(t)+u(t,z)-{\tilde{Q}}_p=0,\\ \alpha {\partial }_{tt}\eta (t,z)+\tilde{k}{\partial }_t\eta (t,z)+\kappa \eta (t,z)-AP(t,z)+A\tilde{P}=0,\\ \rho {\partial }_{t}u(t,z)+\rho u(t,z){\partial }_{z}u(t,z)+{\partial }_{z}P(t,z)+\beta u(t,z)=0,\end{array}$$

$$C^0(\bar{\Omega })=\{f:\bar{\Omega }\to \mathbb{R}|f\text{mbox}iscontinuous\},$$

$$C^s(\bar{\Omega })=\{f:\bar{\Omega }\to \mathbb{R}|f\in C^0(\bar{\Omega })\text{mbox}and{\partial }^{\alpha }f\in C^0(\bar{\Omega })\text{mbox}forany\alpha \text{mbox}with\alpha \le s\},$$

$${\mathcal{H}}^s(\Omega )=\{v\in {\mathcal{D}}^{\prime }(\Omega )|{\partial }^{\alpha }v\in L^2(\Omega )\text{mbox}forall\alpha \le s\},$$

$${\Vert w\Vert }_{{\mathcal{H}}^{r^{\prime }}([0,L])}\le C_r{\Vert w\Vert }_{L^2([0,L])}^{1-\frac{r^{\prime }}{r}}{\Vert w\Vert }_{{\mathcal{H}}^r([0,L])}^{\frac{r^{\prime }}{r}}\text{mbox}foranyw\in {\mathcal{H}}^r([0,L]).$$

$$\{\begin{array}{l}\frac{\partial \{A[h+a(t)+\eta (t,z)]\}}{\partial t}=Q_p-Q_f(t,z),\\ \rho \left[\frac{\partial u(t,z)}{\partial t}+u(t,z)\frac{\partial u(t,z)}{\partial z}\right]+\frac{\partial P(t,z)}{\partial z}=-\frac{8\mu }{r^2}u(t,z),\end{array}$$

$$({\rho }_{w}A\delta )\ddot{\eta }(t,z)+k_d\dot{\eta }(t,z)+k_e\eta (t,z)-A[P(t,z)-\tilde{P}]=0,$$

$$\{\begin{array}{l}{\partial }_t\eta (t,z)+{\partial }_{t}a(t)+u(t,z)-{\tilde{Q}}_p=0,\\ \alpha {\partial }_{tt}\eta (t,z)+\tilde{k}{\partial }_t\eta (t,z)+\kappa \eta (t,z)-AP(t,z)+A\tilde{P}=0,\\ \rho {\partial }_{t}u(t,z)+\rho u(t,z){\partial }_{z}u(t,z)+{\partial }_{z}P(t,z)+\beta u(t,z)=0,\end{array}$$

$$u(0,z)=u_0(z)=f(z)\in {\mathcal{H}}^s([0,L]),$$

$$u(t,0)=u(t,L)=0,$$

$$u(0,0)$$

$$u(0,L)$$

$$f(0)=f(L)=0.$$

$$\eta (0,z)={\eta }_0(z)=g(z)\in {\mathcal{H}}^s([0,L])$$

$$\eta (0,0)$$

$$\eta (0,L)$$

$${\partial }_t\eta (0,0)$$

$${\partial }_{tt}\eta (0,0)$$

$$P(0,0)=s(0)=\tilde{P}-\frac{\bar{\alpha }\alpha {\omega }^2}{A}+\frac{\tilde{k}\bar{\alpha }\omega }{A}+\frac{\tilde{k}{\tilde{Q}}_p}{A}+\frac{\kappa }{A}g(0).$$

$${\partial }_t\eta (0,z)$$

$${\partial }_{t}u(0,z)$$

$${\partial }_{tt}\eta (0,z)$$

$$P(0,z)=s(z),$$

$$s^{\prime }(z)-\frac{A\rho }{\alpha }s(z)=h(z),$$

$$h(z)=\rho \bar{\alpha }{\omega }^2-\rho f(z){\partial }_{z}f(z)$$

$$\displaystyle-\tilde{k}f(z)+\kappa g(z)+A\tilde{P}\bigg{]}.$$

$$P(0,z)=s(z)={\int }_0^z{e}^{\frac{A\rho }{\alpha }(z-\zeta )}h(\zeta )d\zeta +s(0)e^{\frac{A\rho }{\alpha }z}.$$

$${\partial }_t\eta (t,0)=-{\partial }_{t}a(t)+{\tilde{Q}}_p,$$

$$\eta (t,0)=g(0)-a(t)+0,3\bar{\alpha }+{\tilde{Q}}_{p}t.$$

$${\partial }_{tt}\eta (t,0)$$

$$=-\bar{\alpha }\left(-{\omega }^2\sin \left(\omega t-\frac{\pi }{2}\right)+2{\omega }^2\cos \left(2\omega t-\frac{\pi }{2}\right)\right),$$

$$P(t,0)$$

$$-\frac{\bar{\alpha }\tilde{k}\omega }{A}\left(\cos \left(\omega t-\frac{\pi }{2}\right)+\sin \left(2\omega t-\frac{\pi }{2}\right)\right)+\frac{\kappa }{A}g(0)+\frac{\tilde{k}+\kappa t}{A}{\tilde{Q}}_p$$

$$-\frac{\kappa \bar{\alpha }}{A}+\tilde{P}.$$

$$\eta (t,L)$$

Full text

1 Introduction

The main purpose of the present paper is to analyze an intracranic fluid model from a mathematical point of view. One of the main difficulties is related to the complexity of the intracranial dynamics which is origin of many different phenomena: the flow of the cerebrospinal fluid throughout the CSF compartments, the mechanical interaction between the fluid and the brain, the physiology of the brain, coupling with the circulatory-system and production and reabsorption laws.

There is a lack in literature of rigorous mathematical results regarding the analysis of the solutions to the systems of equations which model the different mechanisms in the intracranic pattern. Understanding the mathematical properties of the model is often intimately related to establish its validity in the analysis of the physiological properties. Our aim here is then to provide a comprehensive study of the equations that rule the cerebrospinal fluid (CSF) dynamics.

The physiological process we are interested in originates in the choroid plexus (see Fig. 1), affected by the cardiac cycle, that produces CSF at a constant rate (approximately 500 ml per day): the systole induces an expansion of the choroid that acts like a driving force for the CSF motion. Once secreted the CSF flows through four ventricles linked by different foramina and finally it reaches the subarachnoidal space, that is connected to the fourth ventricle by the foramina of Luschke and Magendie (see Fig. 2). At this point CSF is absorbed across the arachnoid villi into the venous circulation and a significant amount probably also drains into lymphatic vessels around the cranial cavity and spinal canal. The arachnoid villi act as one-way valves between the subarachnoid space and the dural sinuses.

In order to make the cerebral dynamics more comprehensible, detailed mathematical models quantifying forces and their interaction have become fundamental to reveal what medical instruments are not able to do without affecting the data.

In fact since the mid seventies, the mathematical modeling of the intracranial dynamics has gained interest and many researchers have developed and analyzed models of different complexity for each of the three compartments (brain parenchyma, cerebral vasculature and cerebrospinal fluid), either standing alone or coupled. An extensive overview of the models for intracranial dynamics may be found in [2], [12].

More precisely, in order to link the mathematical models for the intracranial dynamics it is possible to identify two classes: models based on partial differential equations and models based on a lumped compartmental description of the intracranial constituents.

Linninger et al. (see [6]) developed a lumped model for describing the CSF pulsatility during the cardiac cycle. This first model included the main constituents present in the cranial vault, however only the CSF compartment was fully treated. In fact the model includes the lateral ventricles, the third ventricle, the fourth ventricle and the subarachnoid space.

The cerebral blood compartment is taken into account only as boundary condition for the CSF system. Indeed, the CSF pulsation is driven by a boundary condition modeling the effects of the choroid plexus which expands and contracts in a prescribed way (modeling the effect of the arterial pressure wave) and the CSF reabsorption is modeled by an equation where the venous pressure is a prescribed constant. The brain parenchyma is also treated in a simplified way. In fact equations to link the CSF compartments pressure and volume are obtained by supposing that each compartment is a cylinder with a thin shell membrane that deforms radially under the pressure difference between the CSF and the parenchyma.

In this paper we adopt a different approach to the study of the cerebrospinal fluid dynamics, in fact we start from the Linninger model by retracing and simplifying the more important steps of the modeling and we provide a comprehensive analysis of the system of equations that represents the pillar of the aforementioned dynamics, in order to fill up the gap between the mathematical theory and modelization.

First of all we clarify the system of equations that describes the CSF hydrodynamics and we obtain the following system

[TABLE]

where $u$ is the CSF velocity flux, $\eta$ is the tissue displacement and $P$ is the pressure. In this first approach we neglect the variation of the cross section $A$ due to the choroid plexus expansion that drives the pulsatile CSF circulation, then we consider it constant in order to study the effects of this statement on the real cerebral physiology. Moreover we choose properly the initial and boundary conditions with the related compatibility conditions that are all treated in detail in Section 2.

Then we proceed by proving the local existence and uniqueness of a solution to the system (1.1). This result is achieved by setting up a proper iterative scheme based on the approximating system strictly related to (1.1). As a first step we prove by induction on the number $n$ of iterations the existence and uniqueness of a solution, $\mathcal{X}^{n}(t,z)=(u^{n},\eta^{n},P^{n})(t,z)$, to the approximating system, then by using higher energy estimates ($\mathcal{H}^{s},s>\frac{9}{2}$) we prove the convergence of the approximating sequence, $\mathcal{X}^{n}(t,z)$, to a local classical solution of (1.1).

This first result allows us to investigate the global existence of a solution for the system (1.1). Since $\eta(t,z)$ and $P(t,z)$ inherit the lifespan of $u(t,z)$, we focus on the third equation of (1.1) that we study first in the homogeneous form in order to obtain a detailed framework. We apply the characteristic method and we get in both cases a Riccati equation that we are able to solve with the standard ODE methods for the homogeneous case and by means of a particular solution that we construct properly for the nonhomogeneous case.

Finally, by using a sharp continuation principle (see [8]) we show that there exists a global solution to (1.1) if and only if it is satisfied a precise condition on the initial datum of $u$.

The global existence of a solution paves the way for another important result due mainly to the coupling of the CSF flux and the cardiac cycle. It is well known that systole and diastole produce a periodic cardiac movement that holds a big role in the intracranic dynamics and this suggests the existence of particular solutions for the system (1.1) if we consider it independent of the spatial variable $z$. Indeed we are able to prove the existence and uniqueness of periodic solutions that are strongly influenced by the forcing term $a(t)$ and whose stability is ensured by the global existence time.

The present paper is organized as follows. In section 2 we begin a detailed discussion with a description of the Linninger model and, once we have reformulated it in a more simple way, we fix the initial and boundary conditions. In section 3 we state our main results on the existence and uniquess of a local solution and on the existence and uniqueness of a global solution to the system of equations (1.1). In section 4, we set up an iterative process by means of an approximating system associated to the (1.1). In section 5, in order to discuss the convergenge of the iterative scheme we derive higher order energy estimates and we prove the convergence for each one of the sequences that appear in the approximating system. Then we conclude the proof of the local existence and uniqueness of a solution to (1.1) in section 6. Section 7 is devoted to the problem of the global existence of solutions. Since the third equation of the system (1.1) has a “Burgers” like behaviour we focus first on that equation and then we deduce the global existence for the remaining equations. Hence we analyze in detail the global existence time in order to adopt a continuation principle that is needed for the proof of the second theorem stated in section 3. In section 8 we analyze the periodic system associated to (1.1) and we study the stability of its solutions. Finally, the last section of this paper is devoted to present numerical simulations for the analyzed cerebrospinal model. In particular, in order to assess the reliability of the results stated in section 3, we carry out the numerical simulations in two different cases: in the numerical method we fix first initial data which satisfy the conditions for the global existence of a solution to the system (1.1), then we choose initial data that violate them. In the first case no blow up of the solutions is displayed and we can observe a parabolic behavior for the velocity profile and, as a consequence, a smooth evolution for the pressure and the tissue displacement. While, in the second case, simulations achieve blow up after a single iteration. These results allow us to validate the statements obtained in the present paper.

Notations.  We need to define the functional spaces that are necessary to study our problem.

Let $\Omega\subset\mathbb{R}$ be an open bounded set. For every $s\in\mathbb{N}$ we define:

[TABLE]

[TABLE]

[TABLE]

where $\mathcal{D}^{\prime}(\Omega)$ is the space of distributions in $\Omega$.

We will denote by $\left\|\cdot\right\|_{L^{\infty}_{t}\mathcal{H}^{s}_{z}}$ the norm $\left\|\cdot\right\|_{L^{\infty}((0,T))\mathcal{H}^{s}([0,L])}$ and by $\left\|\cdot\right\|_{2}$ the norm $\left\|\cdot\right\|_{L^{2}([0,L])}$

Moreover, we recall the Sobolev interpolation theorem by which, if $0\leq r^{\prime}\leq r$, there exists a costant $C_{r}$ such that

[TABLE]

2 Mathematical model and governing equations

In this paper we analyze the model introduced in [6] which simplifies the intracranic dynamics by representing the fluid-structure interactions that originate in the craniospinal pattern and by quantifying the cerebrospinal fluid (CSF) motion such as the variations of the pressure due to the cerebral tissue compression. In order to do that, the ventricles are discretized into cylindrical finite volume with perfect axial dispersion and radial  expandability, the foramina that link the different compartments are treated as elastic tubes and the CSF flow is considered basically laminar.

The continuity equation of the cerebrospinal fluid flow in the ventricles and the axial momentum equation are given by

[TABLE]

where

• •

$A$ and $h$ are the cross section and the height of the ventricular or subarachnoid section respectively;

• •

$r$ is the radius of the foramina and aqueduct;

• •

$\mu$ is the fluid viscosity;

• •

$a(t)=\bar{\alpha}\left(1.3+\sin\left(\omega t-\frac{\pi}{2}\right)-\frac{1}{2}\cos\left(2\omega t-\frac{\pi}{2}\right)\right)$ is the forcing function;

• •

$\eta(t,z)$ denotes the tissue displacement in a section;

• •

$Q_{p}$ is the CSF production rate in the choroid plexus;

• •

$Q_{f}(t,z)=Au(t,z)$ is the CSF flow rate leaving ventricle, i.e., flow in foramina and aqueduct;

• •

$u(t,z)$ is the axial CSF flow velocity.

The dynamics of parenchyma stresses, strains and displacements can be described with the laws of elastodynamics,

[TABLE]

where

• •

$k_{d}\dot{\eta}(t,z)$ is the dissipative force;

• •

$k_{e}\eta(t,z)$ is the elastic force;

• •

$P(t,z)-\tilde{P}$ is the difference between fluid and parenchyma pressure.

Since the deformation of the tissue directly affects the space available to the fluid, the two systems are fully coupled.

This model allows us to compute the pressure $P(t,z)$ in the cerebral ventricles, the velocity $u(t,z)$ of the fluid along the foramina and the tissue displacements $\eta(t,z)$ at each point.

The simplified equations (2.1), (2.2) of motion for the hydrodynamics of the 1D intracranial dynamics can be written together in the following system:

[TABLE]

where $t\in[0,T]$, $z\in[0,L]$, $a(t)\in C^{\infty}([0,T])$, $\tilde{Q}_{p}=\displaystyle{\frac{Q_{p}}{A}}$, $\rho,\beta,\tilde{k},\kappa,\alpha,A,\tilde{P}\in\mathbb{R}$ are physical constants and known values.

First of all, for the unknown $u$ we impose the following initial condition

[TABLE]

where $\displaystyle{s>\frac{9}{2}}$.

Since this paper is a first attempt to perform a rigorous mathematical  analysis on the system (2.3) we start considering simplified boundary conditions neglecting some assumptions imposed by the realistic model. Hence for the single compartment of length $L$, we consider the following boundary conditions for $u(t,z)$:

[TABLE]

and the following compatibility conditions

[TABLE]

moreover we assume for the sake of simplicity

[TABLE]

Similarly for the unknown $\eta$ of the system (2.3), we assume

[TABLE]

and the following compatibility conditions

[TABLE]

Now we want to compute the initial and boundary conditions for the pressure, $P$. With the previous assumptions and by using the first two equations of our system we find that

[TABLE]

hence

[TABLE]

In the same way we get

[TABLE]

By using the notation

[TABLE]

and by replacing everything in the second equation of (2.3) we get

[TABLE]

where

[TABLE]

The solution of the ODE (2.14) with the initial condition (2.11) is given by

[TABLE]

To establish the correct boundary conditions for the unknowns $\eta$ and $P$, we proceed as before with the assumptions (2.5), (2) and the forcing function $a(t)$. Then we get

[TABLE]

and from the previous one we obtain the first boundary condition for $\eta$

[TABLE]

Moreover

[TABLE]

and the first boundary condition for the unknown $P$ is given by

[TABLE]

Similarly we obtain the following boundary condition for $z=L$

[TABLE]

Finally, we observe that further compatibility conditions are given by

[TABLE]

3 Main results

Now we are ready to state the main results of this paper. The local existence and uniqueness of a solution to the system (2.3) with the initial conditions (2.4), (2.8), (2.16) and boundary conditions (2.5), (2.17), (2), is given by the following theorem.

Theorem 3.1**.**

Let us consider the system

[TABLE]

where $t\in[0,T_{0}]$, $z\in[0,L]$, $a(t)\in C^{\infty}([0,T_{0}])$, $\tilde{Q}_{p}=\displaystyle{\frac{Q_{p}}{A}}$, $\rho,\beta,\tilde{k},\kappa,\alpha,A,\tilde{P}\in\mathbb{R}$ are constants, with initial conditions

[TABLE]

$\displaystyle{s>\frac{9}{2}}$, where

[TABLE]

with

[TABLE]

Consider the following boundary conditions for the system (3.1)

[TABLE]

and assume that the following compatibility conditions are satisfied

[TABLE]

Then, there exists a time $T<T_{0}$ such that the problem (3.1), (3.1), (3.1) has a local unique solution

[TABLE]

such that

[TABLE]

Remark 1*.*

By using together the regularity assumptions (3.1) for $u_{0}(z)$ and $\eta_{0}(z)$ and (3.3), (3.1), we have $P(0,z)=s(z)\in\mathcal{H}^{s}([0,L])$.

A natural question which arises now is under which conditions it is possible to get global solutions. In what follows we will show that the answer is yes provided a proper choice of the initial conditions, as is stated in the following theorem on the global existence and uniqueness of a solution to the initial boundary value problem (3.1), (3.1), (3.1).

Theorem 3.2**.**

Let us consider the system

[TABLE]

where $t\in[0,T]$, $z\in[0,L]$, $a(t)\in C^{\infty}([0,T])$, $\tilde{Q}_{p}=\displaystyle{\frac{Q_{p}}{A}}$, $\rho,\beta,\tilde{k},\kappa,\alpha,A,\tilde{P}\in\mathbb{R}$ are constants, with initial conditions (3.1) and boundary conditions (3.1). Assume that the compatibility conditions (3.1) are satisfied and that

[TABLE]

*where $f(z)=u(0,z)$.

Then there exists a global unique solution*

[TABLE]

to the problem (3.8),(3.1),(3.1) such that

[TABLE]

for every $T\geq 0$.

The next sections of the paper are devoted to the proof of the Theorems 3.1 and 3.2.

Moreover the last Section 8 deals with a further development of the global existence theorem. We analyze the particular framework in which we assume initial data very close to periodic conditions that keep the same period of the forcing term $a(t)$ and we state a stability theorem for periodic solutions of the system (1.1).

4 Iterative scheme

In this section we start the proof of the Theorem 3.1. We divide the proof in two steps. In order to find a solution to the system (3.1) we define an approximating system by setting up an iterative process and we proceed by induction. In the second step we prove that the iterative scheme converges to a solution of the problem (3.1) with the initial conditions (3.1).

We define the approximation of the system (3.1) as follows; the sequence $\mathcal{X}^{n+1}=(u^{n+1},\eta^{n+1},P^{n+1})$ satisfies the following system

[TABLE]

where $z\in[0,L]$, $t\in(0,T_{0})$, with $T_{0}\geq T>0$. The initial conditions are given by

[TABLE]

and the boundary conditions are the following

[TABLE]

For the initial step we fix

[TABLE]

In order to prove the existence of a solution to the approximating system (4.1) we proceed by induction on the iteration number $n$.

Basic step: $n=0$. If we fix $n=0$ the system (4.1) becomes

[TABLE]

with the following initial conditions

[TABLE]

Now, we prove the existence of a solution to the problem (4.5), (4). From the first equation of the system (4.5) we obtain a unique solution $\eta^{1}$ of the form

[TABLE]

Then, by using the regularity conditions (4) we have

[TABLE]

Now we replace (4.7) into the second equation of (4.5) and we obtain

[TABLE]

Since $a(t)\in C^{\infty}([0,T_{0}])$ and by using the assumptions (4), we get that

[TABLE]

Moreover, since

[TABLE]

we observe that

[TABLE]

Finally we need to find $u^{1}$. The sequence $u^{1}$ satisfies the following first order hyperbolic system

[TABLE]

In order to find a solution of (4.13), we need to recall the following existence theorem for a first order symmetric hyperbolic system ([10]).

Theorem 4.1**.**

Let $\mathbb{A}_{j}(x,t),\,j=0,...,N$ be symmetric $m\times m$ matrices in $M\times(0,T)$ with $M\subset\mathbb{R}^{N}$, $T>0$, $\mathbf{f}(x,t)$ and $\mathbf{v}_{0}(x)$, m-dimensional vector functions defined in $M\times(0,T)$ and $M$, respectively. Let us consider the problem

[TABLE]

Let

[TABLE]

$\mathbb{A}_{0}(x,t)$* invertible, $\inf_{x,t}\left\|\mathbb{A}_{0}(x,t)\right\|_{M\times M}>0$, $\mathbb{B}\in C((0,T);\mathcal{H}^{s-1}(M))^{m\times m}$, $\mathbf{f}\in C((0,T);\mathcal{H}^{s}(M))^{m}$, $\mathbf{v}_{0}\in\mathcal{H}^{s}(M)^{m}$, where $s>\frac{N}{2}+1$ is an integer. Then there exists a unique solution to (4.14a), (4.14b), i.e. a function*

[TABLE]

satisfying (4.14a) and (4.14b) pointwise (i.e. in the classical sense).

We apply the Theorem 4.1 to the problem (4.13), with $M=[0,L]$, $\mathbf{f}(x,t)=\partial_{z}P^{1}(t,z)$ and we observe that we fullfill the hypothesis of the Theorem 4.1. Then, there exists a unique solution $u^{1}$ such that

[TABLE]

Finally, we have shown that there exists a unique solution

[TABLE]

to the system (4.5) with initial data (4) such that

[TABLE]

Now, we assume that for the n-th step there exists a unique solution

[TABLE]

such that

[TABLE]

and we want to prove the existence of the $(n+1)$-th step.

(n+1)-th step. The $n+1$ iteration corresponds to the system

[TABLE]

with the initial data (4).

Similarly to the basic step we can write

[TABLE]

and we get

[TABLE]

We replace (4.20) in the second equation of the system (4.19) and we obtain

[TABLE]

From the $n$-th iteration we know there exists a unique solution

[TABLE]

by which we get

[TABLE]

After that, we observe that the last equation of the system (4.19)

[TABLE]

satisfies the hypothesis of the Theorem 4.1. In fact

[TABLE]

and by using (4.22) we get

[TABLE]

So by applying the Theorem 4.1, we can assert that there exists a unique solution

[TABLE]

Finally we can summarize the existence of a solution to the approximating system in the following proposition.

Proposition 1**.**

There exists a unique triple of solutions

[TABLE]

to the problem (4.19),(4), where

[TABLE]

5 Convergence of the iterative scheme

In this section we will prove that the iterative scheme previously constructed converges to a solution of the problem (3.1), (3.1), (3.1).

First of all we will recover suitable energy estimate for our problem by analyzing the equation $\eqref{ar50}_{3}$ that is the most tricky one because of its nonlinearity but at the same time it is the equation that allows us to obtain results about the whole system (3.1).

We prove that $\{u^{n}\}$ is a Cauchy sequence and so we get the existence of a classical solution to the equation $\eqref{ar24}_{3}$. Hence, we will be able to prove the strong convergence for the other sequences of unknowns and then the existence of classical solutions to the equations $\eqref{ar24}_{1}$ and $\eqref{ar24}_{2}$.

5.1 Higher order energy estimates for the sequence $\boldsymbol{\{u^{n}\}}$

We proved that

[TABLE]

is a solution to the equation

[TABLE]

Henceforth we will not denote explicitly the space and time dependence and in order to simplify the notation, we set

[TABLE]

[TABLE]

for any $\alpha<s$.

First of all we analyze the third equation of the system (3.1) that, with the new notation, becomes

[TABLE]

If we apply to the equation (5.4a) the operator $\displaystyle{\frac{\partial^{\alpha}}{\partial z^{\alpha}}}$, then we get

[TABLE]

Let’s define

[TABLE]

Applying the time derivative to (5.6) and by using the equation (5.5a) we get

[TABLE]

Now, we estimate each one of the terms of the right hand side of (5.1).

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

In a similar way for the last terms we get

[TABLE]

Then, by summing up (5.8)-(5.1) we get

[TABLE]

where

[TABLE]

[TABLE]

Now, after a detailed analysis of the terms (5.14) and (5.1) we have

[TABLE]

where $C_{1}=C_{1}\left(\left\|u\right\|^{2}_{L^{\infty}_{t}\mathcal{H}^{s-3}_{z}}\right)$ is a locally bounded function of its argument. Summing up (5.16) over $\alpha<s-3$ and integrating over $(0,t)$, where $t\in(0,T)$, $T\leq T_{0}$, we find

[TABLE]

Taking the $\sup$ over $t\in(0,T)$ in (5.17) we get

[TABLE]

then

[TABLE]

By using the notation (5.2) this means that

[TABLE]

We take

[TABLE]

and

[TABLE]

We can show that for all $n\geq 0$

[TABLE]

So, assuming (5.22) true for a given $n\geq 0$, by (5.20) we have

[TABLE]

by induction (5.23) follows for all $n$. By using (5.24) and the equation $\eqref{ar50}_{3}$ we obtain also that

[TABLE]

where

[TABLE]

is the constant inherited by the $n$ iterations.

5.2 Convergence of the sequence $\boldsymbol{\{u^{n}\}}$

Now we are almost ready to prove the convergence of $\{u^{n}\}_{n=0}^{\infty}$. Subtracting two subsequent equations in (5.1) we get

[TABLE]

where

[TABLE]

Following exactly the same line of arguments adopted for the energy estimates of the previous section, we get

[TABLE]

By (5.22), (5.23), (5.28) we get

[TABLE]

Taking into account that

[TABLE]

we have

[TABLE]

Hence, by using (5.29) and (5.32) we obtain

[TABLE]

Assuming $T$ so small that

[TABLE]

we find that

[TABLE]

which implies that there exists $u\in L^{\infty}\left((0,T);L^{2}([0,L])\right)$ such that

[TABLE]

Then by (5.23), (5.36) and the interpolation inequality (1.2), for $0<s^{\prime}<s-3$ and $l,m\in\mathbb{N}$ we get

[TABLE]

From (5.2) it follows that

[TABLE]

and by the Sobolev embedding theorem we get

[TABLE]

Now, in order to increase the time regularity of $u$ we analyze the following difference

[TABLE]

We know that

[TABLE]

Hence, taking the sup over $t\in(0,T)$ and by using the previous computations we get

[TABLE]

Sobolev embedding theorem and (5.38) yield

[TABLE]

By considering together (5.39) and (5.42) we can say that $u\in C^{1}((0,T]\times[0,L])$ is a classical solution of the third equation of our problem.

5.3 Convergence of the sequences $\boldsymbol{\{\eta^{n}\}}$ and $\boldsymbol{\{P^{n}\}}$

In order to study the convergence of $\{\eta^{n}\}$ e $\{P^{n}\}$ we proceed exactly as in Sections 5.1, 5.2. We begin subtracting two subsequent equations in (4.20):

[TABLE]

by which we get

[TABLE]

where $a$ is defined as in (5.34). By using (5.35) we find

[TABLE]

Hence it follows that

[TABLE]

By the interpolation inequality (1.2) we obtain

[TABLE]

We also know that

[TABLE]

and that

[TABLE]

then, taking the sup over $t\in(0,T)$, we obtain

[TABLE]

We can conclude that

[TABLE]

Hence we get

[TABLE]

By using once more (1.2) and taking into account (5.3), we obtain that

[TABLE]

6 Existence and uniqueness

By using the results of the previous section now we are ready to prove the existence of solutions to the system (3.1) with initial data (3.1) and boundary conditions (3.1).

In fact, by considering (5.38), (5.42), (5.3), (5.3) we can pass into the limit in the system (3.1) and we get that $(u,\eta,P)\in[C^{1}([0,L]\times(0,T])]^{3}$ is a classical solution of (3.1).

The last step is to prove the uniqueness. To this aim we consider two solutions $u,v$ of the equation

[TABLE]

with initial condition given in $\eqref{ar25}_{1}$. Then $w=u-v$ satisfies

[TABLE]

Multiplying (6) by $w$ and by integrating in space we have

[TABLE]

where

[TABLE]

By Gronwall inequality and (6) we get $w=0$, hence $u=v$.

For $\eta$ and $P$ we can proceed in the same way and we obtain the uniqueness result.

This concludes the proof of the Theorem 3.1.

7 Global existence of solutions

In the previous sections we proved that the existence and the uniqueness of a solution to the problem (3.1), (3.1), (3.1) is guaranteed in a time interval $[0,T]$ where $T$ is the local existence time defined in Section (5.1) as

[TABLE]

We want to know if it possible to extend the local existence time interval and in particular how much we can push forward the maximal time interval in order to find a global solution to our system. The core of the problem which is related to the study of the global existence of a solution is concentrated on the third equation of (3.1),

[TABLE]

In fact it owns the same behaviour of the Burgers’ equation that has a blow up of the solutions at a certain “breaking” time.

Moreover, by looking at the structure of the system (3.1), we observe that both $\eta$ and $P$ lifespan can be obtained by analyzing (7.2), in other words, as a first step, we can focus only on the lifespan of the velocity flux $u$ that allows us to have the main framework of our problem.

7.1 Homogeneous equation

First of all we analyze the global existence for the homogeneous equation,

[TABLE]

with the initial condition

[TABLE]

with $\displaystyle{s>\frac{9}{2}}$.

To this purpose we define the characteristic curve $\Gamma_{\lambda}$ as

[TABLE]

By using (7.3),(7.4) we find that along the curve $\Gamma_{\lambda}$

[TABLE]

Hence on $\Gamma_{\lambda}$, by using (7.5) we have

[TABLE]

We consider now the space derivative of $u$, namely we set

[TABLE]

which by (7.3), (7.5), (7.7) satisfies along $\Gamma_{\lambda}$ the following ordinary nonlinear differential equation

[TABLE]

with initial data

[TABLE]

By classical ODE method we solve the Riccati Cauchy problem (7.9), (7.10) and we get that on $\Gamma_{\lambda}$

[TABLE]

Hence $u$ is a global solution whenever

[TABLE]

otherwise we have blow-up at the time $T$ given by

[TABLE]

By analyzing (7.13) we have that the solution blows up at the finite time $T>0$ if and only if

[TABLE]

Then we can conclude that the solution to (7.3), (7.4) exists globally if and only if

[TABLE]

7.2 Nonhomogeneous equation

Now we study the nonhomogeneous equation

[TABLE]

with initial condition

[TABLE]

where $\displaystyle{s>\frac{9}{2}}$.

We define the characteristic curve $\Gamma_{\lambda}$ as in (7.5) and we find that along the characteristic $u$ verifies

[TABLE]

As before we denote by

[TABLE]

which by (7.16), (7.18), (7.19) satisfies along $\Gamma_{\lambda}$ the ordinary nonlinear differential Cauchy problem

[TABLE]

In order to study the solutions of the nonhomogeneous Riccati equation (7.20), we need to find a particular solution that allows us to apply the standard methods for ordinary differential equations.

To this aim we look for a particular solution $\bar{\omega}(t)$ of the form

[TABLE]

Hence $a(t)$ and $b(t)$ satisfy the following equation

[TABLE]

In order to solve (7.22) and find $a(t),b(t)$ we set the following linear system

[TABLE]

that reduces the order of the equation (7.22), and we choose the following initial conditions

[TABLE]

The system (7.23) admits a unique solution if and only every element of the system coefficient matrix is bounded, in particular we need that the pressure $P$ is such that

[TABLE]

By the estimates for the local existence of solutions (see Section 5.1) we know that

[TABLE]

and

[TABLE]

where $\mathcal{C}$ is defined in (5.26).

Therefore from the equation (3.1) we obtain

[TABLE]

and, by using the relations (7.26) and (7.2), it yields

[TABLE]

where $\tilde{T}\in(0,T)$ and $T$ is the local existence time defined in (7.1). Hence we can conclude that the condition (7.25) is satisfied.

In order to find the solution of (7.23), first of all we compute the coefficient matrix eigenvalues of the linear system and we get

[TABLE]

Since the model we are dealing with is associated to a physiological flux, we look for real solutions $\omega(t)$. Hence, in order to obtain real eigenvalues (7.30) we need that $P_{zz}$ verifies the following costraint

[TABLE]

The condition (7.31) is fullfilled if we are able to show that

[TABLE]

which is achieved by using again the estimate (7.2).

At this point our goal is to prove that there exists a time $\tilde{T}$ such that the condition (7.32) is satisfied on the interval $[0,\tilde{T}]$. We can show that this is true if and only if the following condition is satisfied,

[TABLE]

In fact, by applying the condition (7.33) to (7.2) we get

[TABLE]

which actually fullfills (7.32) for a certain time $\tilde{T}\leq T$.

With this last result we have proved the existence of a real particular solution of the system (7.20), that is defined in the following way

[TABLE]

if $\lambda_{1}\neq\lambda_{2}$, where

[TABLE]

while

[TABLE]

if $\lambda_{1}=\lambda_{2}$, where

[TABLE]

Therefore we can say that the initial value problem (7.20) admits the  following unique solution

[TABLE]

where

[TABLE]

that is well defined for every $t\in[0,\tilde{T}]$.

7.3 Global existence time

Summarizing, we have proved that under the condition (7.32) there exists a solution $u(t)$ defined for any time $t\in[0,T^{\ast}]$, such that

[TABLE]

with

[TABLE]

where $T$ is the local existence time defined in (7.1).

At this point the main question is wheter it is possible to extend the time $T^{\ast}$ in order to obtain the global existence of the solution $u(t,z)$. For this purpose we need the following Sharp Continuation Principle (see Majda ([8]), Thm 2.2, p. 46).

Theorem 7.1** (A Sharp Continuation Principle).**

Let be $u_{0}\in\mathcal{H}^{s}$ for some $s>\frac{9}{2}$ and $T^{\prime}>0$ some given time. Assume that for any interval of classical existence $[0,T^{\ast}]$, $T^{\ast}\leq T^{\prime}$ for u(t) solution of the equation (7.16), the following a priori estimate is satisfied:

[TABLE]

*where $M$ is a fixed constant independent of $T^{\ast}$.

Then the classical solution $u(t)$ exists on the interval $[0,T^{\prime}]$, with $u(t)$ in $C\left([0,T^{\prime}],\mathcal{H}^{s-3}([0,L])\right)\cap C^{1}\left([0,T^{\prime}],\mathcal{H}^{s-4}([0,L])\right)$. Furthermore, $u(t)$ satisfies the a priori estimate*

[TABLE]

for $T^{\ast}$ with $0\leq T^{\ast}\leq T^{\prime}$ and the constants $C,M$ depend on the space interval $[0,L]$ and on some physics quantities.

By taking into account that

[TABLE]

we want to show that the solution (7.39) satisfies the condition (7.42) of the previous theorem. We observe that

[TABLE]

Since $\lambda_{1}\geq\frac{\beta}{\rho}$ and by applying the condition (7.32) to (7.45) we obtain

[TABLE]

Hence, we can conclude that the hypothesis of the Theorem 7.1 are totally fullfilled and the global existence of the solution $u(t,z)$ is proved, as well, provided the condition (7.33) holds on the initial datum.

Moreover, since $\eta(t,z)$ and $P(t,z)$ are obtained by plugging in the first two equations of the system (3.1) the solution $u(t,z)$ and since they inherit the life span of the latter, the Theorem 7.1 proves actually the global existence of a unique triple solution $\mathcal{X}(t,z)=(\eta,P,u)(t,z)$ to the problem (3.1), (3.1), (3.1).

8 Existence and stability of a periodic solution

This last section is devoted to prove the stability for some particular periodic solution for our system. We denote by $\bar{\eta},\bar{u},\bar{P}$ the periodic solutions, independent of the spatial variable $z$, namely they satisfy the following system

[TABLE]

We recall that the forcing function $a(t)$ is given by

[TABLE]

Moreover the existence of the periodic solutions is proved for a particular choice of the cerebrospinal fluid production rate $Q_{p}$.

We assume that

[TABLE]

where $F(t)$ is a periodic function. For the sake of simplicity and without loss of generality we choose for $F(t)$ the same period of the forcing function $a(t)$, $\bar{T}=\frac{2\pi}{\omega}$. This assumption is allowed by the dynamic we are analyzing, in fact the production rate is well approximated by a periodic function that ensures a constant rate in a proper range of time.

8.1 Existence of a periodic solution

In order to find a periodic solution $\bar{\chi}(t)=(\bar{\eta},\bar{u},\bar{P})(t)$ to the problem (8.1), we define the following initial conditions

[TABLE]

Hence we can compute explicitly the solution to the problem (8.1), (8.1),

[TABLE]

where both $\bar{\eta}(t),\bar{P}(t)$ are $C^{\infty}$ functions on $[0,T]$. The solution (8.1) is the unique periodic solution whose period is influenced by the forcing function $a(t)$ and it corresponds to $\bar{T}=\frac{2\pi}{\omega}$.

8.2 Stability of the periodic solution

Now we study the behaviour of the solution $\mathcal{X}(t,z)$ to the problem (3.1), (3.1) when we assume initial data very close to the periodic initial conditions (8.1). In order to analyze the stability of the periodic solution it is important to observe that we can push the time $T$ up to the period $\bar{T}=\frac{2\pi}{\omega}$ of the solution $\bar{\mathcal{X}}(t,z)=\left(\bar{\eta},\bar{P},\bar{u}\right)(t)$ and we can apply the following theorem on every time interval $[h\bar{T},(h+1)\bar{T}]$ with $h\in\mathbb{N}$. This explains why the estimates in the proof below are bounded no matter how big we choose the time, in fact by using Thm. 7.1  and the periodicity of the system (8.1) we always obtain the following main result.

Theorem 8.1**.**

Let $T>0$ and $s>\frac{9}{2}$. Let $(\eta,P,u)(t,z)$ and $(\bar{\eta},\bar{P},\bar{u})(t)$ be the solutions of the problem (3.1),(3.1) and (8.1),(8.1) respectively. There exist constants $\delta\in(0,1),K(\delta)$ such that, if

[TABLE]

for all $\delta$, then

[TABLE]

for all $\delta$.

Proof.

We start by setting

[TABLE]

which by standard computations satisfy the following system

[TABLE]

Moreover, the previous system is complemented with the following initial data

[TABLE]

and the following boundary conditions

[TABLE]

From now on for the sake of simplicity we will not denote explicitly the time and space dependence.

We observe that the assumptions (8.5) totally fullfills the condition (3.9), which means that we can apply all the classical solutions properties in what follows.

In order to prove our statement we perform, as before, higher order energy estimates of $(N,\Pi,U)$. From now on, for any $\alpha<s-4$ we will use the following notation $N_{\alpha}=\frac{\partial^{\alpha}}{\partial z^{\alpha}}N$, $\Pi_{\alpha}=\frac{\partial^{\alpha}}{\partial z^{\alpha}}\Pi$ and $U_{\alpha}=\frac{\partial^{\alpha}}{\partial z^{\alpha}}U$.

We begin by taking the order $\alpha$ derivative of the third equation of (8.8), that becomes

[TABLE]

in which we plug

[TABLE]

obtained by combining the first two equations of (8.8).

We define the energy as

[TABLE]

and applying the time derivative to (8.13) we get the following estimate

[TABLE]

By estimating and summing up over $\alpha<s-4$ every term in (8.2), we get the following estimate

[TABLE]

where

[TABLE]

We apply the Gronwall lemma to (8.15), and by using the assumption (8.5) we obtain

[TABLE]

Taking into account that $Q_{p}$, the production rate (8), is a totally bounded function with period $\bar{T}$, we can easily observe that (8.2) yields

[TABLE]

and

[TABLE]

where

[TABLE]

Hence, by using (8.2), (8.2) and (8.2) we finally get

[TABLE]

which concludes the proof.

∎

9 Numerical simulations

This section is devoted to the validation of the main results obtained in this paper by performing numerical simulations. The computational parameters are based on magnetic resonance measurements. At locations where measured data are not available, lengths and diameters are estimated by combining literature data with measured and computed flows. All the physical constants that appears in our cerebrospinal fluid model are shown in Table 1 but, first of all, we have to remember that, for the sake of simplicity, we collected them in the following way:

[TABLE]

Moreover we know that the choroid plexus produces CSF at a rate, $Q_{p}$, of $0,32$ cm3/min [1] and from medical literature [3], [1] and published data about ventricular pulsation [4], [9] we deduce that the amplitude of choroid expansion, $a(t)$, is in the range of $1,29-1,55$ cm3/min at each cardiac cycle.

In the inputs list for the simulations an important role is also played by the brain tissue pressure, $\tilde{P}$, that is assumed to be at venous pressure levels ($\sim 10$ mmHg).

As we anticipated in section 2, in this model we assume that the cross sectional area, $A$, is affected only negligibly by the pressure variation and represents for us a constant that we choose in the range of $\sim 3-4$ mm2 according to the clinical data and the compartment section.

Since the problem has an axisymmetric structure, in the numerical approach the model geometry is created in two dimensions. A first order upwind scheme with centered differences is employed for the spatial discretization of the governing equations and for the temporal discretization we adopt a forward Euler scheme with a time step size of $5\times 10^{-3}$.

The boundary conditions involved in numerical simulations are selected by taking into account the conditions (3.1) stated in Theorem 3.1.

For the initial data we proceed in the following way: the pressure at the initial time, $P_{0}(z)$, is implemented exactly as required by the conditions (3.3) and (3.1), while, for the velocity flux, $u$, and the tissue displacement, $\eta$, we choose

[TABLE]

which satisfy the condition condition (3.9).

In Figure 4 we can observe in a section of length $L=1$ the profiles of the pressure (red), the tissue displacement (green) and the velocity flux (blue) at three different instants up to the final time $t=1$.

Remark 2*.*

Because of the no-slip condition, the fluid particles in the layer in contact with the surface of the compartment come to a complete stop. This layer also causes the fluid particles in the adjacent layers to slow down gradually as a result of friction. In general, to make up for this velocity reduction, the velocity of the fluid at the core of the compartment has to increase to keep the mass flow rate through the compartment constant. As a result, a velocity gradient should develop along the compartment. In our simulations we can observe, instead, a reduction in the velocity gradient at the midsection of the compartment modelized as a pipe. This is due to the fact that we are neglecting the compartment cross section variation which, otherwise, would concur to smooth the velocity and pressure profiles and to increase the velocity gradient without loss of energy along the compartment.

Now we want to show what happens when the global existence condition of the solution for the system (3.1) in Theorem 3.2 are violated. In order to do that we fix the following initial data

[TABLE]

and for the pressure we keep the initial condition given by (3.3) and (3.1) in Theorem 3.1 because it automatically inherits the properties of the other initial data.

The Figure 5 shows a velocity shock wave after a single iteration and consequently a blow up in pressure and displacement at the second iteration. This result is exactly what we expected according to the previous detailed analysis of the model.

Finally we can conclude that the simulations are in good agreement with the proof of a global solution for the system (3.1), as well as with the real cerebral phenomena modelized.