Philippe Angot
email: angot@cmi.univ-mrs.fr-calta@enscpb.fr-fabrie@math.u-bordeaux1.fr

Jean-Paul Caltagirone

Pierre Fabrie

-Jean-Paul Caltagirone

Vector Penalty-Projection Methods for the Solution of Unsteady Incompressible Flows

Keywords: Vector penalty-projection methods, projection methods, artificial compressibility, Navier-Stokes equations, incompressible and variable density flows, cell-centered MAC scheme

come

Vector penalty-projection methods for the solution of unsteady incompressible flows

Introduction: overview of methods dealing with the divergence constraint

Let us consider the unsteady Navier-Stokes problem in a bounded domain Ω ⊂ R d (d ≤ 3) with Dirichlet boundary conditions v |Γ = v D on Γ = ∂Ω and f , q (a volumic mass source), Re (Reynolds number), an initial condition v(t = 0) = v 0 given:

∂ t v + (v• ∇)v - 1 Re ∆v + ∇p = f with ∇• v = q in Ω × (0, T ).
[1]

Fully-coupled Navier-Stokes solvers

With usual notations for the semi-discretization in time, the linearly implicit Euler method with v0 = v 0 and a given time step δt > 0 reads: for all n ∈ N such that (n + 1)δt ≤ T , find vn+1 and pn+1 satisfying vn+1

|Γ = v n+1 D such that: vn+1 -vn δt + (v n • ∇)v n+1 - 1 Re ∆v n+1 + ∇p n+1 = f n+1 in Ω [2] ∇• vn+1 = q n+1 in Ω [3]
Here vn , pn are desired to be first-order approximations of the continuous velocity and pressure solutions v(t n ), p(t n ) at time t n = nδt. Of course, higher-order schemes are generally used for practical computations. For the sake of simplicity in the numerical procedure, a semi-implicit scheme where the non-linear term is linearized on the first term is often chosen since it does not suffer from stability conditions. The error analysis of the fully implicit scheme is simpler but it also makes it necessary to use a quasi-Newton algorithm to solve the corresponding non-linear system at each time step in the practical calculations. It is well-known that such a fully-coupled scheme, with any stable space discretization, yields an ill-conditionned and indefinite algebraic system at each time step, especially for small spatial mesh steps h.

Augmented Lagrangian and artificial compressibility methods

Falling into the class of semi-coupled solvers where the velocity components are still gathered, the augmented Lagrangian method which is classical in optimization [FOR 83], with the parameters 0 < λ < 2r + 2/Re and the initial condition p 0 given, reads:

v n+1 -v n δt + (v n • ∇)v n+1 - 1 Re ∆v n+1 -r∇ ∇• v n+1 -q n+1 + ∇p n = f n+1 p n+1 = p n -λ ∇• v n+1 -q n+1 .
This algorithm must in fact be iterated at each time step (not given here for sake of shortness) to produce good approximations of the solution to the reference method, i.e. the coupled implicit Euler scheme [2]-[3] for unsteady flows, see [KHA 00] for details. However, the convergence to a "divergence free" velocity, varying like O(1/r), is all the faster as the augmentation parameter r 1 is larger but the resulting system becomes here also ill-conditioned which greatly increases the solution cost, especially for 3-D problems. On the contrary, the standard Uzawa algorithm with r = 0 converges very slowly or may not converge at all for strongly non-linear problems. With the choice λ = r = δt/ε, the method also yields the implicit artificial (or pseudo) compressibility method studied by Chorin (1967) and Temam (1968) [TEM 86]. Indeed, the scheme then corresponds to the implicit first-order discretization of the following singularly perturbed continuous problem where v ε (0) = v 0 and also p ε (0) are given:

∂ t v ε + (v ε • ∇)v ε - 1 Re ∆v ε + ∇p ε = f in Ω × (0, T ), [4] ε ∂ t p ε + ∇• v ε = q in Ω × (0, T ), [5]
with v ε = v D on Γ. This artificial compressibility problem (for q = 0) can be proved to converge in some sense to the Navier-Stokes system when ε tends to zero, see [TEM 86] and [ANG 08] for new convergence results like Theorem 2.1.

Scalar projection and penalty-projection methods

With a completely different approach, the class of fractional-step or splitting methods, namely the projection methods originally introduced by Chorin (1968) and Temam (1969) [TEM 86], became very popular because they are much cheaper. They produced an important literature with many variants of the type of pressure-correction or velocity-correction schemes and we refer to the complete review of Guermond et al. [GUE 06] and the references therein. Their principle lies in the Helmholtz-Hodge orthogonal decomposition of the space L 2 (Ω) d = H⊕H ⊥ where, see e.g. [TEM 86]:

H = {u ∈ L 2 (Ω) d , ∇• u = 0, u• n = 0 on Γ} and H ⊥ = {∇φ, φ ∈ H 1 (Ω)}.
This allows for a predicted velocity field to be corrected with a "pressure" gradient by the Leray orthogonal projection P H onto the space of solenoidal fields H. As an example, let us consider the so-called penalty-projection method recently proposed and numerically investigated by Jobelin et al. [JOB 06, JOB 07], also theoretically analysed in [ANG 06] -the first one being due to Shen (1992). For r 1 ≥ 0 and p 0 given, the predicted velocity ṽn+1 satisfies in Ω:

ṽn+1 -v n δt +(v n • ∇)ṽ n+1 - 1 Re ∆ṽ n+1 -r 1 ∇ ∇• ṽn+1 -q n+1 +∇p n = f n+1 , [6] with ṽn+1 = v n+1 D on Γ. This system condition number varies as O(r 1 ) [FEV 08].
Then the pressure correction step reads for r 2 ≥ 0:

-∆φ = 1 δt q n+1 -∇• ṽn+1 in Ω with ∇φ• n = 0 on Γ [7] v n+1 -ṽn+1 + δt∇φ = 0 and p n+1 = p n + φ -r 2 ∇• ṽn+1 -q n+1 in Ω. [8]
Under this general form, the standard incremental projection method from Goda (1979) is recovered with r 1 = r 2 = 0 which makes it very simple and attractive. However, the usual projection methods suffer from some drawbacks:

1) First, the original Dirichlet boundary condition on the velocity degenerates into an inconsistent Neumann boundary condition for the pressure increment φ due to the scalar pressure correction step and so the end-of-step velocity only satisfies: 2) When an outflow boundary condition of the "Neumann" type like in Section 3 is imposed on a part of the boundary, it degenerates into a homogeneous Dirichlet condition for φ and the convergence of the standard projection methods is quite poor. Indeed, the l 2 (L 2 ) norm of the velocity and pressure errors converges as O(δt) and 3) In the case of variable density flows, as multiphase flows, the scalar correction step is naturally modified depending now on the density n+1 [GUE 00], as follows:

v n+1 • n = ṽn+1 • n = v n+1 D • n.
O(δt 1 2 ),
-∇• δt n+1 ∇φ = q n+1 -∇• ṽn+1 in Ω with ∇φ• n = 0 on Γ n+1 v n+1 -ṽn+1 + δt∇φ = 0, p n+1 = p n + φ -r 2 ∇• ṽn+1 -q n+1 in Ω.
Thus, it introduces an "artificial" dependence on n+1 of the correction of ṽn+1 : in the case of large density ratios, it generally prevents us from solving φ quasi-exactly which should be to satisfy the divergence constraint, whereas a standard precision in agreement with the discretization errors would be sufficient to get only the pressure! 4) The last drawback, less traditional, addresses the fictitious domain methods using a penalization to impose an immersed boundary condition; see [ANG 99, KHA 00] and the references therein. Here, the L 2 -volumic penalty inside a moving subdomain ω(t) ⊂ Ω with b(t) = χ ω(t) /η, η > 0, modifies the projection method as below:

ṽn+1 -v n δt + (v n • ∇)ṽ n+1 - 1 Re ∆ṽ n+1 + ∇p n + b n+1 ṽn+1 -v n+1 D = f n+1 -∇• δt 1 + δt b n+1 ∇φ = q n+1 -∇• ṽn+1 in Ω with ∇φ• n = 0 on Γ v n+1 -ṽn+1 + δt∇φ + δt b n+1 v n+1 -ṽn+1 = 0, p n+1 = p n + φ in Ω.
Thus, it yields a H 1 -penalty for φ with a diffusion coefficient varying as O(η) when η → 0, i.e. a singularly perturbed problem more difficult to treat numerically. This also gives an additional insight to our first point about the pressure boundary layers.

The vector penalty-projection (VPP r,ε ) methods for incompressible flows

The vector penalty-projection methods are derived to overcome most of the latter drawbacks with a compromise between the best properties of both classes: the augmented Lagrangian (without iterations) and splitting methods under a vector form, i.e. seeking a vector correction v to inherently respect the velocity boundary conditions.

The two-step artificial compressibility (VPP ε ) method

A first two-step VPP scheme with a penalty parameter 0 < ε ≤ δt is devised, for

v 0 ∈ H 1 (Ω) d , p 0 ∈ L 2 0 (Ω) given, as follows: ∀n ∈ N such that (n + 1)δt ≤ T , ṽn+1 -v n δt + (v n • ∇)ṽ n+1 - 1 Re ∆ṽ n+1 + ∇p n = f n+1 ε δt vn+1 δt + (v n • ∇)v n+1 - 1 Re ∆v n+1 -∇ ∇• vn+1 = ∇ ∇• ṽn+1 -q n+1 v n+1 = ṽn+1 + vn+1 , and p n+1 = p n - δt ε ∇• v n+1 -q n+1 with ṽn+1 |Γ = v n+1 D
and vn+1 |Γ = 0 to satisfy the boundary condition on v. Summing the two steps yields the implicit artificial compressibility scheme written in section 1.2.

The two-parameter vector-penalty projection (VPP r,ε ) methods

The two-parameter vector penalty-projection methods with an augmentation parameter r ≥ 0 and a penalty parameter 0 < ε ≤ 1 are then devised. They read for

