Jo ¨rn Bennewitz
email: jbennewitz@tierzucht.uni-kiel.deoped

Norbert Reinsch

Volker Guiard

Sebastien Fritz

Hauke Thomsen

Christian Looft

Christa Ku

Manfred Schwerin

Christina Weimann

Georg Erhardt

Fritz Reinhardt

Reinhard Reents

Didier Boichard

Ernst Kalm

Multiple Quantitative Trait Loci Mapping With Cofactors and Application of Alternative Variants of the False Discovery Rate in an Enlarged Granddaughter Design

The experimental power of a granddaughter design to detect quantitative trait loci (QTL) in dairy cattle is often limited by the availability of progeny-tested sires, by the ignoring of already identified QTL in the statistical analysis, and by the application of stringent experimentwise significance levels. This study describes an experiment that addressed these points. A large granddaughter design was set up that included sires from two countries (Germany and France), resulting in almost 2000 sires. The animals were genotyped for markers on nine different chromosomes. The QTL analysis was done for six traits separately using a multimarker regression that included putative QTL on other chromosomes as cofactors in the model. Different variants of the false discovery rate (FDR) were applied. Two of them accounted for the proportion of truly null hypotheses, which were estimated to be 0.28 and 0.3, respectively, and were therefore tailored to the experiment. A total of 25 QTL could be mapped when cofactors were included in the model-7 more than without cofactors. Controlling the FDR at 0.05 revealed 31 QTL for the two FDR methods that accounted for the proportion of truly null hypotheses. The relatively high power of this study can be attributed to the size of the experiment, to the QTL analysis with cofactors, and to the application of an appropriate FDR.

M UCH effort has been undertaken to identify quan-experiments used the existing male half-sib structure by setting up the so-called granddaughter design (Weller titative trait loci (QTL) associated with genetic et al. 1990). In this design the pedigree consists of sets variation for traits of economic or scientific interest in of genotyped and progeny-tested male half-sib families. livestock species with the aid of genetic markers. Hayes

The power to detect QTL in a granddaughter design is and Goddard (2001) estimated the number of QTL in largely influenced by the size of the half-sib families, dairy cattle for a trait that undergoes selection at about which has its limit in the availability of progeny-tested 50 to 100 depending on the size of the effective populasons. Interval mapping is one of the most applied statistition. With regard to this, the majority of QTL remained cal methods in the analysis of granddaughter designs. undetected by the experiments undertaken so far. The This method uses information from consecutive informain reason for this is that the distribution of the QTL mative markers simultaneously to trace the inheritance effects follow likely a gamma distribution with many of a putative QTL [START_REF] Knott | ple-marker mapping of quantitative trait loci in half-sib popula-Benjamini[END_REF]. A drawback of QTL of small effects and only a few of large effects this approach is that it does not take into account QTL [START_REF] Cattle | The distribution of the effects of Education, Science, Research and Technology. of genes affecting quantitative traits in livestock[END_REF] and that QTL experiments outside the respective marker interval or even QTL on conducted to date are not powerful enough to detect other chromosomes. This might result in a bias of the the QTL of smaller effects.

QTL parameter estimates and in a reduced power of In QTL mapping in pigs it is common to generate the experiment, because the variance explained by addian F 2 cross of different, ideally divergent selected, breeds.

tional QTL appears in the residual of the applied inter-As generating such a cross in dairy cattle is a timeval-mapping model. To overcome these limitations Janconsuming and costly process, most dairy cattle QTL sen (1993) proposed composite interval mapping. This method considers putative QTL in the respective marker interval as well as QTL outside this interval by monly used in outbred populations like granddaughter designs, mainly due to the variability of the information Description of the pedigree and the distribution of the content and due to different haplotypes across half-sib progeny testing and the genotyping of the sires families [START_REF] Hoeschele | Advances in statistical methods to map quantitative trait loci[END_REF] been included in a previous study (Bennewitz et al. 2003a).

(1998) in QTL mapping. Recently, Storey and Tibshi-Genotypes: As it was not possible to analyze the whole genome during this study, nine chromosomes of special interest rani (2003) developed a FDR test procedure that is (BTA2, were selected. less conservative in comparison to the FDR procedure All German sires and a proportion of the French sires were proposed by Benjamini and Hochberg (1995) and is genotyped in Germany for the German set of markers tailored to genomewise experiments such as DNA mi-(Table 1). The remaining French sires were genotyped in croarray experiments.

France during the French genome-analysis experiment for the French set of markers. The numbers of common markers in The aim of this study was to analyze a large-scale both the German and the French sets were limited. Therefore, granddaughter design that included in total almost 2000 03ف German half-sibs of each of these families were additionsons from the German and the French dairy cattle geally genotyped for the French marker set to increase the accunome analysis projects [START_REF]A male bovine linkage map for the ADR granddaughter design[END_REF]; Boichard racy of the haplotype derivation of the common grandsire. It et al. 2003) for nine chromosomes and six traits. The is important to note that the members of these families were genotyped heterogeneously; i.e., not all members of a family data were analyzed using a model that included multiple were genotyped for a particular marker. In total 127 markers QTL. Two versions of the new FDR procedure of Storey were included, almost all microsatellite markers. To avoid a and Tibshirani (2003) were applied. The results show common standardization of the markers and to ensure the the advantage of the multiple-QTL model and demonanonymity of the sires the genotypes were coded. All genotypes strate that the new FDR approach is better suited and were transferred to a central database [START_REF] Reinsch | A multiple-species, multiple-project database for designs: a simple strategy for QTL confirmation and increasing genotypes at co-dominant loci[END_REF]) and checked for their agreement with the Mendelian laws of inher-less conservative in this experiment compared to the itance. For additional information about the genotype scoring classical FDR test procedure. procedure see [START_REF]A male bovine linkage map for the ADR granddaughter design[END_REF] and Boichard et al. (2003). Multipoint marker maps were calculated with the use of CRIMAP (Green et al. 1990). The estimated genetic maps MATERIALS AND METHODS with additional information are available at http:/ /www. tierzucht.uni-kiel.de/QTL_ADR_INRA.htm. Pedigree: The total pedigree consisted of 1977 Holstein sires distributed over 18 families. Seventeen (13) families in-Phenotypes: The traits milk, fat and protein yield, fat and protein percentage, and somatic cell score were included. cluded progeny-tested sires from Germany (France) and 12 families included progeny-tested sires in both countries. The For the French sires daughter yield deviations were used as provided by the French national computing center. For the total number of German (French) sires was 896 (1081). The German sires were included in the German granddaughter German sires no daughter yield deviations were available. Estimated breeding values were therefore taken from the routine design [START_REF]A male bovine linkage map for the ADR granddaughter design[END_REF] and a proportion of the French Powerful QTL Analysis in Dairy Cattle sire evaluation and were deregressed as described by Thomsen

where yЈ ijk is the adjusted phenotype of progeny i within family j for the chromosome h and the remaining variables are as et al. (2001). Note that [START_REF] Thomsen | Mapping of multiple quantitative trait loci by simple tions and de-regressed proofs within a whole genome scan for regression in half-sib designs[END_REF] showed the almost equivalency of daughter yield deviations and dere-defined in (2). For each candidate chromosome the phenotypes were adjusted separately by setting the regression coeffi-gressed breeding values for the use in QTL-mapping experiments. The daughter yield deviations were multiplied by 2 cient b ik for the QTL mapped on that particular chromosome to zero. For the noncandidate chromosomes the phenotypes to make them comparable to estimated breeding values. All phenotypes were expressed in genetic standard deviations as were adjusted for the full set of cofactors. The QTL analysis was now repeated (step 1, model 1) including the permutation provided by each country. Within each family, the phenotypic means of the two half-sib groups (i.e., the within-family mean test, but now with the corresponding adjusted phenotypes instead of using the original phenotypes. If this reanalysis of the half-sibs progeny tested in Germany and the withinfamily mean of the half-sibs progeny tested in France) were revealed new candidate chromosome(s) (p c Յ 0.05), the second step (estimating of cofactor effects, model 2) and the subtracted from the corresponding original phenotype. Additionally, the variances of the phenotypes of these two half-sib third step (phenotype adjustment, Equation 3) were repeated and step 1 was conducted again. The analysis ended when no groups were standardized. The genetic correlation between breeding values for traits milk yield, fat yield, protein yield, new candidate chromosome(s) were identified after performing step 1. See de [START_REF] Thomsen | Mapping of multiple quantitative trait loci by simple tions and de-regressed proofs within a whole genome scan for regression in half-sib designs[END_REF] for a graphical and somatic cell score estimated in Germany and estimated in France are always Ͼ0.87 (results from interbull evaluations, presentation of this protocol.

A grandsire was assumed to be heterozygous at a significant www-interbull.slu.se), indicating that the traits investigated were the same in both countries.

QTL when the haplotype contrast at the estimated QTL position was significant at P Յ 0.05 (t-test). QTL substitution effects Statistical analysis: The most likely marker haplotype of the grandsires was determined using the genotype information of [START_REF] Falconer | Introduction to Quantita-Weller[END_REF] were calculated as the average of the substitution effect estimates from the QTL heterozygous the progeny. For each offspring the probability of inheriting the father's segment of the first chromosome for each centi-grandsires.

In a previous study we detected a statistical QTL-by-environ-morgan was calculated using the genotype information of two consecutive informative markers [START_REF] Knott | ple-marker mapping of quantitative trait loci in half-sib popula-Benjamini[END_REF], where ment interaction when analyzing two half-sib groups that share the same father but were progeny tested in different countries possible, simultaneously. These probabilities were termed QTL transition probabilities (p ) and were retained for the (Bennewitz et al. 2003a). This interaction is defined as the occurrence of significance of a QTL only in one of the two half-rest of the analysis.

The QTL analyses followed the suggestions of de Koning sib groups and can have many reasons such as an interaction of the QTL with the polygenic background that might be differ-et al. (2001) and were performed across families for each trait separately by performing the following three steps. In a first ent in the two groups, a QTL-by-environment interaction in a strict sense, a type I error, or a type II error. In the present step the chromosomes were scanned to identify putative QTL using the following regression model [START_REF] Knott | ple-marker mapping of quantitative trait loci in half-sib popula-Benjamini[END_REF], study the existence of statistical QTL-by-environment interaction as defined above was tested for the chromosomal positions

