Matteo Aletti

Jean-Frédéric Gerbeau

Damiano Lombardi

A simplified fluid-structure model for arterial flow. Application to retinal hemodynamics. 1

Keywords: fluid-structure interaction, blood flow, fibers

A simplified fluid-structure model for arterial flow.

Application to retinal hemodynamics

Introduction

Fluid-structure interaction plays an important role in the cardiovascular system. In many situations, complex nonlinear models that include large displacements and deformations have to be considered. This is, for example, the case for valve simulation [1,[START_REF] Van Loon | A combined fictitious domain/adaptative meshing method for fluid-structure interaction in heart valves[END_REF][START_REF] Van Loon | A fluid-structure interaction method with solid-rigid contact for heart valve dynamics[END_REF][START_REF] Astorino | Fluid-structure interaction and multi-body contact. Application to aortic valves[END_REF] or in the aorta [START_REF] Bazilevs | Isogeometric fluidstructure interaction analysis with applications to arterial blood flow[END_REF][START_REF] Crosetto | Fluid-structure interaction simulation of aortic blood flow[END_REF][START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF]. It is well-known that these simulations are very demanding, and in spite of the progress achieved in recent years ( [START_REF] Fernández | Fully decoupled timemarching schemes for incompressible fluid/thin-walled structure interaction[END_REF][START_REF] Bukač | A modular, operator-splitting scheme for fluid-structure interaction problems with thick structures[END_REF][START_REF] Burman | Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus robin-robin coupling[END_REF] to name but a few), they remain challenging and the subject of active research. In this paper, we consider those situations where it is assumed that the problem under study can be addressed using simplified approaches. The idea is to radically simplify the solid model in order to replace the full fluidstructure problem by a fluid problem with non-standard boundary conditions at the fluid-structure interface. Various approaches have been recently proposed in this direction [START_REF] Figueroa | A coupled momentum method for modeling blood flow in three-dimensional deformable arteries[END_REF][START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF][START_REF] Colciago | Comparisons between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics[END_REF][START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF].

In [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF], F. Nobile and C. Vergara started from a Koiter linear shell model and neglected the flexural terms. After discretization, the resulting fluidstructure equations are reduced to a fluid problem with Robin boundary conditions. In this approach, the fluid domain was moving. In [START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF], O. Pironneau further simplified this approach by fixing the fluid domain, introducing a zero order transpiration boundary condition, and by assuming that the curvature of the artery was constant. More precisely, whatever the geometry of the vessel is, the stiffness term is always computed as if the vessel were a cylinder. With these simplifications, the authors were able to perform a comprehensive mathematical analysis of the problem [START_REF] Chacon-Rebollo | Analysis of a Simplified Coupled Fluid-Structure Model for Computational Hemodynamics[END_REF]. In [START_REF] Figueroa | A coupled momentum method for modeling blood flow in three-dimensional deformable arteries[END_REF], A. Figueroa et al. also assumed that the computational domain was fixed and used a zero order transpiration boundary condition. The structural model was derived assuming homogeneity throughout the thickness. Compared to the two previous approaches, this one requires adding new degrees of freedom to the fluid problem. This drawback is, however, counterbalanced by the fact that the resulting model is more stable on real geometries featuring variations of curvature, according to [START_REF] Colciago | Comparisons between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics[END_REF] where an extensive comparison was proposed.

The main focus of the present study is the simulation of the autoregulation of blood flow in the retinal arteries. This phenomenon is very important since defective autoregulation may play a role in many retinal diseases, including glaucoma which is the second leading cause of blindness worldwide [START_REF] Arciero | Theoretical analysis of vascular regulatory mechanisms contributing to retinal blood flow autoregulation[END_REF]. Autoregulation consists of an active change of the artery diameter in response to a change in the mean perfusion pressure. This is clearly a fluid-structure interaction problem, but it is typically a case in which a full structural model does not seem necessary, at least to render the basic phenomenon, which is a slight contraction of the vessel wall. The application of the model proposed in this article to the autoregulation of blood flow is presented in more detail in [START_REF] Aletti | Modeling autoregulation in three-dimensional simulations of retinal hemodynamics[END_REF].

In our approach, we choose to keep the computational domain fixed, as in [START_REF] Figueroa | A coupled momentum method for modeling blood flow in three-dimensional deformable arteries[END_REF], and we adopt transpiration boundary conditions. Nevertheless, the phenomenon of autoregulation cannot be addressed with the zero order transpiration formula usually adopted in the literature. Our model will therefore be based on first order transpiration. For the structure, we will start from a Koiter shell model without flexural terms as in [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF]. But, as noted in [START_REF] Colciago | Comparisons between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics[END_REF], this leads to a simplified model that may be unstable in real geometries due to inaccuracies in the curvature obtained from medical images. Even when the curvatures are computed accurately, this model, called inertial-algebraic, may be inaccurate in geometries with a locally flat fluid-structure interface. Note that the model proposed in [START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF] does not suffer from this problem since it assumes everywhere a cylinder-like geometry in the terms involving the curvature. In our model, we will show that the introduction of fibers allows us to overcome this problem thanks to the presence of a surface Laplace operator.

In summary, the main features of our simplified fluid-structure model are the following: it mainly consists of a fluid problem on a fixed domain, with Robin-like boundary conditions, which makes it insensitive to the added-mass effect; it takes into account in a simplified manner the presence of fibers in the solid; it is less sensitive than others to strong variations or inaccuracies in the curvatures and, as a consequence, it remains robust in the presence of flat regions in the surface.

The structure of the article is as follows: in Section 2 the structure model is proposed; in Section 3 the fluid model and the coupling are presented; Section 4 deals with numerical discretization issues and Section 5 contains numerical illustrations, including an initial simulation in a retinal arteries network.

Structure model

In this section, a simplified structure model is introduced to describe the dynamics of the wall. Similarly to other studies presented in the literature, the starting point is the Koiter thin shell model (see [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF][START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF]). The resulting model aims to render the motion of a thin shell with one or several fiber layers. When the kinematics of the fibers is considered, it is relevant to keep second order terms in their deformation because they have an important role in the stability of the model (the gradient terms in equation ( 5) below). As a consequence, to be consistent with the approximation made for the fibers, second order terms will also be kept in the shell model. This leads to the inclusion of non-linear terms in the shell model.

Notation

Let Γ be the reference position of the structure and a smooth mapping φ defining its position:

φ : ω ⊂ R 2 → Γ ⊂ R 3 , φ = φ(ξ 1 , ξ 2 ), ∀(ξ 1 , ξ 2 ) ∈ ω. Let (a 1 , a 2 ) be the local covariant basis given by a α = ∂ α φ = ∂φ ∂ξα , α = 1, 2.
In what follows, Greek letters for the indices take values in {1, 2} and latin letters in {1, 2, 3}. The normal unit vector is defined as a 3 = a 1 ×a 2 |a 1 ×a 2 | . Let A and B be the matrix representations of the first and second fundamental forms associated with the reference configuration Γ and let S = A -1 B be the representation of the shape operator. The entries of A and B are respectively given by a αβ = a βα = a α • a β and b αβ = b βα = a 3 • ∂ α a β . The entries of A -1 are denoted by a αβ , thus a ασ a σβ = δ αβ . Given a tensor M = (m σβ ), the entries of A -1 M are denoted by m α β = a ασ m σβ . The surface covariant derivative of a vector field q : ω → R 3 is denoted by

q α||β = ∂ β q α -Γ σ αβ q σ -b αβ q 3 and q 3||β = ∂ β q 3 + b σ β q σ , (1)
where Γ σ αβ are the Christoffel symbols. The covariant gradient of a scalar field q : ω → R is denoted by ∇ c q = (∂ α q) α=1, [START_REF] Van Loon | A combined fictitious domain/adaptative meshing method for fluid-structure interaction in heart valves[END_REF] . In what follows, we denote by u, v the standard L 2 (Ω) inner product and by u, v ω the L 2 (Γ)-scalar product w uv √ a dξ where a = det(A).

Nonlinear Koiter shell model

The equations for the Koiter shell model are introduced following [18]. The hypotheses are the following:

• the displacement of the structure is parallel to the normal of the reference configuration;

• the bending terms are negligible;

• the material is linear, isotropic and homogeneous.

As a consequence of these assumptions, only the membrane part of the Koiter model is considered, the shell deformation is described by the change of metric tensor G and the stress is linear in the deformation. The tensor G is a function of the displacement field η and reads:

g αβ = 1 2 (η α||β + η β||α ) + a ij η i||α η j||β ,
with a 3β = a α3 = 0, a 33 = 1. The constitutive law for the stress-strain relationships is expressed by means of the elastic tensor E, whose contravariant components read:

E αβστ = 4λ s µ s λ s + 2µ s a αβ a στ + 2µ s a ασ a βτ + 2µ s a ατ a βσ , where λ s , µ s are the Lamé coefficients of the structure.

The equilibrium configuration for the shell (see [18] for an extensive discussion) is the stationary point of the energy functional:

ψ κ (η) = 1 2 ω E αβστ g στ (η)g αβ (η) h κ √ a dξ - ω f • η h κ √ a dξ,
where h κ is the Koiter shell thickness and f are the external forces. Given a test field χ, defined in a suitable functional space (according to the boundary conditions of the structure), the equilibrium equations in weak form are obtained by:

Ψ κ (η, χ) := δ η ψ κ (η), χ ω = 0, (2)
where δ η denotes the Fréchet derivative with respect to η.

The form Ψ κ is specialised for the present case according to the assumptions. By using the hypothesis of a pure normal displacement, i.e. η = ηn, and Eq.(1), the covariant components of the change of metric tensor G reduce to:

g αβ = -b αβ η + 1 2 a στ b σα b τ β η 2 + 1 2 ∂ α η∂ β η.
After some algebra (see the details of the computation in Appendix A), the form reads:

Ψ κ (η, χ) := 2E 1 -ν 2 ω c 1 η -3c 2 η 2 + 2c 3 η 3 χ -2∇χ T (C 1 η + C 2 η 2 )∇η+ -∇ T η [(C 1 + 2C 2 η) χ] ∇η + 1 2 ∇η T A -1 ∇η ∇ T χA -1 ∇η h κ √ a dξ- ω f n • χ h κ √ a dξ, (3)
where E is the Young modulus of the material, ν the Poisson coefficient, the constant tensors (C j ) and the coefficients (c k ) are expressed as a function of the mean and Gauss curvatures (respectively ρ 1 and ρ 2 ) and the Poisson ratio as follows:

c 1 = 4ρ 2 1 -2(1 -ν)ρ 2 , (4)
c 2 = 4ρ 3 1 + (ν -3)ρ 1 ρ 2 , c 3 = 4ρ 4 1 -4ρ 2 1 ρ 2 + 1 2 (1 + ν)ρ 2 2 , C 1 = νρ 1 I + 1 2 (1 -ν)S A -1 , C 2 = νρ 2 1 I + 1 2 (1 -ν)S 2 A -1 .
Several remarks are in order. First in the work presented in [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF][START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF] only the linear term in (3) is kept (the one multiplied by c 1 ). In this case, the shell behaves like a linear spring, whose stiffness constant depends upon the local curvatures. The non-linearities introduce two contributions: a non-linear spring and a non-linear membrane part.

Fiber layer

In this section, the equations for a generic fiber layer are detailed. The main hypotheses are the following:

• the energy of the shell and of the fiber layer sum up;

• from a kinematical point of view, the fibers are perfectly attached to the shell;

• the fiber is characterized by an affine stress-strain constitutive law.

The second hypothesis implies that the deformation of the fibers equals the deformation of the underlying shell structure in the direction of the fibers. Let w ∈ T x (Γ) be a unitary vector belonging to the tangent space of Γ defined in the point x ∈ Γ. The deformation of the fiber in the w direction can thus be written as:

ε 1D = w T Gw = -d 1 η + d 2 2 η 2 + 1 2 ∇η T P w ∇η, (5)
where the scalar coefficients d j and the projector P w are defined as d 1 = w T Bw, d 2 = w T BSw, P w = w ⊗ w. Note that the d j may be negative , for instance on a cylinder with outward normal. A constitutive stress-strain relationship for the fibers is assumed of the form:

σ 1D = k 0 + k 1 ε 1D ,
so that k 0 represents the pre-stress of the fiber and k 1 is the linear elastic coefficient.

Let w be the fraction of the total number of fibers aligned with the direction w and h f the thickness of the fibers layer. The elastic energy of the fibers aligned in the direction w is expressed in the form:

ψ w (η) = 1 2 ω w [k 0 + k 1 ε 1D (η)] ε 1D (η) h f √ a dξ + ω r w h f √ a dξ,
where r w represents the potential energy of a force acting on the fibers aligned with the direction w.

The equilibrium equations are introduced in weak form, as the scalar product of the Fréchet derivative of the energy with a test function:

Ψ w (η, χ) = ω w ∇χ T k 0 + k 1 -d 1 + d 2 2 η 2 + k 1 4 W P w ∇η+ w k 0 (-d 1 + d 2 η) + k 1 -d 2 1 η - 3d 1 d 2 2 η 2 + d 2 2 2 η 3 + k 1 2 (-d 1 + d 2 η)W χ+ (δ η r w )χ h f √ a dξ, (6)
where W = ∇η T P w ∇η . Remark that the contribution of the first line is of the membrane type, whereas the second line contains algebraic terms in the test function and hence it renders a non-linear spring-like behavior. When η = 0, the reference configuration is the equilibrium configuration only if the stress exerted by the fibers due to their pre-stress is balanced by the underlying shell. This, in weak form, can be written as:

ω (-w k 0 d 1 + δ η r w )χ h f √ a dξ = 0,
that holds for arbitrary test functions χ, hence

r w = w k 0 d 1 η.
To conclude this section, let us remark that, in general, the fibers are not parallel to only one direction. Consider two linearly independent unitary vectors w, v ∈ T x (Γ) and the associated fiber densities w , v defined in each point of Γ. In such a case we can simply sum the two associated energies ψ w and ψ v .

Equations for the structure dynamics

The equations for the structure dynamics are obtained by adding the inertia terms to the elastic contributions highlighted in the previous sections. In particular the dynamics equations in weak form can be written as:

ω ρ s h s (∂ 2 tt η)χ √ a dξ + Ψ κ (η, χ) + Ψ w (η, χ) + Ψ v (η, χ) = 0,
where Ψ κ , defined in Eq.( 3), represents the contribution of the shell and Ψ w , Ψ v , defined in Eq.( 6), represent the contribution of the fibers aligned in the directions w and v respectively. The total thickness of the structure is denoted by h s and its density is denoted by ρ s . In the following, for the sake of simplicity, the overall contribution of the structure will be denoted as:

Ψ s = Ψ κ + Ψ w + Ψ v . ( 7
)

Fluid-structure coupling

The fluid is governed by the incompressible Navier-Stokes equations:

ρ f (∂ t u + u • ∇u) = ∇ • σ f in Ω t ∇ • u = 0 in Ω t
where u is the velocity, ρ f is the fluid density and σ f = µ f (∇u+(∇u) T )-pI is the fluid stress tensor, where µ f is the dynamic viscosity and p is the pressure. The domain Ω t is, in general, time-dependent, since the wall is an elastic structure which is moving because of the interactions with the fluid. We denote by Ω a fixed reference domain. We normalize both the equations for the structure and the fluid by ρ f and we introduce the kinematic viscosity ν f = µ f ρ f . The quantities E, k 0 , k 1 , f s , ρ s and p are also divided by ρ f , but for the sake of simplicity their notation is not changed.

The boundary ∂Ω t is subdivided into two subsets Γ t , the interface between the fluid and the structure, and Σ t , representing the artificial boundaries of the domain where inlet and outlet boundary conditions are enforced on the normal component of the stress tensor.

Two conditions have to be satisfied on the fluid-structure interface Γ t : the continuity of the velocity and the stress. Since the structure displacement is assumed to be parallel to the normal direction, it holds, for As stated in the introduction, the aim of the present work is to set up a simplified fluid-structure interaction model, that can provide solutions at a moderate computational cost. To this end, the structure equations are treated as a boundary condition for the fluid problem and the problem is discretized on a fixed mesh. In order to render the motion of the wall, a transpiration approach is adopted. A zero-th order transpiration was investigated in [START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF][START_REF] Gostaf | Pressure Boundary Conditions for Blood Flows[END_REF][START_REF] Chacon-Rebollo | Analysis of a Simplified Coupled Fluid-Structure Model for Computational Hemodynamics[END_REF], and proves to be satisfactory to study the propagation of pressure waves. However, in view of the application that motivated this work [START_REF] Aletti | Modeling autoregulation in three-dimensional simulations of retinal hemodynamics[END_REF] (hemodynamic autoregulation), it is important to compute the flow variation induced by the wall dynamics. Thus, a first order transpiration condition is considered.

x ∈ Γ, u(I -n ⊗ n)| x+η(x)n(x) = 0 and u • n| x+η(x)n(x) = ∂ t η.
We denote by y ∈ Γ t a point on the actual boundary and by x ∈ Γ the corresponding point on the reference configuration. The mapping y = Φ(x) is written as y = y(x) = η(x) + x. The displacement is assumed to be given by η(x) = η(x)n(x), where n(x) is the outward normal to the reference domain at the point x. An additional assumption is made: the normal is supposed to remain the same during the evolution, that is n(x) = n(y). This assumption, which was also used in [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF][START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF], can be obtained by assuming small deformations since the difference |n(x) -n(y)| is of order one in ∇η. The first order transpiration conditions are obtained through a first order Taylor expansion around the reference configuration η = 0:

u(y) = u(x) + ∇u(x)(y -x) + O(||y -x|| 2 ), u(y) = u(x) + (η∇un)(x) + O(η 2 ), ( 8
)
where the gradient is taken with respect to the x coordinate. The tangential component of the velocity is computed by multiplying Eq.( 8) by (I -n ⊗ n):

(I -n ⊗ n)u(y) = (I -n ⊗ n)(u(x) + η∇u(x)n(x)) + O(η 2 ).
Using the no-slip boundary conditions and neglecting the high-order terms:

(I -n ⊗ n)u(x) = -η(I -n ⊗ n)(∇u(x)n(x)) on Γ.
From Eq.( 8), the normal component of the velocity can be written as

n ⊗ nu(y) = n ⊗ n(u(x) + η∇u(x)n(x)) + O(η 2 ).
Neglecting high order terms, the continuity of the normal velocity gives:

∂ t η = u • n + η∇un • n.
We make the additional simplifying assumption that the viscous part of

σ f nn | x+η(x)n(x) is negligible compared to the pressure, i.e, σ f nn | x+η(x)n(x) = -p| x+η(x)n(x) . Then, the value of σ f nn | x+η(x)n(x)
is approximated by using a first order Taylor expansion on p

σ f nn | x+η(x)n(x) = -p -η∇p • n
The equations for the coupled system are written in weak form, on a fixed reference frame. Let v, q, χ, w be test functions defined in suitable functional spaces according to the boundary conditions of the problem. In particular let (u, v) ∈ V ⊂ H 2 (Ω), let (p, q) ∈ M ⊂ H 1 (Ω) and (η, χ) ∈ H 1 (Γ), w ∈ H 1 (T(Γ)). Then:

               ∂ t u, v + c(u; u, v) + a(u, v) + b(p, v) = 0 in Ω, t > 0 ∇ • u, q = 0 in Ω, t > 0 ρ s h s ∂ 2 tt η, χ ω + Ψ s (η, χ) + p ref , χ ω = p + η∇p • n, χ ω on Γ ∂ t η, χ ω = u • n + η∇un • n, χ ω on Γ (I -n ⊗ n)(u + η∇un), w ω = 0 on Γ. ( 9
)
The forms a, b, c read:

a : V × V → R, a(u, v) = ν f (∇u + ∇u T , ∇v) Ω ∀(u, v) ∈ V × V b : M × V → R, b(p, v) = -(p, ∇ • v) Ω ∀(p, v) ∈ M × V c(w) : V × V → R, c(w; u, v) = (w • ∇u, v) Ω ∀(u, v) ∈ V × V 4. Numerical discretization

Time discretization

The approximation of the generic quantity f at time t k = k∆t is denoted by f k . The following notation is introduced:

∇ n f = ∇f •n, ∇ nn f = ∇f n•n, f n = f • n.
An implicit-explicit time discretization is adopted in order to avoid the resolution of non-linear problems. As regards the tangential velocity boundary condition, the Taylor contribution accounting for the motion of the wall is taken at the current time step:

(I -n ⊗ n)u k+1 , w ω = -η k (I -n ⊗ n)∇u k n, w ω on Γ. ( 10
)
The normal component of the velocity, which is directly related to the wall displacement η is discretized by adopting a similar strategy:

η k+1 , χ ω = (1 + ∆t∇ nn u k )η k , χ ω + ∆t u k+1 n , χ ω . (11)
The structure equation is discretized as follows:

ρ s h s η k+1 -2η k + η k-1 ∆t 2 , χ ω +Ψ s (η k+1,k , χ)+ p ref , χ ω = p k+1 +η k ∇ n p k , χ ω , (12)
where Ψ s is the energy of the structure defined in Eq. [START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF]. The non-linear terms in Ψ s are treated in a semi-implicit way detailed in Appendix B (from equations (B.1) to (B.3)). After linearization Ψ s is replaced by the sum of a bilinear form of η k+1 and χ and a linear functional of χ. In order to eliminate the current displacement from the the bilinear form, every istance of η k+1 is replaced by Eq.( 11) to obtain:

Ψ s (η k+1,k , χ) = Φ k+1 (u k+1 n , χ) + φ k (χ). ( 13
)
The bilinear form Ψ s (η k+1,k , χ) is now divided into two contributions: the first one Φ k+1 (u k+1 n , χ) is a bilinear form that depends on the current velocity and the second one φ k (χ) is a linear functional of χ where old quantities appear as parameters (see Eq. (B.4) for detailed expressions).

By injecting Eq.( 13) into Eq. [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF] and by collecting all the force terms in one functional, the following is obtained:

ρ s h s ∆t u k+1 n , χ ω + Φ k+1 (u k+1 n , χ) -p k+1 , χ ω = F(χ; η k , p k , u k n ). ( 14
)
The details of the expression for F are reported in Eq. (B.5). Remark that the left-hand side is now made up of the unknowns of the fluid problem, and the right-hand side is computed by using the values of the variables at the previous time step. In this way, the motion of the structure has been implicitly embedded as a boundary condition for the fluid problem. In addition, as the acceleration of the structure is treated implicitly, possible numerical instabilities due to the added-mass effect [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF] are avoided.

The linear part of the Navier-Stokes equations is discretized by an implicit Euler scheme, and the convective term in a semi-implicit way:

c(u k ; u k+1 , v) = u k • ∇u k+1 , v .
Some comments on the stability of the scheme are in order. In the work by Nobile and Vergara [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF], stability was proven for a similar model where the fluid problem was solved on a moving domain and a linear Koiter model was embedded into the fluid boundary conditions. In the work by Pironneau [START_REF] Pironneau | Simplified Fluid-Structure Interactions for Hemodynamics[END_REF] a complete mathematical analysis of a similar method was done on a fixed domain where the coupling was defined by a zero order transpiration. The proof of stability of the formulation proposed in the present work seems to be more complicated because of the first order transpiration that results in a non-linear mixed boundary condition on the surface Γ for the continuous problem.

Application of the boundary conditions

Eq.( 9) is discretized using P1-P1 finite elements with an SUPG stabilization. In this section, several remarks on the imposition of the boundary conditions for the system are presented.

The boundary conditions are separated into a tangential boundary condition of Dirichlet type for the velocity and a generalised Robin boundary condition in the normal direction. In order to impose the conditions via penalization and to avoid cancellation errors, a unique vector-type boundary condition is written. Let z = w + χn be a test function , with w • n = 0. It holds:

-

p k+1 n, z ω + (u k+1 , z) = R(z),
where is a bilinear form and R is a linear functional, deduced from Eq.( 14) and Eq.( 10) (see Appendix B). The implementation of this boundary condition has to be done carefully when working with real geometries. In particular, we observed that with normals constant per element, spurious contributions could appear on the tangential velocity (see Remark below). It is therefore desirable to reinterpolate the normals on the P1 finite elements nodes.

Remark 1. Let us explain the trouble that may come from a piecewise-P0 approximation of the normal when a penalization method on the weak form is adopted. Consider, for instance, the case of a Dirichlet boundary condition u = (I -n ⊗ n)g, where g is a generic vector defined on Γ. When this condition is imposed via penalization the weak formulation reads:

Γ (u -T g) • w ≈ 0
where T = (I -n⊗n) is the projector in the tangent space, and "≈ 0" means "of the order of the inverse of the penalization parameter". If a piecewise-P0 approximation to the normal field is used, the projector operator T el is defined element-wise and is discontinuous. The k-th component of this equation is written by introducing the basis functions used to discretize the problem, namely ϕ, providing:

Γ (u k,i ϕ i -(T el kr g r ) i ϕ i )w k,j ϕ j ≈ 0,
where (T el k,j g j ) i = Γ T el k,r g r φ i . This equation leads to u k,i ≈ (T el k,r g r ) i , and hence it does not guarantee that n el s n el k u k,i ≈ 0, for s = 1, 2, 3. On the other hand, if the piecewise-P1 normal and the corresponding operator T are used, the quantity T g is computed for each node and the following identity holds:

Γ (u k,i ϕ i -T kr,i g r,i ϕ i )ϕ j ≈ 0, l = 1...3
which leads to u k,i ≈ T kr,i g r,i . The velocity satisfies the condition n s n k u k,i = 0, for s = 1, 2, 3. This is particularly relevant for the accuracy in the computation of the normal component of the velocity and thus of the displacement field. ♦

Numerical testcases

In this section, three numerical experiments are presented where the different features of the proposed method are tested.

In the first test case, a pressure wave in a circular cylinder is simulated, in a similar setting to that in [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF] problem is solved in a different geometry, characterized by the presence of a flat region. Finally, the proposed method is applied to an image-based geometry describing a part of the retinal arteriolar network.

Pressure wave in a cylinder

In this test case, the Stokes equations are solved in a cylindrical domain (L = 6cm, R = 0.5cm). Two different configurations are compared, one with fibers and one without. The parameters are reported in Tab.1. The values for the Young modulus and for the fibers pre-stress are chosen such that, after linearization, the spring coefficient (

h k E (1-ν s2 )R 2 + h f k 0 w R 2
) coincides for the two configurations. The remaining structure parameters are: ρ s = 1 g/cm 3 , h s = h κ + h f and ν s = 0.5. For the sake of simplicity, we set k 1 = 0. Fibers are assumed to be aligned with the principal directions of curvature: w is the circumferential direction and v the longitudinal direction. The fibers' fraction is chosen such that fibers are mostly aligned with the circumferential direction ( w = 0.9, v = 0.1).

For such a simple geometry, the fibers' directions and curvatures are known analytically. However, in view of applying the method to realistic geometries, where this information is in general not available, we assume that the directions and the curvatures are unknown, and we estimate them numerically from the computational mesh.

Fig. 1 presents curvature estimations. On the left, the Gaussian curvature (analytically equal to zero) and the mean curvature (analytically equal to one) are plotted against the longitudinal coordinate of the cylinder. Notice the presence of numerical oscillations in the computation of the curvature. This problem may have negative effects on other simplified FSI models, as observed in [START_REF] Colciago | Comparisons between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics[END_REF]. In the proposed model, it is partially overcome by the addition of the fiber layer as shown below. In the right panel of Fig. 1, the estimated principal curvature directions are shown.

The fluid parameters are: ρ f = 1 g/cm 3 and µ f = 0.035 cm 2 /s. The boundary conditions at the inlet and at the outlet are σ(u, p)n = -pn. A pressure equal to zero is assigned at the outlet (p = 0) and at the inlet p = 5000 dyn/cm 2 t ≤ 0.005s 0 t > 0.005s.

In Figures 2 and3 the displacement is shown against the longitudinal coordinate, for two different time instants: t = 0.004 and t = 0.012, respectively. First, there is a difference between the non-linear model and its linear version, both with and without the fiber layer. The pressure wave amplitude is lower for the non-linear model, and the wave velocity is higher. This results from the fact that the non-linear model is characterized by a greater stiffness (due to the non-linear spring contribution).

By comparing the models with and without fibers, it is interesting to notice the regularizing effect of the fibers. The difference in the peak amplitude and in the wave velocity is due to the difference in inertia and stiffness.

Pressure wave in a flat cylinder

The effect of adding a fiber layer to the structure model is even more visible in presence of flat regions in the surface. In realistic geometries, a locally flat region may occur for several reasons including a lack of precision 1) for t = 0.004. Red lines refer to the case without fibers, while black ones to the case with fibers. Solid lines refer to the full model and dashed lines to the corresponding linear version. Figure 3: Displacement in the longitudinal direction for different structure parameters (see Table 1) for t = 0.012. Red lines refer to the case without fibers, while black ones to the case with fibers. Solid lines refer to the full model and dashed lines to the corresponding linear version. in the segmentation process. This is why the fluid-structure interaction in a flattened cylinder (see Figure 4) is investigated.