v 0 ∈ H 1 (Ω) d and p 0 ∈ L 2 0 (Ω) given: for all n ∈ N such that (n + 1)δt ≤ T , ṽn+1 -v n δt + (v n • ∇)ṽ n+1 - 1 Re ∆ṽ n+1 -r∇ ∇• ṽn+1 -q n+1 + ∇p n = f n+1 pn+1 = p n -r ∇• ṽn+1 -q n+1 ε vn+1 δt + (v n • ∇)v n+1 - 1 Re ∆v n+1 -∇ ∇• vn+1 = ∇ ∇• ṽn+1 -q n+1 v n+1 = ṽn+1 + vn+1 , and p n+1 = pn+1 - 1 ε ∇• v n+1 -q n+1 with ṽn+1 |Γ = v n+1 D and vn+1 |Γ = 0.
Notice that (VPP 0,ε ) is different from (VPP ε ) unless δt = 1. Indeed, summing now the two steps gives the discrete problem below which provides an additional dissipation of the velocity divergence for δt < 1:

v n+1 -v n δt + (v n • ∇)v n+1 - 1 Re ∆v n+1 + ∇p n+1 = f n+1 [9] ε p n+1 -p n δt + 1 δt ∇• v n+1 -q n+1 + rε δt ∇• ṽn+1 -q n+1 = 0. [10]
This results in a well-posed generalization of the method early proposed in [CAL 99] considering a singular vector correction step, with ε = 0 in fact, which has a unique solution only with an additional constraint for v, such that e.g.:

∇ × v = 0.
Defining e n = v n -v(t n ) and π n = p n -p(t n ) as the whole velocity and pressure errors with the Navier-Stokes true solution, respectively, we prove in [ANG 08] the following result of first-order convergence in time for the time semi-discrete setting. THEOREM 2.1 (ERROR ESTIMATES FOR (VPP ε -VPP r,ε ) WITH v D = 0 AND q = 0.) Assume (v, p) the solution of [1] smooth enough in time and space, well-prepared initial conditions and for (VPP r,ε ), small enough parameters such that 4r(Re + ε) ≤ 1 and 4c(Ω) √ Re rε ≤ √ δt, c(Ω) being the Poincaré constant, then there exists C = C(Ω, T, Re, f , v 0 , e 0 , π 0 ) > 0 such that we have for all n ∈ N with (n + 1)δt ≤ T ,

(i-VPP r,ε ) ||e n+1 || 2 0 + εδt ||π n+1 || 2 0 + n k=0 δt Re ||∇e k+1 || 2 0 ≤ C δt 2 + ε 2 δt 3 2 , (ii-VPP r,ε ) n k=0 δt ||π k+1 || 2 0 ≤ C δt 2 + ε 2 δt , ||∇e n+1 || 2 0 ≤ C Re 2 δt + ε 2 , (iii-VPP r,ε ) n k=0 δt ||∇• v k+1 || 2 0 = n k=0 δt ||∇• e k+1 || 2 0 ≤ C (δt + ε) εδt 2 . (iii-VPP ε ) n k=0 δt ||∇• v k+1 || 2 0 = n k=0 δt ||∇• e k+1 || 2 0 ≤ C 1 + ε 2 δt 5 2 + ε δt εδt.

Generalizations of (VPP r,ε ) methods for outflow BC or variable density flows

It is straightforward to write the above (VPP) methods with a higher-order scheme, e.g. the usual second-order Gear scheme (BDF2) as in [KHA 00, JOB 06].

1) (VPP r,ε ) methods for incompressible and variable density flows

By omitting here the convection term and with q = 0 for sake of shortness, the vector prediction and correction steps in section 2.2 are naturally modified as follows:

n+1 -n δt + ∇• ( n+1 v n ) = 0 n+1 ṽn+1 -v n δt -∇• µ n+1 ∇ṽ n+1 + (∇ṽ n+1 ) T -r∇ ∇• ṽn+1 + ∇p n = f n+1 ε n+1 vn+1 δt -∇• µ n+1 ∇v n+1 + (∇v n+1 ) T -∇ ∇• vn+1 = ∇ ∇• ṽn+1 .
It is thus clear that the velocity correction v becomes quasi-independent on the density and viscosity µ as ε → 0, see section 4, which is not the case for scalar projection methods. These terms can be then dropped in practical algorithms for small enough ε.

ṽn+1 = v n+1 D on Γ D , -p n n + µ n+1 ∇ṽ n+1 + (∇ṽ n+1 ) T • n = g n+1 on Γ N vn+1 = 0 on Γ D , -(p n+1 -p n ) n + µ n+1 ∇v n+1 + (∇v n+1 ) T • n = 0 on Γ N .

Numerical experiments with the finite volume MAC scheme

The (VPP r,ε ) methods are implemented with a Navier-Stokes finite volumes solver on the staggered MAC mesh of size h issued from previous works; see [KHA 00]. The first test case is the unsteady Green-Taylor vortex such that the mean steady velocity field is of order 1 at Re = 100 like in [CAL 99]. The scheme is O(δt) in time for velocity and pressure with r ≥ 10 -3 , i.e. for rRe = O(1), see Figure 1, whereas it is also O(h 2 ) in space as in [CAL 99]. Moreover, we find that the L 2 -norm of the velocity divergence varies like O(εδt), as expected by the theory; see Theorem 2.1. However, the pressure does not seem to converge for r < 10 -3 , unless performing inner iterations as for the augmented Lagrangian which is not our goal. This is somewhat deceptive since we should expect convergence for smaller values of r in order to degrade the conditioning of the prediction step at the very least. The second benchmark problem is the Rayleigh-Bénard natural convection inside a square differentially heated vertical cavity at Ra = 10 5 , the vertical walls being isothermal and the horizontal walls insulating. Here, we study the convergence properties of the penalty-correction step for this sharp test case at t = 2δt with δt = 1. Again, we get the convergence of the velocity divergence as O(εδt), whatever the viscosity term included in the penaltycorrection step and also for µ = 0. We can reach the machine precision of 10 -15 for double precision floating point computations. Besides, the solution of the penaltycorrection step proves to be all the cheaper as εδt tends to zero: typically one or two iterations of a preconditioned gradient algorithm, e.g. the ILU-BiCGStab solver, are sufficient for εδt ≤ 10 -6 and the condition number of the operator (preconditioned or not) also becomes more and more quasi-independent on the spatial mesh size h. Hence, the central issue for the proposed scheme is that the linear system associated with the projection step can be solved very easily because of the specific form of the right-hand side. Indeed, this latter lies more and more in the image of the left-hand side operator as ε is taken smaller and smaller; see the projection step in section 2.2.

Perspectives

The generalization of the (VPP r,ε ) methods for low Mach number dilatable flows in the same spirit of what was done in [JOB 07] is the subject of a future work. Moreover, their extension to more general vector problems for fictitious domain methods with immersed moving and/or deformable interfaces, e.g. multiphase flows, particulate flows or fluid-structure interaction problems, is an undergoing project.

  This creates artificial pressure boundary layers which spoil the pressure numerical solution. They can be reduced with the rotational projection method from Timmermans et al.(1996) corresponding to the previous algorithm with r 1 = 0, r 2 = 1/Re [GUE 04] or with the standard penalty-projection (r 1 = r 2 = r) or the rotational version (r 1 = r, r 2 = r + 1/Re) for which they are nearly suppressed with moderate values of r between 1 and 10 [JOB 06, FEV 08].

  2) (VPP r,ε ) methods for open boundary conditions on a part Γ N of Γ For a given stress vector (σ(v, p)• n) |ΓN ≡ -p n + µ ∇v + (∇v) T • n = g, we now get to satisfy both the Neumann and Dirichlet velocity boundary conditions:

Figure 1 .

 1 Figure 1. LEFT: (VPP r,ε ) velocity and pressure convergence in time for the Green-Taylor vortex at Re = 100, t = 10h = 1/128, ε = 10 -4 with ||∇• v n || L 2 = O(ε). RIGHT: Convergence of the velocity correction step -divergence L 2 -norm δ as O(ε) -for the natural convection at Ra = 10 5 and t = 2δt with δt = 1, h = 1/256µ = 0 or 1.85 10 -5 (idem) and µ = 1.85 10 -1 .

  respectively for a second-order time scheme, see [GUE 06]. This is clearly improved with the rotational projection [GUE 06] or with the penalty-projection methods if r takes moderate values, typically 1 ≤ r ≤ 10, see [JOB 06, FEV 08].


M A Aguirre

J G Grande

A Calvo

L A Pugnaloni

J.-C Géminard

Pressure independence of granular flow through an aperture

Keywords: numbers: 45.70.-n, 45.70.Mg, 47.80.Jk

We experimentally demonstrate that the flow rate of granular material through an aperture is controlled by the exit velocity imposed to the particles and not by the pressure at the base, contrary to what is often assumed in previous works. This result is achieved by studying the discharge process of a dense packing of monosized disks through an orifice. The flow is driven by a conveyor belt. This two-dimensional horizontal setup allows to uncouple pressure and velocity and, therefore, to independently control the velocity at which the disks escape the horizontal silo and the pressure in the vicinity of the aperture. The flow rate is found to be directly proportional to the belt velocity, independent of the amount of disks in the container and, thus, independent of the pressure in the outlet region. In addition, this specific experimental configuration makes it possible to get information on the system dynamics from a single image of the disks that rest on the conveyor belt after the discharge.

The flow of granular media through an orifice presents some interesting features that have been intensely studied in the last 50 years [1][2][3][4][5][START_REF] Duran | Sands, Powders and Grains[END_REF][START_REF] Ristow | Pattern Formation in Granular Materials[END_REF][START_REF] Nedderman | [END_REF][9][10][11]. What mainly differentiates the discharge of a container filled with granular matter from one filled with a viscous liquid is that the mass flowrate does not depend on the height h of material above the outlet. The explanation most frequently used for this independence is based on the Janssen's effect: the distribution of the weight of the material onto the silo walls, due to the friction forces, leads to a saturation of the pressure at the bottom, which results in a constant flow-rate [12]. In such a reasoning, the pressure is thus implicitly assumed to govern the flow rate.