y ij ϭ gs i ϩ b ik p ijk ϩ e ijk ,
(1) that harbor putative significant QTL by applying the following model, where y ij is the original trait value of sire j of the grandsire i, gs i is the fixed effect of the grandsire i, b ik is the regression y hijk ϭ gs i ϩ E h ϩ b ik ϫ p ijk ϩ ihk ϫ p ijk ϩ e hijk , (4) coefficient of the grandsire i at the chromosomal location k and represents the QTL allele substitution effect (Falconer where E h is the environment h (progeny tested in either Gerand Mackay 1996), and p ijk is the QTL transition probability many or France), ihk is the interaction of the h environment as defined above for sire j in the grandsire family i at chromoand the QTL transition probability, and the remaining varisomal position k. The null hypothesis was that no QTL segreables are as defined above. Note that this was no systematic gates on this chromosome for the trait under consideration, search for the presence of a statistical QTL-by-environment the alternative hypothesis was that one QTL segregates on interaction. this chromosome. The test statistic was an F-ratio defined as For the final set of identified QTL the genomewise error the mean square deviation of regression divided by the mean probabilities (p g ) were calculated using the Bonferroni correcsquare deviation of residuals pooled across families. Chromotion assuming 30 chromosomes, p g ϭ 1 Ϫ (1 Ϫ p c ) 30 . Additionsomewise test statistical critical values were obtained by the ally, a Bonferroni correction approach was applied to calculate use of the permutation test (Churchill and Doerge 1994), experimentwise error probabilities (p e ) assuming nine indeperforming 10,000 permutations. Following the suggestions pendent tests (9 chromosomes investigated) that were conof de [START_REF] Thomsen | Mapping of multiple quantitative trait loci by simple tions and de-regressed proofs within a whole genome scan for regression in half-sib designs[END_REF] a chromosome was declared as a ducted six times (six traits), i.e., p e ϭ 6 ϫ (1 Ϫ (1 Ϫ p c ) 9 ). candidate for carrying a putative QTL when the chromo-Note that this is a rough estimation of the experimentwise somewise error probability (p c ) was Յ0.05. In a second step error rate, because the dependence structure of the tests was the QTL transition probabilities at the position with the highnot quantified. Hence, p e can be Ͼ1 in some cases. Noncentral est test statistic on the candidate chromosomes were included confidence intervals for the estimated QTL position were calas cofactors in the following model to estimate the effect of culated by permutation bootstrapping (Bennewitz et al. all cofactors simultaneously, 2002), performing 250 bootstrap samples. This bootstrap method corrects for the marker impact on the distribution of the QTL position estimates along the chromosome from the

y ij ϭ gs i ϩ ͚ n kϭ1 b ik p ijk ϩ e ijk ,
(2) evaluated bootstrap samples, taking the results of the permutation test into account. where n is the number of total identified candidate chromo-Computing the FDR: In the following a brief description somes, k is here the chromosomal position of the cofactor(s), of the FDR as calculated in this study is given, based on the and the remaining variables are as defined in (1). In a third studies of Benjamini and Hochberg (1995), Storey and step the original phenotypes were adjusted for the estimated Tibshirani (2003), and [START_REF] Storey | Implementation of marker assisted selection consistency of false discovery rates: a unified approach[END_REF]. In multiplecofactor effects using the formula testing procedures the possible outcomes of m tests are the number of hypotheses declared as significant (S ), where F is