For this geometrical setting the inertial-algebraic model proposed in [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF] cannot be used. The main reason is that this model reduces to a structure that behaves pointwise like a spring, whose stiffness constant depends on the curvatures. When both the curvatures vanish, the spring coefficient is zero. Therefore, the displacement becomes very large, leading to unphysical solutions. Adding the non-linear part of the Koiter shell model mitigates this behavior.

In Fig. 5 the displacement is shown in the longitudinal direction in the center of the flat region for t = 0.012. Notice that, in the flat region, the structure without fibers is characterized by a significantly large displacement. This is due to the fact that, in this region, the structure behaves like a membrane, whose stiffness is low. In Fig. 6 the same curves are reported for the displacement in the non-flat region. In this portion, the behavior of the structure is similar to that found in the first test case.

Application

In the last test case the model is applied to a realistic geometry. The network represents the inferior temporal arteriole network in the human retina. The original image was taken in the Drive dataset [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF]. The vasculature has been segmented by applying the imaging methods presented in [START_REF] Al-Diri | An active contour model for segmenting and measuring retinal vessels, Medical Imaging[END_REF][START_REF] Caliva | A new tool to connect blood vessels in fundus retinal images[END_REF]. The 2D data have then been expanded into a 3D-network by assuming a circular section and projecting the results onto a sphere representing the eye. The typical diameter of this network varies between 70µm and 160µm.