Here we show, by using an experimental setup in which the exit velocity of the grains is decoupled from the bottom pressure, that the above argument is improper. Indeed, different flow-rates can be achieved for the same bottom pressure; moreover, the flow-rate remains constant even if the pressure decreases during the discharge.

In general, the discharge of a silo through an orifice can present three regimes: a continuous flow, an intermittent flow, or a complete blockage of the flow due to arching [13][14][15]. In the continuous flow regime, the mass flow rate W is described by the so called Beverloo's law [1,[START_REF] Brown | Principles of Powder Mechanics[END_REF]:

W = Cρ 3D √ g(A -k D) 5/2
where A is the size of the opening, g the acceleration due to gravity, ρ 3D the bulk density and D the diameter of the granules whereas k and C are two empirical dimensionless constants. In a two-dimensional (2D) setup-or similarly in slit shaped apertures-one expects W = Cρ 2D √ g(A -kD) 3/2 [START_REF] Brown | Principles of Powder Mechanics[END_REF]. In the jamming regime, the jamming probability is controlled by the ratio A/D of the aperture size to the grain diameter [14,15,[START_REF] To | [END_REF][18][19][20].

The discharge of particulate systems is not only found under the influence of gravity. In many industrial appli- cations, the grains are horizontally transported at constant velocity (as on a conveyor belt [21] or floating on the surface of a flowing liquid [START_REF] Guariguata | Proceedings of the 6 th International Conference on Micromechanics of Granular Media[END_REF]). Here, we analyze the discharge through an orifice of a 2D pile of disks driven at a constant velocity by a horizontal conveyor belt.

The experimental setup (Fig. 1) consists of a conveyor belt (width 40 cm, length 1 m) above which a confining Plexiglas frame (width 26 cm, length 54 cm) is maintained at a fixed position in the frame of the laboratory. A motor drives the belt at a constant velocity V which ranges from 0.4 cm/s to 4 cm/s. Downstream, the Plexiglas frame exhibits, at the center, a sharp aperture (width A ≤ 10 cm), much wider than the thickness of the walls (0.6 cm). The granular material consists of N = 450 Plexiglas disks [thickness e = (3.1 ± 0.1) mm, and diameter D = (1.04 ± 0.01) cm]. Thus, in our experimental conditions, the ratio A/D ≤ 10. The dimensionless slope Q * ≡ QD/V obtained fitting data reported in the main panel (open squares) and/or from the probability distributions P (Q * ) (crosses) reported in Fig. 4 (errors are ±5%). The flow rate is expected to vanish for the critical flow rate, Ac.

Before the flow is started, the initial state of the system is obtained by depositing, in a disordered manner, N disks on the conveyor belt inside the confining frame. The belt is then moved at a low constant velocity until all grains are packed against the downstream wall, whose aperture is still kept closed (A = 0). The packing fraction of this initial configuration is Φ ≈ 0.82. Then, the aperture is opened to the desired width A = 0 and the belt is moved at constant velocity V . The belt is stopped whenever the Plexiglas frame is emptied or the flow ceases due to the formation of a jamming arch.

A video camera (Pixelink, PL-A741) is used to image the system from above. We record the evolution of the discharge process in either of the following ways: Images of disks that remain inside the confining frame (upstream) are acquired every second during the discharge; One single image of the jet of disks that are at rest on the belt outside the confining box (downstream) is acquired after the discharge.

An intensity threshold is used to convert the pictures into binary images. The number N i of disks that remain inside the confining frame at time t can be calculated from the number of white pixels in the images acquired upstream [Area of a disk (733 ± 8) pixels]. The instantaneous disk flow-rate (averaged over 1 s, because of the acquisition frequency) is defined as Q ≡ -dNi dt . When the aperture is large (A/D ≥ 6), the flow-rate is rather continuous throughout the discharge, N i depends linearly on time and, for any given A, the flow rate is proportional to the belt velocity (Fig. 2). This indicates that, in our experimental conditions, the granular bed rearranges in a time which is smaller or compares with the time necessary for the grains to escape the outlet region. Indeed, if a significant relaxation time is associated with the rearrangement of the granular bed, above a critical velocity, a change in the drive is barely followed by the system and a qualitative change in Q(V ) is observed [21]. In the inset in Fig. 2, we report the dimensionless slope Q * ≡ QD/V as a function of the dimensionless aperture A/D and get, for A/D ≥ 6, the empirical law:

Q = C V D A/D -k (1)
where k = 2.1 ± 0.2 and C = 1.04 ± 0.02. The flow rate is thus expected to vanish at a value A c ≡ kD of the aperture. The value of the constants k and C will be discussed below (see comment on Ref. [21]). Note that the flow rate is not well defined and cannot be extracted directly from N i (t) for A/D < 6, the system being then likely to jam.

In order to estimate the flow rate for A/D < 6, it is particularly interesting to consider the second measurement method in which we take a single picture of the jet at rest on the belt after the experiment is stopped (Fig. 3a). After being released from the confining frame, the grains are simply advected with velocity V and the belt is long enough so that all disks remain on it after the discharge. Thus, the distribution of disks on the belt provides information on the outflow of grains ever since its onset. Therefore, such an image provides the full history of the flow rate: the spatial coordinate along the belt longitudinal axis can be translated into time. Moreover, since the flow rate Q is proportional to V (Fig. 2), one can limit the study to a single value of the velocity.

Let us first comment on the intensity profile in the direction perpendicular to the jet (Fig. 3b): the fraction of white pixels in a row (one pixel wide and about 1500 pixels long) corresponds to the average density Φ of the disks in this row. The profile is smooth and, even for an aperture as large as seven disk diameters, there is no well defined plateau. Such boundary effects, which are present in a large part of the profile, are responsible for the effective aperture width, (A -A c ), pointed out by the finite value of k in Eq. (1).

We now consider the intensity profile along the longitudinal axis of the jet. In order to define an instantaneous flow rate, we arbitrarily integrate over a 30-pixelwide window (a disk diameter) around the considered position. The dimensionless flow rate Q * is obtained by counting the number of white pixels in a column of width D along the jet and dividing it by the surface area of one disk. Other window widths were tested: fluctuations are smaller for wider windows but the average value of Q * remains the same to within 0.2%. The example reported in Fig. 3c shows that, even for a large aperture, Q * is subjected to significant fluctuations (standard deviation around 15% and difference between extreme values of about 60% of the average). The probability distribution P (Q * ) is shown in Fig. 4. On the one hand, P is almost Gaussian when the aperture is large (A/D > 6). The probability of small Q * is negligible and the flow is continuous. On the other hand, when the size of the aperture is reduced, the system might jam (Q * = 0) transitorily or permanently (in agreement with similar experiments carried out in a 2D vertical silo [15,[START_REF] Janda | [END_REF]). In accordance, the probability of small Q * increases.

From the probability distribution, it is possible to define a mean flow rate even when the system presents intermittent jams. The average flow rate during the effective flow events can be estimated by defining Q * ≡ Q * P (Q * )dQ * (Note that jamming events correspond to a peak at Q * = 0 that we disregard in our analysis in order to account for the flow events only). The estimations of Q * obtained in this way for large values of A/D are in full agreement with measurements obtained from the upstream technique (Fig. 2,inset). In addition, the data obtained for smaller apertures also obey to Eq. ( 1), which indicates that there is no drastic change in the empirical law in presence of intermittency. However, we mention that, in 3D configurations, a deviation from the 5/2 Beverloo's scaling for very small apertures has been observed [13]. Thus, a deviation from the linear relation (1) may be present for A/D < 4.5.

We emphasize that, in our experimental conditions, the effective pressure exerted by the granular bed on the downstream wall depends on the dynamical friction force between the disks and the belt. So, at a constant velocity V , the bottom pressure continuously decreases as the number of disks inside the system decreases during the discharge. Thus, during the whole discharge, a constant flow-rate Q is achieved while the pressure continuously varies, which indicates that the granular flow rate is not controlled by it. The horizontal configuration makes it possible to impose the exit velocity and, moreover, to tune independently the overall pressure by changing the weight of the disks or their initial number inside the frame. In order to measure the flow-rate in a different pressure range for the same output velocity, we added an individual extra weight on top of each disk. We thus increased the pressure by increasing the friction force acting on the disks without changing any other property of the system. We did not observe any change in the flow-rate in spite of the change in the pressure, which is again in agreement with a flow-rate only controlled by the escape velocity.

At this point, it is particularly interesting to discuss the physical ingredients leading to Beverloo's law. In our experimental conditions, we obtain a simple linear law relating the flow-rate to the aperture width A, whereas the corresponding relation in a 2D gravity driven configuration corresponds to Q ∝ A 3 2 . Both relations can be derived from a dimensional analysis, taking into account the estimation of the velocity v at which the grains escape the system. In the case of the 2D vertical silo, the free falling arch assumption [START_REF] Nedderman | Statics and kinematics of granular material[END_REF] leads to v ∝ √ gA. Hence, considering the mass flow rate through an aperture of size A, W ≡ ρ 2D Av leads to W ∝ ρ 2D √ gA 3 2 . Boundary effects at the aperture edges lead to two boundary layers of thickness ∼ D (the so-called empty annulus [START_REF] Brown | [END_REF]). Considering the effective aperture size A -kD instead of A, one gets Q ∝ (A/D -k) 3 2 . same arguments lead to the Beverloo's law in 3D. In the case of the conveyor belt v = V and the account of the empty annulus leads to Q ∝ V (A/D -k) in agreement with Eq. (1). It is then worth mentioning that, in the limit of infinitely large apertures, the correction due to the boundary layers can be neglected and Q → ρAV (ρ is the number of disks per unit area) so that C = ρD 2 and, hence, Q = ρV D(A/D -k). Taking into account Φ = ρ πD 2 4 (Φ is the area fraction), we obtain that C = 4 π Φ. From the experimental value C = 1.13 ± 0.02, we estimate Φ = 0.82±0.02, which is indeed compatible with the area fraction initially measured inside the confining frame and which does not noticeably vary during the discharge. We remark that the dimensional analysis could be based on the pressure since it also exhibits a time-scale. In the case of a 3D silo, denoting P the pressure at the outlet, one gets W ∝ √ ρ 3D P A 2 . Be the pressure limited by Janssen's effect or not, W scales with A 2 , which contradicts the experimental A 5/2 scaling.

In most previous works [12], the independence of the flow rate on the silo height h is explained by means of two postulates. The first one is: if the height is larger than twice the diameter of the silo, the pressure at the bottom is constant (Janssen's effect). The second one is: if during the discharge the pressure at the bottom is constant, the flow rate is, as a consequence, also constant. We have proven that this second postulate is not valid. Indeed, experimentally, for the same pressure in the outlet region (same amount of disks inside the system), the flow rate is proportional to the belt velocity. Moreover, the same flow-rate is obtained with different values of the pressure in the outlet region (different amount of disks inside the system or additional weights on the disks). In addition, theoretically, the physical ingredients used to establish Beverloo's law, which is correctly satisfied experimentally, do not involve the local pressure but a characteristic velocity. Thus, Beverloo's law already predicts that the flow-rate should be independent of h even when the Jannsen's effect is not at stake.

In conclusion, using an experimental setup which enables to control the particle velocity independently from the other parameters of the flow, we have shown that the granular flow rate through an orifice is not controlled by the local pressure and that invoking the Jannsen's effect is not pertinent to explain the constant flow-rate mea-sured during a silo discharge. The flow rate is controlled by the mechanism driving the grains out of the system and by the geometry of the outlet, features which govern the exit velocity, and not by the pressure upstream. For applications, our findings thus suggest that granular flow rates can be increased, in gravity driven systems, by locally increasing the particle velocities at the outlet, which might be achieved with devices having a low energetic cost. From the fundamental point of view, they show that more effort must certainly be done to verify some misleading assertions widely stated in material science of dissipative particles during the last five decades.

FIG. 1 .

 1 FIG. 1. Sketch of the experimental setup-Disks, placed on a conveyor belt which moves at a constant velocity V (white arrow), are forced to flow through an aperture in the confining frame. Observation of the grains-either remaining in the confining frame or forming the jet outside-makes it possible to study the dynamical properties of the discharge.

FIG. 2 .

 2 FIG.2. Dependence of the flow rate Q on the belt velocity V for different apertures A/D (Average of 4 runs). Inset: The dimensionless slope Q * ≡ QD/V obtained fitting data reported in the main panel (open squares) and/or from the probability distributions P (Q * ) (crosses) reported in Fig.4(errors are ±5%). The flow rate is expected to vanish for the critical flow rate, Ac.

FIG. 3 .

 3 FIG. 3. Image of the jet and the corresponding density profiles. (a) Image of the jet after the discharge of the container for A/D = 7.0. (b) Area fraction profile Φ in the perpendicular direction. Each data point is obtained as the fraction of white pixels in the corresponding row after thresholding [the vertical dashed line corresponds to the maximum attainable density π/(2 √ 3)]. (c) Instantaneous dimensionless flow rate Q * against the longitudinal coordinate x. The flow rate is obtained from the fraction of white pixels in the corresponding column 1D wide.

FIG. 4 .

 4 FIG.4. Probability distribution of the dimensionless flow rate Q * for different values of A/D. For each aperture A/D, the probability is calculated from all the values of the flow-rate measured during a set of 5 to 10 experiments. For the largest aperture reported (A/D = 6.0), the probability distribution is close to Gaussian, whereas a relatively large probability of small Q * values appears when A is decreased. The continuous lines are only a guide to the eye.

ACKNOWLEDGMENTS

We thank J.P. Hulin for valuables comments. This work has been supported by the program PICT 32888 (ANPCyT) and the International Cooperation Program CONICET-CNRS. J.G.G. acknowledges support from UBA, L.A.P. and M.A.A. from CONICET.


Ph.D Jose L Cantero

Mercedes Atienza

German Gomez-Herrero

Abel Cruz-Vadell

Eulogio Gil-Neciga

Rafael Rodriguez-Romero

Garcia-Solis David

Functional integrity of thalamocortical circuits differentiates normal aging from mild cognitive impairment

Keywords: mild cognitive impairment, cerebral aging, functional connectivity, thalamocortical circuits, Alzheimer's disease, EEG alpha rhythm mild cognitive impairment, cerebral aging, functional connectivity, thalamocortical circuits, Alzheimer's disease, EEG alpha rhythm

Resonance in thalamocortical networks is critically involved in sculpting oscillatory behavior in large ensembles of neocortical cells. Neocortical oscillations provide critical information about the integrity of thalamocortical circuits and functional connectivity of cortical networks, which seem to be significantly disrupted by the neuronal death and synapse loss characterizing Alzheimer's disease (AD). By applying a novel analysis methodology to overcome volume conduction effects between scalp electroencephalographic (EEG) measurements we were able to estimate the temporal activation of EEG-alpha sources in the thalamus and parieto-occipital regions of the cortex. We found that synaptic flow underlying the lower alpha band (7.5-10 Hz) was abnormally facilitated in patients with mild cognitive impairment (MCI) as compared to healthy elderly individuals, particularly from thalamus to cortex (~38% higher). In addition, the thalamic generator of lower alpha oscillations was also abnormally activated in MCI patients. Regarding the upper alpha subdivision (10.1-12.5 Hz), both controls and MCI patients showed a bidirectional decrease of thalamocortical synaptic transmission which was age-dependent only in the control group. Altogether, our results suggest that functional dynamics of thalamocortical networks differentiate individuals at high risk of developing AD from healthy elderly subjects, supporting the hypothesis that neurodegeneration mechanisms are active years before the patient is clinically diagnosed with dementia.

Introduction

MCI is characterized by subtle but clinically significant memory loss beyond what is expected from normal aging. This clinical condition seems to represent a transitional stage between normal aging and some prevalent neurodegenerative disorders such as Alzheimer's disease (AD) [START_REF] Petersen | Mild cognitive impairment. Clinical characterization and outcome[END_REF]. From a neuropathological perspective, the density of neuritic plaques [START_REF] Sabbagh | Pathologic and nicotinic receptor binding differences between mild cognitive impairment, Alzheimer disease, and normal aging[END_REF][START_REF] Markesbery | Neuropathologic substrate of mild cognitive impairment[END_REF][START_REF] Price | Tangles and plaques in nondemented aging and "preclinical" Alzheimer's disease[END_REF] Morris 1999] and neurofibrilary tangles [START_REF] Guillozet | Neurofibrillary tangles, amyloid, and memory in aging and mild cognitive impairment[END_REF] characterizing AD also become apparent in the neocortex of MCI patients. Unfortunately, these instigators of neuronal damage can only be empirically identified in postmortem studies, which make them worthless in clinical practice as early markers of AD.

Cognitive symptoms present in AD are mostly caused by massive cell death and synaptic loss mainly affecting association areas and the pyramidal cells that supply the long projections among distant neocortical regions [START_REF] Morrison | Life and death of neurons in the aging brain[END_REF]. This scenario results in a global disruption of local and long-range neural circuits that are revealed by changes in brain oscillations recorded with the electroencephalogram (EEG). EEG measurements emerge from summated excitatory and inhibitory postsynaptic potentials primarily generated by synaptic activity of large pyramidal neurons [Lopes da [START_REF] Da | Biophysical aspects of EEG and magnetoencephalogram generation[END_REF]]. Among EEG oscillations, the alpha rhythm has been reported to be typically affected in AD patients [START_REF] Huang | Discrimination of Alzheimer's disease and mild cognitive impairment by equivalent EEG sources: a crosssectional and longitudinal study[END_REF][START_REF] Osipova | Altered generation of spontaneous oscillations in Alzheimer's disease[END_REF][START_REF] Babiloni | Sources of cortical rhythms change as a function of cognitive impairment in pathological aging: a multicenter study[END_REF]. Alpha waves are intrinsically generated by synaptic networks of neocortical pyramidal neurons in layer 5 by means of the tonic activation of N-methyl-D-aspartate receptors with endogenous glutamate [START_REF] Silva | Intrinsic oscillations of neocortex generated by layer 5 pyramidal neurons[END_REF]. The neuronal damage evident in diverse neocortical regions of MCI patients [START_REF] Sabbagh | Pathologic and nicotinic receptor binding differences between mild cognitive impairment, Alzheimer disease, and normal aging[END_REF] [START_REF] Gómez-Herrero | Measuring directional coupling between EEG sources[END_REF] was able to reliably determine the temporal activation of the EEG alpha sources and their approximate intracerebral locations. Causal relationships between the neural sources of EEG-alpha rhythm were estimated by using measurements of functional directionality based on Granger causality [START_REF] Kaminski | Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance[END_REF], and the significance of directional strength evaluated with surrogated data. Since previous evidence has shown that the presence of the allele ε4 in the apoliprotein E (ApoE) not only increases the risk of developing AD [START_REF] Caselli | Cognitive domain decline in healthy apolipoprotein E epsilon4 homozygotes before the diagnosis of mild cognitive impairment[END_REF][START_REF] Wang | Olfactory identification and apolipoprotein E epsilon 4 allele in mild cognitive impairment[END_REF][START_REF] Corder | Gene dose of apolipoprotein E type 4 allele and the risk of Alzheimer's disease in late onset families[END_REF]] but also correlates with changes in the cortical dynamics of alpha rhythm [START_REF] Babiloni | Sources of cortical rhythms change as a function of cognitive impairment in pathological aging: a multicenter study[END_REF], both probability of source activation and direction of synaptic flow between EEG-alpha sources were compared between carriers and non-carriers of the allele ε4 in controls and MCI groups.

Methods

Subjects

Twenty MCI patients (9 females, mean age: 66.8 ± 4.7 yr) and 20 cognitively normal volunteers (10 females, mean age: 68.4 ± 6.1 yr) were recruited from the local community and the Dementia Unit of the Neurology Service at the University Hospital Virgen del Rocio, respectively. MCI patients and control subjects were matched in educational years and handedness. Demographics and neuropsychological profile of healthy elderly and The diagnosis of MCI was based on consensus criteria [START_REF] Petersen | Mild cognitive impairment. Clinical characterization and outcome[END_REF]: (i) subjective memory complaints confirmed by the informant (ii) objective memory decline on neuropsychological tests evidenced by scores ≥ 1.5 standard deviations below the ageappropriate mean, (iii) clinical dementia rating (CDR) global score of 0.5 (questionable dementia), (iv) normal independence function both judged clinically and by means of a standardized scale for the activities of daily living, and (v) not meeting DSM-IV criteria for dementia. Depression illness was excluded by clinical interview and the Geriatric Depression Scale (GDS) of Yesavage (shorter form). The cutoff to be included in the study was 0-5. The diagnosis of MCI was finally based on a clinical consensus following evaluation in the Dementia Unit by a senior neurologist and a clinical neuropsychologist.

Inclusion criteria for the healthy elderly group were (i) no subjective memory complaints corroborated by neuropsychological exploration, (ii) CDR global score of 0 (no dementia), and (iii) normal independent function both judged clinically and by means of a standardized scale for the activities of daily living. None of them had a history of neurological, psychiatric disorders and/or major medical illness.

To avoid interference on EEG recordings and neuropsychological performance, the uptake of pharmacological compounds known to significantly affect any cognitive domain was considered cause of exclusion, both in healthy controls and MCI patients.

Apolipoprotein E genotyping

Genomic DNA was extracted from peripheral blood using short proteinase K digestion.

ApoE genotype was determined by PCR (polymerase chain reaction) amplification of the polymorphic fragment of the ApoE gene, and digested by the restriction enzyme CfoI using protocols previously described [START_REF] Wenham | Apolipoprotein E genotyping by one-stage PCR[END_REF].

Electroencephalographic recordings and pre-processing

Continuous EEG recordings were obtained in all participants in relaxed wakefulness with eyes closed between 9-10 AM. Vigilance level was constantly controlled to avoid intermittent alpha oscillatory behavior caused by the drowsiness state. EEG was referenced to linked mastoids from 59 scalp locations (Fp1, Fp2, AF7, AF3, AFz, AF4, AF8, F7, F5, F3, F1, Fz, F2, F4, F6, F8, FT7, FC5, FC3, FC1, FCz, FC2, FC4, FC6, FT8, T7, C5, C3, C1, Cz, C2, C4, C6, T8, TP7, CP5, CP3, CP1, CPz, CP2, CP4, CP6, TP8, P7, P5, P3, P1, Pz, P2, P4, P6, P8, PO7, PO3, POz, PO4, PO8, O1, and O2) according to the International 10-20 system. Vertical and horizontal ocular movements were also simultaneously recorded. Electrode-scalp impedance was always kept below 5 KΩ. All electrophysiological variables were amplified (BrainAmp MR, Brain Vision ® ), filtered (0.1-100 Hz bandpass), digitized (250 Hz, 16-bit resolution), and stored in digital format for off-line analysis. were manually identified and eliminated. A total of 150 s of artifact-free EEG containing alpha rhythm were then available for each participant. The selected epochs were filtered within 6-13 Hz using a real and phase linear bandpass filter.

Determining the temporal dynamics of EEG-alpha sources

Synaptic flows of macroscopic activity between neuronal EEG sources (in our case EEGalpha sources) were modeled using a multivariate autoregressive (MVAR) model [START_REF] Astolfi | Estimate of causality between independent cortical spatial patterns during movement volition in spinal cord injured patients[END_REF][START_REF] Astolfi | Assessing cortical functional connectivity by linear inverse estimation and directed transfer function: simulations and application to real data[END_REF][START_REF] Supp | Directed cortical information flow during human object recognition: Analyzing induced EEG gamma-band responses in brain's source space[END_REF]. This model implies that the timevarying neural current density responsible for the EEG scalp potentials is due to a set of K intracranial signal generators whose mutual dynamics

[ ] T K t s t s t ) ( ),..., ( ) ( 1 = s
can be approximated by:

) ( ) ( ) ( ) ( 1 t t t p n s B s s + - = ∑ = τ τ τ (1)
where p is the order of the MVAR model and ) (τ s

B

is the K x K coefficient matrix corresponding to time-lag τ . The components of the multivariate residual process

[ ] T K t n t n t ) ( ),..., ( ) ( 1 = n
are commonly assumed to be uncorrelated Gaussian processes.

We imposed two additional assumptions on this EEG model. Specifically, we assume that

(1) the components of the residual process are not perfectly Gaussian and (2) the individual components of the multivariate residuals are mutually independent, i.e. they have zero second-order cross-correlations and zero higher-order dependencies.

Conceptually, this is equivalent to assuming that the residuals of the model are the only this approximation implicitely assumes that zero-lag neural flows of information are negligible in comparison to time-lagged flows. This is plausible considering the presence of axonal propagation delays [START_REF] Freeman | Characteristics of the synchronization of brain activity imposed by finite conduction velocities of axons[END_REF]. Within this framework, knowledge of

) (τ s B
is sufficient to characterize the directionality of synaptic flow between brain sources using e.g. the directed transfer function (DTF) [START_REF] Kaminski | Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance[END_REF].

A fundamental problem when studying functional relationships between brain areas is that connectivity between scalp EEG signals is not equivalent to connectivity between their underlying neural sources. This occurs because scalp EEG potentials do not exclusively reveal averaged postsynaptic activity from localized cortical regions beneath one electrode but the superposition of all active coherent neural sources located anywhere in the brain, due to conduction in the head volume [START_REF] Nunez | Electric fields of the brain: The neurophysics of EEG[END_REF][START_REF] Hoechstetter | BESA source coherence: A new method to study cortical oscillatory coupling[END_REF]. Based on the quasi-static approximation of electrical conduction in the head [START_REF] Malmivuo | Bioelectromagnetism: Principles and applications of bioelectric and biomagnetic fields[END_REF] and using Eq. (1), we can assume that the scalp EEG measurements, at M electrodes and at time instant t, denoted by x(t) =

[x 1 (t),…,x M (t)] T , can be modeled as:

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 1 1 t t t t t t p p v x B Φn x Φ ΦB Φs x x s + - = + - = = ∑ ∑ = + = τ τ τ τ τ τ (2)
where Φ is an M x K leadfield matrix containing the coefficients that project activity from each brain source to the scalp electrodes and + denotes Moore-Penrose pseudoinversion. Eq. ( 2) clearly shows that although the scalp EEG measurements are an . Many previous studies have assumed exactly the opposite [START_REF] Bernasconi | On the directionality of cortical interactions studied by structural analysis of electrophysiological recordings[END_REF][START_REF] Baccalá | Partial directed coherence: A new concept in neural structure determination[END_REF]Kus et al., 2004], which makes their estimated functional connectivity diagrams questionable. Therefore, to obtain a valid characterization of the functional connectivity underlying EEG-alpha oscillations we need to estimate: (1) the MVAR coefficients driving the EEG-alpha sources, i.e.

) (τ s B
and (2) the leadfield matrix Φ . The former is needed to compute the DTF between EEG-alpha sources whereas the latter is required for estimating the intracranial location of the corresponding EEG-alpha generators.

To overcome volume conduction effects and to obtain the true underlying functional connectivity and the location of the brain areas involved in the generation of EEG alpha we used MVAR-EfICA, a novel non-invasive approach which has shown its superiority with respect to several other competing approaches [START_REF] Gómez-Herrero | Measuring directional coupling between EEG sources[END_REF]. A schematic representation of the analysis methology used in this study is shown in Figure 1, and consists of the following steps:

Principal component analysis (PCA)

PCA was applied to all healthy elderly subjects in order to reduce the effects of measurement noise and to remove second-order instantaneous cross-correlations caused by volume conduction [START_REF] Richards | Recovering dipole sources from scalp-recorded event-relatedpotentials using component analysis: principal component analysis and independent component analysis[END_REF]

= = = = -+ = -+ ∑ ∑ -1 s x Cx CΦB (CΦ) x CΦn B x r (3)
where C is a K x M matrix implementing the PCA transformation. We used as many PCs as necessary to reconstruct most (99%) of the variance explained by the original EEG signals.

Multivariate Autoregressive (MVAR) modeling

In all healthy elderly subjects, an MVAR model was fitted to the PCs obtained in the previous step using the algorithm ARfit [START_REF] Schneider | Algorithm 808: ARFIT -a matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models[END_REF].

) ( ) ( ˆτ τ PCA PCA B B
≈ denotes the estimated model coefficients and ) ( ) ( ˆt t r r ≈ denotes the estimated model residuals. The model order was 7 in all cases and was automatically selected using Swartz's Bayesian Criterion (SBC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. Namely, we picked the order that guaranteed that greater model orders did not significantly reduce the SBC, i.e. the order for which the reduction in the SBC had reached 90% of the maximum reduction achievable within the tested range of model orders (model orders tested in the present study ranged from 2 to 30). This modified SBC criterion was motivated by the fact that, for our EEG datasets, both SBC and Akaike's information criterion dropped 12 monotonically as the model order increased, which is in close correspondence with previous EEG studies [START_REF] Supp | Directed cortical information flow during human object recognition: Analyzing induced EEG gamma-band responses in brain's source space[END_REF][START_REF] Brovelli | Beta oscillations in a large-scale sensorimotor cortical network: Directional influences revealed by granger causality[END_REF].

Independent component analysis (ICA)

As mentioned above, we assumed that the components of the multivariate residual process driving the EEG-alpha sources in Eq. ( 1) (denoted by ) (t n ) are mutually independent. Moreover, Eq. ( 3) shows that the residual process found in the previous step

( ( ) ) ( ) ( t t n CΦ r =
in Eq. ( 3)) is just a linear mixture of the components of ) (t n . Therefore, we applied ICA [Koldovsky et al., 2006] to invert that mixing by estimating a K x K matrix

( ) 1 ˆ-≈ CΦ W
that minimized the mutual dependencies between the components of the multivariate residual process ) (t r . Then, from the relationship

( ) 1 ˆ-≈ CΦ W
we can easily obtain an estimate of the leadfield matrix:

( ) Φ C W Φ ≈ = + ˆ.
Finally, the activation patterns of the underlying EEG-alpha sources were obtained by spatially filtering the scalp EEG measurements:

) ( ) ( ) ( ) ( ˆt t t t s Φs Φ x Φ s ≈ = = + +
. Each of the K rows of matrix + Φ ˆ is a spatial filter retrieving an individual EEG-alpha source. In the following we will show + = Φ F ˆ. Due to the scale indeterminacy inherent to ICA methods, only the waveform (but not the scale) of the EEG-alpha sources can be recovered. That is, the norm of the columns of the estimated leadfield matrix is arbitrary. Note, however, that this scale indeterminacy is irrelevant for either localizing the EEG sources or computing their mutual directional interactions. et al., 2006] involve stochastic optimization, which raises concerns regarding reliability when applying them to real data [START_REF] Särelä | Overlearning in marginal distribution-based ICA: Analysis and solutions[END_REF]. To overcome this issue, we identified clusters of ICA estimates (EEG-alpha sources) that were consistently found across random initializations of the ICA algorithm, across random bootstrap surrogates of the input data [START_REF] Himberg | Validating the independent components of neuroimaging time-series via clustering and visualization[END_REF], and (across) healthy elderly subjects [START_REF] Gómez-Herrero | Measuring directional coupling between EEG sources[END_REF]. The validity of the clustering results was assessed using the R-index [START_REF] Himberg | Validating the independent components of neuroimaging time-series via clustering and visualization[END_REF].

