Table Of ContentEvolutionary Model of Species Body Mass Diversification
A. Clauset1,∗ and S. Redner2,1,†
1Santa Fe Institute, 1399 Hyde Park Road, Santa Fe NM, 87501, USA
2Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215, USA
We present a quantitative model for thebiological evolution of species body masses within large
groups of related species, e.g., terrestrial mammals, in which body mass M evolves according to
branching(speciation),multiplicativediffusion,andanextinctionprobabilitythatincreaseslogarith-
mically with mass. We describe this evolution in terms of a convection-diffusion-reaction equation
for lnM. The steady-state behavior is in good agreement with empirical data on recent terres-
trial mammals, and thetime-dependent behavioralso agrees with data on extinct mammal species
between 95 – 50 Myr ago.
9
0 PACSnumbers: 87.23.Kg,02.50.-r
0
2
n Animals—both extant and extinct—exhibit an enor- in principle, be applied to any group of related species.
a mously wide range of body sizes. Among extant ter- Letc(x,t)denote the number (density)ofspecieshav-
J
restrial mammals, the largest is the African savannah ing logarithmic mass x=lnM at a time t; we use x as
2 elephant (Loxodonta africana africana) with a mass of thebasicvariableinkeepingwithwidespreadusageinthe
2 107g, while the smallest is Remy’s pygmy shrew (Sun- field [8]. Our model incorporates three fundamental and
cus remyi) at a diminutive 1.8g. Yet the most probable empirically-supported features of biological evolution.
]
E mass is 40g, roughly the size of the common Pacific rat
P (Rattus exulans), is only a little larger than the small- 1. Branching multiplicative diffusion [1, 8]: each
o. est mass. More generally,empirical surveys suggest that species of mass M produces descendant species
i sucha broadbut asymmetric distribution in the number (cladogenesis) with masses λM, where λ is a ran-
b
of species with adult body mass M typifies many ani- dom variable, and the sign of the average hlnλi
-
q malclasses[1,2,3, 4,5],including mammals,birds,fish, denotes bias toward larger or smaller descendants.
[ insects, lizards and possibly dinosaurs. Empirical evidence [8] suggests that hlnλi > 0
(knownasCope’srule[1])forterrestrialmammals.
2 What mechanisms cause species mass distributions to
v assume such a shape? A satisfactory answer would have
4 2. Species become extinct independently, and with a
wide implications for the evolution and distribution of
1 probability that increases weakly with mass [12].
themanyotherspeciescharacteristicsthatcorrelatewith
0
body mass, including life span, metabolic rate, and ex-
4 3. No species can be smaller than a mass M , due,
min
8. tinction risk [6, 7]. Previous explanations for the species for example, to metabolic constraints [13].
massdistributionfocusedondetailedecological,environ-
0
8 mentalandspecies-interactionassumptions[3]. However, The production of descendant species corresponds to
0 empirical data present confusing and often inconsistent growth in the number of species at a rate k that is pro-
v: support for these theories, and none explicitly address portional to the density c itself. Similarly, the prob-
i how species body mass distributions diversify in time. ability p(x) that a species of logarithmic mass x be-
X
InthisLetter,weconstructaphysics-basedconvection- comes extinct may also be represented by a loss term
r
diffusion-reaction model to account for the evolution of that is proportional to c, but with a weak mass depen-
a
the species mass distribution. The steady-state behav- dence [8]. We make the simple choice of linear depen-
ior of this model was recently solved to explain the dence p(x) = A+Bx (but see below). With these in-
species mass distribution for recent terrestrial mammals gredients, c(x,t) obeys the convection-diffusion-reaction
andbirds[8,9],whererecentisconventionallydefinedas equation in the continuum limit
within the past 50,000 years [10]. Here we substantially
extend this approach to give predictions on mammalian ∂c ∂c ∂2c
+v =D +(k−A−Bx)c , (1)
bodymassevolutionthatareingoodagreementwithfos- ∂t ∂x ∂x2
sil data. This model can further be used to estimate the
with bias velocity v = hlnλi and diffusion coefficient
historical rates of body mass diversification from fossil
D =h(lnλ)2i, and where k−A sets the absolute scaleof
data, which are otherwise estimated using ad hoc tech-
species body mass frequencies.
niques. To illustrate this application, we estimate body
To solve Eq. (1), we substitute the eigenfunction ex-
mass diversificationrates from our fossil data, which are
in good agreement with estimates of genetic diversifica- pansion c(x,t)= nAnCn(x)e−γnt, which yields
P
tion from molecular clock methods [11]. Although our γ
focusisonterrestrialmammalevolution,thismodelcan, − DnCn+µCn′ =Cn′′+(α−βx)Cn , (2)
2
−1 6
10 10
KP boundary
5
10
10−2 g)104
y s (
nsit mas103
De y
d
10−3 Bo102
1
10
Modern terrestrial mammals
Model (steady state)
10−4 100
100 101 102 103 104 105 106 107 108 100 90 80 70 60 50
Species body mass (g) Millions of years ago
FIG. 1: (color online) Steady-state solution of the model FIG. 2: (color online) Data on species body mass for 569
for the species body mass distribution and suitably binned North American terrestrial mammals [16] from 95 – 50 Myr
empiricaldatafor4002recentterrestrialmammals(shownas ago. Each horizontal segment represents a species, and end-
anormalizedhistogram with50logarithmically-spaced bins). pointsdenoteitsfirstandlastappearanceinthefossilrecord.
The superimposed curve shows the average of lnM for these
species (smoothed with an exponential kernel).
where the prime denotes differentiation with respect
to x, µ = v/D, α = (k − A)/D, and β = B/D.
result into the definition of z, we can eliminate the scale
We eliminate the first derivative term by introducing
ψ = e−µx/2C and then we use the scaled variable parameter α and write z =zn+β1/3(x−xmin). Thus
n n
z =β1/3x−β−2/3 α− µ2 + γn to transform Eq. (2) each eigenfunction has the form
4 D
(cid:16) (cid:17)
into Airy’s differential equation, ψ′′ −zψ = 0 [14], for C (x)∝eµx/2Ai z +β1/3(x−x ) , (4)
n n min
each eigenfunction ψ . The general solution is ψ(z) = h i
n
a1Ai(z)+a2Bi(z), where Ai(z) and Bi(z) are the Airy in which Ai(z) ∼ e−2z3/2/3 for large z. The competition
functions; here the prime now denotes differentiation betweenthis decayandthe prefactoreµx/2 in C (x) con-
n
with respect to z. Since there can be no species with in-
tributestothebroadnessofthespeciesmassdistribution
finite mass and Bi(z) diverges as z →∞, we set a =0,
2 andthelocationofthemostprobablemass(Fig.1). Par-
while a is determined by the initial condition. (We
1 enthetically, the asymptotic decay of the eigenfunctions
could also incorporate a fixed maximum species body
depends only weakly on the form of the extinction prob-
massMmax,buttheanalysisismorecomplicatedwithout ability p(x). For instance, if we choose p(x)=A+Bxδ,
revealing additional insights.) then as z →∞, the eigenfunctions decay as e−z1+δ/2.
Since the species density vanishes at the minimum
Supposethatagivengroupofanimalsbeganitsevolu-
mass M , the argument of ψ must equal
min n tionary history with a single species of mass M (initial
0
µ2 γ γ condition c(x,t = 0) = δ(x − x0), with x0 = lnM0),
zn=β1/3xmin−β−2/3(cid:18)α− 4 +Dn(cid:19)=z0−Dβn2/3 . after which speciation occurs according to the dynamics
of Eq. (1). We use the fact that the Airy differential
(3)
equation is a Sturm-Liouville problem [15] so that the
The functions ψn = Ai(zn) form the complete set of {Cn} form a complete and orthonormal set. Following
states for the eigenfunction expansion[15]. The first few the standard prescription to determine the coefficients
zerosz areat(roughly)−2.3381,−4.0879,−5.5205and of the eigenfunction expansion, the full time-dependent
n
−6.7867for n=0,1,2,3[14]. (For computing the distri- solution is
bution, we tabulated the first million zeros numerically
usingstandardmathematicalsoftware.) Thecorrespond- c(x,t)= Cn(x0)Cn(x)e−γnt , (5)
ingdecayratesγ arethengivenbyγ =Dβ2/3(z −z ), Xn
n n 0 n
with γ = 0 to give the steady state solution [9]. These where each C (x) is given by Eq. (4). In the long-time
0 n
rates form an increasing sequence so that the higher limit, all the decaying eigenmodes with n > 0 become
termsinthe eigenfunctionexpansiondecaymorequickly negligible and the species mass distribution reduces to
in time. Finally, solving Eq. (3) for α and plugging the c(x,t→∞)∝C (x).
0
3
22
95−90 Mya, t=4.92 90−85 Mya, t=6.72 85−80 Mya, t=6.72 KP boundary
20
y
sit 18
n
e
D 16
14
80−75 Mya, t=7.22 75−70 Mya, t=8.63 70−65 Mya, t=11.74 me
y el ti12
sit od10
n m
De 8
6
65−60 Mya, t=15.45 60−55 Mya, t=17.96 55−50 Mya, t=20.47 4
sity 2
n
e 0
D 100 90 80 70 60 50
Millions of years ago
log mass log mass log mass
10 10 10
FIG.3: (coloronline)Acomparisonofourmodelpredictions FIG. 4: Estimated relation between model time to histori-
(dashed) from Eq. (5) and the mass distributions of extinct cal time, indicating that the broadening of the mammalian
North American mammal species (solid) in nine consecutive speciesmassdistributionbeganapproximately10Myrbefore
timeranges. (Toget reasonable results withsparse empirical theKP boundarywhen non-avian dinosaurs became extinct.
data, we smooth thedistributions with a Gaussian kernel.)
from Eq. (5) and empirical data from Fig. 2. To sim-
We test our predictions for the species body mass dis- plify the comparison,wedividedthe periodfrom95–50
tributionbycomparingwithavailableempiricaldata. In Myragointonine binsofequaldurations,andtabulated
the steady state, our model is characterized by three the distribution of species extant during each of these
parameters: Mmin, the mass of the smallest animal, µ, bins. In each of the nine panels of Fig. 3, we give both
which controls the tendency of descendant species to be the historical time period and the corresponding model
largeror smaller than their ancestors,and β, which con- time that yields a good qualitative match to the data.
trolsthe dependence ofthe extinctionrateonmass. The (The fitting can be made more objective using standard
former two can be estimated directly from fossil data, techniques, but the results are largely the same.)
while the latter is typically estimated by matching the Therelationbetweenhistoricaltimeandmodeltimeis
steady-state solution to the recent data. For terrestrial itself interesting (Fig. 4). During the first 20 Myr (95 –
mammals, we previously found µ≈0.2, β ≈ 0.08, while 75 Myr ago),the species mass distribution is almoststa-
Mmin ≈2g[9]. Usingtheseparametervaluesinthelong- tionary, with model time advancing by only ∆t = 2.30;
time limit, weobtaina goodagreementbetweenthe pre- curiously, between 90 – 80 Myr ago, model time does
dictions of the model and the species mass distribution not advance at all. This period of near stasis may indi-
of recent terrestrial mammals [10] (Fig. 1). cate a lull in the evolutionary dynamics, perhaps due to
Our model also makes predictions about the way the implicit competition from larger species, e.g., dinosaurs.
body mass distribution changes over time, which can be Over the next 20 Myr (75 – 55 Myr ago), however, the
tested with fossil data. Drawing on data from the best distribution broadens considerably (model time advanc-
available source on the evolution of mammalian body ing by ∆t = 10.74), and comes to closely resemble the
masses [16], we plot in Fig. 2 the durations and body recent distribution (Fig. 1). This correspondence sug-
masses of 569 extinct mammal species from 95 – 50 Myr geststhatthe diversificationofmammalianbody masses
ago. ThisperiodincludestheCretaceous-Paleogene(KP) intotheirrecentstatebeganatleast75Myrago,roughly
boundary of 65.5 Myr ago that marks a mass extinc- 10 Myr before the KP boundary and the extinction of
tion event during which more than 50% of then extant the non-avian dinosaurs. This estimate of the timing
species became extinct, including non-avian dinosaurs— of body mass diversification for mammals agrees closely
the dominant fauna for the preceding 160 Myr—and is with some estimates of the timing of mammalian ge-
the subject of many studies regarding the diversification neticdiversificationfromstudiesofmolecularclocks[11],
of mammals. (The Cretaceous period is conventionally andsupportsthe notionthatmammalswerediversifying
abbreviated “K,” after the German translation Kreide.) priortoandindependentlyoftheKPboundaryitself[18].
Usingthesamemodelparametersasabove,andsetting Whether these two forms of diversification are causally
M = 2g, the estimated size of the first mammal [17], linked, however, is unknown.
0
Fig. 3 shows good agreement between model predictions Herewehaveheldallmodelparametersconstantwhile
4
adjusting model time to fit the data. In principle, how- total metabolic flux of a taxonomic group.
ever, model time could advance steadily while varying We thank D. H. Erwin, D. J. Schwab, Z.-X. Luo and J.
some model parameters, perhaps to reflect large-scale Okie for helpful comments, and J. Alroy, A. Boyer and
changes or trends in the selection pressures on species F.Smithforkindlysharingdata. SRgratefullyacknowl-
body size. Empirical evidence supports a stable value of edges support from NSF grantDMR0535503. This work
M (Fig. 2), but little is currently known about how was also supported in part by the Santa Fe Institute.
min
or why µ and β may have varied.
As is typically the case with historical inference using
the fossil record, a few caveats are in order. Our fossil
data are derived from the well-studied North American
∗ Electronic address: [email protected]
region using modern dental techniques, which are less
† Electronic address: [email protected]
prone to biases than older techniques. Still, some biases
[1] S. M. Stanley,Evolution 27, 1 (1973).
and sampling gaps likely persist, and these may explain
[2] J. Kozl owski and A. T. Gawelczyk, Functional Ecology
the slight overabundance of large species, and under- 16, 419 (2002).
abundanceofsmallspecies,inmorerecenttimes(Fig.3). [3] C. R. Allen, A. S. Garmestani, T. D. Havlicek, P. A.
More significantly, recent fossil discoveries suggest that, Marquet, G. D. Peterson, C. Restrepo, C. A. Stow, and
since mammals originated195Myr ago[17], mammalian B. E. Weeks, Ecology Letters 9, 630 (2006).
[4] M.T.Carrano,inAmniotePaleobiology,editedbyM.T.
diversification has proceeded in several waves, and the
Carrano, T. J. Gaudin, R. W. Blob, and J. R. Wible
vast majority of species groups from the earlier waves
(University of Chicago Press, 2006), pp.225–268.
arenowextinct. Ourdatacoveronlythe mostrecentdi- [5] S. Meiri, Global Ecology and Biogeography 17, 724
versification, in which therian (placental and marsupial) (2008).
mammalslargelyreplacedthethendominantnon-therian [6] W. A. Calder III, Size, Function, and Life History
mammal groups [19]. Unfortunately, suitable data on (Dover, Mineola, NY,1984).
these waves of diversification is not currently available, [7] M.Cardillo,G.M.Mace,K.E.Jones,J.Bielby,O.R.P.
Bininda-Emonds, W. Sechrest, C. D. L. Orme, and
and the data we do have is sparse in its coverageof non-
A. Purvis, Science 309, 1239 (2005).
therians. Thus, the application of our model to infer
[8] A. Clauset and D. H.Erwin, Science 321, 399 (2008).
diversification rates should be considered as a proof-of-
[9] A.Clauset,D.J.Schwab,andS.Redner,AmericanNat-
concept, illustrating that a physics-style model can shed uralist 173, 256 (2009).
considerable light on evolutionary dynamics by placing [10] F. A. Smith, S. K. Lyons, S. K. M. Ernest, K. E. Jones,
the fossil record within a theoretical framework. D.M.Kaufman,T.Dayan,P.A.Marquet,J.H.Brown,
In summary, the broad distribution of body masses andJ.P.Haskell,Ecology84,3403(2003),MOMVersion
3.6.1; This data covers 4002 terrestrial mammal species
for mammals appears to be well described by a simple
from the last 50,000 years, and includes extinct species
convection-diffusion-reaction model that incorporates a
to correct anthropogenic extinctions.
small number of evolutionary features and constraints.
[11] M. S. Springer, W. J. Murphy, E. Eizirik, and S. J.
Indeed, our model does not account for many canoni- O’Brien,Proc.Natl.Acad.Sci.(USA)100,1056(2003).
cal ecological and microevolutionary factors, such as in- [12] L. H. Liow, M. Fortelius, E. Bingham, K. Lintulaakso,
terspecific competition, geography, predation, and pop- H. Mannila, L. Flynn, and N. C. Stenseth, Proc. Natl.
ulation dynamics [3]. The fact that our model agrees Acad. Sci. (USA)105, 6097 (2008).
[13] G. B. West, W. H. Woodruff, and J. H. Brown, Proc.
withspeciesmassdatasuggeststhatthecontributionsof
Natl. Acad.Sci. (USA)99, 2473 (2002).
the above-mentionedprocessesto the globalcharacterof
[14] M. A.Abramowitz andI.E. Stegun,Handbook of Math-
body mass distributions can be compactly summarized
ematical Functions with Formulas, Graphs, and Mathe-
by the parameters µ and β in our model. Despite the matical Tables (Dover, 1972), 9th ed.
crudeness of the model, the agreement between its pre- [15] R. Courant and D. Hilbert, Methods of Mathematical
dictions and the available empirical data is satisfying. Physics (Interscience Publishers, 1989).
Our model opens up intriguing directions for theo- [16] J. Alroy (2008), Paleobiology Database Online System-
atics Archive3, http://paleodb.org/.
retical descriptions of evolutionary dynamics. For in-
[17] Z.-X.Luo,A.W.Crompton,andA.-L.Sun,Science292,
stance, our model ignores the populations of individuals
1535 (2001).
within each species; estimating the sizes of these pop- [18] J. Alroy,Systematic Biology 48, 107 (1999).
ulations from body mass vs. population density scaling [19] Z.-X. Luo, Nature450, 1011 (2007).
relationships [20] and body mass vs. home-range size re- [20] E. P. White, S. K. Morgan Ernest, A. J. Kerkhoff, and
lationships [21] may allow us to calculate both the to- B. J. Enquist, Trends in Ecology and Evolution 22, 323
tal biomass contained in a group of related species, and (2007).
[21] J.H.Brown,G.C.Stevens,andD.M.Kaufman,Annual
its temporal dynamics during diversification. Similarly,
Review of Ecology and Systematics 27, 597 (1996).
when paired with scaling relations between body mass
and metabolism [13], we may be able to calculate the