Several physical parameters appear in the equations of the structure and in that of the fluid. The parameters are chosen to be comparable to those found in [START_REF] Arciero | Theoretical analysis of vascular regulatory mechanisms contributing to retinal blood flow autoregulation[END_REF][START_REF] Guidoboni | Intraocular pressure, blood pressure, and retinal blood flow autoregulation: A mathematical model to clarify their relationship and clinical relevance[END_REF][START_REF] Pournaras | Regulation of retinal blood flow in health and disease[END_REF] and they are a realistic and representative set of parameters for retinal arterioles. The parameters used for the arteriole wall are ρ s = 1g/cm 3 , E = 0.05M P a, ν = 0.5, h κ = 5µm, h f = 20µm, k 0 = 0.4M P a, k 1 = 0, w = v = 0.5. For the sake of simplicity, the blood is assumed to be Newtonian, even if for this kind of vessel a non-Newtonian model would be more appropriate. We used ρ f = 1g/cm 3 , µ f = 0.03cPa. Autoregulation is not included in this test. The reader interested in test cases including autoregulation is referred to [START_REF] Aletti | Modeling autoregulation in three-dimensional simulations of retinal hemodynamics[END_REF] Pressure conditions are applied at the inlet: The mean value over time of the incoming pressure is 50mmHg. At the outlets of the domain, the downstream circulation is connected to a venous pressure of 20mmHg via Windkessel compartments. More precisely, we use an RCR model where the parameters are R prox = 6. 10 8 Pa s cm -3 , R dist = 6.10 9 Pa s cm -3 and the capacitance is 1.67 10 -8 cm 3 s -1 Pa-1. On the lateral surface an external pressure, modeling the intra-ocular pressure, of 15mmHg is also imposed.