Assessing the significance of ICA estimates

Each cluster was uniquely represented by a single centrotype EEG-alpha source. The centrotype was defined as the ICA-estimate showing the maximum sum of similarities to other points in one cluster. The similarity between two ICA estimates was assessed using the cross-correlation coefficient between their spatial distribution of scalp potentials.

Only centrotypes of significant clusters were considered as valid EEG-alpha sources. A cluster was considered as significant if it contained ICA-estimates from at least 70% of the analyzed healthy elderly subjects (high inter-subject repeatability) and from at least 90% of the ICA-runs corresponding to those subjects (high intra-subject reliability). Each centrotype EEG-alpha source was characterized by two spatial features: (1) the spatial filter (a row vector of M coefficients f i =[f i,1 ,…, f i,M ]) retrieving the temporal activation of the EEG-alpha source from the scalp EEG measurements, i.e. the corresponding row of the pseudoinverse of the leadfield matrix ( + Φ ˆ) associated with the centrotype ICA estimate, and (2) the scalp distribution of potentials (a column vector of M coefficients

g i =[g i,1 ,…, g i,M
] T ) generated by the centrotype EEG-alpha source, i.e. the corresponding column of the leadfield matrix Φ ˆ associated with the centrotype ICA estimate. Figure 2 shows the analysis steps followed to assess the significance of ICA estimates used to 15 alpha: 10.1-12.5 Hz) for each independent alpha-EEG source, both in healthy elderly subjects and MCI patients. This subdivision of the alpha band was performed considering previous evidence that highlighted the differential effects of normal aging and mild dementia on different portions of the classic alpha band [START_REF] Moretti | Individual analysis of EEG frequency and band power in mild Alzheimer's disease[END_REF].