yЈ hijk ϭ y ij Ϫ ͚ n kϭ1 b ik p ijk , (3)
the number of false positives and T is the number of true [START_REF] Weller | A new approach to the problem of multiple comparisons Montpellier, France. in the genetic dissection of complex traits. Genetics 150[END_REF] emphasized that q(i ) could sometimes TABLE 2 decrease after increasing i. This is, however, incorrect. He Possible outcomes from multiple testing erroneously used the formula q(i ) ϭ m ϫ p(i )/i without the minimization shown in (8). However, this does not correspond Declared as Declared as to the FDR method of Benjamini and Hochberg (1995 (2004) proposed an estimation method for 0 and showed that after applying this method FDR will not be greater than Adapted in slightly modified form from Benjamini and q. Using the method of Storey and Tibshirani (2003) and Hochberg (1995) and Storey and Tibshirani (2003). [START_REF] Storey | Implementation of marker assisted selection consistency of false discovery rates: a unified approach[END_REF], 0 is estimated by taking into account that the P-values of true null hypotheses are distributed uniformly in the interval [0, 1] and the P-values of true alternative positives, and the number of hypotheses declared as nonsighypotheses will be closer to zero. nificant (m Ϫ S ; see also Table 2). Let m 0 denote the number In this study three different values of 0 for the 54 chromoof true null hypotheses. F, T, and S are random variables and somewise error probabilities were used in ( 7) and ( 8). The only S is known and depends on the multiple-test procedure first method for 0 estimation is described by Storey and to be described. The FDR is the expected proportion of false Tibshirani (2003); the estimate is termed ˆ0,S (resulting in positives out of the tests that are declared as significant FDR S ). The second method was introduced by Storey et al. multiplied by the probability Pr(S Ͼ 0) and can be written as ( 2004) and minimizes the mean squared error (MSE) of ˆ0 using a bootstrap approach (resulting in ˆ0,B and FDR B , respec-FDR ϭ E[F/S |S Ͼ 0]Pr(S Ͼ 0)

(5) tively). Both methods use the information provided by the (Benjamini and Hochberg 1995). Storey (2002) denotes distribution of the p c -values. We used the software offered by the conditional expectation E[F/S |S Ͼ 0] as positive false Storey and Tibshirani (2003) for the estimation of ˆ0,S and discovery rate (pFDR). This value seems to be more interesting ˆ0,B , which is available at http:/ /faculty.washington.edu/ as a significance criterion, but it is difficult to control the pFDR فjstorey/qvalue. Additionally, as proposed by Benjamini and in situations when the proportion of true null hypotheses is Hochberg (1995), a value of 1 was used for 0 (resulting in high. When all null hypotheses are true then every discovery 0,BH and FDR BH , respectively). is a false discovery. Hence, for m 0 ϭ m a true discovery is not possible and the pFDR is always equal to one. In this case the user should be interested in avoiding every discovery. This RESULTS means that Pr(S Ͼ 0) should be small. The FDR significance criterion takes into consideration this interest of the user. In

The results from the statistical analysis with the final the case m 0 ϭ m we have FDR ϭ Pr(S Ͼ 0). Therefore we set of cofactors are shown in Tables 3 and4. In total prefer to use the FDR as a significance criterion. Note that we found 25 chromosomewise significant QTL (p c Ͻ both criteria, pFDR and FDR, are almost identical in situations 0.05) distributed over all chromosomes and all traits where it is known that Pr(S Ͼ 0) Ϸ 1.

analyzed, 12 of which were genomewise significant (p g Ͻ Using the FDR as a significance criterion with a chosen FDR level q according to the test procedure of Benjamini and 0.05, Table 4), and 9 of which were experimentwise [START_REF] Jansen | Interval mapping of multiple quantitative trait[END_REF], it is necessary that the m tests are ordered significant (p e Ͻ 0.05, Table 4). A highly significant QTL by their P-values as p(i ) Յ . . . , p(m ) for i ϭ 1 . . . , m (i.e., i for all five milk production traits was found on BTA14.

is a ranking number of the m tests based on their P-values). Winter et al. (2002) and Grisart et al. (2002) found Let i ˆbe the greatest i fulfilling a nonconservative mutation in a strong candidate gene to be most likely responsible for the genetic variance

p(i ) i • m Յ q , (6)
attributable to this QTL. The effect of this mutation was highly significant for all milk production traits in and then the all hypotheses i with i Յ i ˆwill be rejected. This procedure guarantees that FDR Յ (m 0 /m )q (Benjamini and German Holsteins [START_REF] Thaller | Effects of DGAT1 variants on milk production traits in Boichard[END_REF], with a substan-Hochberg 1995). Let 0 be the proportion of true null tial substitution effect for these traits. The estimated hypotheses among all tested hypotheses, i.e., 0 ϭ (m 0 /m ). If number of QTL heterozygous grandsires was between 0 is small then this procedure is very conservative. If 0 is 2 and 5 out of 18 averaged over all QTL. In general, known, then in (6) we could apply q* ϭ q/ 0 instead of q , grandsires of large families were more frequently estigetting FDR Յ q , which increases the power.

If the threshold q of the FDR is not given in advance, then mated to be QTL heterozygous compared to grandsires for each hypothesis tested a measurement, q(i ), of significance of small families. For example, grandsire F04 was in terms of the FDR is questioned. Here q(i ) is defined as the deemed to be heterozygous for 10 of the QTL mapped smallest FDR threshold q for which the hypothesis i will be whereas grandsire F15 was heterozygous only for 1 (not declared as significant. Assume 0 is known. The calculation shown). The estimates of the substitution effects varied of the q(i ) follows directly from the procedure of Benjamini and Hochberg (1995). It begins with q(m ): between 0.36 and 0.9 genetic standard deviations. The estimated widths of the confidence intervals (Ta-

q(m ) ϭ 0 ϫ p(m ). ( 7
)
ble 3) were larger than expected given the size of the Now, for i ϭ m Ϫ 1, m Ϫ 2, . . . , 1 calculate

experiment. An explanation for these contradictory results is that the families were genotyped heterogeneously;

q(i ) ϭ min 0 ϫ m ϫ p(i ) i
, q(i ϩ 1) .

(8)

i.e., not all members of a family were genotyped for e Estimated QTL position in centimorgans from the start of the chromosome. For flanking markers see estimated maps available at http:/ /www.tierzucht.uni-kiel.de/QTL_ADR_INRA.htm.

f 95% confidence interval.

g Average substitution effect. h Number of heterozygous grandsires. i BTA14 harbored genomewise significant QTL for all milk traits. The causal mutation of this QTL is known (Grisart et al. 2002;Winter et al. 2002); see [START_REF] Thaller | Effects of DGAT1 variants on milk production traits in Boichard[END_REF] for the effects of the mutation in the German Holstein.

the same set of markers (Bennewitz et al. 2003a). A only in the analysis without cofactors. The F-value dropped from 2.2 to 1.9 with the full set of cofactors. In general, significant statistical QTL-by-environment interaction was not found in the across-family analysis. Single-family the estimated QTL positions did not change significantly when cofactors were included in the analysis and analyses revealed significant interaction effects for some of the QTL mapped. However, when these results were the size of the estimated confidence intervals tended to be slightly smaller (not shown). A remarkable outcome corrected for multiple testing using the Bonferroni method assuming 12 grandsires with progeny in both is that the threshold values for the chromosomewise error probabilities decreased in general by 1.0ف units countries (Table 1) the significant interaction effects disappeared.

when permuting the adjusted phenotypes during the analysis with cofactors compared to the simple analysis Without cofactor analysis only 18 chromosomewise (6 genomewise) QTL were found. The number of cofac-(not shown). This emphasized the need to apply the permutation test in every round of cofactor selection tors included for each trait was equal to the number of chromosomewise significant QTL (Table 3). It required rather than permuting only the original phenotypes and using the corresponding threshold values in the one, one, three, two, three, and two round(s) of cofactor selection for the traits milk yield, fat yield, protein yield, subsequent analysis.

The q-values estimated by FDR BH , FDR S , and FDR B are fat percentage, protein percentage, and somatic cell score, respectively. In general the F-values were larger shown in Table 4 for 35 hypotheses with the lowest p c -values. Additionally, the experimentwise error proba-in the cofactor analysis. This was most extreme for fat percentage on BTA20 (Table 3). On the other hand, a bilities (p e ) are presented in Table 4. As expected, the FDR BH produced the highest q-values, followed by the QTL for somatic cell score on BTA02 was significant corresponding experimentwise error probability, except for the q-values that were equal to zero (Table 4).

Results from false discovery analysis with the three different

q-values (q S , q B , and q BH )

DISCUSSION i a Trait b BTA p c c p e d q S q B q BH
In this study a QTL mapping was conducted in a large 1 MY 14 0.0000 0.0000 0.0000 0.0000 0.0000 granddaughter design that consisted of almost 2000 2 FY 14 0.0000 0.0000 0.0000 0.0000 0.0000 sires. Compared to other studies the number of mapped 3 FP 14 0.0000 0.0000 0.0000 0.0000 0.0000 QTL, and hence, the statistical power of the experiment, 4 PP 14 0.0000 0.0000 0.0000 0.0000 0.0000 was high, which is due to the large average half-sib family and for the families with observations in both environ-27 PY 19 0.0799 Ͼ0.9 0.0450 0.0412 0.1500 ments the structure is not very balanced. Additionally, 28 FY 02 0.0800 Ͼ0.9 0.0450 0.0412 0.1500 in contrast to the previous study (Bennewitz et al. 29 PP 23 0.0807 Ͼ0.9 0.0450 0.0412 0.1500 30 PP 26 0.0856 Ͼ0.9 0.0450 0.0412 0.1500 2003a), the test was performed only at chromosomal 31 PY 18 0.0861 Ͼ0.9 0.0450 0.0412 0.1500 positions of significant QTL. If an interaction might 32 FY 06 0.1181 Ͼ0.9 0.0598 0.0558 0.1993 have occurred for a QTL the power of the experiment 33 PY 02 0.1748 Ͼ0.9 0.0858 0.0800 0.2860 to map this particular QTL would be reduced and hence 34 FP 18 0.1840 Ͼ0.9 0.0877 0.0818 0.2922 this QTL would probably be missed.

35 MY 20 0.1997 Ͼ0.9 0.0924 0.0863 0.3081

The general findings of the comparison of the analysis a Ranking number of tests.

with and without cofactors in the model (Table 3) are b For trait abbreviations see Table 3.

in good agreement with de Koning et al. (2001). The c Chromosomewise error probability.

increase of the test statistic, and hence the statistical d Experimentwise error probability.

power when including cofactors, is a result of the reduced residual variance. Additionally, the plots of the test statistic along the chromosomes showed in general a more pronounced maximum for significant QTL when q-values of FDR S and FDR B , which were nearly identical. This is also visualized in the plot of the q-values against cofactors were included in the model (not shown). This led, together with the elevated test statistic, to the the ranking number i of the hypotheses tested (Figure 1). For low i (i Ͻ 20) the q-values were on a similar slightly reduced width of the confidence intervals. The potential benefit of the cofactor analysis would even be level, but with increasing i (i Ͼ 20) the differences between the q-values became substantial. For example, greater if the set of cofactors was complete. In this study this set was incomplete because only nine chromosomes when thresholding the q-values at q Յ 0.05, 31 QTL were declared as significant for both the FDR B and the were included and it is reasonable to assume that some chromosomes not included harbor significant QTL. FDR S , but only 20 for the FDR BH (Table 4, Figure 1). The reason for this is the different estimates for the No QTL was declared as chromosomewise (genomewise) significant with a substitution effect of Ͻ0.35 proportion of the truly null hypothesis 0 . For FDR S this estimate was ˆ0,S ϭ 0.30 and for FDR B it was ˆ0,B ϭ 0.28, (0.65) additive genetic standard deviation (Table 3), bearing in mind that the effects are likely overestimated bearing in mind that for FDR BH it was 0,BH ϭ 1 by definition. Every q estimate was significantly below the (e.g., Go ¨ring et al. 2001). Taking the findings of Hayes -Plot of the q-values against the number of hypotheses declared as significant. The hypotheses were ordered decreasing by their P-values. , q BH ; ᮀ, q S ; and ᭡, q B . and Goddard (2001) into account, various QTL of varies in outbred populations, making the definition of well-chosen marker intervals cumbersome. The lowest smaller effects were missed. Two reasons come into question-either they are all located on the chromo-P-values from consecutive marker intervals would show a strong dependence in those families that are not very somes not included in this study or, more likely, some of them are located on the chromosomes included but informative for markers between the intervals. Indeed this is the main reason why it is difficult to apply compos-were not declared as significant. Hence the relatively high power of the experiment mentioned above is only ite interval mapping in outbred populations and why the method of de [START_REF] Thomsen | Mapping of multiple quantitative trait loci by simple tions and de-regressed proofs within a whole genome scan for regression in half-sib designs[END_REF] rather than that for the detection of larger QTL when using the chromosomewise or even genomewise or experimentwise sig-of composite interval mapping was used to account for multiple QTL. nificance level at P Յ 0.05. Thus, the choice of the right significance criterion for controlling the type I error

The FDR was calculated across traits, going against the recommendation of et al. (2002). These authors rate can be seen as a balance between the experimental power and the probability of making a type I error.

showed that adding low heritability traits to the analysis reduces the power when calculating the FDR across As pointed out by [START_REF] Lee | The control of the false loci interval mapping with multiple traits[END_REF], the concept of controlling the type I error rate on the null hypotheses traits. On the other hand, the across-traits FDR calculation is very attractive because it is not necessary to ac-of no QTL might not be appropriate when the trait under consideration is heritable in the population, be-count for multiple-FDR procedures that would arise when the FDR was calculated for each trait separately. cause in these cases there must be QTL on the genome that are responsible for the genetic variance. The au-

The traits included are all heritable and from other studies it is known that QTL segregate for all traits in thors suggested that in these cases it would make sense to test chromosome regions with QTL against those with dairy cattle populations. Additionally, the estimates of 0 were ,3.0ف meaning that 83ف of the 54 tests represent no QTL, which leads automatically to the concept of the FDR [START_REF] Lee | The control of the false loci interval mapping with multiple traits[END_REF]. The fundamental difference true effects. We therefore supposed that the power is not reduced by the FDR calculation across traits. The between the concept of controlling the type I error rate and controlling the FDR is that the type I error rate is applied alternative across-traits significance criterion was the experimentwise error probability. However, as the proportion of false positives among the true null hypotheses whereas the FDR is the proportion of false shown in Table 4, this criterion is overwhelmingly stringent, resulting in a low experimental power. Consider-positives among all accepted alternative hypotheses.

In this study FDR methods were applied under the ing the across-trait FDR calculation, if the p c -values of the tests corresponding to true null hypotheses are posi-general assumption of a maximum one QTL for each trait on a certain chromosome and subsequently used tively correlated and m 0 is known, then the FDR methods applying (6) or ( 7) and ( 8) are valid and conservative the chromosomewise error probabilities as P-values. Note that no hint for two QTL on a single chromosome (Benjamini and Yekutieli 1997). Since the tests used are two sided, positive or negative correlations between for a particular trait could be found in the data (not shown). Alternatively, [START_REF] Lee | The control of the false loci interval mapping with multiple traits[END_REF] recommended the traits would yield positive correlations between the p c , at least for small values of the p c . Therefore we assume the use of the lowest comparisonwise error probabilities within a marker interval as P-values for the FDR calcula-that also for estimated values of m 0 the FDR methods applied in this study are approximately valid. tion in F 2 crosses. However, the informativeness of markers for a certain marker interval as used in this study Three different FDR methods were applied in this more, in the present study it also proved to be less higher than would be expected if all null hypotheses conservative compared to the classical threshold setting. were true. On the other hand, p c Ͼ 0.1 are almost uni-For example, when using the FDR B instead of the chroformly distributed with a slight decrease with higher p cmosomewise or even experimentwise threshold levels values. Additionally it is known from previous work that more QTL will be declared as significant (Table 4). five QTL are real effects in the design that is the QTL Hence, replacing the classical threshold setting based for the five milk production traits on BTA14 (Thaller on P-values by the presented FDR approach brings the et al. 2003). Therefore, although m is not very large, a number of declared QTL closer to the predicted numvalue of 1 for 0 as suggested by Benjamini and Hochber of Hayes and Goddard (2001). berg (1995) is inappropriate and extremely conserva-In summary, a large granddaughter design with progtive.

eny-tested sires from Germany and France was analyzed. Both the FDR S and the FDR B methods accounted for

The comparatively high number of QTL found emphathe mixture distribution of the p c -values, and both prosized the high power of the experiment. This high power duced nearly identical results. Both methods estimated was due to the large design, the inclusion of multiple values for 0 ( ˆ0,S ϭ 0.30 and ˆ0,B ϭ 0.28) that are QTL as cofactors in the statistical model, and the appliapproximately in agreement with the density of the p ccation of the false discovery rate, which accounted for values for p c Ͼ 0.6 in the histogram in Figure 2. This the proportion of truly null hypotheses. As the analysis density reflects the proportion of true null hypotheses revealed no significant QTL-by-environment interaction (Storey and Tibshirani 2003). In his software manual the QTL found can be taken as candidates for the Storey wrote that the FDR S method often works better marker-assisted selection programs currently implethan the FDR B method, but that it can backfire for a mented in both countries, in Germany (Bennewitz et al. small number of tests or in pathological situations. In 2003b) as well as in France [START_REF] Storey | Implementation of marker assisted selection consistency of false discovery rates: a unified approach[END_REF]. general, the authors recommended the use of the FDR B method if the number of tests is small. Similar to these

Figure 1 .

 1 Figure1.-Plot of the q-values against the number of hypotheses declared as significant. The hypotheses were ordered decreasing by their P-values. , q BH ; ᮀ, q S ; and ᭡, q B .

Figure 2 .

 2 Figure 2.-Histogram of chromosomewise error probabilities (n ϭ 54, interval size 0.1). __ indicates the expectation of the density under the assumption of all null hypotheses being true. --and ••• reflect the estimates of the proportion of the true null hypotheses estimated with FDR S and FDR B , respectively.

TABLE 1

 1

  . Recently, de Koning

			across Germany and France
	et al. (2001) proposed a recursive strategy for interval
	mapping in outbred populations that allows the inclu-		Progeny tested in	Genotyped in
	sion of multiple QTL. Their method included the QTL transition probabilities at the estimated position of QTL	Family	Germany	France	Germany	France	Total
	that are outside the marker interval under consideration	F01	42	102	42	102	144
	in the model rather than single markers outside this interval. In chromosome or genome scans multiple tests that are not independent are usually performed. It is a com-mon practice to control the chromosome-or genome-	F02 F03 F04 F05 F06 F07	50 127 126 126 33 23	83 95 168 236 42 27	50 127 126 126 75 50	83 95 168 236	133 222 294 362 75 50
	wise type I error rate of these multiple tests by estimating	F08	19	74	93		93
	the corresponding threshold levels by a permutation	F09	20	26	46		46
	test (Churchill and Doerge 1994). These threshold levels are trait specific and it is to date still not clear how to estimate threshold levels across traits by use of the permutation test. Further, the application of a stringent significance criterion reduces the type I error	F10 F11 F12 F13 F14 F15	25 57 58 28 31 29	62 10 98	87 67 156 28 31 29		87 67 156 28 31 29
	rates, but increases the type II error rate defined as	F16	42		42		42
	the probability that a present QTL will be missed and	F17	60		60		60
		F18		58	58		58

consequently reduces the power of the experiment. The application of the false discovery rate (FDR; Benjamini In total 1977 sires are included. and Hochberg 1995) instead of the type I error rate is a useful statistical tool to overcome these two disadvansires were included in the French granddaughter design tages. Roughly speaking, the FDR is defined as the ex-

(Boichard et al. 2003)

. In Table 1 the pedigree structure is pected proportion of false positives among all rejected summarized. The first 5 families listed in Table 1 had already null hypotheses. It was introduced by Weller et al.

  ). significant nonsignificant Total Unfortunately 0 is unknown in practice. Therefore Benjamini and Hochberg (1995) proposed the adoption of 0 ϭ Ϫ T m Ϫ m 0 mini and Hochberg 2000; Mosig et al. 2001; Fernando et al. true 2002). Storey and Tibshirani (2003) and Storey et al.

	Null hypothesis true Alternative hypothesis mϪ m 0 Total F m 0 Ϫ F T S mϪ S	m 0 m	1. There are some intuitive proposals for estimating 0 (Benja-

TABLE 3 Results from interval mapping (IM) without cofactors (simple IM) and from interval mapping with cofactors (IM with cofactors)

 3

			Simple IM				IM with cofactors
	BTA	Trait a	F b	p c c	F	p c	p g d	Pos e	C.I.95 f	␣ g	No. GS h
	05	MY	2.4	0.0050	2.8	0.0004	0.011	127	[19, 144]	0.69	3
		PY	1.7	0.1630	2.0	0.0340	0.646	117	[21, 148]	0.75	4
		FP	2.1	0.8100	2.5	0.0130	0.334	127	[17, 147]	0.84	3
		PP	2.2	0.0747	2.6	0.0071	0.192	142	[39, 151]	0.69	5
	06	SCS	1.8	0.1028	2.0	0.0440	0.744	99	[16, 135]	0.45	4
		PP	2.0	0.1278	2.6	0.0082	0.219	52	[21, 85]	0.89	2
	14 i	SCS	2.3	0.0074	2.7	0.0012	0.035	149	[126, 149]	0.71	2
	18	FY	2.0	0.0096	2.1	0.0054	0.150	139	[77, 146]	0.74	2
	19	FY	2.2	0.0143	2.5	0.0007	0.021	57	[20, 79]	0.66	3
		FP	2.8	0.0039	2.5	0.0134	0.333	77	[5, 115]	0.62	3
		PP	2.1	0.0762	2.3	0.0275	0.567	77	[7, 126]	0.52	3
		SCS	2.1	0.0236	2.2	0.0095	0.249	58	[7, 130]	0.45	4
	20	FP	1.4	0.5629	2.4	0.0135	0.335	61	[40, 82]	0.90	3
		PP	1.7	0.1889	2.4	0.0055	0.153	45	[15, 73]	0.66	4
	23	MY	2.1	0.0186	2.0	0.0266	0.555	64	[10, 76]	0.50	3
		FY	1.8	0.0379	2.4	0.0004	0.012	64	[1, 68]	0.49	3
		PY	1.8	0.0644	1.8	0.0445	0.745	64	[7, 76]	0.36	3
	26	MY	2.2	0.0045	2.4	0.0015	0.044	30	[7, 37]	0.65	2
		FY	2.4	0.0005	2.6	0.0001	0.003	30	[15, 38]	0.91	4
		PY	2.1	0.0117	2.2	0.0026	0.075	31	[12, 43]	0.67	2

a Trait abbreviations are MY, milk yield; FY, fat yield; PY, protein yield; FP, fat percentage; PP, protein percentage; and SCS, somatic cell score.

b F-test statistic. c Chromosomewise error probability. d Genomewise error probability.

TABLE 4

 4

We thank Mike E. Goddard for helpful discussions and for carefully reading this manuscript. This article has benefited from the critical two methods, Mosig et al. (2001) presented an iterative


Michaël Deschamps
email: michael.deschamps@cnrs-orleans.fr

Dominique Massiot

Two and Three-dimensional Through-Bond Homonuclear-Heteronuclear Correlation Experiments for Quadrupolar Nuclei in Solid-State NMR Applied to 27 Al-O-31 P-O-27 Al networks

Keywords: 3D, J-coupling, quadrupolar nuclei, aluminium, phosphate

    Two and Three-dimensional Through-Bond Homonuclear-Heteronuclear Correlation Experiments for Quadrupolar Nuclei in Solid-State NMR Applied to 27Al-O-31P-O-27Al networks M. Deschamps, D.

Introduction

3D J-mediated correlation experiments have been introduced a long time ago in liquid state NMR experiments [START_REF] Cavanagh | Protein NMR spectroscopy: principles and practice[END_REF]. They characterize triplets of atoms which are connected together, and, as such, are very powerful to characterize the network of chemical bonds in a solution or in a solid-state sample as well. High resolution solid-state experiments can make efficient use of the scalar coupling (J) in MAS experiments on organic and inorganic solids [START_REF] Sakellariou | Through-Bond Experiments in Solids: The Encyclopedia of NMR[END_REF][START_REF] Massiot | Chemical bonding differences evidenced from J coupling in solid state NMR experiments involving quadrupolar nuclei[END_REF], despite the heterogeneous line broadening mechanisms that usually mask the direct spectral expression of the scalar couplings, although benefiting from the increased sensitivity and resolution obtained at high MAS spinning rates and high magnetic fields [START_REF] Gan | Seeking Higher Resolution and Sensitivity for NMR of Quadrupolar nuclei at ultrahigh magnetic fields[END_REF]. Dynamic Angle Spinning (DAS) also enhance the resolution for quadrupolar nuclei [START_REF] Grandinetti | Dynamic-Angle Spinning and Applications[END_REF]. Among mesoporous materials made from AlPO 4 , AlPO 4 -14 has already been studied by X-ray crystallography [START_REF] Helliwell | Structure determination from small crystals of two aluminophosphates CrAPO-14 and SAPO-43[END_REF][START_REF] Fyfe | 1D and 2D solid state NMR investigations of the framework structure of As-synthesized AlPO4-14[END_REF] and NMR [START_REF] Fyfe | 1D and 2D solid state NMR investigations of the framework structure of As-synthesized AlPO4-14[END_REF](8)[START_REF] Wiench | Probing through bond connectivities with MQMAS NMR[END_REF] and will be used for demonstration. Its structure is well known and similar to those of zeolithes, however, the framework is entirely built up with AlO x (x = 4, 5 or 6) and PO 4 polyhedrons. Four, five and six-fold coordination states are observed for the Al atoms and Al are connected together by either bridging oxygen atoms or O-P-O chains (Figure 1). Four different sites can be distinguished for Al and also for P atoms [START_REF] Fyfe | 1D and 2D solid state NMR investigations of the framework structure of As-synthesized AlPO4-14[END_REF] In contrast with liquid state experiments in which the short relaxation times observed for quadrupolar nuclei often precludes their proper manipulation by NMR, several 2D dipolar mediated or J-mediated homonuclear and heteronuclear correlation methods have been introduced for solid-state experiments involving one or two quadrupolar nuclei like 27 Al (S=5/2) (7-14), 71 Ga (S=3/2) (15) or 27 Al and 17 O (S=5/2) [START_REF] Iuga | NMR Heteronuclear Correlation between Quadrupolar Nuclei in Solids[END_REF]. The efficiency of 27 Al / 31 P (S=1/2) experiments benefits from their 100% natural abundance and their good receptivity, with an increased sensitivity at high field thanks to a lower second order quadrupolar broadening for the 27 Al spin [START_REF] Gan | Seeking Higher Resolution and Sensitivity for NMR of Quadrupolar nuclei at ultrahigh magnetic fields[END_REF]. Moreover, the magnetization of the 27 Al central transition (which is only affected by the second-order quadrupolar broadening) can be efficiently enhanced by transferring the populations from the satellite transitions to the central transition using methods like DFS or RAPT [START_REF] Kentgens | Advantages of double frequency sweep in static, MAS and MQMAS NMR of spin I=3/2 nuclei[END_REF][START_REF] Kwak | Enhancing Sensitivity of Quadrupolar Nuclei in Solid-State NMR with Multiple Rotor Assisted Population Transfers[END_REF].

In solid-state, the 2D or 3D HHSQC experiments allow to characterize structural motifs involving three different spins (2 27 Al and 1 31 P in a Al-O-P-O-Al linkage) in a 3 dimensional experiment, correlating the 1D spectra of the three bonded nuclei (2× 27 Al and 1× 31 P). The 2D HHSQC only gives the homonuclear { 27 Al, 27 Al} two dimensional correlation map selectively for the Al-O-P-O-Al motifs. Because the scalar coupling J 2 (X-O-Y) is used, the 2D/3D H-HSQC spectra allow identification of molecular entities extending over 4 chemical bonds and typically over 0.4 to 0.5 nm in size. Despite its low sensitivity, this experiment provides a unique way for obtaining unambiguously this type of information.

Theory

The 2D/3D homonuclear/heteronuclear correlation experiment (3D H-HSQC) makes use of the sole scalar coupling (J) to transfer the magnetization (see Figure 2). This results in a robust and rather simple pulse sequence that only involves manipulation of the spin system with soft π/2 and π pulses, which ensures manipulation of the 27 Al central transition as a fictitious spin 1/2. The transfer occurs from the first aluminium atom to the neighbouring phosphorus and then to a second aluminium through the isotropic part J 2 (Al-O-P) of the scalar coupling in an Al-O-P-O-Al group and provides experimental evidence for through-bond correlations between the Al 1 , P and 1), ( 2) and (3) correspond to n'=2, 3 and 4, respectively, with an initial density operator σ(τ 2 =0) = 2Al 1 z,CT P y :