P in (t
The displacement of the arterial wall is shown in Fig. 7 at the systolic peak and in Fig. 8 at the end of the diastolic phase.

The fact that the displacement is smaller towards the end of the networkis mainly due to two reasons: the pressure at the end of the the network is significantly lower (the pressure drop is between 5 and 10 mmHg) and the diameter of the vessels is smaller. Indeed, the spring contribution to the structure is characterized by a coefficient that, at the first order approximation in η, depends on the 0.0003 0.0006 0.0009 0.0012 0.0000 0.0015

Normal Displacement inverse of the vessel radius squared (as in the Laplace law). This can be quantitavely observed in Fig. 9 by comparing Point 2 and Point 6 (orange and black curves) or by comparing Point 3 and Point 4. In these two cases, we observe similar pressure curves (dashed lines) and different displacement (solid lines).

In the same figure we also observe the propagation of the pressure wave through the vessel network: by comparing Point 1 with Point 5 (red and green curves) we see that the position of the pressure (and displacement) peak is delayed in Point 5 with respect to Point 1. On the other hand the amplitude of the wave is reduced, which is compatible with the fact that no pulsation is observed at the capillary level.

From a practical viewpoint, it is interesting to note that these fluidstructure results have been obtained at a computational cost similar to a fluid problem. It would be interesting to compare the results with those provided by more complex approaches, as was done in [START_REF] Colciago | Comparisons between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics[END_REF] for other simplified models. This will be the object of future work.

Aknowledgements

This research was made possible by a Marie Curie grant from the European Commission in the framework of the REVAMMAD ITN (Initial Training Research network), Project number 316990.

[1] T. B. Le, F. Sotiropoulos, Fluidstructure interaction of an aortic heart valve prosthesis driven by an animated anatomic left the following is obtained:

ψ κ (η) = 1 4 Γ h κ E αβστ g στ g αβ √ a dξ = 1 4 Γ h κ 4λ s µ s λ s + 2µ s tr(A -1 G) 2 + 4µ s tr((A -1 G) 2 ) √ a dξ = 1 4 Γ h κ 2Eν 1 -ν 2 tr(A -1 G) 2 + 2E 1 + ν tr((A -1 G) 2 ) √ a dξ = 1 2 Γ h κ E 1 -ν 2 ν tr(A -1 G) 2 + (1 -ν) tr((A -1 G) 2 ) √ a dξ, (A.3)
where the relationships between the Lamè coefficients µ s , λ s and the Young modulus and the Poisson ratio are used, namely:

λ s = Eν (1 + ν)(1 -2ν) , µ s = E 2(1 + ν) .
The expression of the change of metric tensor is rewritten, by considering that the displacement field is, by hypothesis, aligned with the outward normal:

g αβ = -b αβ η + 1 2 s τ α b τ β η 2 + 1 2 ∂ α η∂ β η,
where s τ α denotes the components of the matrix representation of S.

In order to expand the terms in Eq.(A.3) in terms of η, we write the components of the tensor A -1 G

g δ β = (A -1 G) δβ = a δα g αβ = -a δα b αβ η + 1 2 a δα a στ b σα b τ β η 2 + 1 2 a δα ∂ α η∂ β η.
By introducing the notation ∂ α = a ασ ∂ σ , g δ β it simplifies to

g δ β = -s δβ η + 1 2 s δσ s σβ η 2 + 1 2 ∂ δ η∂ β η.
After some algebraic calculations the quantities tr(A -1 G) and tr((A -1 G) 2 ) can be written in terms of traces of powers of S, which are, in turn, directly related to the curvatures of the surfaces.

It holds:

tr(A -1 G) = -tr(S)η + 1 2 tr(S 2 )η 2 + 1 2 ∇η T A -1 ∇η tr((A -1 G) 2 ) = tr(S 2 )η 2 -tr(S 3 )η 3 -η∇η T SA -1 ∇η + 1 4 η 4 tr(S 4 ) + 1 2 η 2 ∇ηS 2 A -1 ∇η + 1 4 (∇η T A -1 ∇η) 2 .
(A.4) consider only fibers in the w direction.

Ψ s (η, χ) = 2E 1 -ν 2 ω c 1 η -3c 2 η 2 + 2c 3 η 3 χ -2∇χ T (C 1 η + C 2 η 2 )∇η+ (B.1) -∇ T η [(C 1 + 2C 2 η) χ] ∇η + 1 2 ∇η T A -1 ∇η ∇ T χA -1 ∇η h κ √ a dξ- ω f n • χ h κ √ a dξ + ω w ∇χ T k 0 + k 1 -d 1 + d 2 2 η 2 + k 1 4 W P w ∇η+ w k 0 d 2 η + k 1 -d 2 1 η - 3d 1 d 2 2 η 2 + d 2 2 2 η 3 + k 1 2 (-d 1 + d 2 η)W χ h f √ a dξ.
We use a semi-implicit approach to discretize this form in time. The scheme is the following

Ψ s (η k+1 , χ; η k ) = 2E 1 -ν 2 ω c 1 -3c 2 η k + 2c 3 η k 2 η k+1 χ -2∇χ T (C 1 η k + C 2 η k 2 )∇η k+1 + -∇ T η k C 1 + 2C 2 η k+1 χ ∇η k + 1 2 ∇η k T A -1 ∇η k ∇ T χA -1 ∇η k+1 h κ √ a dξ- ω f k+1 n • χ h κ √ a dξ + ω w ∇χ T k 0 + k 1 -d 1 + d 2 2 η k 2 + k 1 4 W k P w ∇η k+1 + w k 0 d 2 η k+1 + k 1 -d 2 1 - 3d 1 d 2 2 η k + d 2 2 2 η k 2 η k+1 + k 1 2 (-d 1 + d 2 η k+1 )W k χ h f √ a dξ.
This form is linear with respect to χ and affine with respect to η k+1 . The coefficients of the form depend on the current displacement η k . It is useful to split the form into two sub-contributions in order to highlight the bilinear part Ψ s (η k+1 , χ; η k ) = B(η k+1 , χ; η k ) + G(χ; η k ). (B.

2)

The bilinear form B(η k+1 , χ; η k ) has the following expression:

B(η k+1 , χ; η k ) = 2E 1 -ν 2 ω c 1 -3c 2 η k + 2c 3 η k 2 -2(∇ T η k C 2 ∇η k ) η k+1 χ -2∇χ T (C 1 η k + C 2 η k 2 )∇η k+1 + 1 2 ∇η k T A -1 ∇η k ∇ T χA -1 ∇η k+1 h κ √ a dξ+ + ω w ∇χ T k 0 + k 1 -d 1 + d 2 2 η k 2 + k 1 4 W k P w ∇η k+1 + w k 0 d 2 + k 1 -d 2 1 - 3d 1 d 2 2 η k + d 2 2 2 η k 2 + k 1 2 d 2 W k η k+1 χ h f √ a dξ.
On the other hand, the functional G(χ; η k ) has the form

G(χ; η k ) = - 2E 1 -ν 2 ω (∇ T η k C 1 ∇η k )χ h κ √ a dξ- (B.3) ω f k+1 n • χ h κ √ a dξ - ω w k 1 2 d 1 W k χ h f √ a dξ.
We use the strong formulation of Eq. [START_REF] Figueroa | A coupled momentum method for modeling blood flow in three-dimensional deformable arteries[END_REF],

η k+1 = (1 + ∆t∇ nn u k )η k + ∆tu k+1 n ,
to replace η k+1 in Eq.(B.2) obtaining:

Ψ s (η k+1 , χ; η k ) = ∆tB(u k+1 n , χ; η k ) + B((1 + ∆t∇ nn u k )η k , χ; η k ) + G(χ; η k ), which is in the same form as Eq.( 13). We, finally, define:

Φ k+1 (u k+1 n , χ) = ∆tB(u k+1 n , χ; η k ) φ k (χ) = B((1 + ∆t∇ nn u k )η k , χ; η k ) + G(χ; η k ).

(B.4)

The expression for F in equation ( 14) is obtained by combining Eq.(B.4) and Eq. [START_REF] Nobile | An effective fluid-structure interaction formulation for vascular dynamics by generalized robin conditions[END_REF].

F = η k ∇ n p k -p ref - ρ s h s ∆t 2 ((∆t∇ nn u k -1)η k + η k-1 ), χ ω -φ k (χ). (B.5)
On the interface surface Γ those two boundary conditions (Eq.( 10), Eq.( 14)) hold:

(I -n ⊗ n)u k+1 , w ω = -η k (I -n ⊗ n)∇u k n, w ω ρ s h s ∆t u k+1 n , χ ω + Φ k+1 (u k+1 n , χ) -p k+1 , χ ω = F(χ; η k , p k , u k n ).

Let z : Γ → R 3 be a vector-valued function of H 1 (Γ), such that z • n = χ and (I -n⊗n)z = w. We also recall that w •n = 0. By introducing the notation T = (I -n⊗n) and by reorganizing the terms ρ s h s ∆t u k+1 n , χ ω +Φ k+1 (u k+1 n , χ) we get T u k+1 , w ω = -η k T ∇u k n, w ω α k u k+1 n , χ ω + Φ k+1 (u k+1 n , χ) -p k+1 , χ ω = F(χ; η k , p k , u k n ),

where

α k = ρ s h s ∆t + 2E∆th κ 1 -ν 2 (c 1 -3c 2 η k + 2c 3 η k 2 -2(∇ T η k C 2 ∇η k ) + w h f k 0 d 2 + k 1 -d 2 1 - 3d 1 d 2 2 η k + d 2 2 2 η k 2 + k 1 2 d 2 W k and Φ k+1 = 2E∆t 1 -ν 2 ω -2∇χ T (C 1 η k + C 2 η k 2 )∇u k+1 n + 1 2 ∇η k T A -1 ∇η k ∇ T χA -1 ∇u k+1 n h κ √ a dξ+ + ω w ∆t∇χ T k 0 + k 1 -d 1 + d 2 2 η k 2 + k 1 4 W k P w ∇u k+1 n h f √ a dξ.
Finally, we replace w by T z and χ by z • n, we multiply the equation for the tangential component by α k and we sum the result

α k u k+1 , z ω + Φ k+1 (u k+1 n , z • n) (u k+1 ,z) -p k+1 , z•n ω = F(z • n; η k , p k , u k n ) -η k T ∇u k n, T z ω R(z)
.

1

  A final version of this manuscript can be found in Computer Methods in Applied Mechanics and Engineering (2016) Preprint submitted to Computer Methods in Applied Mechanics and EngineeringApril 1, 2016

  The continuity of the normal component of the normal stress gives σ f nn | x+η(x)n(x) = -f s -p ref , where p ref denotes a given external pressure acting on the structure.

Figure 1 :

 1 Figure1: On the left panel the mean (brown) and the Gaussian (blue) curvature estimated from the computational mesh are compared to the exact values: one and zero, respectively. On the right panel the principal directions of curvature are depicted: red arrows represent the local direction of maximum curvature, the blue ones refer to the minimum. On the surface the mean curvature is displayed.

Figure 2 :

 2 Figure2: Displacement in the longitudinal direction for different structure parameters (see Table1) for t = 0.004. Red lines refer to the case without fibers, while black ones to the case with fibers. Solid lines refer to the full model and dashed lines to the corresponding linear version.

Figure 4 :Figure 5 :

 45 Figure 4: The domain of the second test case. One side of the cylinder is flat, i.e. all the curvatures vanish.

Figure 6 :

 6 Figure 6: Displacement in the rounded region. Comparison between the model with fibers (black dashed line linear, black solid line full-model) and the model without (red solid line).

  )[mmHg] = 25.12 sin(πt/0.25) + 45 t ∈ [0, 0.25] 45 t ∈ [0.25, 0.8].

Figure 7 :

 7 Figure 7: Displacement field on the retinal vasculature during the systolic peak.

Figure 8 :Figure 9 :

 89 Figure 8: Displacement field on the retinal vasculature at the end of diastolic phase.

Table 1 :

 1 . In the second numerical experiment, a similar ] h κ [cm] k 0 [dyn/cm 2 ] h f [cm] Structure parameters for the two simulations.

		Koiter layer		Fibers Layer
	E[dyn/cm 2 Without fibers 7.5 • 10 5	0.1	0	0
	With fibers	3.75 • 10 5	0.1	2.78 • 10 5	0.1

Appendix A. Derivation of the nonlinear elastic energy

In this Appendix, the expression of the elastic energy for the non-linear Koiter model is derived, when the simplifying hypotheses for the problem (see Section 2.2) are taken into account.

We refer to Section 2.1 for the differential geometry notation.

The expression of the energy functional, when f = 0,

is written in terms of the elasticity tensor defined on the surface and of the change of metric tensor. The elasticity tensor, E αβστ , depends on the Lamé coefficients λ s and µ s and on inverse of the first fundamental form, A -1 , of the surface. In order to simplify the calculations the term E αβστ g στ (η)g αβ (η) h κ is split into two different contributions:

In the following the dependence of G on η is dropped for sake of compactness in the notation. Since A,B,G are symmetric, the two contributions can be further simplified and written in a more compact form as functions of the tensor A -1 G:

By injecting Eq.(A.1) and (A.2) into the expression of the energy functional,

The traces of powers of S can be written in terms of the curvatures as:

where ρ 1 is the mean curvature 1/2(λ 1 +λ 2 ) and ρ 2 is the Gaussian curvature λ 1 λ 2 . By using Eq.(A.4) we derive the final expression for the elastic energy functional

where the coefficients c j depend on the Poisson ratio and on the curvatures (they are reported in Eq.( 4)).

In order to compute the form Ψ κ in Eq.( 2), the first variation with respect to η of the energy functional has to be computed and tested, to get the weak form, against a function χ belonging to a suitable functional space. The procedure to derive the first variation (Fréchet derivative) involves several integrations by parts on the surface that produce extra terms involving the curvatures. These terms, however, disappear when the first variation is tested against the function χ. The final result is reported in Eq.(3).

Appendix B. Details on the time discretization of the boundary condition

The nonlinear form Ψ s (η, χ) depends on the forms associated to the nonlinear Koiter shell and the fibers (see Eq.( 7)). For the sake of simplicity, we


Geoffrey Bruckert

Denis Vivien

Fabian Docagne

Dr Benoit D Roussel
email: broussel@cyceron.fr

Normalization of Reverse Transcription Quantitative PCR Data During Ageing in Distinct Cerebral Structures

Keywords: RT-qPCR, ageing, normalization, GeNorm, Normfinder

  L'archive ouverte pluridisciplinaire

Introduction

Reverse transcription quantitative PCR (RT-qPCR) is widely used in gene expression analysis because of its simple, reproducible and high-throughput features.

It gives a rapid means of understanding gene variation in various organisms/organs by measuring the expression of target genes compared to reference genes across different samples. On top of that, RT-qPCR has become a low cost method allowing a generalization of its use in laboratories. However, to analyse RT-qPCR data, several factors such as sample amount, RNA integrity or cDNA quality, must be controlled because they affect the quantitative measurement of gene expression [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF][START_REF] Bustin | The need for transparency and good practices in the qPCR literature[END_REF]. Thus, in order to obtain a robust analysis of gene expression, reference genes are used as an internal control to normalize and monitor the variations of expression between samples and reactions [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF][START_REF] Bustin | The need for transparency and good practices in the qPCR literature[END_REF][START_REF] Huggett | Real-time RT-PCR normalisation; strategies and considerations[END_REF]. The expression of these reference genes should not vary in the tissues or cells under investigation, or in response to experimental treatment [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF]. As shown by two surveys realized in 2013, there is still many bias in RT-qPCR studies: the stability of reference genes is not systematically checked, and normalization is generally achieved by the use of a single reference gene [START_REF] Bustin | The need for transparency and good practices in the qPCR literature[END_REF], which is in discordance with the MIQE guidelines published by Bustin and collaborators to avoid misinterpretations of RT-qPCR results [START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF].

Ageing is a physiological process associated with deficits in learning and memory tasks, explained by deficits in key brain-functions, like motor coordination [START_REF] Seidler | Motor control and aging: links to age-related brain structural, functional, and biochemical effects[END_REF], and sensory and cognitive functions [START_REF] Baltes | Emergence of a powerful connection between sensory and cognitive functions across the adult life span: a new window to the study of cognitive aging?[END_REF]. Ageing is also a risk factor for many neurological diseases, including cerebral ischaemia [START_REF] Ganesh | Genetics and genomics for the prevention and treatment of cardiovascular disease: update: a scientific statement from the American Heart Association[END_REF], Alzheimer's disease [START_REF] Reitz | Alzheimer disease: epidemiology, diagnostic criteria, risk factors and biomarkers[END_REF],

Parkinson's disease [START_REF] Hindle | Ageing, neurodegeneration and Parkinson's disease[END_REF], and others. To study gene expression during ageing in neurodegenerative pathologies or in normal brain, RT-qPCR has been widely used.

It has been suggested that RT-qPCR data should be normalized to at least two reference genes which expression has been demonstrated to be stable in the conditions studied [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF][START_REF] Bustin | Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): trends and problems[END_REF]. The most frequently used genes for normalization are glyceraldehyde-3-phosphate dehydrogenase (Gapdh) and beta-actin (ActinB). However, they can display variable expressions in many paradigms. For example, Actinb has been shown to vary in several studies, like in ischaemic/hypoxic conditions [START_REF] Harrison | The use of quantitative RT-PCR to measure mRNA expression in a rat model of focal ischemia-caspase-3 as a case study[END_REF][START_REF] Foldager | Validation of suitable house keeping genes for hypoxia-cultured human chondrocytes[END_REF] and Gapdh is not stable in a model of apoptosis induced by age on neurones [START_REF] Ishitani | Evidence that glyceraldehyde-3-phosphate dehydrogenase is involved in age-induced apoptosis in mature cerebellar neurons in culture[END_REF]. Also, their sequences overlap with highly homologous segmental duplications sequences [START_REF] Ghani | Segmental duplications in genome-wide significant loci and housekeeping genes; warning for GAPDH and ACTB[END_REF].