Computing directional connectivity between EEG-alpha sources

The obtained centrotype spatial filters

[ ] T T K T control f f F ,..., ˆ1 =
were employed to determine the temporal activation patterns of each EEG-alpha source in individual subjects (both in controls and MCI). For an arbitrary subject (control or MCI) the activation of the i th centrotype was obtained as s i (t)=f i x(t), where x(t) is a column vector with the scalp EEG measurements of that subject at time instant t. Subsequently, an MVAR model was fitted to these centrotype activations using the same estimation algorithm and model order selection method as described above. Then, the normalized DTF [START_REF] Kaminski | Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance[END_REF] was applied to determine directional causal relationships between these temporal activation patterns. Normalized values of the DTF from the j th to the i th EEG-alpha source at a certain frequency f ranged from 0 to 1. Values close to 1 indicated that the neural dynamics of the source i are mostly explained by a directed flow of synaptic activity from the source j. DTF values close to 0 denoted that there was no significant synaptic flow from source j to source i at that EEG frequency. The DTF between EEG-alpha sources was computed for both healthy elderly subjects and MCI patients.

Computing the alpha peak frequency The alpha peak frequency was computed for the temporal activation of each EEG-alpha source. Spectral power was estimated using the Welch's averaged modified periodogram (4 s segments, 1 Hz resolution, 50% overlapping, and Hanning windowing) as implemented in MATLAB ® v. 7.4 (The MathWorks, Inc.). Alpha peak frequency was identified as the maximum spectral power value within the defined alpha band (7.5-12.5