( ) [ ] ( ) [ ] 2 2 z 2 1 z 4 cos 1 2 / 1 l and 4 cos 1 2 / 1 l τ π τ π J P A I J P A I y y + - - = + - = (1) ( ) [ ] ( ) [ ] 2 2 2 z 2 2 1 z 2 cos 6 cos 4 / 1 l and 2 cos 3 6 cos 4 / 1 l τ π τ π τ π τ π J J P A I J J P A I y y - = + = (2) ( ) [ ] ( ) [ ] 2
+ - = + + = (3)
The two expressions in Eq. ( 3) are plotted in Figure 3 as a function of τ 2 in 1/J units. The 31 P chemical shift evolves during the delay t 2 , and the 31 P scalar couplings to the neighbouring 27 Al spins evolve in a constant time manner during 2τ 2 (see Figure 2).

Discussion

For the 2D experiment, the dataset is Fourier transformed against the two evolution times and the resulting spectrum is shown in 1) which has been established previously (8,[START_REF] Wiench | Probing through bond connectivities with MQMAS NMR[END_REF]. A similar conclusion can be inferred from Figure 5D, where P 1 , P 3 and P 4 connect to both the Al 1 and Al 3 atoms.

Experimental

Weak (selective, i.e. ν 1 ≈ 4 kHz) pulses have been used on the 27 Al central transition (described as a fictitious spin 1/2), which is prepared by enhancement from the satellite transitions by DFS (17) transfer using a 2 ms pulse withν 1 ≈ 8 kHz (RAPT (18) can also be used).

Two Z-filters (300 rotor periods or 21 ms) have been used in the pulse sequence to simplify the phase cycle and suppress unwanted contributions of undetermined origin which were observed without the Z-filter. Unfortunately such long Z-filter delays are approaching the order of the T 1 value, and some signal is lost during the Z-filter. A cleaning B 0 field gradient instead would have the desired effect but was unfortunately unavailable on our MAS probes.

The times τ and τ 2 were set to 40 rotor periods (2.86 ms) and 120 rotor periods (8.57 ms) respectively, and the J 2 ( 27 Al, 31 P) values (unobservable in this case) are estimated to be in the 10-30Hz range. To optimize τ and τ 2 , τ 2 was first set to zero and τ has been optimized to get the largest positive signal. Then τ 2 is increased and optimized to obtain the largest negative signal, as shown in Figure 3. Quantification of the peak intensities requires similar transfer efficiencies:

Different T 2 and J-coupling values would lead to different transfer efficiencies, precluding a proper interpretation of the peak intensities.

The MAS speed was set to 14 kHz and a magnetic field B 0 =17. wide bore Bruker spectrometer. The optimum recycling delay, which depends upon the relaxation behaviour of the central and satellite transitions (because DFS is used), was determined by comparing the signal obtained in the permanent regime for various recycling delays, and 250 ms was considered as optimal, as it allowed a 95% recovery.

Indirect quadrature detection was achieved using the States method, whereas DQD was used during acquisition. A spectral width of 14 kHz was used in the three dimensions, corresponding exactly to the MAS spinning speed, in order to sum any remaining spinning sidebands [START_REF] Massiot | Sensitivity and lineshape improvements of MQ-MAS by rotor synchronized data acquisition[END_REF]. The data have been processed with nmrPipe (20), the spectrum is phased in order to have negative diagonal peaks and positive cross-peaks, and an exponential line broadening of 200 Hz has been applied in the three dimensions (the natural linewidth in the 31 P dimension -without 1 H and 27 Al decoupling-is about 200 Hz, in the 27 Al dimension, second order quadrupolar broadening is larger than the exponential line broadening factor).

Conclusions

2D and 3D experiments correlating Al 1 , P, Al 2 triplets in Al 1 -O-P-O-Al 2 motives are feasible and confirm the connectivity graph of the Al/P network of a complex material such as AlPO 4 -14. The significantly longer experimental time required for the 3D experiment (nearly 6 ½ days), as compared to the 2D Al, Al H-HSQC (1 day), is the price to be paid for the complete edition of connected Al x -O-P-O-Al y motives. However, future hardware improvements such as higher magnetic fields, faster MAS rate, the use of magnetic field gradients which will supersede Zfilters (especially in this experiment) and the eventual inclusion of MQMAS and/or 1 H decoupling and 27 Al satellite transitions refocusing during t 1 might increase the S/N ratio and the available resolution. The principle of this 3D H-HSQC could be extended to any system of three or more coupled spins and applied to many inorganic materials. It provides a detailed map of the coupled network, and a full characterization of chemically bonded molecular motives. When this information remains of little use for crystalline compounds of known structure, recent results

show that it becomes decisive information for the characterization of structural motifs at the nanometer scale in amorphous solids or glasses [START_REF] Iuga | NMR Heteronuclear Correlation between Quadrupolar Nuclei in Solids[END_REF]. Figure 4 H-HSQC spectrum of an AlPO 4 -14 sample, featuring the Al-Al pairs giving rise to negative cross-peaks (in black) and the positive diagonal signals (in red). The spectrum was obtained in the experimental conditions mentioned in the Figure 2 caption. The pulse sequence was run using the following parameters: the recycling delay d 1 was set to 250ms, the τ period for the INEPT transfer was set to 40 rotor periods τ r (optimized for the best overall signal-to-noise ratio, as the different 27 Al-31 P pairs feature distinct scalar couplings and dephasing times), the τ 2 period was set to 120τ r (in order to maximize the negative cross peaks intensities), and a Z-filter of 300τ r was used on a 750MHz wide bore Bruker spectrometer. 16 dummy scans and 1152 scans were acquired for each of the 128 t 1 increments and exponential broadening (LB=100Hz) was applied in both dimensions. Table 1 Al,P connectivity table for AlPO 4 -14, as established in refs. (7, 8).

  Al 1 is five fold coordinated trigonal-bipyramidal (bonded to 4 -O-P-(O-Al) 4 and 1 -O-Al), Al 2 and Al 3 are four fold coordinated tetrahedra (bonded to 4 -O-P-(O-Al) 4 ) and Al 4 is a six fold coordinated octahedra (bonded to 4 -O-P-(O-Al) 4 and 2 -O-Al).

Al 2

 2 nuclei in Al 1 -O-P-O-Al 2 motives with spectral identification of the resonances of these three nuclei in the different projections. Because the experiment uses the scalar coupling it does not require any reintroduction of interaction averaged out by MAS and it undergoes no limitation with the spinning rate which is usually required to achieve better resolution. It only requires long enough transverse relaxation times T 2 ' (typically 15 ms in our case) and a sufficient sensitivity for the acquisition of the 2D or 3D datasets within a reasonable experimental time. The experiment consists in two INEPT transfers separated by a constant time period achieving the magnetization transfer and containing the 31 P encoding time t 2 . The first INEPT transfer leads to the creation of a density operator proportional to 2Al 1 z,CT P y (where the aluminium central transition is described as a fictitious spin 1/2) after an evolution under the scalar coupling during a delay 2τ (1/4J). The second constant time echo period (2 τ 2 ) allows for the conversion of Al 1 z,CT P y into Al 2 z,CT P y , which gives rise to a negative cross-peak in the 2D or 3D H-HSQC spectrum. The positive diagonal signal corresponds to the remaining Al 1 z,CT P y term. This transfer depends upon n', the number of 27 Al spins coupled to the 31 P nucleus. The functions describing the intensities of the diagonal signal (Al 1 z,CT P y ) and of the cross-peaks (stemming from the conversion of Al 1 z,CT P y into Al 2 z,CT P y ) are given below. Equations (

Figure 4 .

 4 The different spectral contributions are indicated on the spectrum. The third evolution time in the 3D experiment allows to separate the (Al x , Al y ) crosspeaks with respect to the 31 P chemical shift in the Al x -O-P z -O-Al y motif. In the 3D experiment, the resulting three dimensional hypercomplex dataset is Fourier transformed against the three evolution times and Figure5B-D shows one 27 Al/ 27 Al plane (with the x and y axis corresponding to the direct and indirect 27 Al dimensions) extracted from the 3D matrix (as shown in Figure5A) and corresponding to the phosphorus chemical shift of the P 2 site δ( 31 P) = -8.03 ppm. Hence, the extracted 2D spectrum will eventually feature the residual negative diagonal peaks of some Al x -O-P 2 group and the positive cross-peaks of any Al x -O-P 2 -O-Al y groups. Three different strips along the 31 P axis, perpendicular to the 27 Al/ 27 Al plane and intersecting with it in the cross (indicating an Al x ,Al y peak), are superimposed on the 2D spectrum. If Al x and Al y are also connected via O-P n -O, the P n peak will appear in the superimposed strip. In (B), the cross is located on the negative (hence diagonal) Al 1 , Al 1 peak in the plane corresponding to P 2 , indicating that a diagonal peak is observed for the triplet Al 1 , P 2 , Al 1 , and the corresponding strip (intersecting with the Al, Al plane in (δ(Al 1 ),δ(Al 1 )) as indicated by the cross) along the 31 P axis features three additional peaks corresponding to P 1 , P 3 and P 4 indicating that similar Al 1 , Al 1 diagonal correlations are observed in the planes corresponding to the chemical shifts of P

Figure 1

 1 Figure 1 Five schemes of the molecular motives involved are shown below in the 27 Al/ 31 P case, featuring: a) Al IV , b) Al V , c) Al VI , d) Al-O-P-O-Al and e) Al-O-Al.