Based on that, Sieber and collaborators published a useful study in 2010 that demonstrated the need to check the stability of normalization genes in ageing and ischaemia for brain studies [START_REF] Sieber | Inter-age variability of bona fide unvaried transcripts Normalization of quantitative PCR data in ischemic stroke[END_REF]. However, the authors did not discriminate different structures/areas but rather used the whole brain.

It is well known that neurodegenerative diseases/models affect specific brain regions, and gene expression also varies in different structures. By using the GeNorm and the Normfinder methods, we studied the stability of common reference genes, such as ActinB (beta actin), Polr2a (polymerase (RNA) II (DNA directed) polypeptide A), Sdha (succinate dehydrogenase complex, subunit A), Ppox (protoporphyrinogen oxidase),

Hprt (hypoxanthine guanine phosphoribosyl transferase), Rpl13a (ribosomal protein L13a), Hmbs (hydroxymethylbilane synthase), Ppib (peptidylprolyl isomerase B, or cyclophilin B), and GAPDH (glyceraldehyde-3-phosphate dehydrogenase) in four mouse brain structures (cortex, hippocampus, striatum and cerebellum) during ageing. We found that the stability of each gene depends on the studied structure, and that reference genes are specific for each structure in ageing. We thus propose, for each of these brain structures, a set of reference genes usable for ageing studies in mice, and we verified our references genes by measuring the expression of the age-related gene tissue plasminogen activator (tPA) in the hippocampus during ageing.

