Éric Fédoroff

Jean-François Ponge
email: jean-francois.ponge@wanadoo.fr

Florence Dubs

Federico Fernández-González

Patrick Lavelle

Small-scale response of plant species to land-use intensification

Keywords: Vegetation, Biodiversity, Land-use diversity, Humus form, Land-use intensification

Plant communities are affected by land-use and landscape heterogeneity and can be used as indicators of environmental change. At small-scale, species composition and species richness of plant communities are influenced by local environment and by diaspores from the surroundings. Thus they reflect the influence of both land-use type and land-use diversity. Plant community composition was studied along a gradient of agricultural disturbance in the Morvan Regional Natural Park (Burgundy, France). Six landscape units of 1 km 2 were selected along a range of increasing land-use intensity.

Sixteen 0.2 m² sampling plots per unit were selected according to a grid-based design to estimate the percent cover of all plant species. Pattern analysis showed that local species richness increased from woodland to crop to grassland, and also increased with land-use diversity. Local plant biodiversity was maximized under intermediate disturbance intensity and minimized at low (woodland) and high (crop) disturbance levels.

Introduction

Plant species are good indicators of soil conditions [START_REF] Dahl | On the vegetation types of Norwegian conifer forests in relation to the chemical properties of the humus layer[END_REF][START_REF] Ellenberg | Zeigerwerte von Pflanzen in Mitteleuropa[END_REF][START_REF] Lee | The calcicole-calcifuge problem revisited[END_REF]. Natural or man-induced disturbance in soil conditions is reflected in the species composition and richness of most plant communities [START_REF] Falkengren-Grerup | Soil acidification and vegetation changes in deciduous forest in southern Sweden[END_REF][START_REF] Pettit | Grazing in remnant woodland vegetation: changes in species composition and life form groups[END_REF][START_REF] Roem | Soil acidity and nutrient supply ratio as possible factors determining changes in plant species diversity in grassland and heathland communities[END_REF]. Some plant species are also able to modify their environment between establishment and maturity [START_REF] Miles | The pedogenic effects of different species and vegetation types and the implications of succession[END_REF][START_REF] Muys | Effects of grassland afforestation with different tree species on earthworm communities, litter decomposition and nutrient status[END_REF][START_REF] Bernier | Altitudinal changes in humus form dynamics in a spruce forest at the montane level[END_REF] and interfere with other species through competition, biochemical control and facilitation [START_REF] Facelli | Interactions after death: plant litter controls priority effects in a successional plant community[END_REF][START_REF] Wardle | An ecosystem-level perspective of allelopathy[END_REF][START_REF] Chou | Roles of allelopathy in plant biodiversity and sustainable agriculture[END_REF]. Land-use history may influence soil and vegetation long after conditions for the original establishment of plant species have disappeared [START_REF] Koerner | Influence of past land-use on the vegetation and soils of present day forest in the Vosges mountains, France[END_REF] due to feed-back processes [START_REF] Ponge | Humus forms in terrestrial ecosystems: a framework to biodiversity[END_REF], and persistence of plant organs and diaspores [START_REF] Zobel | Survival of prolonged burial by subalpine forest understory plants[END_REF][START_REF] Thompson | Seeds in soil and worm casts from a neutral grassland[END_REF].

The present paper aims to discern whether local plant biodiversity expressed in terms of species richness is influenced by land-use intensification and/or landscape diversity, and whether above-and below-ground diversity follow the same trend along a gradient of land-use intensification.

Land-use intensification was defined subjectively by the increasing impact of man on the landscape according to the following scale:

1. old-growth forest 2. managed forest 3. mixed landscape dominated by woodland 4. mixed landscape not dominated by a single land-use 5. mixed landscape dominated by pasture

mixed landscape dominated by arable crops

