Marc Bonnet

A REGULARIZED GALERKIN SYMMETRIC BIE FORMULATION FOR MIXED ELASTIC BOUNDARY VALUE PROBLEMS

HAL is

Although variational Galerkin-type boundary integral formulations have received a lot of attention, mostly in the applied mathematics field (see e.g. Refs 2, 8, [START_REF] Nedelec | Integral equations with non integrable kemels[END_REF][START_REF] Wendland | Mathematical properties and asymptotic error estimates for elliptic boundary element methods[END_REF] they have led to relatively few developments and applications in the field of solid mechanics (see Refs 1,[START_REF] Nishimura | A boundary integral equation method for an inverse problem related to crack detection[END_REF][START_REF] Parreira | On the implementation of the Galerkin approach in the boundary element method[END_REF][START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF], compared to the very popular collocation boundary element methods (CBEM). Usually variational methods involve double surface integrals (see however Ref. [START_REF] Bui | On the variational boundary integral equations in elastodynamics with the use of conjugate functions[END_REF], where a symmetric BIE formulation with a single surface integration is discussed). Variational approaches have nevertheless attractive conceptual features over usual collocation methods: known convergence results are more extensivei3 and actually better than those of CBEM, and they lead to symmetric stiffness-like BEM matrices. Moreover, as it will be seen below, the variational traction boundary integral equation (TBIE) requiresC o,a. interpolation ofthe densities instead ofC l,a. at the collocation point for the collocation TBIE. The usual C 0 conformal shape functions may then be used. On the other hand, the implementation of this approach for general engineering problems is somewhat more involved than using collocation, which explains perhaps its relative lack of popularity.

In the present paper, the author investigates the statement of a Galerkin weighted residual symmetric BIE formulation for the mixed boundary value problem oflinear elastiticy. More precisely, the objective is to present an application of previous works on indirect regularization of displacement CBIEs4•5•6 and of traction BIEs14 to the derivation of a Galerkin formulation for mixed elastic boundary value problems in which either the inner or the outer integration is at most weakly singular while the other integration is regular. As will be seen below, this may improve existing Galerkin formulations for elastic problems found in the recent literature.

We consider the static response of a 3D isotropic, either bounded or unbounded, elastic body n (Lame constants 'A,Jl), subjected to given tractions F on a portion S T of its boundary an while a displacement u d is prescribed over the complementary partS:,= an-Srofthe boundary (body forces are not considered here). The displacement field u(y) (y e Q) is then governed by the homogeneous Navier equation together with Hooke's law: [START_REF] Hamdi | Fonnulation variationnelle par equations integrates pour le calcul de champs acoustiques lineaires proches et lointains[END_REF] and the above-indicated boundary conditions and, ifn is unbounded, suitable decay conditions at infinity. The comma, in eqn( 1) and throughout the paper, indicates partial differentiation with respect to the components of y, and Einstein summation convention is used. Assume that each point of an has an unique tangent plane, the unit normal vectorn being directed towards the exterior of the material body.

Letx e R3, and U/'(xJ') be the i-th component of the (singular) displacement field aty due to the unit point force acting on x along ek direction in an infinite elastic medium (Kelvin fundamental displacement), andi:7...C xJ') the is-component ofthe Kelvin fundamental stress tensor, given by:

(2)

In eqn [START_REF] Bendali | Prob leme aux limites exterieur et interieur pour le systeme de Maxwell en regime harmonique[END_REF] and throughout the present paper, the comma used with two point kernels denotes differentiation with respectto the second argument y. The Kelvin solution possesses the following well-known symmetry properties:

k i k Ui (x, y) = Uk(y, x) = Ui (x, y) a k 8 k ) k (3)
-a Ui (x, y) = --a Ui (x, y = -Ui,r(x, y)

Xr Yr and singularity properties (r = IJx-Yll denoting the euclidean distance betweenx andy):

Ul(x, y) = 0(1/r) � 7s(x, y) = 0(1/r2)

Let z be a fixed point, either interior or exterior ton, (ie. not on the boundary an). The integral boundary representation formulas for displacement u and stress cr read:

KUk(z) + { Ui(y)ns(Y)� fs (z, Y) dS y lan (5)
f ti(y)Ul(z, y) dS y = 0 lan KO'ab(z)-Cabkr r ui(y)ns(Y)�fs,r(z, y) dS y lan �)

+ f ti(y)E� b(z, y) dS y = 0 lan whereK= 1 (zinterior ton) orK=O(zexterior ton). In eqn(6), use has been made of eqn(3), which allows the exchange of differentiations with respect to the source point z and the integration pointy.

Let X E an be a fixed boundary point with an unique outward unit normal n(x) . The limiting case z 4X in eqn [START_REF] Bonnet | Regularized boundary integral equations for three dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF] and in the dot-product of eqn [START_REF] Bonnet | Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods[END_REF] with n(x) yield respectively the so-called displacement and traction collocation boundary integral equations, which are the basis of various collocation BEM formulations. As a matter of fact, due to the well-known singularity properties of the Kelvin fundamental solutions and their derivatives, eqn( 4 ), these boundary integral equations contain integrals which are classically defined as the result of specific limiting processes, namely Cauchy principal value for the r•2 kernels and Hadamard finite part for ther•3 kernels. More recently, it has been shown (see Ref.6), that these singular integrals, suitably rewritten using so called 'indirect regularization' methods, are actually weakly singular integrals provided certain regularity requirements foran, u and t are met, (see also Ref.l l ).

An alternative approach is to consider eqns( 5) and ( 6) in a weighted residual sense, taking the work-like inner products of eqn( 5) with 1iz) and eqn [START_REF] Bonnet | Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods[END_REF]with ua(z)nb(z), where t]z) and uJz) denote fictitious traction and displacement fields on an, and then integrate the result over an with respect to z. This is the basic well-known idea for obtaining Galerkin formulations, which are the subject of the present paper. By performing this manipulation, one can in principle arrive to a symmetric �eighted-residuals formulation if additional requirements onuJz),

t k (z) .
However, difficulties arise from the strongly singular and hypersingularnature ofthe Kelvin kernels, so thatGalerkin formulations cannot be obtained in a straightforward way using the procedure outlined_above. In Ref.l8, the authors con�ider an auxiliary closed surface S such that an lies in the interior of S, and define the fictitious fields �(z) and U0(Z) for Z E s: thus obtaining a formulatign in which the inner integral is taken on an and the outer integral onS. At this stage, both surfaces and all surface fields are discretized in a usual boundary element fashio11:. Finally, in order to tackle the singularity problem, they take the limit S � an at the element integral level. Their numerical implementation of the method hence relies upon an analytical treatment of the limiting process in the singular double element integrals, which they actually perform for 2D situations, straight constant or linear elements, using complex variable techniques. While fully valid in a fundamental viewpoint, this treatment may induce severe practical limitations regarding the choice of elements and interpolation functions.

The purpose of the present paper is the statement of a regularized symmetric Galerkin BIE formulation, identical to the one of Sirtoriet al. in its essence but applicable to any choice of boundary element interpolation. Our formulation relies upon suitable treatments, including regularization, of the displacement and traction boundary integral equations, which are made before any discretization. This avoids the subsequent cumbersome limiting process at the element integral level.

2.Regularized displacement BIE

This step relies upon the generalized rigid-body identity6 which reads:

Vz E R3-{)fl n unbounded ( �>:-'Y)bki + { ns(Y) lan n bounded (?)
'Efs (z, y) dS y = 0

Now eqn( 7) is multiplied byu/z) and subtracted from eqn(5), yielding:

{Uk(z) + { (tli(Y)-ui(z)) 'Efs(z, y) dS y lan -{ ti(y)U i k(z, y) dS y = 0 lan (8)
where the coefficienty assumes the same meaning as in eqn [START_REF] Bui | On the variational boundary integral equations in elastodynamics with the use of conjugate functions[END_REF] ( { ti(y)ua(z)nb(z)'E�b(z,y)dS y dS z ls lan l n t;(y) {is (iia(z)-u.(y)) n,(z)E�,(z, y) dS z + iia(Y) [is n,(z )E�,(z, y) dSz]} dS y (1\) l n t;(y) { ii;(YHt-�<)+ is (Ua(z) -U.( y)) n,(z)E�,(z, y) dS z } dSy where identity' eqn [START_REF] Bui | On the variational boundary integral equations in elastodynamics with the use of conjugate functions[END_REF] has been used and integrations over § and an exchanged. At this point, one can consider the limiting caseS� an in eqn( 11 ). Since the inner integral over Sbecomes only weakly singular, while the subsequent integral over an is regular, the limiting <;_,ase of eqn( 11) is simply the same formula wherean is substituted toS andx E antoz E S:

{ Biqks(z, y)elhsnl(y)uk,h(Y) dS y lan Equation ( 14) shows that the kernel B; "-'.(x,y) has a llxy ll"1 singularity.

Hence, the first (either inner or outer) s'ii rface integration in eqn( 16) is of aweaklysingularnature, while the otheris regular. As a consequence, the limiting case s � an of eqn( 16),j_lere again, is simply_the same formula where an is substituted to sand X E an to X E S:

uk(y)dSydSz = { ejfqnj(z)iii,J(x) lan

(17)
{ Biqks(z, y)erhsnr(y)uk,h(Y) dS y dSx lan -Finally, collecting eqns( 12), ( 17), the limiting form of (IO) whenS � an is given by: E�, ( z, y) dS y } dS., -l. n e; 1, n; ( z )U;,J( z) (18) { Biqks(z, y)erhsnr(y)uk,h(Y) dS y dS x 0 lan where the coefficient a assumes the same meaning as in eqn [START_REF] Bui | On the variational boundary integral equations in elastodynamics with the use of conjugate functions[END_REF].

One can note that, for the result ( 18) to be valid, the surface fields t and u must respectively be piecewise continuous and QJ.a. These regularity requirements are identical to those for eqn [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF]. Likewise, the surface fictitious displacement u must also be C o.a. In contrast, one can note that the collocation traction BIE eqn( 6) requires that u be C 1 •a, see Refs 6, 9 and 11. The latter fact is recognized as serious as far as numerical implementation is considered because the development of general C 1 boundary elements for 3D problems is problematic.

4.Galerkin boundary integral formulation for the mixed boundary-value problem

Now let us apply the previous results to the statement of a variational BIE formulation for the mixed elastic boundary-value problem. More precisely, eqns( 9) and ( 18) are considered with u E V 11 , t E vT [START_REF] Wendland | Mathematical properties and asymptotic error estimates for elliptic boundary element methods[END_REF] V T = fill piecewise continuous on an and t(x) = 0 V X S T } (20) Splitting an into its components S i t and S 7 , eqns [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], and ( 18) then lead to the following variational equations: { Buu(u, u) + Btu(t, u) But(u, t) + Btt(t, t)

.Cuu(u) + .Ctu(u) .Cut(t) + .Ctt(t)

where But(u, t)

Btt(u,t) -{ { ti(y)tk Jsu Jsu (x)Ui k ( x , y) dS y dS x Btu(t, u) Buu(u,u) .Cu ( U) -L t f ( Z ) { ')'ii; (X ) (21) (22) (23) (25)
+ L T (U.(z)-u.(y)) n,(x)E�,(x, y) dSy } dS z

5.Numerical implementation of the regularized variational BIE

The symmetry properties (28) allows a reduction of computational effort for both the matrix building and the linear system solution steps, which have been extensively discussed elsewhere. We rather wish to focus the discussion on the treatment of the singular surface integrals ({ E Eo)

Then, for (a, P)= 1,2, the natural basis (a.), metric tensor (gcx�) and unit normal n onE are given by:

({ E Eo)

The shape functions Nq are also used for the interpolation of the fields u, t, ii, i. Moreover, it can be shown that: and the regular quantitiesNq( p,a;11), r(p, a;11), Uf' (p,a;rt), L.! (p, a;11 ), � i q ks, elp ,a ;T'I ) are defined by: (33) 1 ns(Y)�fs(x, y) ( tli(Y)-Ui( X)) dSy The variational formulation (21) involves double surface integrals, but at most one of the two surface integration is singular, the other being always regular. In such case, the singular integration should be performed first.

VyE E y/g([jn( {)eabcnb(Y( {) )f, c(Y( { )) at = ecxf33 O �cx a/3 (

6.Concluding comments

In this paper, a symmetric Galerkin regularized BIE formulation is established. Its main feature is that, before any discretization, one of the successive surface integrations is at most weakly singular while the other is regular. Thus, the present formulation, which may be understood as an improvement of the one developed in Ref. IS, provides a basis for the numerical treatment of general 3D situations. Conventional boundary element interpolations of any degree for the geometry as well as the unknowns can be implemented in a straightforward manner, using conventional techniques for numerical integration together with a now well-established singularity cancellation method at the shape function level.

The present approach can be readily applied to other similar situations, such as potential problems or elastodynamic problems. Moreover, similar formulations are already known in the literature, notably about variational formulations for crack problems (see references quoted herein).

0

  Finally, let x e an and assume u. e Ql•a. at x, where C"'• .... denotes the set of functions m times contin�ously differentiable such that3(a,C) > 0, I u ; (x) -u ; (y) I�IJx-yJJc:t. Then all integrals ineqn(8)areweakly singular in the limft z � x. As a consequence, eqn(8) taken for x e an is the regularized DBIE.4Now, one can take the inner product of eqn[START_REF] Cortey-Dumont | Simulation numerique de probU::m esde diffraction d'ondes par une fissure[END_REF] written forz=x E an with a fictitious surface traction fieldfcx) and integrate the result over an again. This leads to:'Y { ik (z) {uk(z)+ { (ui(Y) lan lan-u;(x))Ef,(x,y)dS y } dS x(9)-{ { ti(y)tk(z)U i k(x,y)dS y dS x lan lanIt is worth noting that, for the result of eqn(9) to be valid, the surface fields t and u must respectively be piecewise continuous and C o,a (the latter condition ensures that the factor u(y)u(x) actually weakens the singularity of the inner integral). These are indeed the very conditions under which the integrals in the collocation DBIE, eqn[START_REF] Bonnet | Regularized boundary integral equations for three dimensional bounded or unbounded elastic bodies containing curved cracks of arbitrary shape under dynamic loading[END_REF], are convergent, see Ref.l2.

3 .

 3 Regularized traction BIEHere, following the approach of Sirtori et al., I& an auxiliary regular surface S, which completely surrounds an, is inJroduced, such that there exists a one-to-one mapping X E an � Z E S. !_:et US introduce a fictitious displacement field u(z) defined for z E S. Now th� inner product of eqn[START_REF] Bonnet | Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods[END_REF] with ii(z) ® n(z) is taken and integrated overS, which gives:"is u.(z) { O"ab(z)n,(z) + l n l;( y)n,(z) E�,( z, y) dS y } dS z -Cabh is n,(z)Ua(z) (10) l Ui( y)n �(y)'Efs,r(z , y) dSy dS z 0 lan Let us investigate the two double surface integrals in eqn(1 0). One has:

_lim(

  { ti(y)ua(z)nb(z)'E�b(z, y) dS y dS z ZES .... . zEanls lan f ti(Y) {ui(Y)(!-K)+ f (ua(z)[START_REF]Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoe1asticity[END_REF] lan lan -u.(y)) n,(x )E�,( z, y) dSx } dSyIn eqn[START_REF]Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoe1asticity[END_REF], the outer integral is weakly singular integrals while the inner is regular (they have to be considered in this order, because the inner integral, considered alone, is of the strongly singular type).Next, let us consider the integral in eqn(10) containing the hypersingular kernel Cabkr :L�. v• In this case, we will use two successive integrations by parts. Following8ecache1 and Nishimura& Kobayashi,14 there exists a fourth-order tensor B(ZJ' ), associated to the Kelvin stress tensor, such as, for any z -:f. y: the components of the permutation tensor. Moreover, the Stokes formula for a closed surfacer reads: [ eabcna(Y)f,b(Y) dS y = 0 (15) for any continuous and piecewise continuously differentiable function fly) (if fly) is only piecewise continuous, identity ( 15) would also present contour integrals over lines accross which/has a jump). It is worth emphasizing that njy) fb(y) can be expressed in terms of tangential derivatives off Using eqn(13) and applying eqn(15)jwice (once for each of the two integration variables y e an and z e S), one has: ( u i (y)ni(z)Cijab { ��1 ,b(z , y)nr(y)uk(Y) dS y dS z ls lan -( iii(y)ejfqnj(z) { Biqks,t(z, y)erhsnr(Y) ls lan uk,h(Y) dS y = is ejfqnj(z)iii,t(Y)

  given displacements u 0 and tractions tD are incorporated in the linear forms� " and� , . From eqns(21) to (2 7) and taking the symmetry properties of eqn(3) into account, the following symmetry properties, which are expected as classical features ofGalerkin BIE methods like e.g. in Ref.l8, are apparent:Vu, u E Vu, Vt , t E VT Buu(u, u) = Buu(u, u ) Btt(t, t) = Btt(t, t) (28) But(1t, t) = Btu(t, u )The variational BIE formulation (21) is the main result of the present paper. Indeed it must be understood as a reformulation of the variational symmetric BIE formulation developed in Ref.[START_REF] Sirtori | A Galerkin symmetric boundary element method in elasticity: formulation and implementation[END_REF]. The difference between the two formulations, which constitutes the main motivation of this paper, is that our result (21) has undergone a preliminary regularization treatment, using singularity exclusion methods developed in our previous works (e.g. in Ref.[START_REF] Bonnet | Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods[END_REF]) together with a stress function and integration by parts method which has been previously used, e.g. in Refs 1 and 14 for the BIE modelling of crack problems. As a result, our variational formulation can be expected to be more flexible in aB EM implementation point of view, since boundary element geometrical and field interpolation of any degree can be incorporated into eqn(21) without difficulty, as will be briefly shown in the next section.

  which appear in eqn(21 ), in the simplest case ofisoparametric conformal BEM interpolation. The surface on is divided into boundary elements, which are mapped onto areforence element E0, which is generally the square�= (�1 ' � 2 ) e [-1,1] 2 or the triangle 0 � � 1 + �2 � 1. Let the discussion be restricted to the consideration of a single element£. The location of a pointy onE is expressed in terms of n shape functions A* and n geometrical nodesAk (k = 1 ..... n) located on the boundary of£:

  31) Let rt = ( rt I ' rt 2 ) denote the antecedent of x on £ 0 • Foil owing a common practice in BEM (see e.g. Ref.l7), set �1 = 111 + pcos a, � 2 = 112 + psin a. Then: dSy = y/g([jd6d6 = jll(fj pdpda = J(e)pdpda (32)

  f(p,a;11) * 0 and U(p,a;11), L.(p,a;11) and Aiqh.eg(p,a;11) are regular at p = 0. Hence the singular surface integrals wliich occur in eqn(21) can be written:

E ( 34 )

 34 ul L n,( y )f:7,(p, a; •l)M•(p, a; '1)J(e)dpda (35) L uhx, y)ti(Y) dSy tlj Ul(p, a; 1J ) l\1q ( { ) J ( { )dpda E (36) Equations (34), (35) and (36) take full advantage of the regularization. Their numerical evaluation of can be performed with standard product Gaussian quadrature formulas, provided a further coordinate change (p,a) � (vl'v2) is made in eqn(34) in order to recover an integral over the square [-1, 1 f.4


Keywords: Diversity, Estimation, Extinction, Maximum likelihood, Phylogeny, Simulation, Speciation

There is considerable interest in the possibility of using molecular phylogenies to estimate extinction rates. The present study aims at assessing the statistical performance of the birth-death model fitting approach to estimate speciation and extinction rates by comparison to the approach considering fossil data. A simulation-based approach was used. The diversification of a large number of lineages was simulated under a wide range of speciation and extinction rate values. The estimators obtained with fossils performed better than those without fossils. In the absence of fossils (e.g., with a molecular phylogeny), the speciation rate was correctly estimated in a wide range of situations; the bias of the corresponding estimator was close to zero for the largest trees. However, this estimator was substantially biased when the simulated extinction rate was high. On the other hand the estimator of extinction rate was biased in a wide range of situations. Surprisingly, this bias was lesser with medium-sized trees. Some recommendations for interpreting results from a diversification analysis are given.

Introduction

The tempo and mode of evolution has been one of the fundamental questions in evolutionary biology [START_REF] Simpson | The major features of evolution[END_REF]. Recent advances in phylogenetics have given a fresh look at this issue [START_REF] Barraclough | Phylogenetics and speciation[END_REF]. The reconstruction of the relationships among species allows one to test whether the shape of the reconstructed phylogeny agrees with a given model of diversification (e.g., [START_REF] Kirkpatrick | Searching for evolutionary patterns in the shape of a phylogenetic tree[END_REF][START_REF] Slowinski | Testing whether certain traits have caused amplified diversification: an improved method based on a model of random speciation and extinction[END_REF][START_REF] Mckenzie | Distributions of cherries for two models of trees[END_REF]. In the situation where the branches of the phylogeny have estimated lengths, it is possible to fit models in order to estimate the rates of diversification of the studied lineage (Nee et al., 1994b;[START_REF] Paradis | Assessing temporal variations in diversification rates from phylogenies: estimation and hypothesis testing[END_REF].

Conceptually, the diversification of a lineage may be seen as a series of speciation and extinction events through time. Speciation events give birth to new species, and extinction events result in the death of species. This representation is practical for modelling the diversification of lineages since, with some further assumptions, this agrees with the birth-and-death processes which have been extensively studied in the past (Kendall, 1948a(Kendall, ,b, 1949;;[START_REF] Darwin | The behaviour of an estimator for a simple birth and death process[END_REF][START_REF] Keiding | Maximum likelihood estimation in the birth-and-death process[END_REF]. Nee et al. (1994b) proposed a method to estimate both speciation and extinction rates of a lineage using the reconstructed phylogenetic relationships of the living species, for instance using molecular data. They developed a likelihood-based approach to estimate both speciation and extinction rates. It follows that the extinction rate of a lineage could be estimated even in the absence of fossils, and thus without observing any event of extinction [START_REF] Nee | Inferring population history from molecular phylogenies[END_REF]. This conjecture may seem counterintuitive since the information on extinction comes from extinct species: for instance, when a phylogeny with fossils is analysed, the estimation of extinction rate is done with the ratio of the number of extinction events on the cumulative numbers of species. Clearly, the number of extinction events cannot be observed without fossils. [START_REF] Kubo | Inferring the rates of branching and extinction from molecular phylogenies[END_REF] developed a method close to Nee et al.'s (1994b): they used the same birth-death model but they fitted this model with a polytope algorithm instead of maximum likelihood. [START_REF] Kubo | Inferring the rates of branching and extinction from molecular phylogenies[END_REF] then showed, using simulations, that the estimator of extinction rates has a too large variance to be reliable.

It is important to assess the precision of the method introduced by Nee et al. (1994b) since it could have many potential applications in the study of organisms poorly represented in the fossil record, such as soft-bodied organisms. This method may be used also to study the dynamics of viral populations [START_REF] Nee | Inferring population history from molecular phylogenies[END_REF]. The goal of the present study is to assess the statistical performance of this method in a wide range of parameter values. Nee et al.'s (1994b) method performance was compared to the performance of the method where fossils, and thus past extinction events, are observed.

Methods

A large number of phylogenies were simulated with different values of speciation and extinction probabilities. The approach adopted was to simulate the speciation and extinction events, rather than generate sets of branching times from theoretical distributions, in order to mimic as close as possible the evolutionary process.

Simulations were started with a single species. At each time step, each species living in the lineage had a probability (denoted µ) to die. If it survived, each species had then a probability (denoted λ) to generate two daughter-species, otherwise it simply survived to the next time step. This process was simulated during 1000 time steps; all speciation and extinction events were recorded.

The probability of speciation λ varied between 0.0001 and 0.005 with a step of 0.0001, and for each of these values, µ varied between 0 and λ -0.0001 with the same step of 0.0001. The constraint λ > µ was imposed so that the risk of extinction of the entire lineage was minimized. For each combination of these two parameters, the simulation was replicated 1000 times. Thus, a total of 2,500,000 simulations were started (50 values of λ × 50 values of µ × 1000 replicates). Since the simulations were fully stochastic, so were the number of species at any time step. The lineages with 0, 1, or 2 living species at the end of the simulation were discarded. For the other lineages, two phylogenies were recorded: the first one with all speciation and extinction events so that the extinct species are included, and the second one with only the species living at the end of the simulation.

The phylogenies with only the living species were analysed with Nee et al.'s (1994b) method. This method fits a birth-death process (Kendall, 1948b) to the set of branching times calculated from the phylogeny. A re-parametrization is needed to allow the model fitting so that the parameters under consideration are r = λ -µ, and a = µ/λ (see Nee et al., 1994b). The model is then fitted by maximum likelihood to give the estimates of both parameters, denoted r and â, respectively.

The phylogenies with the living and the extinct species were analysed with the method described by [START_REF] Keiding | Maximum likelihood estimation in the birth-and-death process[END_REF] who gave the following two maximum likelihood estimators of the speciation and extinction rates:

b = B T T 0 X t dt , d = D T T 0 X t dt , (1)
where B T is the number of speciation events between time 0 and T , D T is the number of extinction events during the same time, and X t is the number of species living at time t.

The integral at the denominators of eqn (1) was calculated as the sum of the branch lengths of the phylogeny.

The analyses of the simulated data were all done with R [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF]R Development Core Team, 2003) using a package specially developed for phylogenetic analyses [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]. For all simulated lineages, the estimates r, â, b, and d

were stored together with the number of tips in both phylogenies, and the depth of the phylogeny without fossils (age of the youngest common ancestor of all living species).

Note that the depth of the phylogeny with fossils was always 1000 time steps.

The relative error of each estimated parameter was calculated using:

r -(λ -µ) λ -µ , â -µ/λ µ/λ , b -λ λ , d -µ µ . (2)
Using relative errors avoids having large errors due to the fact that the values of the parameter are large as well. In the simulations with µ = 0, the relative error of â could not be calculated and was replaced by its absolute error: â -µ/λ. On the other hand, when µ = 0 the relative error of d was not considered because there was no extinction event (D T = 0), and thus d = 0 in all replications.

Giving the very large amount of simulated data, summary statistics of the results were computed. The variations in the relative errors of the estimated parameters were investigated with respect to four variables: the true values of λ and µ, the number of tips of the phylogenies, and its depth. To give a picture of how the mean tendency and the dispersion of the relative errors varied, the results were summarized using box-and-whiskers plots. The relative errors were first dispatched in different categories with respect to: (i) each value of λ, (ii) each value of µ, (iii) an interval of the number of tips, or (iv) an interval of the depth of the tree. For the number of tips and the depth of the tree, the intervals were created with an algorithm in R in order to make categories with approximately similar numbers. The numbers of phylogenies analysed in each category are given in the Appendix.

In a second step, for each category a box-and-whiskers was drawn where the box represents the lower and upper quartiles (thus 50 % of the data are within the box), the median being indicated by a horizontal line inside the box. The upper and lower limits of the whiskers extend up to the quantiles Q 0.05 and Q 0.95 , respectively (so that 90 % of the data are between the limits of the whiskers. Giving that a few points in these plots are extreme outliers, particularly towards the positive values (data not shown), showing the extreme of these distributions with the whiskers (as usually done with such plots, [START_REF] Venables | Modern applied statistics with S[END_REF] would increase the range of the y-axis by several orders of magnitude, and would result in tiny unmeaningful boxes. These summary statistics were computed for each category of λ, of µ, of the phylogeny depth, and of the number of tips of the phylogeny without fossils for the relative error of r and the absolute error of â.

They were also computed for each category of λ, of µ, and of the number of tips of the phylogeny with fossils for the relative errors of b and of d.

To allow a more direct comparison of both kinds of estimators, estimates of the speciation and extinction rates without fossils (denoted bo and do , respectively) were computed by back transformation of r and â:

do = r â 1 - â , bo = r + do . (3)
The relative errors of bo and do were calculated in the same way than for b and d.

To give a global picture of the bias in bo and do , the relative errors of these estimates were plotted with respect to the number of tips of the analysed tree for different values of λ and of µ. Five values of λ were selected (0.003, 0.0035, 0.004, 0.0045, and 0.005), and six values of µ (0, 0.0009, 0.0014, 0.0019, 0.0024, and 0.0029) for this analysis. A local polynomial fit (Venables and Ripley, 2002, p. 230) was performed on each plot.

Results

The 2,500,000 simulations yielded 604,165 lineages that were amenable to further analyses (that is they had at least three species living after 1000 times steps). Fig. 1 shows an example of such a simulated lineage with both derived phylogenies. Among the 604,165 phylogenies considering only living species, estimation of r and a was not possible in 112,083 cases (there was no convergence of the fitting algorithm). This failure was clearly related to the number of tips in the tree: the mean number of tips among these 112,083 trees was 5.73 (sd = 2.94, median = 5, maximum = 32), whereas it was 31.63 (sd = 46.67, median = 16, maximum = 1140) among the 492,082 other trees.

The number of tips varied between 3 and 1140 for the phylogenies without fossils, and between 3 and 1161 for the phylogenies with fossils. In both cases, the distribution had its maximum at the smallest value, and was highly skewed towards the largest ones.

Fig. 2 shows the results for r. There was a positive bias in this estimator for the smallest values of λ, and the median of the relative error of r converged to zero for increasing values of λ to be close to zero at λ ≈ 0.002. The dispersion of the relative error of r also decreased with increasing values of λ. The opposite results were observed with respect to the values of µ: the median of the relative error of r was very large for the largest values of µ with a median of about 15 for µ = 0.0049. The median of the relative error of r was somewhat unsensitive to the number of tips with just a small negative bias for the phylogenies with three tips. However, the dispersion of this error was greatly influenced by the number of tips: this dispersion decreased continuously for all categories of number of tips. A positive bias in r was observed when the depth of the tree was less than about 400 time units. This bias progressively decreased to zero when the depth of the tree increased. The dispersion of relative error of r also decreased when the depth of the tree increased to stabilize at a depth of about 700 time units.

Fig. 3 shows the results for â. With respect to λ, the relative error of â was large for the small values of λ and slowly converged to zero when λ increased. The relative error of â was greatly affected by the value of µ: there was a systematic negative bias in â meaning that it underestimated the actual value of a in most situations. On the other hand, the dispersion of the relative error of â decreased greatly with increasing values of µ. The effect of the number of tips of the tree on â was complex: there was a negative bias in â for the smallest numbers of tips, and this bias converged to zero with increasing numbers. The median of the relative error of â was 0 for trees with 16 tips, but a negative bias was observed for trees with more tips. With respect to tree depth, the median relative error of â converged progressively to zero with increasing values of tree depth.

The dispersion of the relative error of â was large for all values of tree depth.

Fig. 4 shows the results for b. The relative error of b was very large for very small values of λ, and its median quickly converged to zero with increasing λ. It is remarkable that the effect of λ on b was very similar to that on r with two differences, however: the dispersion of the relative error of r was greater (for λ = 0.005, the first and third quartiles were -0.118, 0.067 for b, and -0.264, 0.160 for r), and the median of the relative error of b converged to zero more quickly though the difference was slight but systematic (for λ = 0.0016, the median relative errors of b and r were 0.099 and 0.203, respectively).

The relative error in b was only slightly affected by the value of the extinction rate with a small positive bias for the large values of µ. In all cases, the dispersion of b was small with 50 % of the values being between -0.2 and 0.2. As previously for r and â, a negative bias in b for the smallest numbers of tips was observed. However, the dispersion of the relative error of b was much smaller than for r. As for the other estimators, the relative error in d was affected by the number of tips in the tree, but the dispersion was greater than for b.

The estimates of speciation rate with ( b) and without ( bo ) fossils were compared for the different values of λ and µ (Fig. 6). With respect to λ, the performance of both estimators were very close except for the very small values of λ where the estimates with fossils were somehow better than without, though a positive bias was also observed for the former. In terms of dispersion (as measured by the inter-quartile range), the estimates with fossils were slightly better than without fossils, particularly for the largest values of The same comparison was done between d and do (Fig. 6). With respect to λ, there was always a negative bias for both estimators, and they both converged to zero with increasing values of λ but the convergence was much quicker for d than for do . The dispersion of these estimates showed opposite trends, though the results were somewhat more complicated. With fossils, the dispersion of d increased with increasing values of λ, but then decreased for λ > 0.0023. On the other hand, the dispersion of do increased for all values of λ, though it reached a plateau for 0.0003 < λ < 0.0023. With respect to µ, both estimators showed parallel results but the performances of d were much better than for do , in terms of bias (the median of the relative error of d was much closer to zero than the one of do ) as well as in terms of dispersion of the estimates.

The analysis of the relative error of bo with respect to the number of tips showed very similar results for the different values of λ and µ (Fig. 7). In all cases the fitted curve was close to zero for a number of tips between 10 and 20 or more. The remarkable difference among the different plots was the maximum number of tips which was clearly related to the difference λ -µ. All plots also showed similar results for the dispersion of the relative error of the estimates with a range of values quite similar, and a decrease in dispersion with an increasing number of tips.

The same analysis for do showed more complex results (Fig. 8). When µ = 0 the dispersion of the relative error of do (as shown by the range of the y-axes) was very low, and the fitted curves converged to zero with growing numbers of tips. On the other hand, when µ > 0 the fitted curves were almost always negative, increased with growing numbers of tips up to 10-20 tips, and then decreased afterwards. There was a continuous variation though from the smallest to the largest values of µ: positive values of the fitted curves were observed for µ = 0.0009, whereas the fitted curves were well below zero for µ = 0.0029. The dispersion of the relative error of do decreased with increasing values of both λ and µ.

Discussion

There is undoubtedly considerable interest in the possibility to use molecular phylogenies (or any phylogeny inferred from extant species) to estimate extinction rates since this approach does not require to collect data through time. However, it is necessary to assess the biases and limits of this approach which is the goal of this paper.

The assessment of the bias of the estimators with fossils was done mainly for comparison: since the lineages that went extinct before 1000 time steps were not considered, the biases of these estimators were not correctly estimated.

Overall, the estimators obtained with fossils ( b and d) performed better than those without fossils (r and â): the former had generally smaller bias and smaller variance than the former. All four estimators performed better with increasing values of speciation rate which is obviously due to the fact that high values of λ result in more speciation events in the simulated trees eihter with or without fossils. On the other hand, when the extinction rate was high, all estimators were biased since most trees went extinct and thus only those with a low realized value of µ were effectively analysed.

All estimators behaved generally well with respect to the number of tips: the median errors were close to zero for 15 or more tips. However, the dispersion of all four estimators was continuously influenced by the number of tips: it was lowest for the largest number of tips. An exception to this pattern was â which showed a negative bias for the lowest and highest numbers of tips. Clearly, more tips result in more data to analyse, and thus it is expected that the estimators are more accurate (median close to zero, and low dispersion). In the case of â, the negative bias for the highest numbers of tips may be due to the fact that most large trees were simulated with a low value of µ and a high value of λ. Consequently, these trees had a very low ratio λ/µ (= a) which critically affected the relative error of â. It is noteworthy that this negative bias for the highest numbers of tips was also observed when considering the absolute error of â, but was much slighter than for the relative error (results not shown).

In the trees without fossils, the most recent common ancestor to all living species (i.e. the root of the tree) varied randomly depending on the simulation. The median relative error of r was close to zero for tree depth values of about 450 or more, whereas a value of about 900 or more was required to obtain a median error of â close to zero. The depth of a tree is clearly related to extinction rate: the highest the value of µ, the lowest the probability of a species appearing early in the simulation to survive until present.

It should be noted that all the results in this study were obtained with trees, either with or without fossils, which were known without error in terms of unlabelled topology and branch lengths. This is unlikely to be always true in real situations since there are many sources of error in estimating phylogenetic trees as clearly illustrated by various phylogenetic studies (see Whelan et al., 2001, for a review). If these errors in estimating trees are uniformly distributed along the tree, it should be expected that the median errors of the estimators studied here are not affected, though their dispersion are likely to be increased since a supplementary source of variation is added. On the other hand, if a systematic bias in estimating the trees exists, this will add a bias in the estimators of diversification.

Another assumption of the present study was that all species of the lineage are included in the reconstructed phylogenies. This is likely to be untrue in real situations.

With fossils, many species are likely to have not been fossilized and thus cannot be included in a possible phylogenetic study. Without fossils, it is rare to have all living species of a clade to be included in a phylogenetic reconstruction. In the context of testing for temporal variation in diversification, it was shown that missing taxa in phylogenies induce a bias resulting in a substantial increase in the type I error rate (Nee et al., 1994a;[START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF]. There has been no assessment of the possible bias of the estimators r and â when a phylogeny is incomplete. However, it could be speculated from the present study that missing taxa may have an effect on the estimation of r and a. Removing randomly some species from a clade can be seen as similar to analysing a smaller clade with the same parameters λ and µ. Thus the effects of missing species can be predicted from the effects of the number of tips observed on Figs. 7 and8.

Consequently, missing taxa in phylogenies are likely to bring about a negative bias in r (and bo ), whereas the effect on â (and do ) would be more complex. In both cases, the bias is likely to be slight if the number of missing species is low.

It appears that accurate estimation of both speciation and extinction rates can be achieved only with fossil data. In the absence of fossils (e.g., with a molecular phylogeny), only the difference between these rates (r) can be estimated with some accuracy in a wide range of situations, notably when the speciation rate is relatively large compared to the extinction rate. On the other hand, the estimation of a was inaccurate in a wide range of situation. This result was anticipated by Darwin (1956, p. 30) who stated that "If N 1 , . . . , N k [the cumulative numbers of species] are the only observed quantities, estimation of µ/λ is likely to be very inaccurate since the range of values of µ/λ giving the same set N 1 , . . . , N k with a reasonable probability is very large." The present study gives some empirical numerical support to Darwin's logical argument. Remarkably, accurate estimation of a was achieved when the extinction rate was close to zero and the tree depth was close to the actual age of the lineage, suggesting that a may be correctly estimated only when the tree without fossils is close to the tree with fossils (i.e. the 'real' tree of the lineage).

An unexpected result comes from the fact that the error in do increases with increasing number of tips in the tree . By contrats to â, this cannot be explained by the fact that larger trees are produced by larger values of λ (see above). It seems rather that large trees produced with a moderate value of µ do not show a typical distribution of their branch lengths. This is further evidenced by the fact that trees with a moderate number of tips (≈ 20) gave good results in terms of relative error of do . Using the formula for the expected mean number of species in a clade after a time T , e (λ-µ)T (Kendall, 1948b), it can be found that a clade with an expectation of 20 species after T = 1000 is characterized by λ -µ ≈ 0.003.

Conclusions and Recommendations

The results from the present study can be used to define recommendations when interpreting an analysis of diversification. It is generally not possible to choose between both situations considered here, with or without fossils. Lineages with fossil data are often extinct, and those which are studied with a molecular phylogenetic approach usually have no or a poor fossil record. Indeed, the comparison between both kinds of estimators was not intended to define guidelines but to give a comparison for the estimators without fossils.

Phylogenies should be analysed as complete as possible: missing species are likely to introduce a negative bias in the estimation of λ, and an undetermined bias for µ. A clade with at least 15 species is appropriate to estimate λ. This parameter is likely to be underestimated with smaller trees. Surprinsingly, µ is likely to be well estimated with medium-sized trees (with 10-20 species). However, µ is likely to be underestimated in most situations except if µ ≈ 0 where an overestimation is expected. In the latter case, the bias will be very small though. In all cases, it seems better to consider do as a lower bound of µ.

If some informations are available on the age of the studied clade and the age of the most recent common ancestor of the species included in the phylogeny (called tree depth in the present paper), this may be used in interpreting the estimates of speciation and extinction rates. If the ratio of the latter on the former is less than 0.5, then the estimate of λ is likely to be positively biased. If this ratio is less than 0.9, then the estimate of µ is likely to be negatively biased.

The present study is the first extensive analysis of the statistical performance of the birth-death estimators as applied to phylogenetic data. [START_REF] Nee | Inferring speciation rates from phylogenies[END_REF] used simulations to compare the statistical properties of several estimators but he considered only the birth-only model (also called Yule model). Some further studies are clearly needed, particularly to assess the robustness of the birth-death estimators when rates vary through time or across lineages since such situations are likely to be more biologically plausible than the homogeneous rates case considered here. F igure 3: Results for â. On each plot, the effect is considered alone and the trees are combined for all the different values of the other effects (e.g., the trees with the same value of λ may have different number of tips).

Fig. 5

 5 Fig.5shows the results for d. The relative error of d was greatly influenced by λ with a strong negative bias for the smallest values of speciation rate and a progressive convergence of the median with increasing values. With respect to µ, there was a strong negative bias for the smallest values of extinction rate with a quick convergence to zero with increasing µ. However, a negative bias was observed for the largest values of µ. The dispersion of the relative error in d continuously decreased with increasing values of µ.

λ.

  With respect to µ, the performances were more contrasted. With fossils, the bias in b was close to zero for all values of µ but showed a positive trend with increasing values of µ. Without fossils, the bias in bo showed a negative trend, and was much stronger than for b for the largest values of µ. The contrast was even stronger when considering the dispersion of the estimates: the dispersion of bo was always stronger than the dispersion of b and increased with increasing values of µ, whereas the dispersion of bo decreased with increasing values of µ.

Figure 1 :

 1 Figure 1: Illustration of the two phylogenies derived from the simulation of a lineage with speciation rate λ = 0.005 and extinction rate µ = 0.001. (A) Complete phylogeny with the extinct species. (B) Phylogeny with only the species living at the end of the simulation. Note that the three most ancient species have disappeared from the data in (B). The parameter estimates for these data are: (A) b = 0.0058, d = 0.0015, and (B) r = 0.0052, â = 0.

  for b. On each plot, the effect is considered alone and the trees are combined for all the different values of the other effects (e.g., the trees with the same value of λ may have different number of tips).

  for d. On each plot, the effect is considered alone and the trees are combined for all the different values of the other effects (e.g., the trees with the same value of λ may have different number of tips).

Figure 7 :

 7 Figure 7: Relative error of the estimates of speciation rate without fossils ( bo ) with respect to the number of tips of the analysed tree. The parameters λ and µ are the speciation and extinction rates of the simulated trees, and are arranged as rows and columns, respectively. Each grey point is a replication; the lines are local polynomial fits. 29

Figure 8 :

 8 Figure8: Relative error of the estimates of extinction rate without fossils ( do ) with respect to the number of tips of the analysed tree. The parameters λ and µ are the speciation and extinction rates of the simulated trees, and are arranged as rows and columns, respectively. Each grey point is a replication; the lines are local polynomial fits.

Table A2 .

 A2 Number of simulated lineages with respect to each value of extinction rate µ.

	µ	Number of trees µ	Number of trees
	0	37360 0.0025	8196
	0.0001	35528 0.0026	7571
	0.0002	34013 0.0027	6840
	0.0003	32422 0.0028	6448
	0.0004	30931 0.0029	6006
	0.0005	29406 0.003	5446
	0.0006	27971 0.0031	4845
	0.0007	26425 0.0032	4587
	0.0008	25073 0.0033	4162
	0.0009	23634 0.0034	3727
	0.001	22234 0.0035	3346
	0.0011	20566 0.0036	3040
	0.0012	19411 0.0037	2625
	0.0013	18449 0.0038	2461
	0.0014	17213 0.0039	2109
	0.0015	16106 0.004	1796
	0.0016	15516 0.0041	1571
	0.0017	14419 0.0042	1342
	0.0018	13450 0.0043	1117
	0.0019	12629 0.0044	918
	0.002	11751 0.0045	762
	0.0021	10774 0.0046	542
	0.0022	10137 0.0047	382
	0.0023	9492 0.0048	266
	0.0024	9030 0.0049	120

Table A3 .

 A3 Intervals defined to make groups with respect to the depth of the tree for the phylogenies without fossils.

	Depth of the tree Number of trees Depth of the tree Number of trees
	2-148	9725	783-794	9745
	149-220	9849	795-806	10010
	221-277	9950	807-817	9636
	278-324	9751	818-828	10059
	325-367	9907	829-838	9319
	368-404	9751	839-848	9766
	405-438	9912	849-858	9931
	439-469	9882	859-868	10358
	470-498	9783	869-877	9322
	499-524	9700	878-886	9544
	525-549	9963	887-895	10048
	550-571	9639	896-904	10072
	572-592	9679	905-912	9432
	593-612	9779	913-921	10816
	613-632	10067	922-929	9754
	633-650	9757	930-936	8740
	651-668	10159	937-944	10378
	669-684	9564	945-952	10558
	685-700	9999	953-959	9557
	701-715	9724	960-966	9814
	716-729	9484	967-973	9670
	730-744	10433	974-980	10065
	745-757	9531	981-986	8786
	758-770	10090	987-993	10436
	771-782	9568	994-1000	10620

Table A5 .

 A5 Intervals defined to make groups with respect to the number of tips for the phylogenies with fossils.

		Number of tips Number of trees Number of tips Number of trees
		3	19538	24-25	16913
		4	24095	26	7670
		5	25499	27-28	14293
		6	25061	29-30	12877
	(A)	7	24074	31-32	11635
		8	23054	33-34	10696
		9	21616	35-37	14303
		10	20508	38-40	12640
		11	19086	41-43	11157
		12	18043	44-47	13336
		13	16724	48-51	11397
		14	15725	52-56	12161
		15	14767	57-62	12340
		16	13505	63-70	13300
	(B)	17	12804	71-79	11831
		18	12076	80-91	12212
		19	11457	92-107	11679
		20	10518	108-133	12472
		21-22	19609	134-181	12108
		23	9154	182-1166	12232

Acknowledgements

I am grateful to two anonymous referees for their constructive comments on a previous of this paper. This research was financially supported by the Institut Franc ¸ais de la Biodiversité and the Centre National de la Recherche Scientifique. This is publication 2004-009 of the Institut des Sciences de l' Évolution (Unité Mixte de Recherche 5554 du

Depth of the tree (lower bound of the interval)

Relative error of r ^ F igure 2: Results for r. On each plot, the effect is considered alone and the trees are combined for all the different values of the other effects (e.g., the trees with the same value of λ may have different number of tips).


Emmanuel Paradis
email: paradis@isem.univ-montp2.fr

Analysis of diversification: combining phylogenetic and taxonomic data

Keywords: diversity, estimation, extinction, maximum likelihood, phylogeny, speciation

The estimation of diversification rates using phylogenetic data has attracted a lot of attention in the last decade. In this context, the analysis of incomplete phylogenies (e.g., phylogenies resolved at the family-level but unresolved at the species-level) has remained difficult. I present here a likelihood-based method to combine partially resolved phylogenies with taxonomic (species-richness) data in order to estimate speciation and extinction rates. This method is based on fitting a birth-and-death model to both phylogenetic and taxonomic data. Some examples of the method are presented with data on birds and on mammals. The method is compared with existing approaches that deal with incomplete phylogenies. Some perspectives and generalizations of the approach introduced in this paper are further discussed.

INTRODUCTION

Evolutionary biology deals with changes in species diversity through time and their mechanisms. Such changes are brought about by speciation and extinction events, and quantifying the rates of these events contributes to a better understanding of the biological mechanisms of evolution. During the past decade, there has been a lot of interest in using phylogenetic data to estimate speciation and extinction rates. This follows from the increasing number of phylogenetic studies using molecular techniques which has led to resolve the relationships among more and more taxa. Most of the methods to study diversification with phylogenetic data consider phylogenies at the species level, and thus require that most (if not all) species of the studied lineages be included in the reconstructed phylogenies [START_REF] Hey | Using phylogenetic trees to study speciation and extinction[END_REF]Nee et al. 1994a;[START_REF] Sanderson | Reconstructing shifts in diversification rates on phylogenetic trees[END_REF][START_REF] Paradis | Assessing temporal variations in diversification rates from phylogenies: estimation and hypothesis testing[END_REF]Paradis , 1998a;;[START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF]. However, quite often phylogenetic data are available at a higher level than the species, for instance giving the relationships among families whereas the relationships among species within the families are unknown though the numbers of species in these groups are known. Examples of such situation are provided by the phylogeny of birds where almost all 144 families are included but about 1000 species out of 9646 are represented [START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF], or the phylogeny of Insects where the relationships among the 30 orders have been inferred but the relationships among the more than a million species cannot be drawn [START_REF] Mayhew | Shifts in hexapod diversification and what Haldane could have said[END_REF].

In this paper, I present a method to estimate speciation and extinction rates which combines phylogenetic and species richness data. Some examples using birds and mammals are presented.

METHODS

The approach used here is an extension of the birth-death model developed by Nee et al. (1994a). They used a birth-death process to model phylogenetic diversification and estimate speciation and extinction rates with the branching times from an observed phylogenetic tree. Basically, this method requires that all species belonging to the studied clade are included in the reconstructed tree (Nee et al. 1994a,b).

The present approach assumes that two kinds of information are available.

First, the phylogenetic relationships (tree topology and branch lengths) of a subset of all species belonging to the studied clade are known; this is referred to as the phylogenetic data in this paper. It is further assumed that the reconstructed phylogenetic tree is ultrametric, that is the branch lengths are linearly related to time. Second, the number of species belonging to each tip of this reconstructed phylogeny is known as well, and this is called the taxonomic data. This is likely to be a very common situation. For instance, it is usual in a phylogenetic study to resolve the relationships among higher level taxa (e.g., among orders within a class, or among families within an order) and that the number of species in each of these taxa is known. It is assumed below that both the phylogenetic and the taxonomic data are known without error.

Denote N the number of tips in the reconstructed phylogeny, and n i the number of species belonging to the ith tip (i = 1, . . . N ).

Let us assume that a living species has an instantaneous probability of splitting in two species called the rate of speciation and denoted λ, and an instantaneous probability of disappearing called the rate of extinction and denoted µ. We further assume that both rates are constant through time and across lineages. Under this "birth-and-death" model, we can formulate the probabilities of the observed data using formulae in [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF] and in Nee et al. (1994a).

The phylogenetic data are comprised by a set of N lineages that originated from the nodes of the reconstructed phylogeny and survived until present. The probability for a lineage originated at time t of not going extinct at time T is (Nee et al. 1994a, eq. 2): -t) .

P (t, T ) = λ -µ λ -µe -(λ-µ)(T
(1)

The probability density of the branching events in the observed tree is proportional to the number of extant lineages multiplied by λ. It is thus possible to compute the likelihood of the phylogenetic data by multiplying together the probabilities of the observed branching events and the probabilities that the lineages survived until present. This likelihood is:

L P = (N -1)!λ N -2 N i=3 P 2 (t i , T )P 2 (t 2 , T ). (2)
This is actually a modified likelihood from Nee et al. (1994a, their eq. 20) to take into account the fact that the lineages survived with one or more species untill present.

Let us rewrite equation 1 as:

P (t, T ) = 1 -a 1 -ae -rx , (3)
where a = µ/λ, r = λ -µ, and x = T -t. If T is the present, x is in fact the branching time. Some algebra lead us to rewrite equation 2 as:

L P = (N -1)!r N -2 (1 -a) N exp 2r N i=2 x i N i=2 (e rx i -a) -2 . ( 4
)
A logarithmic transformation of the above equation leads to an expression easier to manipulate:

ln L P = ln(N -1)! + (N -2) ln r + N ln(1 -a) + 2r N i=2 x i -2 N i=2
ln(e rx i -a).

(5)

The taxonomic data are represented by a set of monophyletic groups of species which numbers and dates of origin (the lengths of the terminal branches of the reconstructed tree) are both known. The probability for a lineage originated from a single species to have n species after a time t is (Kendall 1948, eq. 16):

µ(e (λ-µ)t -1) λ(e (λ-µ)t -µ) if n t = 0, (6) 1 - µ(e (λ-µ)t -1) λ(e (λ-µ)t -µ) (1 -ζ t )ζ n-1 t if n t > 0, (7)
where

ζ t = λ(e (λ-µ)t -1) λe (λ-µ)t -µ . (8)
The same reparametrization used above for the phylogenetic data (a = µ/λ, r = λ -µ) and an algebraic development lead us to rewrite equation 7 as:

P (n t = n) = (1 -a) 2 e rt (e rt -a) n+1 (e rt -1) n-1 for n ≥ 1. ( 9
)
Multiplying together the corresponding probabilities for the N observed lineages gives the likelihood of the taxonomic data:

L T = (1 -a) 2N N i=1 e rt i (e rt i -1) n i -1 (e rt i -a) -(n i +1) , (10)
and after logarithmic transformation:

ln L T = 2N ln(1 -a) + r N i=1 t i + N i=1 (n i -1) ln(e rt i -1) - N i=1 (n i + 1) ln(e rt i -a). ( 11
)
The log-likelihood of the combined phylogenetic and taxonomic data (denoted ln L) is the sum ln L P + ln L T . Given some data, finding the maximum of ln L over a and r gives their respective maximum likelihood estimates (MLEs) denoted â and r. Under the assumption that ln L is approximately normally distributed around its maximum, it is possible to estimate the standard-errors of the MLEs using its second derivatives.

Practically, maximizing ln L would require us to find its first partial derivatives with respect to a and r; however these expression are too complex to be solved analytically. I used instead a numerical minimization method. The likelihood was transformed as the deviance (equal to -2 ln L); finding the maximum likelihood was equivalent to minimizing the deviance. This was done with the nonlinear minimization function of R [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF] which uses a Newton-type algorithm for unconstrained minimization [START_REF] Schnabel | A modular system of algorithms for unconstrained minimization[END_REF]. This method allows one to find the minimum of a function even when its partial and second derivatives are unknown. However, these derivatives can be computed numerically, allowing to calculate the standard-errors of the MLEs.

If the normal approximation of the likelihood surface is not satisfied (e.g., due to asymmetries), the standard-errors of the MLEs based on the second derivatives of the log-likelihood may be incorrect. An alternative approach is to use the profile likelihood [START_REF] Hudson | Interval estimation from the likelihood function[END_REF]): this consists of looking at the likelihood surface around its maximum. Since a difference between two log-likelihoods less than or equal to 1.92 is not significantly different with P > 0.05 [START_REF] Edwards | Likelihood[END_REF]), all likelihood values greater than or equal to the maximum likelihood minus 1.92 define a 95 % confidence region for both parameters a and r. This is best visualized by plotting the likelihood surface for different values of the two parameters. The confidence intervals of a and r can then be transformed in confidence intervals of λ and µ by developing the inequalities defining these intervals (see APPENDIX 1).

The programs implementing the method developed here are freely distributed in the package APE (analysis of phylogenetics and evolution, [START_REF] Paradis | APE: an R package for analyses of phylogenetics and evolution[END_REF].

EXAMPLES

Birds

The extensive study of [START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF] provided an estimate of the phylogenetic relationships for more than 1000 species among the 9646 living bird species based on DNA/DNA hybridization experiments (their famous "tapestry").

This covered all avian orders, and almost all families. [START_REF] Sibley | Distribution and taxonomy of birds of the world[END_REF] made a review of the taxonomy of birds, and listed the number of species for all taxonomic levels, thus providing the numbers of species for the 23 orders and 144

families of living birds.

The goal of the present analysis was to estimate diversification rates for birds at the order and family levels. The data for 137 families were considered. The seven families not included here were groups with uncertain relationships [START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF]) and contained very few species: Brachypteraciidae (5 species), Raphidae (3), Mesitornithidae (3), Philepittidae (4), Callaeatidae (3), Picathartidae (4), and Hypocoliidae (1). The resulting tree was not fully dichotomous and contained one trichotomy (see Sibley & Ahlquist 1990, Fig. 356).

This was resolved by assuming a null branch length below this node, so that the corresponding branching time was counted twice in the analysis. The original units of Sibley & Ahlquist (∆T 50 H, a measure of the distance among the tips of the tree) were used here for the branch lengths.

Twenty three orders were considered: the resulting phylogeny was fully dichotomous (Fig. 1). As for families, the number of species were from Sibley & Monroe 1990.

The parameter estimates for the bird family data were â = 0, and r = 0.401 (s.e. = 0.006). The maximum log-likelihood was -427.45. A plot of the log-likelihood around its maximum shows a 95 % confidence interval of [0, 0.008] for a, and [0.385, 0.418] for r (Fig. 2a). The 95 % confidence intervals for speciation and extinction rates (in ∆T 50 H time units) were 0.385 ≤ λ ≤ 0.421, and 0 ≤ µ ≤ 0.003.

The same estimates for the analysis at the order-level were â = 0, and r = 0.287 (s.e. = 0.007); the maximum log-likelihood was -144.58. The log-likelihood surface indicates a 95 % confidence interval of [0, 0.040] for a, and [0.268, 0.308] for r (Fig. 2b). The 95 % confidence intervals for speciation and extinction rates (in ∆T 50 H time units) were 0.268 ≤ λ ≤ 0.321, and 0 ≤ µ ≤ 0.013.

It is noteworthy that the confidence intervals for λ estimated from family and order data do not overlap. This is discussed below.

Eutherian Mammals

Wilson & Reeder (1993) listed 26 orders of mammals together with a complete taxonomic listing at the species level. On the other hand, there is, for the moment, no complete phylogeny with branch lengths of the orders of mammals. Douzery et al. (in press) reconstructed a phylogeny with dated nodes for 18 orders of mammals based on molecular data for 42 species. They used a nonparametric method which does not assume a molecular clock [START_REF] Sanderson | A nonparametric approach to estimating divergence times in the absence of rate constancy[END_REF] to estimate the divergence dates of the tree. I discarded the data for the order Diprotodontia (represented by the kangaroo Macropus) which was the only non-Eutherian order in this study. I used data for 17 orders among which the species richness of Cetacea (78 species)

was pooled with that of Artiodactyla (220 species), thus leaving 16 orders. I considered only Eutherian mammals since the divergence between this clade and the other mammals is not yet well dated (Douzery pers. comm.).

The tree was not fully dichotomous and contained one trichotomy (Fig. 3).

This was resolved in the same way than for the phylogeny of avian families. The unit of branch lengths was in million years ago (Ma) calculated from the ages of the nodes reported by Douzery et al. (in press).

The parameter estimates for this mammalian phylogeny were â = 0, and r = 0.080 (s.e. = 0.003). The maximum log-likelihood was -102.21. A plot of the log-likelihood around its maximum shows a 95 % confidence interval of [0, 0.050] for a, and [0.074, 0.088] for r (Fig. 4). The 95 % confidence intervals for speciation and extinction rates are 0.074 ≤ λ ≤ 0.093, and 0 ≤ µ ≤ 0.005.

DISCUSSION

There is great interest in using phylogenetic data to study macroevolutionary processes as shown by recent reviews [START_REF] Sanderson | Reconstructing shifts in diversification rates on phylogenetic trees[END_REF][START_REF] Barraclough | Phylogenetics and speciation[END_REF]). Significant progress has been accomplished when the available phylogeny is complete (i.e., has been reconstructed with all species; [START_REF] Slowinski | Testing whether certain traits have caused amplified diversification: an improved method based on a model of random speciation and extinction[END_REF][START_REF] Harvey | Phylogenies without fossils[END_REF]. However, the development of approaches for analyzing incomplete phylogenies has been more difficult [START_REF] Paradis | Assessing temporal variations in diversification rates from phylogenies: estimation and hypothesis testing[END_REF](Paradis , 1998b;;[START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF]. The method presented in this paper is an attempt to combine phylogenetic and taxonomic data to estimate speciation and extinction rates, and thus can be used to analyse incomplete phylogenies provided some data on species richness are available.

The present method can be viewed as an extension of the method introduced by Nee et al. (1994a). Both methods are based on the birth-and-death process studied by [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF]. The main distinction between both methods, particularly apparent in the computation of their respective likelihoods, is that Nee et al. (1994a) consider that all surviving lineages have one species living at present, whereas the present method simply assumes that these lineages survive (with one or more species).

It is interesting to point out that both approaches aim to estimate parameters rather than to test hypotheses. Nee et al.'s (1994a) method can be extended to more general birth-death models (e.g., assuming time-dependence in both parameters, see Nee et al. 1994a) but, as far as I know, this does not seem to have been used in practice. Whether the approach presented here could be extended to more complex birth-death models requires more work. It would be interesting in the future to compare extensively the estimates from Nee et al.'s (1994a) method and the present one.

In a previous work, I developed a method to take into account missing species in phylogenetic data when estimating diversification rates [START_REF] Paradis | Assessing temporal variations in diversification rates from phylogenies: estimation and hypothesis testing[END_REF](Paradis , 1998b)). This is based on making an analogy between branching times in a phylogeny and survival events: thus missing branching times can be treated as censored data (i.e.

when death is not known precisely but it is known that it occurred after some time). This approach has the advantage of making available likelihood-based fitting of alternative models, and thus testing alternative hypotheses on variation in diversification rates [START_REF] Paradis | Assessing temporal variations in diversification rates from phylogenies: estimation and hypothesis testing[END_REF](Paradis , 1998b)). However, this has the inconvenient property that when the proportion of missing data is too large there is not enough information to allow efficient estimation. Thus, the approach based on survival models is not applicable to the analysis of the kind of data analysed in the above examples.

Though the approach presented here aims to analyse incomplete phylogenies, the data considered (branch lengths, tree topology, and number of species) are assumed to be known without error. This assumption is likely to be untrue in some situations since the tree is usually an estimate and species numbers in some groups may be inacurrate particularly in speciose or poorly-known groups. The maximum likelihood approach developed here could be, at least in theory, extended to take into account such uncertainties but this would certainly require complex mathematical developments. Furthermore, it is likely that in most situations the uncertainty in the data reflects conflicts between a limited number of scenarios (such as the attribuation of a sub-group to a particular sub-clade) and thus it is possible to consider alternatively all of them. [START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF] proposed a method to assess the impact of incompleteness in phylogenetic data on testing for constant diversification. This is based on simulating phylogenies with constant rates then sampling randomly a subset of the species "living" at the end of the simulation. This Monte Carlo approach makes it possible to estimate the distribution of their test statistic (the γ-statistic) in the presence of missing data. Nee et al. (1994b) used simulations to assess the effect of missing data on analysis of diversification with a graphical method: the lineages-through-time plot [START_REF] Harvey | Phylogenies without fossils[END_REF][START_REF] Nee | Inferring population history from molecular phylogenies[END_REF]. They showed that the analysis of incomplete phylogenies may indicate a spurious decline through time of diversification (Nee et al. 1994b). There is still an obvious need for methods that would tease apart the effects of incomplete sampling from those of temporal variation in phylogenetic diversification rates.

The analyses presented in this paper have mainly an illustrative purpose, but the results obviously call for a few comments. We will first convert the rates estimated for birds in meaningful units (Ma -1 ). The depth of the bird tree of [START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF] is 28 in ∆T 50 H units. [START_REF] Paton | Complete mitochondrial DNA genome sequences show that modern birds are not descended from transitional shorebirds[END_REF], using various methods, dated the origin of the modern bird orders between 110 and 130 Ma depending on the method. We need thus, roughly, to multiply [START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF] time-scale by 4 to obtain a time-scale in Ma. Therefore, the corresponding rates converted in Ma -1 (after dividing the above speciation and extinction rates by 4) are 0.096 ≤ λ ≤ 0.105, and 0 ≤ µ ≤ 0.001 for the family-based analysis, and 0.067 ≤ λ ≤ 0.080, and 0 ≤ µ ≤ 0.003 for the order-based one.

Concerning birds, the analysis at the level of orders yielded smaller estimates and larger confidence intervals than the same analysis at the level of families.

There could be two explanations for this discrepancy. First, the family-based data can be thought as more accurate than the order-based ones. Consequently, the analysis of the former may give more accurate (with narrower confidence intervals)

estimates than for the latter. Second, heterogeneity in rates may affect the present estimates. The present method assumes that speciation and extinction rates are constant through time and the same for all lineages. This is certainly not true for the present data since species richness is extremely uneven among bird families, with some Passerines being extremely diverse (e.g., Fringillidae). It could be that at the level of orders, the effect of heterogeneity was lessened since the heterogeneity among orders may be smaller than among families. In other words, the family-based estimates may have been "pushed-up" by the many diversified families of Passerines.

The estimates for mammals are close to those for birds, particularly for the order-based analysis of the latter group. This is in agreement with the slightly younger origin of the modern orders of Eutherian mammals (ca. 90 Ma) compared to birds (ca. 110 Ma). The Eutherian mammals have 4260 species (Mammals total 4628 species according to Wilson & Reeder 1993) against 9623 for the birds analysed here. However, we need to keep in mind that rates are certainly heterogeneous in both groups: rodents, like passerines for birds, have surely diversified at a higher rate than other lineages of mammals (see [START_REF] Purvis | Macroevolutionary inferences from primate phylogeny[END_REF]Paradis 1998b for analyses of heterogeneity in diversification rates among lineages of primates; and Mooers & Heard 1997 for a review). However, the present analyses suggest that, on average, birds and mammals have diversified at similar rates.

The three analyses reported in this paper yielded an estimate of a equal to zero. It is obvious that some extinction events occurred during the history of mammals and birds as evidenced by the fossil record. However, the estimation of extinction rates with phylogenies without fossils is still an open issue. For instance, [START_REF] Purvis | Macroevolutionary inferences from primate phylogeny[END_REF] reported estimates of a equal to zero in primate lineages in two cases out of four, and Nee et al. (1994b) reported similar instances in other groups. The method presented in this paper can, in theory, yields estimates of a greater than zero (APPENDIX 2). One referee of the present paper pointed out that there is currently some confusion in the literature with respect to the issue of estimating extinction rates with phylogenies without fossils (e.g., [START_REF] Sims | Nonstochastic variation of species-level diversification rates within angiosperms[END_REF]. In the absence of an extensive theoretical treatment of this issue, it seems sound to take the estimates of extinction rates inferred from phylogenies without fossils (complete or not) with great care. An interesting point would be to assess whether the confidence interval for the extinction rate inferred with the profile likelihood correctly covers the true value of µ.

The method presented in this paper has surely a wide range of potential applications. Future research will need to address the issues of temporal variation and heterogeneity in rates. Ultimately, linking the methods dealing with phylogenetic diversification and the methods dealing with the evolution of characters, such as the various phylogenetic comparative methods, should be a significant progress in evolutionary analysis.

APPENDIX 1

This appendix shows how the confidence intervals of the parameters a and r are transformed to give confidence intervals of the speciation and extinction rates (λ and µ, respectively).

Suppose the profile likelihood of the model gives the two 95 % confidence intervals:

a 1 ≤ a ≤ a 2 , (12)
r 1 ≤ r ≤ r 2 . ( 13
)
Since we have a = µ/λ, and r = λ -µ, and under the condition that λ is positive nonnull, we can re-write both inequalities as:

a 1 λ ≤ µ ≤ a 2 λ, (14)
r 1 + µ ≤ λ ≤ r 2 + µ. (15)
Thus we have for λ the inequalities:

r 1 + a 1 λ ≤ λ ≤ r 2 + a 2 λ, (16)
which can be written separately as:

r 1 + a 1 λ -λ ≤ 0, 0 ≤ r 2 + a 2 λ -λ. (17)
Developing algebraically both inequalities gives (under the condition that both a 1 and a 2 are less than one):

r 1 1 -a 1 ≤ λ, r 2 1 -a 2 ≥ λ, (18)
giving the 95 % confidence interval for λ:

r 1 1 -a 1 ≤ λ ≤ r 2 1 -a 2 . ( 19
)
Substituting in the above inequalities of µ gives its 95 % confidence interval:

a 1 r 1 1 -a 1 ≤ µ ≤ a 2 r 2 1 -a 2 . ( 20

I

  am grateful to three anonymous referees for their constructive comments on an earlier version of this paper, and to E. Douzery for sending me his results on divergence dates on mammals. This research was financially supported by the Institut Français de la Biodiversité and the Centre National de la Recherche Scientifique. This is publication 2003-052 of the Institut des Sciences de l' Évolution (Unité Mixte de Recherche 5554 du Centre National de la Recherche Scientifique).

Figure 1 :Figure 2 :

 12 Figure1: The phylogenetic relationships among the orders of birds according to[START_REF] Sibley | Phylogeny and classification of birds: a study in molecular evolution[END_REF]. The numbers on the right-hand side are the number of species in each order according to[START_REF] Sibley | Distribution and taxonomy of birds of the world[END_REF].

Figure 3 :Figure 4 :

 34 Figure 3: The phylogenetic relationships among the orders of birds according to Douzery et al. (in press). The numbers on the right-hand side are the number of species in each order according toWilson & Reeder (1993).

APPENDIX 2

This appendix shows that the method presented in this paper can yield estimates of a greater than zero. This can be shown by examining the first partial derivative of the log-likelihood function ln L. Since ln L = ln L P + ln L T , we have:

Deriving equations 5 and 11 with respect to a is straightforward:

giving for the log-likelihood function:

If the partial derivative of ln L can be positive at a = 0, then ln L can have a maximum for positive values of a. Fixing a = 0 in the above equation gives:

The first two terms above are always positive, whereas the third one is always negative. Thus, the above expression can be positive depending on the values of r,

x i , t i , and N . Consequently, the log-likelihood can potentially achieve a maximum for positive values of a.