Results

RT-qPCR efficiency

Quantification cycle Cq was correlated to the starting concentration of the dilution series (from 1:10 to 1:10000; five orders of magnitude). The coefficient of correlation (R²) ranged from 0.977 to 0.998 in the cortex, 0.971 to 0.986 in the hippocampus, 0.977 to 0.993 in the striatum, and 0.975 to 0.994 in the cerebellum.

The efficiency was then calculated as described in the material and methods section and ranged from 82.45% to 126.625% as shown in the table1. Non-normalized Cq values in our experiments were all within the range of the dilution series.

Expression of reference genes in the brain

We measured the raw expression (here called starting quantity) of our nine genes in four cerebral structures: cortex, hippocampus, striatum and cerebellum.

Among the tested genes we observed variations during ageing in the different brain structures. For example, in the cortex, hmbs, Actinb, Ppox, Polr2a, and GAPDH do not vary between 8 weeks, 7 months, 17 months and 22 months in contrast to Ppib, Rpl13a, Hprt and Sdha (Table 2). In the hippocampus, only Ppib, Hprt and Polr2a

show no statistical differences of expression during ageing (Table 2). In the striatum, Ppib, Rpl13a, Actinb, Ppox, Sdha, Polr2a and GAPDH do not show statistical differences of expression during ageing, whereas in the cerebellum only Ppib, Rpl13a, and GAPDH show no statistical differences (Table 2).