Hz).

Statistical analysis

In order to assess the statistical significance of the normalized DTF values obtained, the surrogates-based procedure proposed by [START_REF] Theiler | Testing for nonlinearity in time series: The method of surrogate data[END_REF] was used. DTF values below the confidence threshold (p < 0.01) were considered not significant and removed from the analysis.

Alpha peak frequency, probability of source activation and directionality of the synaptic flow were evaluated with two types of mixed ANOVAs: i) without including covariates;

and ii) with scores provided by the Mini Mental State Examination (MMSE) and age as covariates (if they were significant, otherwise this ANOVA was not performed). Each covariate was tested independently for significance, and only significant covariates (p < 0.05) were included in the final ANOVA. Mauchly's W was computed to check for violations of the sphericity assumption. When Mauchly's W test was significant, the Greenhouse-Geisser correction was applied to all repeated-measures analyses (original degrees of freedom and the corresponding epsilon value were reported; p value reflects the epsilon correction). Homogeneity of variance was evaluated with Levene's test. The In order to confirm previous reports of alpha slowing in MCI patients, alpha peak frequency was evaluated with a mixed ANOVA including, as the within-subject factor, the electrodes showing the maximum amplitude within the alpha band for each EEGalpha centrotype as determined by the MVAR-EfICA procedure, and group (healthy elderly and MCI) as the between-subject factor.

Activation probabilities of different EEG-alpha sources were analyzed for the lower and upper alpha band separately in controls and MCI patients. The ANOVA included EEGalpha source (see below the three different sources that achieved statistical significance)

as the within-subject factor and group (healthy elderly and MCI) as the between-subject factor. A similar statistical design was applied to evaluate the directionality of the synaptic flow between EEG-alpha sources.

In the present study, a higher number of MCI patients showed the allele ε4 in the ApoE (p < 0.04) when compared to healthy elderly subjects. As the presence of this allele has been associated with an increased risk of conversion to AD [START_REF] Caselli | Cognitive domain decline in healthy apolipoprotein E epsilon4 homozygotes before the diagnosis of mild cognitive impairment[END_REF][START_REF] Wang | Olfactory identification and apolipoprotein E epsilon 4 allele in mild cognitive impairment[END_REF][START_REF] Corder | Gene dose of apolipoprotein E type 4 allele and the risk of Alzheimer's disease in late onset families[END_REF], a new set of mixed ANOVAs was performed in MCI patients to evaluate its effect on each dependent variable separately. In this particular case, genotype (presence or absence of ε4) instead of group was introduced in the different ANOVAs as the between-subject factor.

Results

Spatial patterns of EEG-alpha rhythm

The R-index statistic indicated that 10 clusters described best the observed data among all the possible partitions obtained by the hierarchical clustering algorithm. Out of those 10 clusters, only three (clusters 7, 9 and 10) satisfied the established statistical criteria for being considered EEG-alpha sources (Figure 2). First, clusters 7, 9 and 10 were by far the largest with 1125, 1050 and 1050 ICA-estimates, respectively. Second, they contained estimates from 75%, 70%, and 70% of the healthy elderly subjects, respectively, indicating a high inter-subject repeatability.

Intracranial location of EEG-alpha centrotypes

Coordinates of the Talairach-Tournoux atlas were employed to determine the anatomical location of the single electric dipole best explaining the scalp distribution of each centrotype ICA-estimate for the three significant clusters mentioned above. The electric dipole corresponding to cluster 7 was located in caudal regions of the thalamus (x=9, y=-25, z=9), cluster 9 showed its maximum activation in the precuneus (x=2, y=-60, z=28)

and cluster 10 was located in the middle occipital gyrus, within the limits of the cuneus (x=11, y=-97, z=13). Figure 3 displays the brain sources of the alpha rhythm based on the selected cluster.

Alpha peak frequency in healthy aging and MCI

The alpha peak frequency corresponding to the thalamic, precuneus and cuneus EEGalpha sources was obtained for controls and MCI patients. As there were no significant On the average, the alpha peak frequency in MCI patients was 0.6 Hz (SE = 0.17) lower than in controls. This abnormal slowing was evident for the three EEG-alpha sources (see Table 2), independent of the presence of the allele ε4.

EEG-alpha source activation in healthy aging and MCI

Directional connectivity between EEG-alpha sources

The two-way ANOVA (with MMSE as covariate) performed to evaluate the directionality of the synaptic flow between thalamus, cuneus and precuneus within the lower alpha band showed a significant group effect and a significant directionality x group interaction (p < 0.04). In agreement with the enhanced activation of the thalamic alpha source observed in the MCI group, the synaptic transfer was also abnormally facilitated in these patients (~24% higher) relative to controls, particularly from thalamus to cortex (~38% higher in MCI; p < 0.008). This effect was similar in MCI ε4-carriers and non-carriers. No differences between controls and MCI patients were obtained for cortico-cortical interactions in the lower alpha band. Analyses only revealed a higher flow of synaptic activity from cuneus to precuneus than in the opposite direction [F(1,38) = 22.61; p < 0.00003]. Figure 5 shows comparisons of DTF values obtained between brain sources participating in the generation of lower alpha oscillations.