Figure 2

 2 Figure 2Pulse sequence for: A) the 2D (t 1 incremented) and B) the 3D (t 1 and t 2 incremented) H-

Figure 3

 3 Figure 3 Plots of the intensities of the diagonal (2Al 1 z,CT P y ) and cross-peak (2Al 2 z,CT P y ) and their sum (overall signal) as a function of τ 2 in 1/J units. The optimum transfer is obtained for the

Figure 5 A

 5 Figure 5 A 2D plane (as shown in (A)) is extracted from the 3D 27 Al, 31 P, 27 Al H-HSQC, corresponding to the phosphorus chemical shift δ( 31 P)= -8.03 ppm (corresponding to P 2 ) and is displayed three times in (B), (C) and (D). Three different strips or slices along the 31 P axis have

Figure 1 Figure 2 Figure 3

 123 Figure 1

Figure 4

 4 Figure 4

Figure 5

 5 Figure 5

Michael

  Deschamps is a Maître de Conférences at the University of Orléans in France. He did his undergraduate work at the Ecole Normale Supérieure in Paris, where he was trained in NMR during his Ph.D. under the supervision of Prof. Geoffrey Bodenhausen. After a post-doctoral EMBO fellowship at the University of Oxford, where he worked with Prof. Iain Campbell and Dr. Jonathan Boyd, he joined the University of Orléans, working on solid-state NMR in the Centre de Recherche sur les Matériaux à Haute Température (CRMHT).Dominique Massiot is a CNRS Directeur de Recherche at the CRMHT in Orléans. He did his Ph.D. in Geochemistry at the Paris VII University in 1983, and has been working ever since at the CRMHT using solid-state and high temperature NMR to tackle a wide variety of problems in inorganic chemistry. He is the author of "dmfit", a widely used software for the interpretation of complex solid-state NMR spectra.