The choice of local plant biodiversity over regional biodiversity [START_REF] Tilman | Diversity and production in European grasslands[END_REF][START_REF] Poiani | Biodiversity conservation at multiple scales: functional sites, landscapes, and networks[END_REF][START_REF] Loreau | Biodiversity and ecosystem functioning: current knowledge and future challenges[END_REF] was due to the need to account for site factors but also for interference and facilitation processes between plant individuals and species (Hester et al., 1991a;[START_REF] Facelli | Interactions after death: plant litter controls priority effects in a successional plant community[END_REF][START_REF] Nilsson | Effects of plant litter species composition and diversity on the boreal forest plant-soil system[END_REF].

Material and methods

The Morvan Regional Natural Park (Burgundy, Central France) is under submontane-atlantic climate with continental influence, a mean annual rainfall of 1000 mm, an temperature of 9°C. The parent rock is granite, with moderate to strong soil acidity, and dominant humus form as mull [START_REF] Perrier | La forêt au Mont Beuvray: phytoscociologie[END_REF][START_REF] Ponge | Collembolan communities as bioindicators of land-use intensification[END_REF]. Forested areas (45% in the region) comprise coniferous plantations [45%, i.e. silver fir, Abies alba Mill., white fir, Abies grandis (Dougl. Ex D. Don) Lindl., Douglas fir, Pseudotsuga menziesii (Mirb.) Franco, spruce, Picea abies (L.) Karst.], and deciduous stands (55%, i.e. common beech, Fagus sylvatica L. and sessile oak, Quercus petraea (Mattusschka) Lieblein).

Agricultural areas (55% in the region) consist of grassland (80%, half of which are permanent pastures and the other half temporary hay meadows) and mostly cereal crops (20%, i.e. wheat, Triticum aestivum L., barley, Hordeum vulgare L., rye, Secale cereale L. and Christmas trees).

Agricultural systems vary from extensive to organic farming whereas rural depopulation and forestry policies influenced the landscape over the last five decades [START_REF] Plaisance | Évolution des paysages forestiers et arborés en Bourgogne et Franche-Comté de l'an mille à nos jours[END_REF]. Many forests have been transformed into coniferous plantations and more recently agricultural land was afforested using national subsidies. In parallel, agriculture shifted from a short (food crop/fallow) to a long rotation system (grassland/cereal crop).

Six representative 1 km 2 landscape units were intuitively selected on the basis of regional knowledge and aerial photographs:

 Unit 1 in a 100-150 year deciduous forest landscape managed by the public sector, with natural regeneration and selection by man.

 Unit 2 in a 50 year coniferous (mainly silver fir) forest landscape managed by the public sector, with artificial regeneration (clear-cut followed by plantation). Unit 2 used to be a deciduous forest similar to Unit 1.

 Unit 3 in a landscape privately afforested 50 years ago with Douglas fir, spruce, and meadows. Remains of the old deciduous forest and some cereal crops were also present.

 Unit 4 in a mosaic of wet meadows, with some 30 years old Douglas fir and spruce.

 Unit 5 in a meadow landscape dominated by organic farming, with a rotation between meadow (6-10 years) and crop (2-4 years), with a few 20-50 year Douglas fir or spruce plantations.

 Unit 6 in an agricultural landscape dominated by conventional cereal crops with some fallow and Christmas tree.

Using aerial photographs a grid of 16 regularly spaced points (200 m mesh size) was projected on each of the six units, and retrieved using a calibrated GPS system. Each sample point was indicated by a central post used as a reference by all teams participating to the research program which involved springtails, lichens, butterflies, birds, soil macroinvertebrates, carabids, and remote sensing. The random selection of sample points was not the best method to detect vegetation patterns [START_REF] Gillison | The use of gradient directed transects of gradsects in natural resource surveys[END_REF], but allowed all scientists to use the same sampling design, and comparisons between countries and taxonomic groups to be made. Vegetation was studied from July to November 2001 across a range of six nested areas (50, 25, 12.5, 5, 1 and 0.2 m 2 ), at fixed distance and orientation from the central post. Small-scale study of vegetation was conducted on 0.5 x 0.4 m (0.2 m 2 ) areas. These areas were located 10 m East from the central post. Percentage cover by the various plant species, vegetation categories (moss, lichen, herb, shrub, tree, total), stone, and litter was estimated visually to the nearest 5%.

The Humus Index was measured at each sampling plot according to [START_REF] Ponge | Humus Index: an integrated tool for the assessment of forest floor and topsoil properties[END_REF] and [START_REF] Ponge | Collembolan communities as bioindicators of land-use intensification[END_REF]. Waterlogging in the top 10 cm was recorded and landscape heterogeneity measured by the diversity of land-use within each unit, using [START_REF] Shannon | A mathematical theory of communications[END_REF] formula.

Correspondence analysis (CA) was used to discern trends in the composition of plant communities [START_REF] Greenacre | Theory and Applications of Correspondence Analysis[END_REF]. The principal variables were species cover, and plant categories.

Additional variables were litter and stone cover, local biodiversity, units and land-use types (Table 1), humus forms, Humus Index, and waterlogging. Discrete variables were coded as 1 or 0. To give the same weight to all parameters all variables were transformed into X = (x-m)/s + 20, x being the original value, m the mean of a given variable, s its standard deviation. The addition of a constant factor of 20 allowed all values to be positive. Following transformation, factorial coordinates of variables were interpreted directly in terms of their contribution to the factorial axes. Variables were doubled to allow for the dual nature of measurement [START_REF] Greenacre | Theory and Applications of Correspondence Analysis[END_REF].

Results

Mixed or herbaceous dominated landscapes (units 3 to 6) had more plant species at the local scale than forested landscapes (Table 1). The curve for local plant biodiversity followed that of landuse diversity, and increased from unit 1 to unit 4 then decreased from unit 4 to unit 6 (Fig. 1). The Spearman rank correlation coefficient between local plant biodiversity and land-use variety was significant at the 0.01 level (r s = 0.89, d.f. = 4).

Strong differences were observed between land-uses when pooling samples (Table 1). At the local scale grassland had more plant species than arable crops, fallows being intermediate, fallows

and crops having more plant species than forests. These differences were of the same order of magnitude as those observed across landscape units, but standard errors of the means were smaller.

The difference between agricultural land and woodland, with a high and a low number of plant species respectively, was reduced when comparing units (Fig. 2). Species richness increased by a factor of 3 in woodland from unit 2 to 4, and decreased by a factor of 2 in agricultural land from unit 3 to 6.

Species richness showed a marked decrease in arable crops along the gradient of land-use intensification, and was positively influenced by landscape heterogeneity in both woodland and agricultural land.

In the CA of plant communities 1 , axis 1 was neglected because it separated only one sample with hydrophilic vegetation from all other samples. The projection of samples in the plane formed by axes 2 and 3 showed three branches, i.e. woodland, grassland and arable crops, fallows being positioned according to previous land-use for the most recent or as successional stages towards woodland.

The different units were positioned in the plane of axes 2 and 3 according to the dominance of land-use types within each of them (Fig. 3). Axis 2 showed that hay meadows and pastures clearly differed from crops by a denser and richer herb vegetation. Axis 3 depicted a gradient of increasing species richness and herb cover, decreasing litter and tree cover, and decreasing Humus Index according to land-use intensification. The significance of axis 3 associated the gradient of land-use intensification from unit 1 to 6 with increasing soil biological activity, as inversely measured by the Humus Index [START_REF] Ponge | Humus Index: an integrated tool for the assessment of forest floor and topsoil properties[END_REF], and increasing local plant biodiversity, as expressed by species richness per 0.2 m 2 . Axis 3 was significantly correlated with Humus Index (r s = -0.61, P < 10 -4

) and local plant biodiversity (r s = 0.57, P < 10 -4

).

Discussion

Grassland exhibited a higher level of plant biodiversity than arable crops (Fig. 2), an effect that can be explained by grazing or repeated mowing which increased nutrient cycling, created gaps, giving more chance to poorly competitive species to reach the canopy, and increased dispersion by animals or wind [START_REF] Watt | The effect of excluding rabbits from acidiphilous grassland in Breckland[END_REF][START_REF] Grime | Floristic diversity in a model system using experimental microcosms[END_REF][START_REF] Silverstone | Gaps in the canopy: the missing dimension in vegetation dynamics[END_REF]. Heavy grazing may, however, exert a detrimental influence on species richness [START_REF] Pettit | Grazing in remnant woodland vegetation: changes in species composition and life form groups[END_REF], whereas conventional farming may contribute to decrease plant biodiversity [START_REF] Edwards | Effects of pesticides on soil invertebrates and plants[END_REF]. The decreasing local plant biodiversity observed in arable crops from unit 4 to 6 (Fig. 2) can be explained by a decrease in the number of landscape components, and by the absence of rotation between meadows and crops in unit 6. On the other hand, germination of dormant seeds after cultivation [START_REF] Brenchley | The weed population of arable soil. I. Numerical estimation of viable seeds and observations on their natural dormancy[END_REF] can contribute to increase plant biodiversity when crops replaced meadows.

The land-use type did not explain all the variability observed in local plant biodiversity, since for a given land-use more heterogeneous landscapes exhibited more plant species per unit surface (Fig. 2). Strong differences in plant species composition occurred in the three habitats woodland, grassland and arable crops. Given the ability of most plant species to disperse, especially through birds and wind, the seed rain may be more diverse provided woodland, grassland and arable crops exchange diaspores at the scale of a unit. Dispersal ability (Hester et al., 1991b) and permeability of the landscape [START_REF] Honnay | Permeability of ancient forest edges for weedy plant species invasion[END_REF] are important components of local plant biodiversity. In the Morvan region spatial heterogeneity reflects also temporal heterogeneity, through rotation between arable crops and meadows, which increased the effect of immediate dispersal [START_REF] Rescia | Ascribing plant diversity to historical changes in landscape: a methodological approach[END_REF]. In unit 4, mostly represented by small thickets and woodland contiguous to agricultural land, local plant biodiversity was high and accounted for strong landscape heterogeneity [START_REF] Collinge | Ecological consequences of habitat fragmentation: implications for landscape architecture and planning[END_REF].

Landscape diversity may have contrasting effects on soil-dwelling animal communities [START_REF] Ponge | Collembolan communities as bioindicators of land-use intensification[END_REF]. Above-and belowground biodiversity may reinforce each other through positive feedbacks [START_REF] Perry | Bootstrapping in ecosystems[END_REF][START_REF] Hooper | Interactions Between Aboveground and Belowground Biodiversity in Terrestrial Ecosystems: Patterns, Mechanisms, and Feedbacks[END_REF][START_REF] Ponge | Humus forms in terrestrial ecosystems: a framework to biodiversity[END_REF]. The observed decrease in springtail species richness when landscape variety increases may be related in Morvan to more rapid changes in land-use in heterogeneous than in homogeneous landscapes [START_REF] Ponge | Collembolan communities as bioindicators of land-use intensification[END_REF]. Plants being much more rapidly dispersed than springtails (Collembola), at least at the scale of the study unit (Hester et al., 1991b), factors which can explain changes in local biodiversity of these two taxonomic groups are not necessarily similar. However, the high nutrient status of the ecosystem could enhance biodiversity of all categories of organisms [START_REF] Ponge | Humus forms in terrestrial ecosystems: a framework to biodiversity[END_REF], as was observed in unit 6. Soil fertility has been made responsive for losses of biodiversity of plants by by favouring fast-growing species [START_REF] Huston | A general hypothesis of species diversity[END_REF][START_REF] Tilman | Resource Competition and Community Structure[END_REF][START_REF] Loreau | Biodiversity and ecosystem functioning: a mechanistic model[END_REF], whereas [START_REF] Brandt | Effects of limestone dust accumulation on composition of a forest community[END_REF]Rhoades (1972), Miller et al. (1977), [START_REF] Roem | Soil acidity and nutrient supply ratio as possible factors determining changes in plant species diversity in grassland and heathland communities[END_REF] demonstrated the contrary. We suggest that more attention should be paid to the plant toxicity of N compounds [START_REF] Knight | The role of earthworms in managed permanent pastures in Devon, England[END_REF][START_REF] Thomas | Plant diversity in managed forests: understory responses to thinning and fertilization[END_REF] before claiming that nutrients are responsible for collapses in biodiversity through active competition.

Land-use type and diversity had some major effects on local plant biodiversity in the present study, although trends only were detected by CA [START_REF] Benzécri | Statistical analysis as a tool to make patterns emerge from data[END_REF][START_REF] Loreau | Biodiversity and ecosystem functioning: a mechanistic model[END_REF]. Land-use intensification increased plant species richness by providing herbs with more light and nutrients. The highest level of disturbance (arable crops) decreased biodiversity against perennial plants and favoured annuals tolerant to agricultural practices. Compared to homogeneous ones, heterogeneous landscapes tend to include more species in the seed rain and create edge effects [START_REF] Harris | Edge effects and conservation of biotic diversity[END_REF]. As a result, the intermediate disturbance level, (hay meadows and pastures) exhibited the highest level of local biodiversity, particularly in units with high land-use diversity.

The Intermediate Disturbance Hypothesis [START_REF] Connell | Diversity in tropical rain forests and coral reefs[END_REF] probably applies to vegetation [START_REF] Molino | Tree diversity in tropical rain forests: a validation of the intermediate disturbance hypothesis[END_REF], and to taxonomic groups with colonization rate exceeding the disturbance rate [START_REF] Sheil | Disturbing hypotheses in tropical forests[END_REF]. The reverse is true for soil microarthropods, which benefit from a more stable environment but are vulnerable to changes in land-use [START_REF] Ponge | Collembolan communities as bioindicators of land-use intensification[END_REF].

Figure captions

Fig. 1 .

 1 Fig. 1. Distribution of mean (± S.E.) plant species richness per 0.2 m 2 plot (local plant biodiversity) and land-use variety (Shannon Index) across six landscape units.

Fig. 2 .

 2 Fig. 2. Distribution of mean (± S.E.) plant species richness per 0.2 m 2 plot across six landscape units, in woodland, grassland and arable crops.

Fig. 3 .

 3 Fig. 3. C.A. projection of additional and some main variables along axes 2 and 3 (bold = higher; lower values).

Species list available upon request from the corresponding author

Acknowledgements

This study was part of the European Community program BioAssess EVK2-CT-1999-00041, coordinated by A. Watt (UK), who is greatly acknowledged for financial support and fruitful exchange of ideas between partners. Y. Capon and F. Choquet are also acknowledged for their skillful contribution to plant sampling. Language and clearity have been improved by anonymous referees, who are greatly acknowledged.


Théodore Papadopoulo

Sylvain Vallaghé

Maureen Clerc

The adjoint method for general EEG and MEG sensor-based lead field equations

Keywords: Electroencephalography, Magnetoencephalography, Forward problem, Adjoint method

The inverse source problem in electroencephalography (EEG) and magnetoencephalography (MEG) generally uses a lead field, which relates any source in the brain to its measurements at the sensors. For complex geometries, there is no analytical formula of the lead field. The common approach is to numerically compute the value of the lead field for a finite number of point sources (dipoles). There are several drawbacks: the model of the source space is fixed (a set of dipoles) and the computation can be expensive for as much as 10 000 dipoles. The common idea to bypass these problems is to compute the lead field from a sensor point of view. We show how the adjoint method can be used to derive general EEG and MEG sensorbased lead field equations. Within a simple framework, we provide a complete review of the explicit lead field equations, and we are able to extend these equations to non-point like sensors.

I. INTRODUCTION

Electroencephalography (EEG) and magnetoencephalography (MEG) inverse problems traditionally optimizes a source configuration by "comparing" actual measurements to quantities simulated using a forward model that predicts the electric potential or magnetic field. The EEG forward problem is governed by the following partial differential equation:

∇ • (σ ∇V ) = ∇ • J p in Ω σ ∇V • n = 0 on ∂ Ω ( 1
)
where V is the electric potential, σ is the conductivity, J p is the primary current vector representing brain electrical activity. Ω represents the head domain (an open bounded region of the ambient space R 3 ). Similarly, the magnetic field is given by the Biot-Savart law:

B(r i ) = B 0 (r i ) - µ 0 4π Ω σ ∇V × ∇ 1 R dr , (2)
with B 0 (r i ) = µ 0 4π Ω J p × ∇1 R dr and R = r ir . It is obvious that both V and B are linear functions of J p . Because of this property, the predictions of the forward problem at the sensor locations (the only information used by inverse problems) can be conveniently summarized by a leadfield, which is a linear operator L that maps source activations to simulated potentials of magnetic fields at sensor locations. Since it is linear, this operator can be written using the Riesz representation theorem:

L(J p ) =< ∂ L ∂ J p , J p > , (3)
where <, > denotes for the inner product 1 of functions in R 3 :

< u, v >= Ω u(r) • v(r)dr .
Traditionally, the forward problem equations are discretized and the leadfield is represented by the matrix ∂ L ∂ J p whose number of lines equals the number of sensors and whose number of columns equals the number of parameters of the dipoles considered as possible sources (around 10 000 in distributed source models). Each column of this matrix gives the values of the electric or magnetic measurements at the sensor locations for a given dipole parameter. A simple way to compute the leadfield matrix is thus to solve as many forward problems as there are source parameters, which can be quite time consuming with detailed realistic geometries for the head. One solution to bypass this problem relies on some matrix manipulations on the discretized problem so as to estimate rows of the leadfield matrix (instead of columns) [START_REF] Weinstein | Lead-field bases for electroencephalography source imaging[END_REF][START_REF] Wolters | Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem Inverse Problems[END_REF]. This is quite interesting as there are generally much less sensors than source positions and as the resulting computation is very similar to the original forward problem. Another solution is to work in the continuous domain, and to use the Helmholtz reciprocity principle to express the leadfields as the solutions of different forward problems [START_REF] Weinstein | Lead-field bases for electroencephalography source imaging[END_REF][START_REF] Nolte | The magnetic lead field theorem in the quasi-static approximation and its use for magnetoencephalography forward calculation in realistic volume conductors[END_REF][START_REF] Schimpf | Application of quasi-static magnetic reciprocity to finite element models of the MEG lead[END_REF].

Both these approaches can been explained using the adjoint method [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF], which is applied here in the continuous domain to both the EEG and MEG leadfields. As the approach is very general, it is quite easy to include within this framework new information such as sensor geometry.

II. THE ADJOINT METHOD

Let p represent a parameter which can be a real scalar or vector function. For any p, we consider the (hopefully unique) function v which satisfies the state equation:

Av = b(p) , (4)
where A is a linear operator of R 3 and b is a differentiable function of R 3 . To give some intuition, in our case the state function v is the electric potential V (r), the parameter p is the primary source current vector J p (r), the operator A is the generalized Laplacian ∇ • (σ ∇.) and b is the divergence operator.

Let us consider a "measurement function" g which is a linear functional of p and v. From the Riesz representation theorem, there exists m and n ∈ E such that g(p, v) =< m, v > + < n, p >. To compute the differential of g with respect to p (knowing that v also depends linearly on p), we first build a Lagrangian L by adding to the measurement equation the inner product of the state equation with a Lagrange multiplier w:

L (p, v, w) =< m, v > + < n, p > + < w, Av -b(p) > . Notice that g(p, v(p)) = L (p, v(p), w) for all w as Av(p) = b(p). Assuming that L is differentiable with respect to both v and p, the Lagrange multiplier w is chosen such that ∂ L ∂ v (p, v(p), w)δ v =< m, δ v > + < w, Aδ v >= 0 for all δ v .
Introducing the adjoint operator H * of a linear operator H, such that < x, Hy >=< H * x, y >, the previous equation can be rewritten as < m + A * w, δ v >= 0 for all δ v, which is equivalent to:

A * w = -m . ( 5
)
This equation is called the adjoint equation. With this choice for w, the differential of g can be written as:

δ g = ∂ L ∂ p (p, v(p), w)δ p =< n - ∂ b ∂ p * w, δ p > . ( 6
)
Since g is linear in p, g(p) =< ∂ g ∂ p , p >. g(p) being a measurement, ∂ g ∂ p is exactly the leadfield L and Eq. ( 6) exactly expresses this operator as a function of the adjoint variable w:

L = n - ∂ b ∂ p * w . (7)

III. APPLICATION TO M/EEG LEADFIELDS

This section briefly details the application of the previous result to M/EEG leadfields. For more detailed calculations, the reader is referred to [START_REF] Vallaghé | The adjoint method for general EEG and MEG sensor-based lead field equations Physics in[END_REF]. Using the previous result first requires the computation of the adjoint operators associated to our state equation (1) (i.e. A * and ∂ b ∂ p * ). In our case, v is the potential V (r), p is the primary source J p (r) and b is the divergence operator, which is linear so that ∂ ∇•J p ∂ J p J p = ∇ • J p . Using the divergence theorem:

< w, ∇•J p >= Ω w∇•J p dr = - Ω J p •∇wdr =< -∇w, J p > ,
because sources are restricted to be strictly inside the head, so the boundary integrals vanish. Thus:

∂ b ∂ p * w = ∂ ∇ • J p ∂ J p * w = -∇w . (8)
For A * , we start with < w, Av >:

< w, Av >= Ω w∇ • (σ ∇V )dr .
Assuming that w is continuous on Ω, but that σ ∇w • n can be discontinuous at tissue interfaces, the divergence theorem applied twice leads to (using the facts that σ ∇V • n = 0 on ∂ Ω and that V and σ ∇V • n are continuous functions over Ω):

< A * w, v > = Ω w∇ • (σ ∇V )dr = Ω V ∇ • (σ ∇w)dr -∑ k S k V [σ ∇w • n] S k ds , (9)
where S k are the tissue interfaces and [.] S k denotes the jump of a function on a given surface S k .

A. EEG leadfield

One line of the EEG leadfield is the potential of an electrode located at point r i with respect to some reference electrode located at r 0 (mean reference can also easily be handled with slightly more complicated formulae). Computing a leadfield for dipoles with fixed orientations can also be done easily following the same principle.

g EEG (p, v) = g EEG (J p ,V ) = V (r i ) -V (r 0 ) = ∂ Ω V (r)(δ (r i ) -δ (r 0 )) dr ,
where the notation δ (r) refers to the Dirac distribution at r ∈ ∂ Ω. The integral is restricted to ∂ Ω since all the EEG sensors are located on this surface. As g(J p ,V ) does not depend on J p , this leads to m EEG = δ (r i ) -δ (r 0 ) and n EEG = 0. Equating all the surfacic and volumic terms in the linear forms shows that [σ k ∇w • n] S k = 0 for all inner interfaces (as there are no other surfacic terms for those), hence only the last surface (∂ Ω) is to be kept in the sum in Eq. (9):

< A * w,V >= Ω V ∇ • (σ ∇w)dr - ∂ Ω V [σ ∇w • n] ∂ Ω ds .
From Eq. ( 5), w satisfies:

∇ • (σ ∇w) = 0 in Ω σ ∇w • n = δ r i -δ r 0 on ∂ Ω (10)
Combining Eq. ( 7) and Eq. ( 8), the EEG leadfield for punctual electrodes is:

L EEG = ∇w (11)

B. MEG leadfield

For MEG, the only difference is in the type of measurements: g MEG (p, v) = g MEG (J p ,V ) = B(r i ) • d i , with B given by the Biot-Savart law (2). This law has two terms, B 0 depending on J p and the other term depending on V . Applying the adjoint state approach requires the computation of the corresponding m MEG and n MEG .

n MEG corresponds to the B 0 term:

B 0 (r i ) • d i = µ 0 4π Ω J p × ∇ 1 R • d i dr = µ 0 4π Ω J p • ∇ 1 R × d i dr .
Thus

n MEG = µ 0 4π ∇ 1 R × d i . ( 12
)
m MEG corresponds to the Ohmic current term:

< m MEG , J p > = - µ 0 4π Ω σ ∇V • ∇ 1 R × d i dr = µ 0 4π Ω V ∇ • σ 1 R × d i dr- µ 0 4π Ω ∇ •V σ 1 R × d i dr = µ 0 4π Ω V ∇ • σ 1 R × d i dr- µ 0 4π ∑ k S k V [σ ] S k 1 R × d i • n ds .
Thus m MEG has both a volumic term and several surfacic terms. Equating all the surfacic and volumic terms in the linear forms yields the following equation for the adjoint:

∇ • (σ ∇w) = -µ 0 4π ∇ • σ ∇ 1 R × d i in Ω k , [σ ∇w • n] S k = -µ 0 4π [σ ] S k ∇ 1 R × d i • n on S k
Note that if the conductivity σ is constant and isotropic in a domain Ω k , then ∇ • σ ∇ 1 R × d i = 0 and w is harmonic in each domain Ω k . This result has been used previously for the numerical computation of the MEG leadfield [START_REF] Nolte | The magnetic lead field theorem in the quasi-static approximation and its use for magnetoencephalography forward calculation in realistic volume conductors[END_REF][START_REF] Schimpf | Application of quasi-static magnetic reciprocity to finite element models of the MEG lead[END_REF].

Combining Eq. ( 7), Eq.( 12) and Eq. ( 8), the MEG leadfield for a punctual magnetometer is:

L MEG = ∇w + µ 0 4π ∇ 1 R × d i .

IV. INCORPORATING SENSOR GEOMETRY

A. MEG: Squids geometry

A magnetometer measures the flux of the magnetic field through a small loop. For a magnetometer i, let M i be the surface enclosed by the loop, and d i the unitary vector normal to M i . The function g can thus be written as g(J p ) = M i d i • B(r ′ )dr ′ . Switching the integrals over M i and Ω and since the only dependence of B(r ′ ) on r ′ is in R, it is quite easy to show that following the lines of section III.B yields:

∇ • (σ ∇w) = -µ 0 4π M i ∇ • σ ∇ 1 R × d i dr ′ in Ω k , [σ ∇w • n] S k = -µ 0 4π [σ ] S k ∇ 1 R × d i • n on S k
As in the point wise case, the term

M i ∇ • σ ∇ 1 R × d i dr ′ vanishes if σ is constant
and isotropic in a domain Ω k . The resulting expression for the leadfield is:

L MEG = ∇w + µ 0 4π M i ∇ 1 R × d i dr ′ .
This formulation is the same as the one given in [START_REF] Wolters | Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem Inverse Problems[END_REF]. Gradiometers can be treated similarly. Generally, MEG manufacturers give a set of positions and weights for each sensor, and the linear combination of the magnetic field at these positions using these weights is meant to recreate the measurement. But this requires as many adjoint problems as there are positions. The main advantage of the proposed approach is to require a single adjoint problem per sensor, reducing notably the needed for the leadfield computation.

B. EEG: Electrode surfaces

Incorporating the surface of EEG electrodes in the model is more complicated as the electrode perturbs the electrical field due to the shunt effect: potential is constant at the surface of the electrodes due to the high conducting electrode (gel and metal). This effect has been modeled for EEG [START_REF] Jorma | Effects of electrode properties on EEG measurements and a related inverse problem[END_REF] and requires a modification of the boundary condition of the state equation [START_REF] Weinstein | Lead-field bases for electroencephalography source imaging[END_REF], which becomes:

   V + z k σ ∇V • n = v k on e k e k σ ∇V • n dr = 0 σ ∇V • n = 0 on S\ ∪ e k (13)
where e k is the k-th electrode, v k is the constant value of the potential on the k-th electrode and z k is the effective contact impedance which models the electrochemical effect at the skin-electrode interface. Following the adjoint approach yields the same leadfield equation as in the pointwise EEG case (Eq. ( 11)). However, the boundary condition of the adjoint state equation for w is transformed into:

             w + z k σ ∇w • n = W k on e k , e k σ ∇w • n dr =      1 k = i -1 k = 0 0 otherwise , σ ∇w • n = 0 on S\ ∪ e k , (14)
where W k is a constant.

V. RESULTS

To illustrate the approach, lead field for magnetometers are computed for a standard three nested spheres model with radii of 0.87, 0.92, 1, meant to represent brain, skull and scalp tissues with respective conductivities of 1, 0.02, 1. While it is known that in a spherical geometry the magnetic field outside the conductor does not depend on the conductivities, we intentionally put different conductivities in our model to test that the numerical solution is actually similar to the case of a homogeneous sphere. Computations were made using a finite element method presented in [START_REF] Théodore | Implicit Meshing for Finite Element Methods using Levelsets in Proceedings of MMBIA 07[END_REF] and results were compared with the ground truth given by the analytical formulation. The magnetometer is oriented in the x direction, positioned at a distance of 0.03 of the outermost sphere, and has a radius of 0.015. The differential equation of section IV A was solved, using our finite element method. The integral term over the magnetometer loop was computed after applying Stockes theorem with a Gauss-Kronod method using 61 points. Figure 1 shows the absolute error is similar and small for three dipole orientations. More results can be found in [START_REF] Vallaghé | The adjoint method for general EEG and MEG sensor-based lead field equations Physics in[END_REF].

VI. CONCLUSION

The adjoint method is a powerful general tool to derive efficient leadfield computation for various sensors types.

Fig. 1 :

 1 Fig. 1: Absolute error between numerical and analytical solution with respect to dipole depth.

For scalar functions the inner product inside the integral becomes a simple product.

ACKNOWLEDGEMENTS

This work is partly supported by contract EADS 2118.


Emmanuel Olivi

Maureen Clerc

Théodore Papadopoulo

Domain decomposition for coupling finite and boundary element methods in EEG

Keywords: EEG, BEM, FEM, Coupling, Domain Decomposition

The forward problem in electroencephalography aims to simulate on the scalp the potential V of an electromagnetic field generated by a simulated source. It must fit precisely with the electromagnetic propagation in the patient head. Yet, the skull anisotropy happens to be highly anisotropic, and must then be modeled. Although boundary element methods cannot deal with anisotropy like finite element methods, the symmetric BEM offers a higher accuracy than FEM wherever the conductivity can be considered as constant (i.e. for the brain and the scalp). A domain decomposition (DD) framework allows to split the global system into several ones with smaller computational domains. Then, one method (BEM or FEM) can be used per volume. This work presents such a coupling formulation of a 3-DD method solving iteratively a BEM for the brain, a FEM for the skull layer, and finally a BEM for the scalp.

I. INTRODUCTION

The problem of finding the sources in the brain that are responsible for a measured EEG signal is an inverse problem of localization whose resolution requires many resolutions of the forward problem. The latter computes the potential V on the scalp of an electromagnetic field generated by a simulated source (in our case a dipole J p ). This propagation through the head is governed by the Maxwell equations in their quasi-static approximation. The forward problem resolution is achieved solving for V the Poisson equation:

   ∇ • (Σ∇V ) = ∇ • J p in Ω (Σ ∇V ) • n = 0 on Γ (1)
,with Σ the head tissues' conductivity and Γ the Scalp/Air interface (outflow current through the neck is neglected). The forward model must fit precisely with the electromagnetic propagation in the patient head. Although the white matter makes the brain anisotropic, the skull anisotropy has a strong influence on the results, and must therefore be taken into account. Indeed, the skull conducts electricity ten to eighty times stronger in tangential directions than in the ra-dial one. Finite Element Methods (FEM:s) can deal with this anisotropy, but offer a poor accuracy when compared to the symmetric Boundary Element Method (sBEM) in homogenous media. Therefore the coupling of both should yield better accuracy ; using the BEM wherever the conductivity can be considered as constant (i.e. for the scalp and the brain when we neglect white matter anisotropy), and the FEM for the skull.

II. THE SYMMETRIC BEM, THE IMPLICIT FEM AND COUPLING FORMULATIONS

A. The symmetric BEM, and the implicit FEM Wherever conductivity is considered as constant in a set of nested regions (see Fig. 1) then the first equation in Eq. 1 becomes σ ∆V = ∇ • J p and can be solved using integral formulations. Instead of the classical formulation that involves double-layer potentials, we use a formulation that involves both double-layer and single-layer potentials, to yield a symmetric system of equations. The symmetric BEM [START_REF] Kybic | A common formalism for the integral formulations of the forward EEG problem[END_REF] involves both the potential and the normal current as unknowns on triangulated surfaces describing the different regions. FEM:s usually require to mesh the entire volume to spread the information from the source through the computational domain. The geometry used to solve the forward problem in EEG comes from anatomical MRI:s. These 3D-images are then segmented to obtain a levelset description of the regions to model. Generation of the (most often tetrahedral) meshes can be achieved but is usually expensive computationally and time consuming. The implicit FEM [START_REF] Papadopoulo | Implicit meshing for finite element methods using levelsets[END_REF] bypasses this mesh generation step, going directly from the levelsets of the interfaces separating the various domains to the matrix associated to the FEM. It works on Cartesian grids with the potential V at the node locations as unknown, described with Q1 finite elements. The sBEM developped in the open-source code OpenMEEG, turns out to be much more accurate than other BEM:s, and also than the implicit FEM (iFEM).

B. A BEM-FEM coupling

Domain decomposition is a widely used numerical technique that allows to split the computational domain into regions of constant conductivities on the one hand and the skull on the other hand. We will then be able to use the appropriate method (BEM or FEM) per volume. Boundary conditions ensure the communication between the sub-problems. Several iterations are needed to solve the global system, and a relaxation at the interfaces is compulsory to ensure convergence. We first propose a Neumann-Dirichlet approach which iteratively solves a BEM for the brain region, and a FEM for the skull and the scalp region. At iteration k, we first solve a BEM with Neumann boundary condition, to obtain V at the brain/skull interface. Secondly, a mixed Dirichlet -homogeneous Neumann problem is solved with the FEM using the newly computed V at the brain/skull interface as boundary condition. This gives V on the scalp:

               σ 1 ∆V k 1 = ∇ • J p in Ω 1 , σ 1 ∂ n V k 1 = λ k on Γ 1 , ∇ • (Σ∇V k 2,3 ) = 0 in Ω 2 ∪ Ω 3 , V k 2,3 = V k 1 on Γ 1 , σ 3 ∂ n V k 2,3 = 0 on Γ 3 , (2)
Finally, λ k+1 is updated by combining the normal current computed by the FEM on the brain/skull interface and the previous λ k with a relaxation parameter ω:

λ k+1 = (1 -ω) λ k + ω(Σ 2 ∇V k 2 ) • n on Γ 1 , ω > 0
We also present a Neumann-Dirichlet-Neumann coupling procedure, restricting the FEM to the skull only (see Fig. 2), and completing the computations with a BEM for the scalp:

                               σ 1 ∆V k 1 = ∇ • J p in Ω 1 , σ 1 ∂ n V k 1 = λ k 1 on Γ 1 , ∇ • (Σ 2 ∇V k 2 ) = 0 in Ω 2 , V k 2 = V k 1 on Γ 1 , V k 2 = V k-1 3 on Γ 2 , σ 3 ∆V k 3 = 0 in Ω 3 , σ 3 ∂ n V k 3 = λ k 2 on Γ 2 , σ 3 ∂ n V k 3 = 0 on Γ 3 , (3)
With relaxations made on the Neumann's parameters λ 1 and λ 2 :

λ k+1 1 = (1 -ω 1 ) λ k 1 + ω 1 (Σ 2 ∇V k 2 ) • n on Γ 1 , ω 1 > 0 λ k+1 2 = (1 -ω 2 ) λ k 2 + ω 2 (Σ 2 ∇V k 2 ) • n on Γ 2 , ω 2 > 0 (4)
Initial condition for this coupling scheme is: V 0 3 = 0, with the relaxation parameters: ω 1 = ω 2 = 1 for k = 1. For k > 1, relaxation parameters are set manually to ω 1 = ω 2 = 0.7 ; a higher relaxation parameter would make the scheme diverge. Special care must be observed for the Neumann problem, indeed no accumulation of current should appear, therefore numerically we ensure that:

Ω 1 ∇ • J p + Γ 1 λ k 1 = 0 , ∀k Γ 2 λ k 2 = 0 , ∀k (5)
Fig. 3: RDM for the isotropic case, for pure and coupled methods.

III. NUMERICAL RESULTS.

In order to validate results of the forward problem, analytical solutions are computed on a three-layer concentric sphere model [START_REF] Zhang | A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres[END_REF], for both isotropic and anisotropic skull layer. Radii of the spheres and conductivities of the different layers are respectively {0.87, 0.92, 1.0} and {1.0, 0.0125, 1.0} for the isotropic case. Considering the anisotropic one, we have set up the conductivity in the tangential direction to ten times the normal one. The sBEM is computed on 642 point mesh per surface, and the iFEM considers a Cartesian grid of 90 points in each direction. The BEM-FEM coupling as well as the BEM-FEM-BEM coupling use the same previous grid sizes. Computations have been done for 5 dipoles oriented in Cartesian coordinates: (1, 1, 0), and locations along the Z-axis: {0.465, 0.615, 0.765, 0.8075, 0.8415}. Accuracies of the numerical solutions are given by the Relative Difference Measure (RDM) of the potential on the scalp:

RDM = V analytic V analytic - V computed V computed (6)
Fig. 3 shows the RDM for both pure and coupled methods in the isotropic case. The iFEM method in blue (dashed) is the less precise due to the poor grid definition. We see that the sBEM is the most accurate with a RDM always below the 2%. The BEM-FEM coupling in red (full with diamonds) shows a slight improvement compared to the accuracy of the FEM, whereas the BEM-FEM-BEM coupling really inherits its accuracy from the BEM. Coupling process have been stopped at 40 iterations. At this point the relative residuals computed as:

(Σ 2 ∇V k+1 2 )•n| Γ i -λ k i λ k i
, for i = 1, 2 were all below 6.10 -5 . These results do not show the full benefit of the coupling, because the coupling takes its sense for anisotropic conductivities. As BEM cannot handle anisotropy, they must model the skull as isotropic and therefore their high accuracy is partly wasted because of the coarse model approximation. On Fig. 4, we have plotted the result of the BEM with an isotropic skull whose conductivity equals 0.0125 compared with the anisotropic analytical solution with conductivities (radial,tangential) = (0.0125, 0.125). The iFEM as well as the BEM-FEM-BEM coupling have been computed with an anisotropic skull. This time, results have been compared on 15 dipoles: the five previous locations for three different orientations: (1, 0, 1) in blue, (1, 1, 0) in green and (0, 0, 1) in red color. One can see that the coupling (full line) has a better precision than the BEM (dashed line with diamonds) for all dipoles expected the two located at (0, 0, 0.8415) with orientations (1, 0, 1) and (0, 0, 1). As it comes closer to the skull, the singularity of the dipole makes the approximation due to the iFEM's grid worse. Using a higher resolution FEM's grid, this error should go diminishing. Results on the blue curve for the coupling are not gratifying; there are some convergence problems with our BEM-FEM-BEM algorithm, for some dipole orientations, that we have not yet solved.

IV. CONCLUSION

The domain decomposition framework really allows to take advantage of both the versatility of the FEM in modeling conductivity and the accuracy of the BEM. It gives hope for providing accurate solutions of the forward problem without resorting to diminishing the mesh size in finite element methods. One could try to find a relaxation parameter changing through iterations to faster convergence. Parallelism is often used in the domain decomposition community, but it was not pertinent in our case since the sBEM has a direct solver, and only requires a matrix vector product to solve a new Neumann problem for each iteration ; most of the time is spent in the iterative solver of the FEM.

Fig. 1 :

 1 Fig. 1: A three-layer head model with nested regions.

Fig. 2 :

 2 Fig. 2: The three-layer head model with the FEM's grid for the skull.

Fig. 4 :

 4 Fig. 4: RDM for the anisotropic case, for pure and coupled methods.


M Porfiri

F Dell'isola

Multimodal beam vibration damping exploiting PZT transducers and passive distributed circuits

Structures are subject to extensive vibrations that can reduce their life and contribute to failure. Piezoelectric transducers in conjunction with appropriate electric networks, can be used as a mechanical energy dissipation device. Most of the previous work has focused on active control in which electric power is supplied to the actuators that exert actions on the host structure to suppress its vibrations. Alternatively, undesired mechanical energy of a host structure could be converted into electrical energy that can be dissipated through a set of resistor. This does not require an external power unit and is a more economical means of controlling vibrations of a structure, but an effective transduction of mechanical energy into electric energy has to be guaranteed. Therefore, in the present work a completely passive distributed electric controller is found to achieve simultaneous beam multimodal vibration suppression. The optimal electric network (interconnecting the piezoelectric transducers) may be synthesized as a finite differences approximation of the derived distributed (infinite dimensional) optimal controller.

INTRODUCTION

In [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF] the possibility of damping beam vibrations by means of a piezoelectric transducer positioned on the structure and shunted with an (LR) impedance is investigated. The passive impedance combined with the inherent capacitance of the piezoelectric transducer creates a damped electromechanical beating. The method proposed in [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF] allows for an efficient single-mode control of structural vibrations, whenever the resonant circuit is tuned to the mechanical mode to be suppressed; nevertheless the efficiency of the electromechanical coupling strongly depends on the position of the transducer over the host structure. Moreover, the technical feasibility of the passive piezoelectric controller proposed in [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF] is limited, since too large inductors are required to produce low-frequency electrical resonance with the small inherent capacitance of the piezoelectric transducer. A generalization of the "piezoelectric shunt damping" technique proposed in [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF] to the multimodal control of vibrations is presented in [START_REF] Hollkamp | Multimodal Passive Vibration Suppression with Piezoelectric Materials and Resonant Shunts[END_REF] and [START_REF] Fleming | Optimization and Implementation of Multimode Piezoelectric Shunt Damping Systems[END_REF], where complicated impedances are placed across the terminals of the piezoelectric transducer. Such a control methodology seems to present severe inconveniences, indeed the used inductances are still very high and the damping efficiency, even for few modes, is drastically reduced.

In order to cure these drawbacks, without forsaking the advantages featured by passive control, in [START_REF] Alessandroni | A revival of electric analogs for vibrating mechanical systems aimed to their efficient control by PZT actuators[END_REF] it is proposed to position an array of piezoelectric elements on the beam and interconnect their electric terminals via a circuit "analog" to the beam (i.e. governed by the same evolution equations) completed with proper resistors. The electric analogs synthesized in [START_REF] Alessandroni | A revival of electric analogs for vibrating mechanical systems aimed to their efficient control by PZT actuators[END_REF] show several inconveniences, e.g. the presences of negative inductors and multiport transformers. In order to overcome these technical problems, in [START_REF] Andreaus | Piezoelectric passive distributed controllers for beam flexural vibrations[END_REF] an electric analog exploiting only inductors and single port transformers is derived.

In both [START_REF] Alessandroni | A revival of electric analogs for vibrating mechanical systems aimed to their efficient control by PZT actuators[END_REF] and [START_REF] Andreaus | Piezoelectric passive distributed controllers for beam flexural vibrations[END_REF], it is a priori assumed that, in order to achieve effective damping of several modes, the electric interconnection should be the beam electric analog. The aim of this paper is to substantiate this assumption by proving that the electric analog represents the optimal passive controller. Indeed, we will consider a rather wide class of passive electric controllers and we will formulate and solve the problem of finding the optimal controller evolution equations in order to maximize the time decay rates of every mechanical mode.

Section 2 is devoted to the derivation of the equations of motion of a Piezo-ElectroMechanical (PEM for brevity) beam, i.e. a beam coupled to a passive electric network via distributed piezoelectric transduction. The beam is modelled as a continuous system, while the electric circuit is regarded as a lumped system. From the mechanical point of view, the piezoelectric elements introduce additional stiffness and mass, and exert electrically driven couples at their ends; while, from the electric point of view they are treated as capacitors in parallel connection with current sources driven by average beam strain rates. Therefore, the resulting system state is ruled by the beam deflection field and the discrete set of circuit independent flux linkages. In previous works, e.g. [START_REF] Maurini | Comparison of piezoelectronics networks acting as distributed vibration absorbers[END_REF], no investigations were made towards unstructured passive electric networks, i.e. circuits with not prescribed topology (in [START_REF] Maurini | Comparison of piezoelectronics networks acting as distributed vibration absorbers[END_REF] several electric interconnections constituted solely by two terminal impedances are analyzed). The resulting governing equations cannot be easily handled without exploiting numerical techniques.

In order to achieve a more maniable model of the overall system, in Section 3 a homogenized model of a PEM beam is deduced from the previously introduced refined model. The state of the homogenized system is determined by the beam deflection field and the electric flux linkage field, and its motion is governed by two gyroscopically coupled partial differential equations. In order to avoid spillover phenomena in between different mechanical modes, the electric operators appearing in the governing equations are assumed to share the mechanical eigenfunctions; namely, their eigenfunctions are assumed to be equal to the mechanical ones and only their eigenvalues are left as unknowns.

Considering the simply supported case, the modal analysis of the PEM beam can be easily performed as it is done in Section 4. In particular, a simple expression for the eigenvalues of the PEM beam is derived.

In Section 5 the optimization problem is formulated and completely solved, by the use of the aforementioned modal analysis. Indeed, the hypothesis made on the eigenfunctions of the electric operators, i.e. to be coincident with the mechanical modal shapes, together with the modal expansion worked in Section 4 lead to the simpler problem of optimizing a sequence of decoupled finite dimensional controllers. The solution of these standard problems enables the determination and spectral representation of the optimal electric operators. In particular it is proven that the optimal interconnection is the "analog" one, completed with resistors represented by a second order spatial derivative.

REFINED MODEL FOR PEM BEAMS

We consider a host beam of length l, width w and thickness h on which an array of N uniformly distributed piezoelectric transducers is positioned as shown in Figure 1 (covering both the beam faces). The length of the transducers is assumed to be equal to l p , while the width is assumed to be equal to that of the beam; d denotes the distance between the adjacent patches.

The piezotransducers are placed along the beam axis and polarized in the transverse direction, in the so called bender configuration (see e.g. [START_REF] Andreaus | Piezoelectric passive distributed controllers for beam flexural vibrations[END_REF]). These transducers will be interconnected by an electrical network which will be synthesized in order to accomplish given optimality conditions on the mechanical vibration decay. Introducing a set of nodes {xi} i=1,...,N representing the geometrical centers of the transducers, defined by:

l w h N l p , d Electrical circuit x l p /2 + dl p d x i i ,...N x i l p d ii, ...N, x x N l M x M x, t E b I b k mm N i l p x -x i u x, t k me N i l p x -x i ψ i t ,
x i = (l p + d)i i = 1, . . . N, (1)
and also including the boundary nodes x 0 = 0 and x N +1 = l, the contact couple M at the section labelled by the abscissa x, over the beam span, can be expressed as the sum of mechanical and piezoelectric contributions as follows:

M (x, t) = E b I b + k mm N i=1 (RECT lp (x -x i )) u ′′ (x, t) + k me N i=1 RECT lp (x -x i ) ψ(l) , (2)
where RECT lp (x -x i ) is a unit step function centered at x i of width l p , ψ i is the flux linkage of the i-th transducer measured with respect to a common reference ground for every transducer, u is the beam deflection field, E b is the Young modulus of the material of the beam, I b is the beam section moment of inertia, k mm is the transducers bending stiffness, k me is the piezoelectric coupling coefficient, t denotes the time variable and superposed dot and prime respectively mean time and space derivative.

The distributed inertia can be accounted for by introducing the following constitutive equation for the applied load (external mechanical forcing is excluded):

b T (x, t) = -ρ b + N i=1 (RECT lp (x -x i ))ρ p ü(x, t). (3)
where ρ b and ρ p are respectively the mass density per unit length of the beam and the transducer layers.

The balance equations for the considered electrically excited vibrating beam yield:

M (x, t) ′′ -b T (x, t) = 0 (4)
when M (x, t) and b T (x, t) are given by Equations ( 2) and ( 3) respectively.

From a purely electrical point of view the i-th piezoelectric bender transducer can be described as a capacitor kee in parallel connection with a "mechanically driven" current source, which injects into the electrical circuit the current J i driven by its mechanical time rate of deformation:

J i = k me ( u′ (x i + l p /2, t) -u′ (x i -l p /2, t)).
The electrical system interconnecting the electrical terminals of the bender transducers is assumed to be a linear, time invariant, reciprocal, passive N terminals network (see e.g. [START_REF] Newcomb | Linear Multiport Synthesis[END_REF]) constituted only of resistors, inductors and transformers. Thus by Kirchhoff balance of currents at node i, the evolution equations for the variables ψ j are given by the following second order system of ODEs:

k ee ψi + N j=1 K ij ψ j + N j=1 D ij ψj -k me ( u′ (x i + l p /2, t) -u′ (x i -l p /2, t)) = 0, i = 1, . . . N. (5)
[K ij ] being the electric stiffness and [D ij ] the electric dissipation.

A CLASS OF HOMOGENIZED MODELS FOR PEM BEAMS

When it can be assumed that the number of piezoelectric transducers is sufficiently large, a continuous flux linkage field ψ(x, t) can be introduced and the governing equation expressed in terms of the deflection u is seen to be (for more details on the needed homogenization techniques see e.g. [START_REF] Kader | Distributed control based on distributed electronic circuits: application to vibration control[END_REF]):

k M u IV + ρü + c f K me ψ′′ = 0, k M = E b I b + cf k mm ρ = ρb + cf ρ p
where the covering factor is defined as:

c f = 1 1 + d/l p .
Thus the dimensionless form of the mechanical evolution equation becomes:

ü + α 4 u IV + β 2 em ψ′′ = 0,      α 4 = K M l A ρω 2 0 β 2 em = c f k me ψ 0 l 2 ρω 0 u 0 , (6)
where a characteristic frequency ω 0 , a characteristic length equal to the beam length l, a characteristic deflection u 0 and a characteristic flux linkage ψ 0 have been introduced. As this cannot cause misunderstanding, we have adopted the same letters to denote both dimensional and dimensionless differential operators and kinematical descriptors.

Furthermore, when the number of transducers can be assumed to be sufficiently large, the electrical system can be described by a sole PDE which, expressed in terms of the dimensionless flux-linkage ψ, reads:

ψ + K[ψ] + D ψ -β 2 me u′′ = 0, β 2 me = k me u 0 l 2 kee lp ω 0 ψ 0 ,
where K and D are linear operators. In the following we will restrict our investigation to simply supported beams PEM beams, i.e. we will assume simply supporting boundary conditions for both the mechanical and electrical systems. Furthermore, in order to avoid undesired spillover phenomena in between different vibration modes, we will look for a dissipative electric network which does not modify the electrical and mechanical modal shapes. Therefore, in the present paper the following restrictive assumption on K and D are made:

K [•] = k (•) IV , D [•] = d (1) IV (7)
k and d being analytic functions which take positive values on the real positive axis a , see e.g. [START_REF] Pedersen | Functional Analysis in Applied Mathematics and Engineering[END_REF].

a The positiveness of k(λ) and d(λ) , when λ is positive assures the passivity of the electric network.

This hypothesis on the structure of the electrical stiffness and damping, respectively K and D, seems to be general enough for our design aims and leads to a simple solution of the synthesis problem of determining an optimal lumped electrical circuit to be piezoelectrically coupled to the given flexible structure (see e.g. [START_REF] Alessandroni | A revival of electric analogs for vibrating mechanical systems aimed to their efficient control by PZT actuators[END_REF] and [START_REF] Andreaus | Piezoelectric passive distributed controllers for beam flexural vibrations[END_REF]). Thus the dimensionless form of the evolution electromechanical equations is:

ü + α 4 u IV + β 2 ψ′′ = 0 ψ + K[ψ] + D ψ -β 2 u′′ = 0 ;      α 4 = K M l 4 ω 2 0 β 2 = k me l 2 ω 0 c f l p k ee ρ (8)
where, in order to preserve the form of a gyroscopic coupling,the characteristic flux linkage and displacement are chosen to satisfy the following relation:

c f K ee l p ρ = u 0 ψ 0 .
The eigenfunctions of the both the mechanical and electrical operators are easily found to be:

v m (x) = √ 2 sin(mπx), m = 1, 2, . . . ; (9)
The eigenvalues of the elastica operator are

(mπ) 4 , m = 1, 2, . . . ;
while those of the electric stiffness and dissipation operators are respectively k(m) and d(m).

MODAL ANALYSIS OF HOMOGENIZED PEM BEAMS

In this subsection a modal analysis for a simply-supported PEM beam, based upon the Galerkin method, is performed.

The set of basis functions chosen to represent the deflection u(x, t) and the flux linkage ρ(x, t) of the PEM beam is constituted by the eigenfunctions v m (x) specified in (9) of the elastica differential operator (•) IV determining the modal forms for both the uncoupled electrical and mechanical motions [START_REF] Pedersen | Functional Analysis in Applied Mathematics and Engineering[END_REF]. Consequently it is useful to consider a Fourier expansion for the solution of (8):

u(x, t) ψ(x, t) = ∞ m=1 p m (t) q m (t) v m (x),
in terms of the mechanical and electrical Fourier coefficients p m (t) and q m (t) of the expansion; achieving with simple algebraic manipulations the following system of ODE s governing their time evolution, in which unitary dimensionless modal masses are made to appear:

1 0 0 1 pm (t) qm (t) + ω 2 p (m) 0 0 ω 2 q (m) p m (t) q m (t) + 0 -c(m) c(m) 2ζ q (m)ω q (m) ṗm (t) qm (t) = 0 0 (10)            ω 2 p (m) = (αmπ) 4 ω 2 q (m) = k(m) ζ q (m) = 1 2 d(m) 2 √ k(m) c(m) = (βmπ) 4 (11)
Let us remark that ω p (m) and ω q (m) respectively the natural angular frequency of the mechanical and electrical m-th modes, while c(m) denotes the m-th gyroscopic coupling coefficient and ζ q (m) the damping factor of the m-th electrical mode. Let us explicitly remark that the chosen boundary conditions and basis functions lead to a set of uncoupled evolution problems for each electromechanical pair of Fourier coefficients, hence it is possible to easily derive an analytical solution for each m-th pair. The poles of the system in (10) can be determined as the four zeros s m 1 , . . . s m 4 of the following polynomial:

P m (s) = s 4 + 2ζ q (m)w q (m)s 3 + (w 2 q (m) + c(m) + w 2 q (m))s 2 + + 2ζ q (m)ω q (m)ω 2 p (m)s + ω 2 p (m)ω 2 q (m).

DESIGN OF THE PASSIVE OPTIMAL CONTROLLER

In this section we will design the optimal passive distributed controller following the pole placement technique applied to the evolution of each pair of electrical and mechanical modes. In particular, we will determine an optimal expression for the stiffness and dissipative functions k(m) and d(m) appearing in the time evolution of the m-th electrical mode; consequently, we will establish the forms of the operators K and D, exploiting their spectral representation.

The chosen optimality condition requires the determination of the values of ω q (m) and ζ q (m) in order to maximize the exponential time decay rate τ m of the solution of the fourth order system (10), governing the time evolution of the m-th pair of electrical and mechanical modes, defined as:

τ m = min i=1,...4 {|Re[s m i ]|}.
In order to maximize the exponential time decay τ m , the optimal values of ω q (m) and ζ q (m) by requiring the four roots of P m (s) to be coincident (pole placement technique, see e.g. [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF]). Hence we enforce the polynomial P m (s) to be factorized as:

P m (s) ≡ (s 2 + 2σ(m)s + (σ(m) 2 + ω(m) 2 )) 2 .
Equating the coefficients of the above mentioned polynomials the following the following matching conditions are imposed on the electrical controller:

ω q (m) = ω p (m) ζ q (m) = √ c(m) ωp(m) . (12)
Condition (12) 1 establishes that in order to maximize the time rate decay of the m-th pair of electrical and mechanical modes, the electrical m-th mode has to be resonant at the mechanical resonance frequencyq ω q (m) . Thus, by use of (11), the following conditions on the values of the functions k and d in correspondence with the mode number m under control are stated:

k(m) = (αmπ) 4 d(m) = 2(βmπ) 2 . ( 13
)
The values of the real and imaginary part of the coincident roots are found to be:

     σ(m) = c(m) 2 , ω(m) = ω p (m) 2 -c(m) 4 (14)
hence, the damping ratio (defined as the sine of the phase of the coincident roots measured from the imaginary axis) becomes:

ζ(m) := σ(m) ω(m) 2 + σ(m) 2 = c(m) 2ω p (m)
.

Introducing the relations (11) into (14), we get the following expressions for the real and imaginary parts of the placed roots:

     σ(m) = (βmπ) 2 2 ω(m) = (mπ) 2 α 4 -β 4 4
Taking into account the definitions (8) of the parameters α and β , we get the following relation in terms of the properties of the piezo-electromechanical beam

ζ(m) = c f k me 2 √ γK M (15)
Conditions (13) establish that the optimal distributed passive circuit to be mechanically fed by the piezotransducers has to be governed by the following PDE:

ψ + K[ψ] + D[ ψ] = 0, (16)
where the stiffness and damping operators are defined by:

K[ψ] = α 4 ψ IV D[ ψ] = -2β 2 ψ′′ (17)
Remark that the self-adjoint operators K and D are both differential operators and satisfy the relationship K ∝ D 2 .

(18)

The found electrical circuit can be thought as the analog circuit of a simply supported vibrating beam, with an internal dissipation proportional to the curvature velocity u′′ . The synthesis of a lumped electrical network governed by a finite difference approximation of ( 16) with (17) has been addressed in [START_REF] Andreaus | Piezoelectric passive distributed controllers for beam flexural vibrations[END_REF] and [START_REF] Alessandroni | A revival of electric analogs for vibrating mechanical systems aimed to their efficient control by PZT actuators[END_REF].

CONCLUSIONS

In the present paper the concept of piezoelectric shunting by means of a single transducer suitably positioned on a structural member is generalized. Indeed we consider a beam hosting an array of piezoelectric transducers and look for a multiterminal electric network shunting them in an optimal way. In the first two sections a maniable model for a PEM beam is found, while in the subsequent sections its electric components are optimized to get the most efficient damping of mechanical vibrations. We have proved that, in principles, a multishape, multifrequency (i.e. a complete multimodal) vibration suppression is possible by means of a completely passive electric controller. This is obtained by synthesizing a distributed electric circuit governed by the same PDE governing the beam flexural vibrations, and completed with a resistive network governed by the second order spatial derivative. Even if we have considered only a simply supported PEM beam, the same reasonings can be successfully applied to differently constrained PEM beams. Nevertheless, in these cases the optimal electric operators may not be purely differential as in (17).

Figure 1 .

 1 Figure 1. Geometry of the problem.

Acknowledgements

The partial support of the Engineering Science and Mechanics Department of Virginia Polytechnic Institute and State University is gratefully acknowledged by the authors. This presented research has been also partially supported by MIUR, Ministero per l'Innovazione, l'Università e la Ricerca Fondi Ricerca PRIN "Sintesi di circuiti piezoelettrici e tecniche di disaccoppiamento per il controllo di vibrazioni meccaniche" (protocollo 2001097882 003) and by the Università di Roma "La Sapienza".


A Vaçulik

C Kounda-Kiki

C Sarthou

J F Ponge

Soil invertebrate activity in biological crusts on tropical inselbergs

Granite inselbergs protrude from forest and savanna in the tropics. They are exposed to harsh climates (alternation of heavy rain and severe drought) and provide little nutrient for plants. Soil animals and humus components were investigated in cyanobacterial crusts close to patches of epilithic vegetation on the surface of the Nouragues inselberg (French Guiana). Three biological crust samples, corresponding to bromeliacean carpets of increasing size (supposed of increasing age), were sampled for faunal and micromorphological studies. Arthropods (mainly mites and insects) were abundant and highly diversified, the more so after enchytraeid worms ate and tranformed the cyanobacterial mass. Below the superficial cyanobacterial crust, humus was made of a loose assemblage of enchytraeid faeces where these animals were present, or of a compact assemblage of cyanobacteria and amorphous organic matter where mites were the dominant animal group. Roots abounded in the humified part of the crust. We conclude that soil invertebrates, in particular enchytraeid worms, are important for the accumulation of organic matter on granite outcrops, and so therefore for the encroachment of plant succession.

Introduction

Inselbergs protrude from rain forests and savannas in the tropics. They support an epilithic vegetation which differs abruptly from that on the surrounding land. Plant successional trajectories have been described in several instances, pointing on difficulties for the establishment of woody vegetation. Vegetation needs to cope with strong erosion, caused by heavy storms and lack of protection, alternation of severe drought and heavy rain and abrupt changes in temperature. Filamentous cyanobacteria with heterocysts are particularly well-adapted to these harsh environments, since they live autotrophically for carbon and nitrogen in the absence of soil, can recover after desiccation and take opportunity of short wet periods for growth and reproduction [START_REF] Belnap | Factors influencing nitrogen fixation and nitrogen release in biological soil crusts[END_REF]. Cyanobacteria are most responsible for mineral weathering through the excretion of organic acids and repiration. Despite the common occurrence of cyanobacterial biofilms on rock surfaces, differences in the development of higher vegetation have been recorded according to countries and even in the same country, as a result of weak exchange of diaspores between isolated outcrops [START_REF] Sarthou | Shrub vegetation on tropical inselbergs (French Guiana)[END_REF].

In French Guiana, the Bromeliaceae Pitcairnia geyskesii L.B. Smith is common on well-drained slopes, while forbs and herbs abound around temporary rain-fed ponds.

This species is the dominant component of low vegetation, locally called 'savane-roche' [START_REF] Sarthou | Mécanismes de colonisation par la végétation d'un inselberg granitique en Guyane Française[END_REF]. At this stage biological crusts, a few centimetres thick, occur at the contact between bromeliacean carpets and cyanobacterial biofilms (Figure 1). These crusts constitute the first step in the development of humus, and they are niches for the establishment by seed of several plant species. In later stages of soil development, organic matter may accumulate to a depth of several decimetres, together with the appearance of a higher, woody vegetation which sheds litter and protects the ground from rapid desiccation.

We may wonder to what extent animals living in cyanobacterial crusts of tropical inselbergs differ from those in other organic soils and biological crusts and whether they affect the genesis of the soil. Do they contribute to the accumulation of a stable structure in these frequently disturbed environments, or are they only subordinate members of a microbial-dominant biocenosis? The answer may also help to throw light on the role of soil animals in primary successions, since previous studies [START_REF] Ponge | The forest regeneration puzzle: biological mechanisms in humus layer and forest vegetation dynamics[END_REF] have shown that soil animals play a decisive role in the regeneration of mountain coniferous forests.

Study site

We studied the crusts on the Nouragues inselberg, which is 100 km SSE of Cayenne (French Guiana) at 4°5'N and 52°42'W, with its peak at 430 m above sea level and protruding from the rain forest. The rock is a tabular outcrop of pinkish Caribbean granite, containing 27% potassium-feldspar, 37% plagioclase, 33% coarse-grained quartz and 2% accessory minerals [START_REF] Grimaldi | Geography and climate[END_REF]. The inselberg is domeshaped and elongated from East to West. It lacks any forest cover on the eastern flank, where rock is denuded or covered by savane-roche vegetation. The rock is covered by cyanobacterial biofilms made of several species of filamentous cyanobacteria, among which the genera Stigonema, Scytonema and Schizothrix are dominant on free-drained slopes. The dominant plant species is the Bromeliaceae Pitcairnia geyskesii, which forms dense, continuous carpets that cover flat and moderately sloping areas. On the periphery of these impenetrable, spiny carpets, are cyanobacterial crusts (Figure 1), especially in the upslope part. Pitcairnia geyskesii and Clusia minor L. establish themselves from seed, and mosses and lichens develop in these spongy, undulating crusts, which can collect water and nutrients [START_REF] Sarthou | Mécanismes de colonisation par la végétation d'un inselberg granitique en Guyane Française[END_REF]. The crusts are thus the main regeneration niche of the savane roche.

Materials and methods

Selection of sampling plots

We sampled the soil fauna and humus profiles at the upslope border of three bromeliacean carpets of varying size, considered as representative of the variation we observed on the inselberg. Pit1 was a micro-carpet (60 cm diameter) made of only a few rosettes of P. geyskesii, Pit2 was a small carpet (3 to 5 m diameter), and Pit3 was a wide carpet (10-15 m diameter). Supposing that carpets of P. geyskesii increase in size in the course of time through the sympodial development of this plant, we judged that our samples represented three stages of the development of bromeliacean carpets in the order Pit1<Pit2<Pit3. In all three cases the cyanobacterial crust was well developed, its thickness increasing in the same order as above (Table 1). Notice that the absence of replication was dictated by the need to avoid destruction of the regeneration niche of the savane roche at the surface of the inselberg. Sampling of plants and animals (invertebrates included) should be minimized, too, in such fragile, surface-restricted, environments. We sampled the humus profile and the soil for arthropods on 1 and 2 May 2002, during a rainy period.

Sampling of fauna

For sampling soil arthropods, a 15-cm diameter and 10-cm deep stainless steel cylinder was forced by hand into the crust, without discarding any debris, until the granite was reached. The whole content of the cylinder was thoroughly excavated with a spoon, then put in a plastic bag before transport to the laboratory (Brunoy, France). The extraction of soil arthropods was done by the dry funnel method [START_REF] Edwards | A comparison of extraction methods for terrestrial arthropods[END_REF].

Animals escaping the desiccating samples were preserved into 95% ethyl alcohol until we could study them. They were separated into morphospecies, on the basis of characters observable under a dissecting microscope, then classified in broad taxa (Table 1). This allows thorough comparisons between samples and calculation of biodiversity indices, using morphospecies in place of true species as unit taxa, in case species names are unknown [START_REF] Oliver | Designing a cost-effective invertebrate survey: a test of methods for rapid assessment of biodiversity[END_REF].

Densities of animals were used for the calculation of total abundance of arthropods and number of morphospecies per sample. We measured the size of all morphospecies (overall length of the body, appendices excluded) to the nearest 0.04 mm, using a calibrated reticle in the eye piece. For each morphospecies we selected the size of fully developed specimens (adults). This allowed us to calculate the mean potential size of morphospecies in a given sample, by averaging that of all individual taxa.

Sampling of humus profiles

Three samples (which we denote Pit1, Pit2, and Pit3) were taken in the immediate vicinity of samples for fauna extraction. A block of humus 25 cm² in area was cut with a sharp knife, with as little disturbance as possible, until the granite was reached (2, 7.5 and 2 cm depth, respectively), and the crust surrounding it was gently excavated, preserving the zone where soil animals were to be sampled the day after. The various layers within the blocks were distinguished by eye in the field. Each layer was fixed immediately in 95% ethanol. The thickness of each layer was measured and annotated according to the nomenclature of [START_REF] Brêthes | Classification of forest humus forms: a French proposal[END_REF], which was adapted for biological crusts. The organic horizons were classified into OL (intact crust) and OH (humified crust), other soil horizons being absent. When several layers were sampled in the same horizon (on the basis of visible differences), sub-samples were numbered successively, for example OH1, OH2 and OH3. Samples of plants growing nearby (P. geyskesii, C. minor, Scleria cyperina Willdenow ex Kunth) were also taken and preserved in ethyl alcohol to allow subsequent identification of plant debris.

In the laboratory, we spread each sub-sample (layer) gently with our fingers in a petri dish, taking care not to break the aggregates; the petri dish was then filled with 95% ethanol. Plant samples from the site were also placed in alcohol. The soil specimens were examined under a dissecting microscope at 50 X magnification with a cross reticle in the eye piece. A transparent film with a 200-point grid was placed above the preparations. At each grid point, using the reticle as an aid for fixing the position, we identified and counted the material beneath it. The relative volume percentage of a given component was estimated by the ratio of the number of points identified to the total number of points inspected.

The various kinds of plant debris were identified visually by comparison with reference plant samples. Litter components (leaves, twigs) were classified according to plant species and decomposition stages on the basis of morphological features. Animal faeces were classified according to animal groups when possible, on the basis of previous experience with temperate (Ponge, 1991a, b;[START_REF] Topoliantz | Earthworm and enchytraeid activity under different arable farming systems, as exemplified by biogenic structures[END_REF] and tropical soils [START_REF] Topoliantz | Réponse Fonctionnelle de la Pédofaune à la Mise en Culture Itinérante et Permanente des Sols du Sud-Ouest de la Guyane Française[END_REF]. When necessary, the identification of humus components was checked at higher magnification. For that purpose, a small piece of a given humus component was collected with scissors then mounted in a drop of chlorallactophenol for examination in a phase contrast microscope at 400 X magnification.

After quantification of humus components was completed, each petri dish was thoroughly inspected under the dissecting microscope, in order to collect all enchytraeids that were present in the corresponding layer. This allowed us to add these terrestrial annelids to morphospecies extracted as mentioned above, by extrapolating their numbers to the same surface area (1.8 dm 2 ). Notice that estimates of enchytraeid densities were not done on the same samples as for arthropods, but on samples taken in the immediate vicinity. Some animals were mounted in chloral-lactophenol for examination of their gut contents by transparency under phase contrast (Ponge, 1991a, b).

Data analysis

In the absence of replication within each biological crust studied, we did not attempt statistical tests on our faunal data. However, to compare the species richness of our three samples on the same volume basis (thus avoiding artefacts due to variable sampling effort), we calculated the theoretical number of morphospecies per 20 cm 3 by a jackknife method [START_REF] Legendre | Numerical Ecology[END_REF]. For that purpose we divided the total number of arthropods found in a given sample (20 cm 2 section, variable depth) by its thickness in cm. This allowed us to estimate the number of arthropods, had our sampling effort been restricted to only one cm of humus. This estimate was used for a random selection of arthropods within the total population, this sub-sample being used for a new calculation of the number of morphospecies. The same process was replicated ten times, allowing to compute an average value. We found that ten replicates were enough to give a confidence interval (standard error of the mean times Student's t) of less than 10% of the mean. This procedure allowed non-biased estimates of biodiversity when comparing samples of varying size (here the volume of soil sampled varied from 0.35 to 0.7 dm 3 ).

Percentages of occurrence of humus components in the various layers of the three humus profiles investigated were subjected to correspondence analysis (CA), a multivariate method using the chi-square distance between individuals and between variables in a symmetrical manner [START_REF] Greenacre | Theory and Applications of Correspondence Analysis[END_REF]. The different classes of humus components were the active (main) variables, coded by their percentage of occurrence by volume, estimated by the corresponding number of points divided by the total number of points counted. The three plots (Pit1, Pit2, and Pit3) were treated as passive (additional) variables, i.e. they were projected on the factorial axes as if they had been involved in the analysis, without contributing to the factorial axes.

In order to give the same weight to all variables (active and passive) they were transformed to mean 20 and unit variance by

X = (x-m)/s + 20,
where x is the original value, m is the mean of a given variable and s is its standard deviation. The addition to each standardized variable of a constant 20 made all values positive, because correspondence analysis deals only with positive numbers, commonly counts. Following this transformation, factorial coordinates of variables can be interpreted directly in terms of their contribution to the factorial axes. The transformation used here gives to correspondence analysis most properties of the wellknown principal components analysis, while keeping the advantage of the simultaneous projection of variables and samples on the same factorial axes and the robustness due to the principle of distributional equivalence.

Results

The composition of humus profiles

We identified 71 humus components in the 10 layers we sampled in the three pits (Table 2). They were pooled into 14 gross classes on the basis of affinities in their composition.

Notice that these gross classes were not mutually exclusive. For instance, all components comprising pieces of intact cyanobacterial crust were included in the gross class 'intact crust', while some of them were also included in other gross classes (Table 2). This allowed an overview of the vertical distribution of humus components in the three profiles (Figure 2).

The two thinnest biological crusts, Pit1 and Pit2, corresponding to the smallest bromeliacean carpets, exhibited strong affinities (Figure 2). Both comprised an intact crust overlying a thin layer made of roots of P. geyskesii and assemblages of humus (strongly transformed organic matter) and cyanobacteria (visible as intact cell files under the dissecting microscope). No animal biogenic structure was visible, and living cyanobacteria were present over the whole humus profile. Fungi (mycelial wefts) were visible in the top first cm, but disappeared below. Plant litter was almost absent. In these two cases plants (here only P. geyskesii) contributed by their root systems, not by their aboveground parts. No trace of animal activity (faecal deposits) was visible. We observed that the humus material was not firmly aggregated, spilling out easily in the petri dishes even after ethyl alcolhol had precipitated colloids.

In contrast, Pit3 contained a large volume of enchytraeid faeces, making up to 35% of the total solid matter (Figure 2). The intact cyanobacterial crust was overlaid by a thin litter, made of flower stems of P. geyskesii and leaf fragments of C. minor. Below the crust, decayed litter was present, making up to 52% of the total solid matter, except in the bottom-most cm where roots dominated the matrix (46%). The percent in volume of amorphous humus without cyanobacteria increased from the top to the bottom of the humus profile. There was a small amount of mite faeces, mainly in the cyanobacterial crust. Few free cyanobacteria were visible, unlike Pit1 and Pit2. More mineral particles were also present in Pit3 than the two other profiles, and we observed that the humus material was firmly aggregated.

Correspondence analysis, made on all humus components (71), revealed that two major factors influenced the composition of humus layers (Figure 3). Axis 1 (29% of total variance) displayed the vertical distribution of humus components, while axis 2 (21% of total variance) displayed differences between the three crusts. The surface layers (OL horizon) were projected on the positive (right) side of Axis 1, while deeper layers (OH horizon) were projected on the negative (left) side of this axis. The distribution of humus layers along axis 1 (Figure 5) showed that more stratification occurred in the composition of the biological crust in Pit3 than in Pit1 and Pit2. Figure 3 (see Table 2 for codes) showed that Pit3 was characterized by several litter components above the cyanobacterial crust (right side of axis 1), which were replaced by enchytraeid faeces and humified organic matter in the OH horizon (left side of axis 1). Axis 2 showed that some components were present in Pit1 and Pit2, but not in Pit3. This concerned mainly the composition of the crust, which included mosses, mycelial wefts, humified crust and some particular microbial colonies, which were absent or poorly represented in Pit2.

The soil fauna

We separated the total arthropod fauna into 54 morphospecies (Table 1). Among arthropods, mites were the most diverse group, amounting half the total number of morphospecies. The group was the most dense per unit surface. Springtails were also dense, although they were much less diverse than the mites. Enchytraeids were present only in Pit3, where they were denser than all other animal groups, amounting 270 000 individuals per square metre, more than 3 times the total number of arthropods.

Pit3 was by far the most populated biological crust, both in individuals and in morphospecies (Table 1). We may wonder whether the richness of species observed in Pit3 compared with Pit1 and Pit2 was due to sampling effort, the volume of soil extracted being nearly 4 times that of Pit1 or Pit2. In other words, suppose we had taken the same volume of humus, we may wonder whether the number of morphospecies would be the same or not in our three samples. The jackknife procedure showed that, even when reweighted to 1 cm of humus (177 cm 3 ), the total number of arthropod morphospecies in Pit3 was still more than twice that in Pit1 or Pit2. Thus the increase in species richness we observed in Pit3 was not an artifact caused by our taking a larger volume of humus.

The mean size of arthropod taxa, averaged over all morphospecies present in a sample, was greater in Pit3 than in Pit1 and Pit2 (Table 1). This was due to the presence in Pit3 of bigger mites such as MOR 35' (1.2 mm) and of macro-arthropod species (arachnids other than mites, millipedes, insect larvae and adults), the size of which often exceeds 1 mm.

Examination of the contents of the enchytraeid guts from Pit2 revealed that these annelid worms ingested a variety of humus components, among which mite faeces, cyanobacteria and fungi were dominant, but a small amount of quartz grains was also observed. The disintegration of mite faeces was evident, suggesting that these animals transformed mite faeces into their own humus (faecal material). Examination of mite faeces within enchytraeid guts revealed a diet only of plant material, with a large amount of dark brownish-red undigested condensed tannins. Humified organic matter (whether included in animal faeces or not) showed, when examined under phase contrast, that bacteria, living and dead, were abundant, and formed the bulk of humus.

Discussion

The examination of cyanobacterial crusts under the light microscope revealed profound differences in the composition of humified organic matter depending on whether enchytraeids were present or not. Arthropods, mostly mites, are abundant and diverse in the three profiles, but their faeces do not accumulate, contrary to what has been observed in temperate organic soils [START_REF] Bal | Morphological investigation in two moder-humus profiles and the role of the soil fauna in their genesis[END_REF]. In the presence of enchytraeids (in Pit3), most mite faeces are probably ingested by these worms and are transformed into enchytraeid faeces, which accumulate below the cyanobacterial crust. The same substrate succession from mites to enchytraeids has been already observed by Ponge (1991a) in temperate Scots pine litter. Aggregates of humus and living cyanobacteria are common in Pit1 and Pit2, but no mite faeces are recognizable, contrary to expectation, given the high density and diversity of these microarthropods. It seems that mite faeces, which are probably produced in large amount in Pit1 and Pit2, disintegrate rapidly. Part of the resulting colloidal matter is agglomerated by cyanobacteria, forming aggregates of humus and cyanobacteria as observed by [START_REF] Rao | The effect of surface growth of blue-green algae and bryophytes on some microbiological, biochemical, and physical soil properties[END_REF]. In our samples, such aggregates formed a compact horizon, while enchytraeid faeces formed loose assemblages. However, in our alcohol-preserved samples, these aggregates lost particulate matter easily when suspended in alcohol, whereas enchytraeid faeces did not.

We suspect that in the presence of intense run-off, crusts without enchytraeids are less stable than those with enchytraeids. Thus, the development of enchytraeid populations creates a structure [START_REF] Didden | Involvment of Enchytraeidae (Oligochaeta) in soil structure evolution in agricultural fields[END_REF], probably resistant to run-off, allowing crusts to increase in thickness insofar as carbon enters the system via photosynthesis and litter input. The cementation of organic and mineral particles by enchytraeid polysaccharides (mucus) has already been reported [START_REF] Chan | Occurrence of enchytraeid worms and some properties of their casts in an Australian soil under cropping[END_REF].

We observed the establishment by seed of P. geyskesii and C. minor in the biological crusts upslope of the bromeliacean carpets. The favourable influence of biological crusts on seed capture, germination and establishment of seedlings has been reported [START_REF] Belnap | Influence of biological crusts on soil environments and vascular plants[END_REF], and we can readily imagine that the pore network created by enchytraeids [START_REF] Didden | Involvment of Enchytraeidae (Oligochaeta) in soil structure evolution in agricultural fields[END_REF] can help the development of plant root systems, and can decrease nutrient losses due to run-off [START_REF] Isichei | Mineral nutrient flow from an inselberg in south-western Nigeria[END_REF]. The disturbance caused by the burrowing of these animals might also contribute to alleviate competition from already established vegetation. In the same way, the destruction of cyanobacteria in enchytraeid guts can locally decrease their toxicity [START_REF] Burja | Marine cyanobacteria: a prolific source of natural products[END_REF], suspected to occur from the census of dominant taxa [START_REF] Sarthou | Cyanophycées de l'inselberg des Nouragues (Guyane française)[END_REF]. Such, still imperfectly known mechanisms could be added to interactions among plants which have been suspected to explain facilitation by early successional plant species on granitic outcrops [START_REF] Houle | Seed availability and biotic interactions in granite outcrop plant communities[END_REF].

We observed fungi and some other microbial colonies in cyanobacterial crusts in the absence of any lichen symbiosis, as did [START_REF] States | Soil fungi as components of biological soil crusts[END_REF]. We also found that these fungi were preferentially consumed by invertebrates, in particular mites and springtails, although fungal spores and mycelia were much less abundant than cyanobacteria. Thus, fungi could serve as food for microarthropods, although they do not seem to influence the structure of cyanobacterial crusts to a great extent.

The large number of arthropod species in Pit3 was probably mediated by enchytraeid activity, as [START_REF] Loranger | Impact of earthworms on the diversity of microarthropods in a vertisol (Martinique)[END_REF] demonstrated to occur with earthworms.

In bromeliacean carpets, which are very poor in plant species, the large number of invertebrate species we found in associated biological crusts is worth to note, even in the absence of enchytraeid activity. Densities of microarthropods and species richness are similar to values reported by [START_REF] Petersen | A comparative analysis of soil fauna populations and their role in decomposition processes[END_REF] for temperate organic soils.

The large arthropod morphospecies in Pit3 can be ascribed to the favourable effect of annelid activity, too. [START_REF] Haarløv | Vertical distribution of mites and Collembola in relation to soil structure[END_REF] showed that the size of arthropods correlated well with pore and particle size. Also, [START_REF] Loranger | Impact of earthworms on the diversity of microarthropods in a vertisol (Martinique)[END_REF] observed that pores and channels created by earthworms favour bigge microarthropods without disfavouring smaller ones, thus increasing biodiversity. Enchytraeids might play the same role as earthworms, by creating space for bigger arthropods in the pore network they create. This, among other processes, might explain why invertebrates of increasing size can be observed in the course of primary successions [START_REF] Rusek | Pedozootische Sukzessionen während der Entwicklung von Ökosystemen[END_REF].

In conclusion, we have identified two distinct steps in the appearance of a pristine soil structure in biological crusts on tropical inselbergs: (i) a thin crust formed by cyanobacteria, mixed with fungi and other bacterial colonies, that retain humified organic matter loosely attached to cyanobacterial filaments (Pit1 and Pit2)

(ii) a thick crust formed by cyanobacteria and plant litter, covering a structure made of organic enchytraeid faeces, into which plant roots proliferate (Pit3).

In that sense, enchytraeids could be considered as ecosystem engineers, as earthworms, termites and ants have been in other kinds of environments.

Figure legends

Figure 1 A bromeliacean carpet of Pitcairnia geyskesii with its upslope biological crust (Pit2)

Figure 2 The vertical distribution of humus components (pooled in bulk categories, see

  Fig. 1

Table 2 )

 2 in the three studied profiles. Depth and organic soil horizons in abscissa and percent volume of solid matter in ordinate (scale indicated at the top of the

	graph)
	Figure 3 Correspondence analysis on micromorphological data. Projection of samples
	and variables in the plane of the first two factorial axes (vertical stratification
	and differences among profiles, respectively). Top: projection of humus layers
	(black dots, linked by lines according to vertical succession). Bottom: projection
	of humus components (meaning of codes in Table 2) and three profiles in the
	plane of the first two factorial axes. Scale intervals 0.01. Symbols for humus
	components and profiles centered on point coordinates

Acknowledgements

We thank several state and regional organizations for financial support and commodities, P. Charles-Dominique and his staff at the Nouragues Field Station and P.

Chatelet and R. Garrouste for technical help and fruitful discussions.

Code Component

Gross category