For the upper alpha subdivision, the three-way ANOVA (with age as covariate) revealed a significant group effect [F(1,37) = 9,32; p < 0.005]. In particular, MCI patients showed a decreased bidirectional synaptic flow of about 32% between thalamic and cortical sources when compared to healthy elderly subjects. Analyses also showed a significant interaction between the factors group and age [F(2,37) = 8.17; p < 0.002]. In agreement with this finding, we found that the bidirectional synaptic flow between thalamocortical sources of the upper alpha band significantly decreased with age in healthy elderly subjects but not in MCI patients. The scatter plot shown in Figure 6 Like the remaining dependent variables, no differences between controls and MCI were found when the influence of genotype on directional connectivity was compared in both subdivisions of the alpha band.

Discussion

In the present study, we have examined whether alterations in the directionality of the synaptic flow between EEG-alpha sources were able to distinguish amnesic MCI patients from cognitively intact elderly individuals. In addition to the typical EEG-alpha slowing associated with neurodegenerative disorders, MCI patients also showed an abnormal enhancement of the thalamic source activation together with a facilitation of synaptic transmission from thalamus to cortex restricted to the lower alpha range (7.5-10 Hz).

These findings reveal, for the first time, that thalamic dysfunctions account for abnormalities of the alpha rhythm in MCI patients, thereby providing new insights into the role of thalamocortical integrity in the preclinical stages of AD.

Thalamocortical generation of alpha oscillations in healthy aging

Animal studies have shown that much of the patterning of the alpha rhythm is subject to controlling influence from synchronizing pacemakers within thalamocortical networks [Lopes da Silva, 1991;Steriade et al., 1990a]. Our results provide additional support for This finding is in agreement with results obtained by Schreckenberger and colleagues [START_REF] Schreckenberger | The thalamus as the generator and modulator of EEG alpha rhythm: a combined PET/EEG study with lorazepam challenge in humans[END_REF] who reported a positive correlation between EEG alpha power and metabolism of the lateral thalamus as well as occipital cortex (cuneus) and adjacent parts of the parietal cortex (precuneus) in humans.

The main obstacle to estimate the sources of electrical activity in the thalamus comes from the spatial arrangement of neurons in thalamic nuclei resulting in a closed field [Lopes da [START_REF] Da | Biophysical aspects of EEG and magnetoencephalogram generation[END_REF]. Indeed, the three-dimensional symmetry of the dendritic organization in thalamic neurons makes it difficult to determine distant neuronal contributions to scalp-recorded electric potentials [START_REF] Druckman | Brain and cognition. Some new technologies[END_REF]].