Table

  1 , P 3 and P 4 respectively. Hence Al 1 is never connected to another Al 1 via an O-P-O group, because only negative diagonal peaks have been observed, and Al 1 is bonded to P 1 , P 2 , P 3 and P 4 . In (C), the cross is located on the Al 2 , Al 4 cross-peak, indicating that Al 2 and Al 4 are connected via an O-P 2 -O linkage (because the extracted 2D spectrum corresponds to P 2 ). Two peaks are observed in the strip for the chemical shifts of P 2 and P 4 , hence the cross-peaks between Al 2 and Al 4 are observed in the two following groups: Al 2 -O-P 2 -O-Al 4 (observable in the plane δ(P 2 ) which is represented in the figure) and Al 2 -O-P 4 -O-Al 4 (observable in the plane δ(P 4 ), not represented here). This complies with the known connectivity table (

  6 T (750 Mhz 1 H Larmor frequency with 27 Al and 31 P Larmor frequency equal to 194.5 and 303.7 MHz respectively) was

	used. 1024 scans, 16 dummy scans, 40 increments in F1 and F2, 512 points in F3 and a recycling
	delay of 250 ms have been used for an overall experimental time of 6 days and 13h on a 750MHz

Acknowledgments

We acknowledge financial support from CNRS UPR4212, FR2950, MIAT and Région Centre and we are grateful to Pr. P.J. Grandinetti, Drs. F. Fayon and V. Montouillout for stimulating discussions.


Julie Groleau

Matthieu Chabanas

Christophe Marecaux

Natacha Payrard

Brice Segaud

Michel Rochette

Pascal Perrier

Yohan Payan

A biomechanical model of the face including muscles for the prediction of deformations during speech production

Keywords: Face models, Muscle modeling, Lip/teeth interaction

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Many biomechanical models of the human face have been proposed in the literature. They were generally developed either in the context of computer graphics animation [START_REF] Lucero | A model of facial biomechanics for speech production[END_REF][START_REF] Sifakis | Simulating Speech with a Physics-Based Facial Muscle Model[END_REF], or of computer-aided maxillofacial surgery [START_REF] Gladilin | Towards a Realistic Simulation of Individual Facial Mimics[END_REF][START_REF] Chabanas | Patient specific Finite Element model of the face soft tissue for computer-assisted maxillofacial surgery[END_REF] or of speech production studies [START_REF] Gomi | A physiologically based model of perioral dynamics for various lip deformations in speech articulation[END_REF]. Most of them propose to model the face with a volumetric mesh defined by an external (the "visible" part of the face) and an internal surface (the part in contact with the skull), with some nodes or layers in between. Mechanics of the tissues (epidermis, dermis, hypodermis, fat and muscles) is then modelled through the relation between displacements and forces of mass points, or through strain/stress relations in the case of finite element models. These studies have raised a number of important issues:

(1) how to model muscles fibres and their action on the 3D mesh; (2) how to account for the subject/patient specific muscular morphology (in terms of fibres insertions and interweaving); (3) how to control the large number of muscles in order to produce a given speech articulation or facial mimics. Both latter points have already been addressed and discussed by our group, respectively for computer-aided craniofacial surgery [START_REF] Chabanas | Patient specific Finite Element model of the face soft tissue for computer-assisted maxillofacial surgery[END_REF] and through a motor control model of tongue muscles activation for speech production [START_REF] Buchaillard | A 3D biomechanical vocal tract model to study speech production control: How to take into account the gravity[END_REF].

This paper deals with the first issue and proposes a method to define from Computer Tomography (CT) images a subset of muscles fibres within a 3D mesh of the face. Contacts between lips and teeth are also handled. First results of facial mimics' simulations are presented.

A BIOMECHANICAL MODEL OF THE FACE INCLUDING MUSCLES FOR THE PREDICTION OF DEFORMATIONS DURING SPEECH PRODUCTION

Julie Groleau 1 , Matthieu Chabanas The starting point of this modeling work is the 3D Finite Element model of the face soft tissues, built out of CT scan of a single patient, which was originally proposed in [START_REF] Chabanas | Patient specific Finite Element model of the face soft tissue for computer-assisted maxillofacial surgery[END_REF]. It relies on a volumetric mesh consisting of hexahedrons and wedges elements (Fig. 1, left). The displacements of several nodes located on the internal surface of the facial mesh are constrained in order to represent attachments of the facial tissues on the skull.

While biological soft tissues are known to behave nonlinearly [START_REF] Gérard | Non linear elastic properties of the lingual and facial tissues assessed by indentation technique. Application to the biomechanics of speech production[END_REF], they were first represented by a homogeneous, isotropic, linear material. This hypothesis was retained in a first stage in order to focus on muscles modelling and contact management, before ongoing with more realistic modelling. Simulations were computed using the ANSYS TM v11 finite element software.

The first part of our study has consisted in building the muscles involved in facial mimics' generation. In order to ensure anatomical and physical reliability, muscles courses and insertions were directly defined from medical images and anatomical charts (Fig. 1, middle), with the help of a maxillofacial surgeon. The locations of points describing the muscle fibres were measured in the different scan slices. These points were then integrated into the mesh to model muscle insertions. They were linked with piece-wise uniaxial tension cable elements to model muscle fibres (Fig. 1, right). The Orbicularis Oris muscle was designed slightly differently: it is represented by two ellipsoid cable elements centred on the mouth opening, and representing the marginal and the peripheral parts of the muscle.

The cable elements based approach allows integrating muscles into the model independently of the mesh itself. Therefore, the mesh can be easily refined or modified, without requiring any change in the muscles structure definition. The fibers cable elements are controlled in tension by their cross section area, their initial strain and an activation parameter. They generate forces that are applied to the soft tissues mesh thanks to the notion of dependencies. In other words, muscular fibres extremities are linked with the facets of the surrounding mesh elements. When a muscle is activated, the corresponding cable elements exert forces on the mesh elements and induce, then, soft tissue deformations.

The second part of our study concerns lips-teeth and upper-lip/lower-lip contacts, which are of primary importance in lips movements and deformations. Teeth are materialized in the model by surfaces extracted from the CT data and interpolated with Spline functions. ANSYS contact elements, which provide collision detection and sliding reaction, are used to mesh lips and teeth surfaces.

a. Simulating smiling lips.

Figure 2 presents the mesh deformations (in colours/grey scale) and the final face shape when the Zygomaticus Major, the Risorius and the Levator Labii Superioris are activated simultaneously. The bottom panels show the results when lip/lip and lip/teeth contacts are taken into account, while the top panels show the results obtained without contact. In presence of contacts, the mimic seems more realistic. This can be particularly well assessed on the side views (right panel), where, in the absence of contacts, an interpenetration of the lips can be noticed. On the contrary,, lips are slightly opened in the other condition, which is in agreement with data on smiling.

b. Orbicularis activation and lips protrusion.

An important issue for speech production concerns the control of protruded lips, such as in the production of /u/ or /y/ in French. It has been suggested in the literature [START_REF] Honda | Physiology of the lips and modeling of lip. gestures[END_REF] that the interaction between an activation of the Orbicularis Oris and lip/lip and lip/teeth contacts could be responsible for this particular lip gesture. This hypothesis was tested with our face model. Two types of simulations were then run with the activation of the Orbicularis Oris, without (Figure 3 top panels) and with (Figure 3 low panels) handling contacts. In the absence of contacts, a rounding of the lips is observed (top panel, left). Rounding is classically associated with protrusion. However, the side view (top panel, right) shows a strong retraction of the lips, which is at the opposite of a protrusion. Including the contact limits the retraction, but it does not generate any protrusion or rounding. The absence of protrusion can be explained by the fact that our model does not separately control the marginal and the peripheral parts of the Orbicularis Oris. Honda et al. [START_REF] Honda | Physiology of the lips and modeling of lip. gestures[END_REF] found namely an EMG activation during protrusion only in the peripheral part of the muscle. Gomi et al. [START_REF] Gomi | A physiologically based model of perioral dynamics for various lip deformations in speech articulation[END_REF] confirmed this finding with their biomechanical lip model. However, the absence of rounding is in agreement with Gomi et al.'s [START_REF] Gomi | A physiologically based model of perioral dynamics for various lip deformations in speech articulation[END_REF] statement (p.130) who suggested that not only the Orbicularis Oris, but also "additional muscles combinations (jaw opening, peripheralis or other perioral muscles) would be required to form rounded lips." "

IV. DISCUSSION AND CONCLUSION

The face model presented in this paper integrates an original representation of muscle fibres and muscle force generation in a 3D mesh based on piece-wise uniaxial tension cable elements. It also models contacts between upper lip and lower lip and between lips and teeth in a realistic way.

Simulations of smiling lips show that the proposed muscle representation is adapted to the generation of lip deformations that are realistic both in amplitudes and in directions. This is an important result since this representation can be implemented independently of the mesh. It will then facilitate the generation of speaker specific mesh using mesh matching algorithm [START_REF] Couteau | The Mesh-Matching algorithm: an automatic 3D mesh generator for finite element structures[END_REF]. It will also increase the efficiency of such a modelling approach to study the impact of face surgery on smiling and on mimics in general. Indeed, the geometrical structure of the mesh can easily be modified to account for different kinds of surgeries, without inducing a careful, difficult and long redefinition of all muscles fibres in the mesh.

Our results show also the importance of contacts modelling, both for lips/teeth interactions, but also for upper lip/lower lip interaction. It is important, not only because it prevents for unrealistic interpenetrations, but also because it allows sliding movements that constrain and guide the movement. This is certainly a phenomenon underlying complex lip shaping such as protrusion and rounding.

On the other hand, the simulations of the consequences of the activation of the Orbicularis Oris, for which no distinction is made in our model between the marginal and the peripheral parts, do not generate rounding. These results are in contradiction with the simulations carried out by Gomi et al's [START_REF] Gomi | A physiologically based model of perioral dynamics for various lip deformations in speech articulation[END_REF] who did make this distinction. This shows that collecting accurate neurophysiological and anatomical data is a major challenge to test hypotheses about the control of speech gestures, once realistic biomechanical models are available.

.

FiberFigure 1

 1 Figure 1 Left: the finite element mesh of the face soft tissues. Middle: interactive segmentation of muscles fibers on CT data. Right: location of eleven muscles involved in facial mimics on the left side of the skull.

Figure 2

 2 Figure 2 Simulation of smiling lips without (top panels) and with (bottom panels) handling of lips-teeth and upper-lower lips contacts. Displacements are in mm.

Figure 3

 3 Figure 3 Simulation of the Orbicularis Oris contraction, without (top panels) and with (bottom panels) handling of lips-teeth and upper-lip/lower-lip contacts.


Nicolas Normand

Convex Structuring Element Decomposition for Single Scan Binary Mathematical Morphology

Keywords:

This paper presents a structuring element decomposition method and a corresponding morphological erosion algorithm able to compute the binary erosion of an image using a single regular pass whatever the size of the convex structuring element. Similarly to classical dilation-based methods, the proposed decomposition is iterative and builds a growing set of structuring elements. The novelty consists in using the set union instead of the Minkowski sum as the elementary structuring element construction operator. At each step of the construction, already-built elements can be joined together in any combination of translations and set unions. There is no restrictions on the shape of the structuring element that can be built. Arbitrary shape decompositions can be obtained with existing genetic algorithms with an homogeneous construction method. This paper, however, addresses the problem of convex shape decomposition with a deterministic method.

Introduction

Mathematical morphology operators are time-consuming with large structuring elements and brute force algorithms. In the past, several methods have been described to reduce the cost of these operators. Two main approaches exist, the first one uses a decomposition of a large structuring element into a set of smaller ones. The result is obtained by a series of operations with small structuring elements. The overall cost is then directly connected to the number of operations and depends on the size of the initial structuring element. The second one consists in binarizing a distance map, which can be computed in a fixed number of image scans. It requires the structuring element to be expressed as a distance disk but then the computational cost is constant whatever the size of the structuring element.

The method proposed here is based on a new generalized distance transform (GDT). The algorithm is quite similar to local distance propagation algorithms but in our case, distance increments are not constant over disks size, which allows for much more flexibility in the disk construction. As an example, we describe an algorithm to decompose any convex 2D polygon in a series of pseudo-distance disks for single scan distance map computation.

The overall cost of the mathematical morphology operators derived from this GDT is constant with the size of the structuring element. Moreover, they can be used in a pipeline fashion and have very low memory requirements. In section 2 existing structuring element decomposition methods will be recalled.

Distances and Structuring Elements

Mathematical Morphology Operators

Let A and B be two sets of points in the discrete grid E with origin O, the neutral element for the symetry in E (p symetric element is denoted as p). The erosion of A by the structuring element B is defined as:

A B = {p| (B) p ⊆ A} (1)
where (B) p is B translated by p:

(B) p = {x + p|x ∈ B}.
The erosion dual operator, the dilation can be defined as:

A ⊕ B = (A c B) c = {p| (B) p ∩ A = ∅}. ( 2
)
The notation ⊕ denotes the Minkowski sum of two sets i.e. the set of sums of two elements, one taken from the first set the other from the second one. These basic operators lead to a great variety of image transformations [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. However, the algorithm directly derived from the fundamental definition given in (eq. 1) is not efficient for a large structuring element B, as it requires the exploration of all translated points of B for each point of the image.

Distance Map

The distance transform associates to any point x the smallest distance to a point outside the set X:

d X (x) = min y ∈X (d(x, y)) (3)
The distance map is linked to the mathematical morphology erosion by the fact that the set of points whose distance map values are at least r is the eroded set of X by D(r):

A Ď(r) = {p|d X (p) ≥ r} (4)
with D(r) = {p|d(O, p) < r}. Eroding a shape from a distance map consists in thresholding distance values, so the erosion cost depends only on the cost of the distance transform. Since algorithms exist to compute a distance map in a fixed number of scans [START_REF] Rosenfeld | Distances functions on digital pictures[END_REF][START_REF] Yokoi | On generalized distance transformation of digitized pictures[END_REF][START_REF] Borgefors | Distance transformations in digital images[END_REF], they can be used to erode with a constant cost whatever the size of the structuring element. For usual distances, each disk is constructed by the dilation of the previous disk with a basic structuring element as illustrated in fig. 1.a. In a sequential distance map computation, the symetric neighborhood is divided in two halves which are passed over the image once, in reverse order scans [START_REF] Rosenfeld | Distances functions on digital pictures[END_REF].

By moving the morphological center of the disks to the last scanned pixel, some one-pass algorithms can be obtained [START_REF] Wang | An algorithm for a generalized distance transformation based on minkowski operations[END_REF]. We can not refer to these disks as distance disks since the symetry property of distances is not verified anymore. Hence the transform is called a generalized distance transform (GDT). However, there is a strong constraint on the shape of the disks since the basic structuring element is unique for the whole set of disks, so this method will only apply to very specific structuring elements.

By mixing different structuring elements, distances like the octagonal distance add some variability: each disk can be built from the previous by a different structuring element (fig. 1.b). Any shape that is decomposed in a series of dilations can be constructed with this method. However, distance transfom algorithms are only known for some specific cases (for instance when two building structuring elements are used periodically [START_REF] Wang | An algorithm for a generalized distance transformation based on minkowski operations[END_REF][START_REF] Wang | Some sequential algorithms for a generalized distance transformation based on minkowski operations[END_REF]).

On the other side, another way of mixing different neighborhoods is used in chamfer distances (fig. 1.c): each disk is built from different-size disks according to local distances described in a neighborood mask.

Structuring Element Decomposition

The structuring element decomposition methods rely on the fact that a series of erosions with a set of structuring elements is equivalent to a single erosion with the Minkowski sum of the structuring elements:

(A B) Č = A ( B ⊕ Č). (5)
Decomposition methods generally use a series of basic structuring elements computable by specific hardware machines in one clock cycle. The shape of the building structuring elements depends on the hardware platform and convex polygon decomposition algorithms were presented for instance for linear shaped building structuring elements [START_REF] Gong | On decomposition of structure element for mathematical morphology[END_REF] and for 4 and 8-neighborhood parallel machines [START_REF] Xu | Decomposition of convex polygonal morphological structuring elements into neighborhood subsets[END_REF]. These decompositions lead to optimal morphological operator implementations for parallel or pipeline architectures, but conversely to distance-based methods, the computational complexity depends on the size of the structuring element. Since some convex polygonal structuring elements can not be decomposed by Minkowski sums, an extra final set union can be needed as displayed in fig. 1.d [START_REF] Ji | Erosion and dilation of binary images by arbitrary structuring elements using interval coding[END_REF]. In this case, the initial decomposition can be obtained from a single scan GDT. However, the complexity of the last step depends on the shape of the structuring element. Other methods use a fixed number of scans, but are still restricted to simple shapes such as lines [START_REF] Van Herk | A fast algorithm for local minimum and maximum filters on rectangular and octogonal kernels[END_REF] or rectangles [START_REF] Van Droogenbroeck | Algorithms for openings of binary and label images with rectangular structuring elements[END_REF] and also need combination for other kinds of elements [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF]. In order to deal with arbitrary shapes, combinatorial and genetic algorithms have been proposed [START_REF] Anelli | Decomposition of arbitrarily shaped binary morphological structuring elements using genetic algorithms[END_REF].

Convex Polygon Decomposition for Single Scan Erosion

The proposed method is the combination of a construction scheme used to recursively build structuring elements (section 3.1), a generalized distance transform (section 3.2) and a decomposition algorithm (section 3.3) which determines how structuring elements have to be assembled to obtain a given convex polygon.

A sample polygon P is shown in fig. 2. It is convex since it is equal to the intersection of all the half-planes supported by its sides. The aim of the method is to obtain structuring element B, the discrete counterpart of P. B is the set of discrete points of the square grid included in the closed polygon P. The construction is directed by a series of increasing polygons {P(i)} i∈ [2..N ] used as templates for the structuring elements assembling. Each structuring element B(i) is the discrete counterpart of its corresponding polygon P(i), defined in the continuous plane.

p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p r

B(1) b r B(2) b b r b B(3) b b b r b b b b B(4) r 'T p1 p2 p

p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p r

I1(i) 0 1 2 3 I2(i) 0 0 1 3

Structuring Element Construction

Like the methods recalled in the previous section, the proposed structuring element construction scheme recursively builds a family of increasing elements. However, each structuring element can be built from different smaller elements (conversely to dilation-based construction) and size increments are not fixed for each neighborhood (conversely to chamfer disks). This method operation can be compared to local distance increment with varying weights. Each structuring element B(i) is the union of smaller structuring elements translated according to a set of neighbors {p k }. For instance, in fig. 2,

B(2) = B(1) ∪ (B(1)) p1 B(3) = B(2) ∪ (B(2)) p1 ∪ (B(1)) p2 B(4) = B(3) ∪ (B(3)) p1 ∪ (B(3)) p2
where B( 1) is the simplest element, only containing the origin {O}.

A general expression is given by introducing I k (i), the index of the element used in neighborhood p k for B(i), B(0) the empty set and neighbor p 0 the origin:

∀i = 2 . . . N B(i) = k∈[0,K] (B(I k (i))) p k ( 6
)
The values of I k (i) are summarized in a construction

(i) = i -1, we have (B(I k (i))) p0 = B(i -1), so B(i -1
) is always a subset of B(i). Without loss of generality, we can assume that each

I k table contains increasing values (I k (i) ≥ I k (i -1)).
Comparison with Other Methods. This construction scheme generalizes the disk or structuring element construction methods previously recalled. Chamfer distances use constant local distance increments which correspond to a fixed difference between a constructed disk size i and the included disk size

I k (i). Chamfer distance d a,b is obtained with I k (i) = i -a or I k (i) = i -b depending on p k .
Dilation series are obtained by taking I k (i) = i -1 for each p k belonging to the structuring element used to build B(i).

Each p k can be any point in the discrete plane. The neighbor set is determined from the shape of the structuring element (section 3.3).

Single Pass Generalized Distance Transform and Erosion

The value of the distance map at point x is the index of the largest structuring element centered in x contained in X. It is built from elements located on x neighbors: {x + p k }. The current element size is the greatest one that contains all the neighbor elements:

d X (x) = max{i|∀k, I k (i) ≤ d X (x + p k )} = min k {max{i|I k (i) ≤ d X (x + p k )}}
In order to speed up the distance transform computation, we introduce M k (j):

M k (j) = max(i|I k (i) ≤ j). ( 7
)
The distance transform is then:

d X (x) = min k (M k (d X (x + p k ))) if x ∈ X 0 otherwise (8)
M k (j) represents the index of the largest element B(i) that can be built with B(j) in the neighborhood k. M k (j) is at least equal to 1 due to column 1 filled with 0 in table 1. M k can be computed once from the construction table I k . table 2 shows the M k values corresponding to the example construction values displayed in table 1 The overall complexity is linear with the number of image pixels like all GDT. Furthermore, if all the neighborhoods are chosen to be causal (all p k precede O in the scan order) then only one image scan is needed. While it is also true for some GDT for few restricted shape classes, this GDT works with any convex polygonal shape as it will be shown in next section.

The erosion of X by B = B(N ) is finally:

p ∈ X B ⇔ d X (p) = N
The causality hypothesis implies that the last vertex in scan order must be equal to the origin O. The single-scan algorithm structure permits to use it in a pipeline chain, with one stage for each morphological operation (for instance, a morphological opening requires two pipeline stages). A first implementation has been realized on a Xilinx Spartan IIE FPGA educational card fed with a PAL video signal. The FPGA handles the input synchronization signal and regenerates it on the output. Due to the low cost of the algorithm, at least 8 morphological pipeline stages with different structuring elements can be handled at video rates without extra resources (only in-chip memory is used). The inputoutput delay is only a fraction of a pixel for each stage and an extra delay can be introduced in the output synchronization signal to compensate the translation of the structuring element center. Table 3. Half-plane location table

i 0 1 2 3 4 A0,-1(P(i)) -∞ 0 0 0 0 A-1,2(P(i)) -∞ 0 1 2 4 A-2,1(P(i)) -∞ 0 1 2 3 A1,0(P(i)) -∞ 0 0 0 0

Convex Structuring Element Decomposition

The proposed structuring element construction evokes the anisotropic growth of a single crystal in which epitaxial layers of atoms are successively deposited on a crystal seed. The shape of the crystal is influenced by the physical properties of atoms which constrain the orientations and by the speed of the deposit which may differ from an orientation to another. The orientation of its sides remain constant during the growth. The shape of the structuring element is controlled by artificial constraints which maintains the direction of its sides. However, as the discrete plane produces orientation artifacts especially for small structuring element sizes, the growth is proceeded on a family of continuous polygons which are then used as templates for the structuring elements. The decomposition method is able to process any convex polygon i.e. any closed shape that can be obtained from the intersection of half-planes. For instance, P(4) shown in fig. 2 is bounded by the following half-planes:

(x, y) ∈ P(4) ⇔        -y ≤ 0 -x + 2y ≤ 4 -2x + y ≤ 3 x ≤ 0 ⇔        A 0,-1 (P(4)) ≤ 0 A -1,2 (P(4)) ≤ 4 A -2,1 (P(4)) ≤ 3 A 1,0 (P(4)) ≤ 0 (9) with: A p,q (X) = max (x,y)∈X (px + qy)
The decomposition algorithm consists in moving the planes from their initial seed position (tangent at the origin with A p l ,q l = 0) to their final position. A series of positions is computed for all half-planes as displayed in table 3 for fig. 2 polygons.

Half-planes locations are set in such a way that the sides of intermediate polygons P(i) have a constant orientation and an increasing length:

∀l, ∀i, ∀j, det[v l+1,i -v l,i ; v l+1,j -v l,j ] = 0 (10) ∀i > 0, ∀j ≥ i, ||v l+1,i -v l,i || ≥ ||v l+1,j -v l,j || (11)
Structuring elements from polygons The index of the structuring element used in neighborhood p k is determined as the largest polygon translated by p k that is included in P(i):

I k (i) = max(i : (P(I k (i))) p k ⊆ P(i) = max(i : ∀l, A p l ,q l P(I k (i)) + A p k ≤ A p l ,q l P(i))
This expression of I k (i) ensures that every structuring element B(i) is a subset of the corresponding polygon P(i): ∀i, B(i) ⊆ P(i).

Polygon set As a result of the polygon side properties (constant orientation and increasing length, eq. 10, 11), the series of polygons can be iteratively constructed by Minkowski sums in the continous plane [START_REF] Ohn | Morphological decomposition of convex polytopes and its application in discrete image space[END_REF]. A direct consequence on the construction is that:

v l,i = v l,j + p k ⇒ I k (i) ≥ j ⇒ (v l,j ∈ B(j) ⇒ v l,i ∈ B(i)) Therefore, if the set of intermediate vertex v l positions {v l,i } i∈[i.
.N ] contains a path from O to v l using neighbor moves (plus extra non discrete positions), then v l is necessarily contained in P. Algorithm 1 takes this point into consideration.

Half-plane positions are guided by the movement of vertices. Each vertex is initially located at the origin and follows a path to its final position using the two neighbors of its influence cone. The algorithm ensures that each position in the path is correctly reached by the half-plane, i.e. that half-plane boundaries meet exactly at the vertex intermediate positions.

Neighbor selection This phase is actually the first in the decomposition process, it must ensure that the obtained structuring elements are convex and that paths to vertices can be obtained with polygons of increasing side length (eq. 11). Each pair of successive neighbors defines an influence cone that have some similarities with chamfer disks geometry [START_REF] Thiel | Chamfer masks: Discrete distance functions, geometrical properties and optimization[END_REF]. In an influence cone, each pixel is reached by a series of moves along the two neighbors. The main difference with chamfer distances is that the vertices of the structuring element do not necessarily belong to boundaries between cones. Therefore the number of needed neighbors is generally less than the number of vertices in P. There are two constraints on the pair of neighbors (p k , p k+1 ): i all pixels from the influence cone must be reachable from the neighbors. A necessary and sufficient condition is that p k and p k+1 form a regular cone (det[p k , p k+1 ] = 0) [START_REF] Thiel | Chamfer masks: Discrete distance functions, geometrical properties and optimization[END_REF]. ii all pathes to a point must be included in the structuring element. We have introduced a unified structuring element construction scheme, the corresponding generalized distance transform algorithm and a convex polygon decomposition method. Eroding an image only requires a single regular scan of the image pixels which differs from the classical chamfer distance transform by table lookups instead of constant local distant increments. The computational properties of these algorithms allow their use in a pipeline manner, optimizing time and memory consumption for series of morphological operations.

Fig. 1 .

 1 Fig. 1. Some disk and structuring element construction examples. a: d4 disks (first row). b: octagonal distance disks (second row). c: chamfer distance d2,3 disks (third row). d: line elements obtained by a GDT are gathered in the last step (last row)

Fig. 2 .

 2 Fig. 2. Series of structuring elements (top), series of polygons (bottom)

Table 1 .

 1 Structuring element construction table (see text concerning column 1).

	i	1	2	3	4

Table 2 .

 2 Generalized distance transform table M k . M1(1) = 2 because disk 2 can be built with disk 1 in neighborhood 1 but disk 3 can not

	j	0	1	2	3	4
	M1(j)	1	2	3	4	4
	M2(j)	2	3	3	4	4

  If p ∈ B(i) is in the cone (p k , p k+1 ) and p = ap k + bp k+1 then all the points in the parallelogram (O, ap k , p, bp k+1 ) must be inB(i). Ap m ,qm ({v l,i }) ≤ AH m then {All half-planes contain v l,i , move it to the next intermediate location} choose neighbor k v l,i+1 ← v l,i + p k else v l,i+1 ← v l,i end if end for {Half-plane shift} for l ← 2 to L -1 do {Reach the closest vertex v l or v l+1 } AH l ← min(Ap l ,q l {v l , v l+1 }) end for i ← i +1 end while Algorithm 2 Determination of neighbors {Selection of the two initial neighbors} p1 = (v2 x / gcd(v2 x , v2 y ), v2 y / gcd(v2 x , v2 y )) p2 = (vL x / gcd(vL x , vL y ), vL y / gcd(vL x , vL y )) {Neighbor insertion for condition i} while k < K do if det(p k , p k+1 ) = 1 then {(p k ,p k+1 ) is not a regular cone} find a et b with extended Euclide's algorithm such that bp k xap k y = 1 end if k ← k + 1 end while {Neighbor insertion for condition ii} for n = 1 to L do {Detection of the cone (p k ,p k+1 ) containing vn} while vn is not atteignable do {division of the cone} insert neighbor p k +p k+1 after k if vn is in the second half-cone p k+1 +p k+2 (indices after insertion) then

	Algorithm 1 Half-planes shift computation
	i ← 0 while ∃l : v i,l = vi do {Update reached vertices intermediate position} for l ← 1 to L do {Test of vertex v l,i } if ∀m, n ≥ bp k+1 x -ap k+1 y
	k ← k + 1 end if
	end while
	end for

p k x p k+1 y -p k y p k+1 x > n -1 insert (a + np k x , b + np k y ) after k (indices above k are shifted)