These results clearly demonstrate variations in many commonly-used reference genes, but also differences between brain structures during normal ageing.

It also demonstrates that one reference gene is usually not enough, and normalization to multiple reference genes needs to be done for each cerebral structure during ageing. To address this, we used two different methods: by using the QBase+™software that contains the GeNorm algorithm, and by using the Normfinder algorithm.

Determination of the correct reference genes in cerebral structures

Cortex

In the cerebral cortex, multiple analyses of the nine reference genes revealed a GeNorm V value under the 0.15 threshold for the use of 2 genes including Actb, Polr2a, Sdha, Ppox, Hprt, Rpl13a, Hmbs, Ppib, and GAPDH (Figure 5A). The GeNorm M value [START_REF] Hellemans | qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data[END_REF] showed that the two best genes to use as reference genes in the cerebral cortex during ageing are Rpl13a and GAPDH (figure 5B): between 8 weeks and 22 months, these two genes display the same profile of expression in all our studied conditions.

To ensure the choice of the reference genes in the cortex, we also performed a Normfinder analysis that calculates a stability value after the estimation of an intraand intergroup variation (see material and methods [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF]). The ranking of the reference genes (figure 5C) showed that the best genes for normalizing RT-qPCR data in the cortex during ageing are Rpl13a and GAPDH.

Hippocampus

The same method was applied in the hippocampus, aged from 8 weeks to 22 months-old. Again, among the nine selected genes, two paired genes were enough to normalize RT-qPCR data in the hippocampus during ageing, as shown by the GeNorm V value (Figure 6A). According to the GeNorm M value, Rpl13a and Hprt is the best pair of reference genes to normalize RT-qPCR data in the hippocampus during ageing (Figure 6B).

The Normfinder analysis gave Rpl13a and Sdha as the most stable genes (figure 6C). There was a slight difference in the choice of genes between Genorm and Normfinder. However, both analyses ranked Rpl13a as the most stable gene in the hippocampus during ageing.

Striatum

In the striatum of C57BL6/J mice aged from 8 weeks-to 22months-old, analysis of the studied reference genes revealed that 2 genes of the list are enough to normalize RT-qPCR data (Figure 7A). The GeNorm M value showed that Rpl13a

and Ppib are the best candidates for the normalization of RT-qPCR data in the striatum during ageing (Figure 7B), and the Normfinder analysis (Figure 7C) gave the same result.

Cerebellum

Regarding the cerebellum in ageing, the analysis of the GeNorm V value showed that two reference genes are also enough to normalize RT-qPCR (Figure 8A). According to the GeNorm M value, these genes are GAPDH and pPIB (Figure 8B). The Norminder analysis gave a slightly different result, as the two best genes were GAPDH and Hprt (Figure 8c).

Evaluation of reference genes in RT-qPCR experiments

To test our hypothesis, we decided to check the expression of a gene during ageing. As it was previously demonstrated that tissue plasminogen activator (tPA)

was impacted during ageing [START_REF] Roussel | Age and albumin D site-binding protein control tissue plasminogen activator levels: neurotoxic impact[END_REF][START_REF] Cacquevel | Ageing and amyloid-beta peptide deposition contribute to an impaired brain tissue plasminogen activator activity by different mechanisms[END_REF][START_REF] Obiang | Enriched housing reverses age-associated impairment of cognitive functions and tPAdependent maturation of BDNF[END_REF], we decided to study its expression in the hippocampus of mice from 8 weeks to 22 months old.

By using Rpl13a and Hprt, the two recommended reference genes of the Genorm analysis, we observed that the expression of tPA was decreased during ageing (Figure 9A, p=0.0079; Kendall's correlation test). However, when using Hmbs and ActinB, commonly considered as "standard" reference genes, and also under the threshold of 0.5 on the Genorm M value, no changes in the expression of tPA was observed during ageing (Figure 9B, p=0.7399; Kendall's correlation test). These results highlight the importance of efficient normalization for reliable results in gene expression studies in ageing.

Discussion

RT-qPCR is a useful and sensitive technique that gives a large amount of data within a short amount of time and is able to detect small differences of expression.

However, because of technical variations (storage condition, RNA isolation and purification, amount of starting material, enzymatic efficiency …), normalization is an obligatory step [START_REF] Bustin | Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): trends and problems[END_REF][START_REF] Bustin | Pitfalls of quantitative real-time reverse-transcription polymerase chain reaction[END_REF][START_REF] Huggett | Real-time RT-PCR normalisation; strategies and considerations[END_REF]. Since RT-qPCR has become a routine technique for molecular biology laboratories, numbers of reference genes have been used, with glyceraldehyde-3-phosphate dehydrogenase (Gapdh) or beta-actin (ActinB) as goldstandards. However, in many cases, these genes have been demonstrated to be unstable, like ActinB during cerebral ischaemia [START_REF] Harrison | The use of quantitative RT-PCR to measure mRNA expression in a rat model of focal ischemia-caspase-3 as a case study[END_REF], or in many other conditions (reviewed in [START_REF] Suzuki | Control selection for RNA quantitation[END_REF]).

In the present study, we show that (i) during ageing, in different cerebral structures, there are variations of expression for nine commonly used reference genes (Actinb, Polr2a, Sdha, Ppox, Hprt, Rpl13a, Hmbs, Ppib and GAPDH), (ii) at least two genes are necessary for a good normalization process in brain structures during ageing, (iii) each cerebral structure studied here (cortex, hippocampus, striatum and cerebellum) requires specific paired-reference genes (Table 3).

It has already been suggested by Sieber and colleagues in 2010 that in physiological brain during ageing, RT-qPCR needs to be normalized to more than one reference gene [START_REF] Sieber | Inter-age variability of bona fide unvaried transcripts Normalization of quantitative PCR data in ischemic stroke[END_REF]. They performed the study on whole brains during physiological ageing (from 8 weeks to 24 months) but also in animals subjected to a model of stroke. It appeared that ischaemia leads to the modification of the reference genes to use for normalization of RT-qPCR data in ageing (Hprt and Hmbs Vs ActinB and Hmbs).

Despite the importance of reference genes, especially for ageing, many studies use only one reference gene, and do not check its stability in experimental conditions as recommended by the MIQE [START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF]. Here, we show the importance of using at least two reference genes to normalize RT-qPCR data, but also the importance of checking their stability of expression. It is confirmed by studying the expression of tPA during ageing in the hippocampus. By using the GeNorm recommended reference genes, we observed a decrease in tPA expression during ageing, as previously reported [START_REF] Roussel | Age and albumin D site-binding protein control tissue plasminogen activator levels: neurotoxic impact[END_REF][START_REF] Cacquevel | Ageing and amyloid-beta peptide deposition contribute to an impaired brain tissue plasminogen activator activity by different mechanisms[END_REF][START_REF] Obiang | Enriched housing reverses age-associated impairment of cognitive functions and tPAdependent maturation of BDNF[END_REF]. However, when we used two other genes commonly considered as "standard" reference genes, we did not see any difference in tPA expression with age.

It is interesting to note that for the hippocampus and the cerebellum, there is a small difference between the two proposed reference genes. Indeed, for the hippocampus for example, the GeNorm analysis recommends Rpl13a and Hprt, while the Normfinder analysis recommends Rpl13a and Sdha. This is probably due to the difference between the two algorithms: GeNorm gives a stability value on the principle that the expression ratio of two ideal internal control genes is identical in all samples [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF]; whereas Normfinder combines the intra-and the intergroup variation into a stability value [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF]. However, when there is a difference between the two methods, the reference genes are always close to the most stable genes, reflecting their ability to be a good reference gene.

Our data also reveal that each structure needs its own reference paired-genes.

This is very important as many neurological diseases are structure-specific. In stroke, depending of the model used, the majority of the infarct lesion is in the cortex and in the striatum [START_REF] Orset | Mouse model of in situ thromboembolic stroke and reperfusion[END_REF]; various models of Parkinson's disease affect the substancia nigra [START_REF] Duty | Animal models of Parkinson's disease: a source of novel treatments and clues to the cause of the disease[END_REF]; in models of Alzheimer's disease the first amyloid deposits are structurespecific depending of the model used [START_REF] Seeman | Alzheimer's disease: beta-amyloid plaque formation in human brain[END_REF]. Moreover, a co-factor of all of these diseases is ageing, and it has been previously shown that ageing plays important roles in stroke [START_REF] Roussel | Age and albumin D site-binding protein control tissue plasminogen activator levels: neurotoxic impact[END_REF] for example. This is particularly important for pathological studies in ageing, as the pathology by itself will change the reference genes needed for a good normalization process [START_REF] Sieber | Inter-age variability of bona fide unvaried transcripts Normalization of quantitative PCR data in ischemic stroke[END_REF], however the same attention should be used for cerebral structure selection, instead of working on whole brain, as we demonstrated here.