In contrast, previous studies have determined unambiguous thalamic sources of scalp EEG phenomena [Gobbelé et al., 2004;[START_REF] Trujillo-Barreto | Bayesian model averaging in EEG/MEG imaging[END_REF][START_REF] Isaichev | Rhythmic alpha-activity generators in the human EEG[END_REF] suggesting that thalamic oscillations can be captured under specific experimental manipulations. Further evidence has shown that the highest degree of cross-correlation in the alpha range between the lateral pulvinar nucleus and different cortical regions was observed in the precentral area [START_REF] Fujita | Cross-correlation analysis of the lateral pulvinar and scalp EEG in man[END_REF] where the electrode that best represented the thalamic centrotype was located in the present study.

Our results also revealed that thalamocortical synaptic transmission remained alike from thalamus to cortex and vice versa. These findings, together with results from human studies employing 3D equivalent dipole modeling [START_REF] Isaichev | Rhythmic alpha-activity generators in the human EEG[END_REF]Başar et al., 1997;[START_REF] Pritchard | Segregation of the thalamic alpha rhythms from cortical alpha activity using the Savit-Green S-statistic and estimated correlation dimension[END_REF][START_REF] Schurmann | Electroencephalogram alpha (8-15 Hz) responses to visual stimuli in cat cortex, thalamus, and hippocampus: a distributed alpha network?[END_REF], support the notion that complex interactions between local and non-local EEG sources, instead of a single or [START_REF] Nunez | Spatial-temporal structures of human alpha rhythms: theory, microcurrent sources, multiscale measurements, and global binding of local networks[END_REF].

Thalamic dysfunction at generating alpha oscillations in MCI

Alpha slowing together with the enhanced activation of the thalamic source for the lower alpha band indicates a dysfunctional state of this subcortical structure in MCI patients.

The alpha slowing has been consistently reported in EEG studies of AD [e.g., [START_REF] Moretti | Individual analysis of EEG frequency and band power in mild Alzheimer's disease[END_REF], but results are contradictory in MCI patients. To date, only one study has shown that the mean frequency peak of alpha rhythm was significantly lower in MCI when compared to healthy elderly subjects, and significantly higher when compared to AD patients [START_REF] Fernandez | MEG spectral profile in Alzheimer's disease and mild cognitive impairment[END_REF]. Other studies, however, either found no differences between MCI and controls [START_REF] Huang | Discrimination of Alzheimer's disease and mild cognitive impairment by equivalent EEG sources: a crosssectional and longitudinal study[END_REF][START_REF] Jelic | Quantitative electroencephalography power and coherence in Alzheimer's disease and mild cognitive impairment[END_REF], or only showed a slight tendency toward a reduction of the alpha peak frequency [START_REF] Osipova | Source estimation of spontaneous MEG oscillations in mild cognitive impairment[END_REF].

A chronically slowed alpha rhythm during resting wakefulness has been suggested as a definitive indicator of underlying neurological pathology [START_REF] Niedermeyer | Alpha rhythms as physiological and abnormal phenomena[END_REF], but mechanisms by which alpha oscillations slow down have not yet been revealed in humans. Evidence from in vitro preparations of the cat lateral geniculate nucleus has shown that a moderate activation of metabotropic glutamate receptors can lead to a deceleration of alpha waves through a slowing of high-threshold bursting neurons interconnected by gap junctions [START_REF] Hughes | Synchronized oscillations at alpha and theta frequencies in the lateral geniculate nucleus[END_REF]. Should this cellular mechanism be affected to some extent by the pathogenesis of AD, the EEG-alpha rhythm would also be affected. In support of this hypothesis, amyloid deposits are not only present in the neocortex but also in deep gray matter structures as the thalamus [START_REF] Rudelli | Morphology and distribution of Alzheimer neuritic (senile) and amyloid plaques in striatum and diencephalon[END_REF]; plaques may disrupt the glutamatergic system by affecting simultaneously different signaling pathways [reviewed by [START_REF] Parameshwaran | Amyloid beta peptides and glutamatergic synaptic dysregulation[END_REF]. Should this also be the case in MCI, inhibition of glutamate receptors might account for the slowing of alpha oscillations and for the abnormal enhanced synchrony between in-phase and anti-phase groups of high-threshold bursting neurons of the thalamus [for a description of electrophysiological properties of the thalamocortical neurons see [START_REF] Hughes | Just a phase they're going through: the complex interaction of intrinsic high-threshold bursting and gap junctions in the generation of thalamic alpha and theta rhythms[END_REF].

Other authors [e.g., [START_REF] Llinás | Thalamocortical dysrhythmia: A neurological and neuropsychiatric syndrome characterized by magnetoencephalography[END_REF] have reported a shift from a normal alpha rhythm to a robust low-frequency theta across a diverse selection of neurological and psychiatric disorders, which strongly correlated with the presence of positive symptoms. Based on these results, it could be postulated that the thalamic dysfunction at generating alpha oscillations might not be specific to AD. However, we would like to bring to the attention of the reader two important factors that could have influenced the results reported by Llinás and colleagues. First, controls and patients were not of the same age (controls ranged between 24 and 45 yrs; patients ranged between 28 and 73 yrs) which might be an important issue of confusion in this study because alpha oscillations become slower with age. And second, results were based on averaging cerebral oscillations of four patients diagnosed with Parkinson's disease, one patient with tinnitus, two patients with neurogenic pain, and two patients with major depression. The latter approach impeded statistical confirmation of the decrease of alpha and increase of theta in each patient population. Although the present study could not determine whether alpha slowing is specific to AD, it describes for the first time specific abnormalities underlying alpha

Thalamocortical dysfunctions at generating alpha oscillations in MCI

Animal studies have shown that the posterior parietal gyri and the occipital lobes receive major inputs from caudal portions of the thalamus, including pulvinar and the lateral geniculate nucleus [START_REF] Steriade | Thalamus[END_REF]. Changes in the behavior of large ensembles of thalamocortical relay cells or alterations in this complex circuitry are reflected in the cortical EEG oscillations [Steriade et al., 1990b], providing an in vivo global measurement of thalamocortical integrity in a variety of cerebral diseases. In agreement with this idea, MCI patients showed abnormal synaptic transmission between thalamus and cortex in generating alpha oscillations. In particular, synaptic flow from thalamus to cortex was facilitated for the lower alpha band in MCI relative to healthy elderly subjects;

whereas for the upper alpha band, transmission between thalamus and cortex decreased in both directions. These changes in thalamocortical functional connectivity may reveal modifications in intrinsic membrane properties of the thalamic neurons or synaptic dysfunctions between excitatory and inhibitory neuronal populations involved in the generation of rhythmic EEG activity. Whereas a progressive decrease of thalamocortical transmission in the upper alpha band is expected with increasing normal aging (see Figure 6), the parallel gradual increase of thalamocortical synchrony in the lower alpha might signal the route to functional decline in MCI, these patients becoming more White matter damage [e.g., [START_REF] Duan | White matter damage of patients with Alzheimer's disease correlated with the decreased cognitive function[END_REF][START_REF] Bozzali | White matter damage in Alzheimer's disease assessed in vivo using diffusion tensor magnetic resonance imaging[END_REF] and aberrant patterns of functional connectivity [e.g., [START_REF] Wang | Altered functional connectivity in early Alzheimer's disease: A resting-state fMRI study[END_REF][START_REF] Rombouts | Altered resting state networks in mild cognitive impairment and mild Alzheimer's disease: an fMRI study[END_REF][START_REF] Greicius | Default-mode network activity distinguishes Alzheimer's disease from healthy aging: evidence from functional MRI[END_REF] between different brain regions have been suggested as accounting for the cognitive deficits observed in AD patients. Thus, studies using resting-state functional magnetic resonance imaging (fMRI) have shown that abnormal functional connectivity patterns in AD not only include decreased correlations between anterior and posterior regions of the brain but also increased positive correlations within the prefrontal lobe, parietal lobe and occipital lobe [START_REF] Wang | Altered functional connectivity in early Alzheimer's disease: A resting-state fMRI study[END_REF]]. Consistent with these results, it has been found that MCI patients show less deactivation of the default mode network (involving anterior frontal, precuneus and posterior cingulate cortex) as compared to healthy elderly subjects, but more deactivation when compared to AD patients [START_REF] Rombouts | Altered resting state networks in mild cognitive impairment and mild Alzheimer's disease: an fMRI study[END_REF]. Taken together, metabolic and electrophysiological results

indicate that abnormal enhancement of functional connectivity could be used as surrogate markers at improving the ability to differentiate patients with high risk of AD from healthy elderly subjects.

Relationship between ApoE4 and thalamocortical dysfunction in MCI

Gene dosage of the ApoE ε4 allele is a major risk factor for familial AD of late onset [START_REF] Corder | Gene dose of apolipoprotein E type 4 allele and the risk of Alzheimer's disease in late onset families[END_REF]. But the interaction between this genotype and the integrity of thalamocortical circuits involved in alpha generation has not been explored to date in either healthy aging or MCI patients. The present study showed that a significantly higher confirming previous findings that the presence of the ε4 allele is more frequent in persons with MCI than in healthy elderly [START_REF] Caselli | Cognitive domain decline in healthy apolipoprotein E epsilon4 homozygotes before the diagnosis of mild cognitive impairment[END_REF][START_REF] Wang | Olfactory identification and apolipoprotein E epsilon 4 allele in mild cognitive impairment[END_REF]. However, no relationship was found between directional changes of thalamocortical synaptic flow during alpha generation and the presence of the ε4 allele in the present study. Despite the increased risk of ε4 carriers to develop AD, MCI population also contains ε4 carriers who may or may not develop a neurodegenerative condition. The intrinsic heterogeneity of this group and the small size of our sample population may account for the absence of significant relationships between changes in brain mechanisms of alpha generation and the presence of the ε4 allele in our MCI sample. Alternatively, it might happen that thalamic facilitation within the lower-alpha band was not specific to MCI but a more global marker of neurological/psychiatric dysfunction, accounting for the lack of correlation with the APOE ε4 genotype that seems to be more specific to AD disease.

Finally, it could be further hypothesized that once changes in thalamocortical dynamics underlying EEG-alpha slowing appear, the carrier status of apoe4 could no longer be relevant. Experimental support for this hypothesis comes from large-community based studies that did not find a significant influence of the APOE ε4 on the rate of AD progression in cognitive or functional domains [START_REF] Kleiman | Apolipoprotein E epsilon4 allele is unrelated to cognitive or functional decline in Alzheimer's disease: retrospective and prospective analysis[END_REF][START_REF] Hoyt | Individual growth curve analysis of APOE epsilon 4-associated cognitive decline in Alzheimer disease[END_REF].

In summary, results from the present study suggest that abnormal functioning of thalamocortical networks accounts for differences in EEG-alpha slowing between MCI patients and healthy elderly subjects. Future follow-up studies will confirm whether or not these changes parallel cognitive decline, thereby providing a reliable in vivo functional marker to predict which patients diagnosed with MCI will convert to AD. Figure 1. Flow chart displaying a scheme of the main analysis steps followed in this study. Note that the MVAR-EfICA methodology (represented by the inside of the dotted square line) is a combined approach using PCA, MVAR and ICA which was initially applied to the EEG datasets of the healthy old subjects to remove effects of volume conduction from multichannel EEG recordings (see Methods section). Next, significant clusters of ICA estimates were represented, identified by independent centrotypes, and located in the 3D brain space of healthy elderly and MCI patients.

F

  . Study protocols were previously approved by the Ethical Committee for Clinical Investigations at the University Hospital Virgen del Rocio, and the Ethical Committee for Human Research at the University Pablo de Olavide.

  conditions that may affect brain structure or function were also excluded.

  prominent ocular, muscular and/or any other type of artifacts

  new" or intrinsic information in each brain generator and that the MVAR coefficient matrices are responsible for synaptic transfer between EEG sources. Indeed,

  the corresponding MVAR coefficients are, in general, quite different from the MVAR coefficients driving the underlying brain sources, i.e.

  and specifically EfICA [Koldovsky

  intracranial location of EEG-alpha sources and the directionality of synaptic flow between them.The features of the K centrotype EEG-alpha sources underlying the healthy elderly population was therefore summarized by a single leadfield matrixDetermining the intracranial location of EEG-alpha centrotypesWe used the global leadfield matrix obtained for the healthy elderly population ( control Φ ˆ)to determine the intracranial locations of the EEG-alpha centrotypes. Each column of control Φ ˆ defines the distribution of scalp potentials generated by each individual EEGalpha source. The brain locations corresponding to these potentials were obtained using a 3D brain space. The activation probability of each grid element in this brain space was computed by the standardized weighted low-resolution brain electromagnetic tomography (swLORETA)[START_REF] Palmero-Soler | swLORETA: a novel approach to robust source localization and synchronization tomography[END_REF]. A realistic head model of three layers (scalp, skull and brain with conductivities of 0.33, 0.0042, and 0.33, respectively) was used for this purpose. Although previous studies have shown that inter-subject variability in the electric resistivities of these layers may affect the error sources in EEG modeling[START_REF] Gonçalves | In vivo measurement of the brain and skull resistivities. Using an EIT-based method and the combined analysis of SEF/SEP data[END_REF][START_REF] Oostendorp | The conductivity of the human skull: Results of in vivo and in vitro measurements[END_REF], we did not have the possibility of measuring them in vivo for each participant in the present study. Source reconstruction solutions were projected onto the 3D MR images of the Collin's brain provided by the Montreal Neurological Institute. Probabilities of source activation based on Fisher's Ftest were obtained for two subdivisions of the alpha band (lower alpha: 7.5-10 Hz; upper

  -Keuls or the Games-Howell procedures were applied for multiple post hoc comparisons depending on population variances.

  the MMSE scores and age on the alpha peak frequency, these two covariates were not included in the ANOVA. As expected, MCI patients showed a significant decrease in the alpha peak frequency as compared to controls [F(1,38) = 7.02; p < 0.02].

Figure 4

 4 Figure4shows source activation probabilities of the lower and upper alpha subdivisions

  depicts the degree of linear relationship between these two factors in both groups as derived from Pearson's correlation. As shown, regression fitness reached statistical significance only for healthy p < 0.01) suggesting that the strength of the thalamocortical synaptic flow within the upper alpha decreases with age as a result of normal aging processes. No significant effects were found for cortico-cortical interactions in the upper alpha band in any of the groups.

  high-risk population for the development of AD. It remains to ascertain whether or not the enhanced facilitation of thalamic and thalamocortical synaptic transmission is common to abnormal alpha generation in other neurological and neuropsychiatric disorders.

  development of AD. Future longitudinal studies in MCI patients who will later convert to AD are necessary to confirm this hypothesis.

  patients had the allele ε4 in the ApoE when compared to controls,

Figure legends Figure 1 .

 legends1 Figure legends

Figure 2 .

 2 Figure 2. Analysis protocol followed to assess the significance of ICA estimates at

Figure 3 .

 3 Figure 3. Localization of electric dipole sources for the three significant EEG-alpha

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Mean activation probability of EEG-alpha sources in healthy elderly subjects

Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 23456 Figure1. Flow chart displaying a scheme of the main analysis steps followed in this study. Note that the MVAR-EfICA methodology (represented by the inside of the dotted square line) is a combined approach using PCA, MVAR and ICA which was initially applied to the EEG datasets of the healthy old subjects to remove effects of volume conduction from multichannel EEG recordings (see Methods section). Next, significant clusters of ICA estimates were represented, identified by independent centrotypes, and located in the 3D brain space of healthy elderly and MCI patients. Finally, the directed transfer function (DTF) was applied to determine directional causal relationships between EEG-alpha sources obtained from temporal activation patterns deduced from centrotype spatial filters. 160x181mm (300 x 300 DPI)

  We found that thalamus together with the cuneus and precuneus were the main generators of alpha oscillations under resting EEG conditions.

	Page 21 of 50	Human Brain Mapping
	this hypothesis in humans.
	F o
	r
	P
	e e r
		R e v i e w
	John Wiley & Sons, Inc. 22

Page 22 of 50 John Wiley & Sons, Inc.

	Human Brain Mapping
	multiple isolated neural generators, are responsible for the genesis of the human alpha
	rhythm in normal aging
	F o
	r
	P
	e e r
	R e v i e w
	23

Table 2 .

 2 Means and standard errors of alpha peak frequency in healthy elderly subjects and MCI patients, carriers and non-carriers of ε4. Note that the alpha peak frequency was measured at significant centrotypes of the thalamus, cuneus, and precuneus.

	Controls	9.47 ± 0.19	9.62 ± 0.18	9.57 ± 0.22
	MCI
	ε4 carriers	8.67 ± 0.22	8.83 ± 0.27	9.22 ± 0.32
	ε4 non-carriers	8.77 ± 0.20	9.09 ± 0.24	9.00 ± 0.29

John Wiley & Sons, Inc.Human Brain Mapping

Acknowledgements

This research was supported by grants to JLC from the European Union (FP6-2005-NEST-Path 043309), Spanish Ministry of Education and Science (SAF2005-00398), and Regional Ministry of Innovation, Science and Enterprise, Junta de Andalucia (CTS-229).