Our study highlights the importance of the brain structure studied, by showing that each structure needs its own specific reference genes. It also shows that a stability expression control is primary, as demonstrated by the expression of tPA within the hippocampus. It comes in addition with the MIQE guidelines [START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF] and other publications which appeal for "good practices" in RT-qPCR studies [START_REF] Bustin | The need for transparency and good practices in the qPCR literature[END_REF][START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF][START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF].

Age is a common factor to many neurological disorders and RT-qPCR made from brain tissue in ageing should be well designed to avoid miss-interpretations of the data. Sieber and colleagues were the first ones to underline this specific point in neurological disorders [START_REF] Sieber | Inter-age variability of bona fide unvaried transcripts Normalization of quantitative PCR data in ischemic stroke[END_REF]. However they performed their study on the whole brain.

Many neurological studies are more focused on a specific cerebral area (hippocampus, cortex, striatum…) depending on the neurological disease studied.

Here, we demonstrate that each cerebral structure needs its own reference genes validation, and we propose the best couple of genes for the cortex, the hippocampus, the striatum and the cerebellum for ageing studies in the C57BL6/J mice. However, for pathological studies in ageing, the stability expression of reference genes needs to be checked.

Material and methods

Aged animals

Our study was performed on male C57BL/6 mice of distinct age and approved by local (French Legislation act no. 87-848, Ministère de l'Agriculture et de la Forêt)

and European (directives from the European Communities Council, November 24, 1986, 86/609/EEC) authorities. C57BL/6 mice have a lifespan of ~26 months [START_REF] Forster | Genotype and age influence the effect of caloric intake on mortality in mice[END_REF].

Four groups were designed: 8 weeks-old, 7 months-old, 17 months-old and 22 months-old, with 5 mice per group. Animals were housed in a temperature-controlled room on a 12-h-light/12-h-dark cycle, with food and water give ad libitum.

Brain sample preparation

Mice were anesthetized with 5% isofluorane in a 70%/30% mixture of NO2/O2 and perfused with cold heparinised NaCl 0.9% to wash out blood cells. Brains were excised and precisely dissected under a binocular microscope to isolate cerebral cortex, hippocampus, cerebellum and striatum. Tissue sections were maintained at -80°C until mRNA extraction.

Total RNA extraction and quantification

Total RNA was isolated from brain structures with TRI-reagent (Sigma; Lyon, France) according to the manufacturer. Briefly, tissues were dissociated in TRIreagent and RNA was isolated by the addition of chloroform. Total RNA was washed by ethanol and treated with TURBO DNase (Ambion; Saint Aubin, France) to avoid DNA contamination. Total RNA was quantified by spectrophotometry (Nanodrop Technologies; Wilmington, USA).

Reverse transcription

We performed the reverse transcription from 1µg of total RNA and the iScript kit (iScript Select cDNA Synthesis Kit, Biorad; Marnes-la-Coquette, France) in a total volume of 20µL with the following cycle conditions: 42°C (90min); 85°C (5min). The cDNA products were then stored at -20°C until their use.

Quantitative PCR

RT-qPCR were performed from 1µL of 1:20 diluted cDNA in 15µL total of a 1X solution of IQ SYBR Green Supermix (Biorad; Marnes-la-Coquette, France) containing 200nM of each primer.

Based on mRNA coding sequences (www.ensembl.org), mouse specific primers were designed by using the Primer3Plus software (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi) [START_REF] Rozen | Primer3 on the WWW for general users and for biologist programmers[END_REF] (sequences are in Table 1). Assays were run in triplicate on the CFX96 Real-Time system c1000 thermal cycler (Biorad; Marnes-la-Coquette, France), with the following cycle conditions: 95°C (3min); [95°C (2sec), 60°C (20sec)]x39; 70°C (30sec).

Data analysis

Each plate contains the whole samples (triplicate of the 5 mice per group) for one gene and one cerebral structure (8 genes, 4 structures: 32 plates of RT-qPCR)

and Cq values were obtained from the Bio-rad CFX manager software (Biorad;

Marnes-la-Coquette, France). Relative gene expression values (Starting quantity)

were calculated with the Bio-rad CFX manager software and corrected with each efficiency (E) inferred from 5-fold standard dilution series, as described by Pfaffl [START_REF] Pfaffl | A new mathematical model for relative quantification in real-time RT-PCR[END_REF].

E= 10-(1/slope) Kruskal Wallis test was performed on the starting quantities to check the stability of each gene in the studied structures.

To investigate the best normalization process to reference genes, we used the GeNorm algorithm associated to QBase+™ software, and the Normfinder algorithm.

GeNorm allows the calculation of stability (M value) and the ranking of selected reference genes, on the basis of non-normalized expression levels [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF][START_REF] Hellemans | qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data[END_REF]. It also suggests the number of internal controls to use in the selected experimental conditions by giving the V value. It is stated that the benefit of using an extra (n+1) reference gene is limited as soon as the V n/n+1 value drops below the 0.15 threshold. Moreover, the stability calculation is based on a geometric mean rather than an arithmetic mean by following these equations [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF]:

For every combination of 2 internal control genes j and k, an array A jk of m elements is calculated which consist of a log 2 -transformated expression ratios a ij /a ik (Equation 1). We define the pairwise variation V jk for the control genes j and k as the standard deviation of the A jk elements (Equation 2). The gene-stability measure M j for control gene j is the mean of all pairwise variations V jk (Equation 3).

Equation1:

Legends:

Table1: list of the nine reference genes used and their forward and reverse primers sequences. Ensembl (CCDS) numbers, length of the RT-qPCR transcripts, and efficiency of each primers are also indicated gene.

Figure 1 :

 1 Figure 1: raw data of RT-qPCRs performed on cortex of mice aged from 8 weeks to 22 months old. Starting quantity was directly obtained from the Cq value using the Bio-rad CFX Manager software for Ppib (A), Hmbs (B), Rpl13a (C), ActinB (D), Hprt (E), Ppox (F), Sdha (G), Polr2a (H), and GAPDH (I). n=5 independent animals in each group.

Figure 2 :

 2 Figure 2: raw data of RT-qPCRs performed on hippocampi of mice aged from 8 weeks to 22 months old. Starting quantity was directly obtained from the Cq value using the Bio-rad CFX Manager software for Ppib (A), Hmbs (B), Rpl13a (C), ActinB (D), Hprt (E), Ppox (F), Sdha (G), Polr2a (H), and GAPDH (I). n=5 independent animals in each group.

Figure 3 :

 3 Figure 3: raw data of RT-qPCRs performed on striatum of mice aged from 8 weeks to 22 months old. Starting quantity was directly obtained from the Cq value using the Bio-rad CFX Manager software for Ppib (A), Hmbs (B), Rpl13a (C), ActinB (D), Hprt (E), Ppox (F), Sdha (G), Polr2a (H), and GAPDH (I). n=5 independent animals in each group.

Figure 4 :

 4 Figure 4: raw data of RT-qPCRs performed on cerebellum of mice aged from 8 weeks to 22 months old. Starting quantity was directly obtained from the Cq value using the Bio-rad CFX Manager software for Ppib (A), Hmbs (B), Rpl13a (C), ActinB

Figure 9 :

 9 Figure9: expression of tPA measured by RT-qPCR according to the GeNormselected reference genes (A), or to random reference genes (B). The correlation between the expression and the age is calculated with the Kendall correlation coefficient.

Figure 1

 1 Figure 1: Cortex

Table 2 :

 2 Kruskal wallis test of the nine genes in the four studied cerebral structures during ageing. The test was performed on non-normalized starting quantities with

	Graph Prism software. *: p<0.05; **: p<0.01; ***: p<0.001.

Table 3 :

 3 choice of reference genes for RT-qPCR normalization after GeNorm analysis.

Acknowledgements:

We would like to thank the Basse-Normandie region for funding Bruckert G and Dr Roussel BD, and its credits for the laboratory.

Then, according to the M value, the best reference genes are selected with the lowest pairwise variability.

Normfinder estimates first variations within the group, and then variations between different groups. Having estimated both the intra-and intergroup variation, the algorithm combines the two into a stability value, which adds the two sources of variation and thus represents a practical measure of the systematic error that will be introduced when using the investigated genes [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF].

(D), Hprt (E), Ppox (F), Sdha (G), Polr2a (H), and GAPDH (I). n=5 independent animals in each group.