S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
M A S S E X T I N C T I O N
Thresholds of catastrophe in the Earth system
Daniel H. Rothman
The history of the Earth system is a story of change. Some changes are gradual and benign, but others, especially
those associated with catastrophic mass extinction, are relatively abrupt and destructive. What sets one group
apart from the other? Here, I hypothesize that perturbations of Earth’s carbon cycle lead to mass extinction if they
exceed either a critical rate at long time scales or a critical size at short time scales. By analyzing 31 carbon isotopic
events during the past 542 million years, I identify the critical rate with a limit imposed by mass conservation.
Identification of the crossover time scale separating fast from slow events then yields the critical size. The modern
critical size for the marine carbon cycle is roughly similar to the mass of carbon that human activities will likely
have added to the oceans by the year 2100.
Copyright © 2017
The Authors, some
rights reserved;
exclusive licensee
American Association
for the Advancement
of Science. No claim to
original U.S. Government
Works. Distributed
under a Creative
Commons Attribution
NonCommercial
License 4.0 (CC BY-NC).
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
INTRODUCTION
Five times in the Phanerozoic [the past 542 million years (My)], more
than three-fourths of marine animal species have vanished in mass
extinctions (1). Each of these events is associated with a significant
change in Earth’s carbon cycle (2, 3). However, there are a great many
disturbances of the carbon cycle (3, 4) whose imprint in the geologic
record is almost exclusively environmental. These other events are
relatively benign, with extinction rates barely, if at all rising, above back-
ground levels (5–7).
Two well-studied examples illustrate these distinctions. The end-
Permian extinction [~252 million years ago (Ma)], the most severe
mass extinction in the Phanerozoic (8), plays out over a period of 104
to 105 years; the extinction interval immediately follows a perturbation
of the carbon cycle of similar duration (9). The Paleocene-Eocene
Thermal Maximum (~55.5 Ma) is a carbon cycle event of roughly simi-
lar time scale, with unambiguous signs of global warming and ocean
acidification (10, 11). It is associated with a significant extinction of ben-
thic foraminifera (12, 13), but there is no evidence of other major ex-
tinctions (11, 13, 14). What makes these two events so different? The
unambiguous co-occurrence of the end-Permian extinction with massive
volcanism (15 ) provides one indication, as does the greater magnitude
of the end-Permian carbon isotopic excursion compared to that of the
Paleocene-Eocene Thermal Maximum (10, 16). However, no general
distinction, applicable throughout the Phanerozoic, exists (5–7).
The idea that mass extinction results from catastrophic environmental
change was largely developed more than two centuries ago by Cuvier
(17). But what defines such an environmental catastrophe? Newell (18)
implicitly provides a practical definition: “Extinction… is not simply
a result of environmental change but is also a consequence of failure
of the evolutionary process to keep pace with changing conditions in the
physical and biological environment” (18). In other words, catastrophe
results when the rate of environmental change exceeds a threshold.
To consider such ideas within the context of the carbon cycle, sup-
pose that the mass m of inorganic carbon in the oceans is increased
from its steady-state value, m*, by an amount Dm over a duration tenv
of time. In terms of the relative change M = Dm/m*, Newell’s
statement hypothesizes that mass extinction results when M/tenv ex-
ceeds a critical rate rc. However, one recognizes immediately that this
condition cannot apply for impulsive change as tenv → 0. In this case,
Lorenz Center, Department of Earth, Atmospheric, and Planetary Sciences, Massa-
chusetts Institute of Technology, Cambridge, MA 02139, USA. Email: dhr@mit.edu
M must instead exceed a finite critical size Mc. Putting it all together,
catastrophe follows when
Mc;
rctenv;
tenv ≪ tx
tenv ≫ tx
ð1Þ
ð1Þ
ð2Þ
M >
where tx is the crossover time scale separating fast from slow change.
Because the asymptotic regions 1 and 2 intersect at tenv = tx, the three
unknowns Mc, rc, and tx can be assumed to satisfy
Mc ¼ rctx
ð3Þ
In principle, Eq. 3 shows how the record of slow environmental
change at geologic time scales informs our understanding of rapid
change at human time scales, even though mechanisms at the different
time scales may differ. Here, I provide an analysis of the geological
record that makes contact with this picture.
RESULTS
I begin by constructing an empirical database of significant global per-
turbations of Earth’s carbon cycle during the Phanerozoic. Each per-
turbation is considered to be an isolated event, or excursion, in the
evolution of the isotopic composition d1 of inorganic (carbonate)
carbon. Except for the most recent event (beginning at the Last Glacial
Maximum), the analysis focuses exclusively on isotopic changes that de-
crease d1. Each isotopic event is parameterized as shown in Fig. 1A.
The database contains each event’s geologic age, the time tenv over
which d1 decreases, and the isotopic compositions d1(0) and d1(tenv)
that bracket the change from time t = 0 to t = tenv. Whenever possible,
estimates of tenv are obtained from isotope geochronology or astro-
chronology. The compilation contains a total of 31 events distributed
throughout the Phanerozoic (Fig. 1B), including each of the so-called
Big Five mass extinctions. The resulting data are contained in Table 1.
Further information is available in Materials and Methods.
Figure 2A displays the isotopic response Dr = d1(tenv) − d1(0) of each
event as a function of its time scale tenv. Two features of the data imme-
diately stand out. First, there is a relative absence of large excursions
at short time scales, suggesting that there may be a characteristic lim-
itation to the rate at which the marine inorganic carbon reservoir can
change its isotopic composition. Second, four of the Big Five mass
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
1 of 12
S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
δ1(0)
τ
δ1( )
T
A
B
–500
–400
–300
Date [My]
–200
–100
0
δ1
r
e
b
m
u
N
7
6
5
4
3
2
1
0
Fig. 1. Parameterization and temporal distribution of carbon isotopic events
in the database. (A) Schematic diagram illustrating how the parameters d1(0) and
d1(tenv) are derived from an isotopic time series d1(t). The time scale tenv repre-
sents the duration of time during which the isotopic composition d1 decreases.
(B) Histogram depicting the distribution of events in time.
extinctions appear to exhibit anomalously fast rates of change. The
fifth mass extinction, in the late Devonian at the Frasnian-Famennian
boundary, is clearly different, lying in the slow region in the lower right-
hand corner (labeled “FF”). McGhee (19) and Bambach et al . (20 ) have
suggested that the Frasnian-Famennian event is not an extinction event
but instead represents a decrease in the rate at which new species orig-
inate. Its relatively small and slow isotopic perturbation provides in-
direct support for this idea.
To understand these observations more deeply, I transform the
geochemical changes to physical fluxes (Materials and Methods). The
transformation assumes that the global carbon cycle is composed of two
reservoirs: inorganic and organic carbon. Carbon is exchanged between
the two reservoirs by photosynthesis and respiration. Inorganic carbon
of isotopic composition di enters the cycle with flux ji. Both inorganic
and organic carbon exit the cycle when buried as rock.
i
=j*
i , where j*
ðtÞ of inorganic carbon with isotopic composition d′
The downturn of d1 is assumed to be driven by a time-dependent
flux j′
i , which is
i
set equal to the average Phanerozoic isotopic composition of organic
carbon [− 27.8‰ (21)]. In other words, j′
i represents an imbalance in
the loop between photosynthesis and respiration on time scales up to
tenv. Its origin is unspecified, but its growth is assumed to be linear.
Taking each event to begin at t = 0 and defining the normalized flux
JðtÞ ¼ j′
i is the unperturbed input flux, one then has
JðtÞ ¼ ft=tenv;
ð4Þ
where f is an unknown, event-specific, dimensionless constant that
specifies the maximum normalized flux J(tenv). An alternative model
could be a constant perturbation J = f. But such a choice would result
in equilibration to a new steady state in the many instances in which
¼ 140 ky(thousand
tenv is greater than the turnover timet0 ¼ m*=j*
years) (Materials and Methods). The analysis here instead assumes
that perturbations grow until t = tenv. The linear model (Eq. 4) can
0 ≤ t ≤ tenv
therefore be viewed as a first-order approximation of a flux function
that satisfies this condition. The normalized mass M in Eqs. 1 and 2
derives from its integral:
Mðf; tenvÞ ¼ 1
t0
∫tenv
0
JðtÞdt ¼ ftenv
2t0
ð5Þ
1
To determine the flux J(t) consistent with each isotopic event, I
calculate the value of f that solves differential equations for the evo-
lution of d1(t), M(t), and J(t) given the initial values d1ð0Þ ¼ d*
1, J(0) = 0,
and M(0) = 0 and the final value d1ðtenvÞ ¼ d*
1, Dr,
and tenv are obtained for each event from the database. The solution of
this boundary value problem for f provides the evolution of d1(t)
that matches its observed initial and final states while remaining
consistent with differential equations describing changes in the isotopic
composition of marine inorganic carbon in response to the flux J(t). It
also provides, via Eqs. 4 and 5, the rate at which J grows and thus the
mass M.
þ Dr. The values of d*
Figure 2B shows M as a function of tenv, rescaled according to Eq. 5,
for each event. The gross features seen in the raw data of Fig. 2A remain,
but now, the appearance of a limiting rate, suggested once again by the
relatively barren upper-left half of the plot, is sharper. Moreover, the
limiting rate now has a physical interpretation. From Eq. 4, one sees
that the normalized flux perturbation reaches its maximum value, J = f,
when t = tenv. The limiting rate therefore represents a critical dimen-
sionless flux fc beyond which the carbon cycle rarely ventures. When
it does, mass extinction often follows.
I hypothesize that fc derives from organic carbon that would nor-
mally be buried as rock but is instead converted to CO2 and remobilized
into the marine carbon cycle. The rate of this remobilization can be no
greater than the rate at which organic carbon is immobilized (that is,
buried) in an unperturbed steady state. During the Phanerozoic, this
limiting rate equals 23 ± 7% of the steady input (or output) flux j*
i
(Materials and Methods). Therefore, fc = 0.23 is the maximum possible
normalized flux perturbation in a mass-conserving carbon cycle with no
anomalous sources or sinks. The straight line in Fig. 2B, which shows
M(fc, tenv) as defined by Eq. 5, represents this hypothesis quantitatively.
This line separates two kinds of carbon cycles. Below the line, per-
turbations can be explained by a partial cessation of organic carbon
burial. Above the line, they cannot. On the line, the small leak repre-
sented by burial is temporarily closed when the perturbation reaches
its maximum strength at t = tenv. Roughly half of the events in the
database fall within the error bars of this boundary. Four of the Big Five
mass extinctions lie above the upper error bar, whereas nearly all other
events fall below it. These observations suggest that noncatastrophic
perturbations represent flux changes along a fixed set of pathways,
whereas catastrophic perturbations result from the introduction of new
pathways acting as anomalous sources or sinks. I refer to the borderline
perturbations as critical events. The preponderance of critical events
supports the idea that large yet benign flux perturbations tend to grow
until their source—the small fraction of organic matter that would nor-
mally be destined for burial—is fully depleted via respiration.
The critical events are likely driven to their limit by the intrinsic
dynamics of the carbon cycle over a time scale set by the characteristic
relaxation time associated with a negative feedback. However, the
variation of tenv over 104 to 106 years would seem to rule out a unique
damping mechanism. This disqualification would necessarily hold if
the carbon cycle did not change over geologic time. But evidently it
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
i
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
2 of 12
S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
Table 1. Event names, event abbreviations, and data used to construct Figs. 2 and 3. Values of d1 and Dr are given in per mil (‰) and derive from the
formula d1 = 1000 × (R − Rstd)/Rstd, where abundance ratios R = 13C/12C are obtained from samples of inorganic (carbonate) carbon and Rstd is a standard ratio.
The variables Dr, tenv, f, and M are defined in the text. In Fig. 2, the six Eocene and four Miocene events are denoted by Eo and Mio, respectively.
Event name
Ediacaran-Cambrian
Nemakit-Daldynian-Tommotian
Cambrian Spice
End-Ordovician
Silurian Mulde
Silurian Lau
Frasnian-Famennian
Tournaisian
Mid-Capitanian
End-Permian
Early Triassic
Triassic-Jurassic
Toarcian
Aptian
Albian-Cenomanian
Mid-Cenomanian
Cenomanian-Turonian
Cretaceous-Paleogene
Early late Paleocene
Paleocene-Eocene Thermal Maximum
Eocene Thermal Maximum 2
Eocene Hyperthermal H2
Eocene Hyperthermal I2
Eocene Hyperthermal I1
Eocene Thermal Maximum 3
Eocene-Oligocene boundary
Miocene Climatic Optimum 1
Miocene Climatic Optimum 2
Miocene Climatic Optimum 3
Miocene Climatic Optimum 4
Last Glacial Maximum–to–Holocene
Abbreviation
Date[Ma]
d1(0)
Edi
NDT
Spice
Ord
Mul
Lau
FF
Tou
Cap
PT
Tr
TJ
Toar
Apt
Al/Cn
mCn
CT
KT
ELPE
PETM
ETM2
H2
I2
I1
ETM3
Oi1
MCO1
MCO2
MCO3
MCO4
LGMH
541.00
524.36
497.00
445.80
428.20
423.50
372.20
351.55
262.45
251.94
251.22
201.64
182.60
120.21
100.50
95.90
94.20
65.50
58.90
55.50
53.70
53.60
53.20
53.10
52.50
33.50
16.90
16.40
16.00
15.60
0.02
2.50
5.00
4.30
4.90
4.00
7.00
3.50
5.60
4.00
4.00
2.00
0.00
2.10
2.95
2.50
3.10
5.30
3.13
3.48
1.90
1.40
1.40
1.56
1.56
0.20
1.60
1.94
1.80
2.33
2.17
0.11
d1(tenv)
−4.00
−3.00
1.30
−0.70
1.00
1.00
2.00
2.10
−1.00
−1.50
−1.50
−2.00
−0.90
1.80
1.80
2.25
3.70
1.98
3.07
−0.80
0.27
0.80
0.95
0.84
−0.60
1.00
1.44
1.04
1.67
1.57
0.45
Dr
−6.50
−8.00
−3.00
−5.60
−3.00
−6.00
−1.50
−3.50
−5.00
−5.50
−3.50
−2.00
−3.00
−1.15
−0.70
−0.85
−1.60
−1.15
−0.42
−2.70
−1.12
−0.60
−0.61
−0.72
−0.80
−0.60
−0.50
−0.76
−0.66
−0.60
0.34
tenv[My]
2.000
0.982
1.000
0.230
0.260
0.500
1.540
1.510
0.500
0.060
0.250
0.050
0.150
0.047
0.110
0.138
0.553
0.026
0.039
0.083
0.045
0.033
0.040
0.040
0.037
0.430
0.028
0.022
0.032
0.029
0.018
f
0.32
0.41
0.12
0.42
0.19
0.30
0.06
0.13
0.27
1.13
0.26
0.49
0.29
0.26
0.08
0.08
0.07
0.44
0.11
0.41
0.28
0.19
0.16
0.19
0.24
0.03
0.18
0.35
0.21
0.21
M
2.302
1.450
0.443
0.344
0.181
0.538
0.311
0.727
0.480
0.243
0.228
0.087
0.157
0.043
0.030
0.038
0.135
0.041
0.015
0.121
0.045
0.023
0.023
0.028
0.032
0.047
0.018
0.028
0.024
0.022
−0.20
−0.013
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
did: Fig. 3A shows that the upper envelope of event size M decreases
with time. To understand this evolution, note that Eq. 5 predicts that
changes in the size of critical events, for which f = fc, derive directly
from changes in tenv. Figure 3B shows that this time scale—which
should be roughly equivalent to the damping time—decreases toward
the present.
The time scale of critical events appears to have declined since the
early Triassic (~250 Ma) by about one order of magnitude, to about 20
to 30 ky since the Miocene (~16 Ma). The smaller time scale is similar
to the ~10-ky time scale of the modern oceans’ homeostatic response
to a change in pH (22–25). The time scale’s decrease since the Triassic
is consistent with the evolutionary expansion of the modern biological
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
3 of 12
S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
M
O
L
I
E
M
O
year
L
M
D
A
B
Fig. 2. Depictions of the size and time scale of carbon isotopic events. (A) Mag-
nitude of the isotopic shift Dr = d1(tenv) − d1(0) as a function of its duration of time,
tenv. (B) Dimensionless mass perturbation M = |Dm|/m* as a function of the di-
mensionless time scale fctenv/2t0 for each of the events depicted in (A). The
straight (identity) line denotes the equality predicted by Eq. 5 when f = fc =
0.23 ± 0.07. Event abbreviations are defined in Table 1. Error bars (Materials
and Methods) are guides for interpretation.
pump beginning approximately 220 Ma (26, 27), corresponding to the
earliest observations of coccolithophores and associated calcareous phy-
toplankton in the fossil record (28, 29). As the biological pump
strengthened, the export of precipitated calcium carbonate from surface
waters to deep ocean sediments increased, leading to growth of the deep-
sea sedimentary carbonate reservoir (26, 27, 30). When pH decreases in
modern oceans, dissolution of sedimentary carbonate minerals is one of
the principal processes via which these disturbances are damped (31, 32).
Because this mechanism strengthened over the last 220 My, critical per-
turbations of the marine carbon cycle would have become shorter in time
and smaller in magnitude, but not necessarily less dangerous, because
their maximum flux would have remained fc. These observations are
each consistent with the suggestion that the carbon cycle became less var-
A
B
100
10–1
10–2
106
105
Mass extinction
Other
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
104
–500
–400
–200
–300
Date [My]
–100
0
Fig. 3. Evolution of size and time scale. (A) Dimensionless mass perturbation M =
|Dm|/m* as a function of the geologic age of the event. The relative absence of smaller
events in the deeper past likely derives from poor data and a lack of interest. The
absence of large events in the more recent past, however, cannot be explained by this
bias. (B) Time scale tenv of critical events (events falling within the error bars of the
straight line in Fig. 2B) as a function of geologic age. Error bars (Materials and Methods)
are guides for interpretation.
iable after the introduction of the deep-sea carbonate sink (26, 27, 30, 33).
They also identify the crossover time scale tx with the homeostat time
scale, implying tx ~ 10 ky (22–25) in modern oceans.
With these results in hand, the modern critical size Mc of an asymp-
totically fast perturbation of the marine carbon cycle now follows. First,
insertion of M(fc, tenv) from Eq. 5 into Eq. 2 provides the critical rate
rc = fc/2t0. Equation 3 then predicts
tx
Mc ¼ f
c
2t0
¼ 0:0082 ± 0:0041
ð6Þ
ð7Þ
where, as a guide, a total uncertainty of ±50% is assumed. The mod-
ern oceans contain about 38,000 Pg C; the critical dimensional mass
is therefore about 310 ± 155 Pg C.
From 1850 to the present, human activities have resulted in the
addition of about 155 ± 25 Pg C to the oceans (34). Projections for
further carbon uptake depend strongly on the trajectory of fossil fuel
emissions and land use, among other factors (34, 35). Figure 4 com-
pares the critical mass to the present accumulation and four projec-
tions to 2100 obtained from coupled climate–carbon cycle models
(34). The strictest emission scenario results in oceanic carbon uptake
whose mean falls just below the critical mass; at the other extreme,
the mean model uptake is about 74% greater than critical. Although the
uncertainty of each prediction in Fig. 4 is considerable, all scenarios
for cumulative uptake at the century’s end either exceed or are com-
mensurate with the threshold for catastrophic change.
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
4 of 12
S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
RCP8.5
RCP6.0
RCP4.5
]
C
g
P
[
e
k
a
t
p
u
e
v
i
t
l
a
u
m
u
C
600
500
400
300
200
100
0
RCP2.6
Present
Fig. 4. Cumulative modern ocean uptake of carbon since 1850, up to the
present (green) and projected to 2100 (blue), compared to the predicted crit-
ical mass of 310 Pg C (solid red line) and an assumed uncertainty of ±50%
(dashed red lines). Projections (34) are given for four representative concentra-
tion pathway scenarios RCPx, where x represents the radiative forcing, in units of
W/m2, deriving from accumulated emissions in the year 2100 (35). At the extreme
ends of the projections, RCP2.6 represents the range of lowest-emission scenarios
in the scientific literature, and RCP8.5 represents the high range of nonclimate
policy scenarios. Of the two intermediate pathways, RCP4.5 corresponds to an
emission pathway resulting from many climate policies found in the literature,
whereas RCP6.0 is representative of most scenarios without limitations on emis-
sions (35). The present cumulative uptake is obtained by adding 6 years of an
annual uptake rate of 2.3 Pg C year−1 (34) to the 2011 total of Ciais et al. (34).
DISCUSSION
The value of the empirical results in Figs. 2 and 3 lies not in any indi-
vidual data points but rather in their collective trends. The following
picture emerges. Carbon isotopic events rarely exceed a maximum iso-
topic shift that grows roughly like the logarithm of their time scale. This
upper bound appears related to the minimum rate—zero—at which or-
ganic carbon can be immobilized as rock. Events outside this limit result
from a fundamental disturbance of the carbon cycle, possibly related to
unstable dynamics, mass extinction, or both.
These conclusions follow from analyzing all isotopic events the
same way. Exceptions are, however, expected. For example, four events
(Ediacaran-Cambrian, Nemakit-Daldynian-Tommotian, Paleocene-
Eocene Thermal Maximum, and Miocene Climatic Optimum 2) un-
accompanied by mass extinction exceed the upper error bar of the
critical rate. If, say, these events were driven by dissociation of methane
hydrates [for example, (36)] rather than respired organic carbon, the
isotopic composition of the perturbative flux would be much lighter
(about −60‰). The estimated size of these events would then drop
by more than 50%, and they would each lie below the critical line. This
observation may help explain the Paleocene-Eocene Thermal Maximum’s
modest biotic impact. The Frasnian-Famennian extinction provides
another exception. Supposing that it is indeed a mass extinction, its
presence well below the critical line illustrates an important point: Mass
extinctions need not be caused by disruptions of the carbon cycle (2).
Modern investigations of mass extinctions often emphasize a plural-
ity of causes. Erwin’s “complex web of causality” (8, 37) addresses
how a combination of volcanism, climate change, marine anoxia, meth-
ane release, and other environmental stressors may have contributed
to the end-Permian extinction. Recent studies of the end-Cretaceous
extinction consider massive volcanism (38) in addition to a bolide impact
(39). Flood basalt eruptions are also clearly associated with the end-
Triassic (40) and end-Permian (15) extinctions, but their contribution
to CO2 levels is ostensibly modest (41). Evidently, the carbon cycle both
indicates and excites Earth system change. These dual roles merge, how-
ever, if external perturbations cause the cycle to respond by magnifying
the initial disturbance. System-wide instability may then follow. Because
the critical rate rc bounds qualitatively different dynamical regimes, per-
turbations that exceed rc (at time scales much greater than tx) suggest
such unstable evolution. The carbon cycle thus becomes one of many
environmental stressors, and an array of causes is naturally implicated.
Multiple causes may also be responsible for a “sixth extinction”
(1, 42). However, the anthropogenic disturbance of the carbon cycle
merits its own appraisal. If an anomalous flux greater than Mc/tx brings
the carbon cycle past a stability threshold at time scales much greater
than tx ~ 104 years, then the instantaneous addition of an anomalous
mass greater than Mc will excite similar behavior. The existence of these
thresholds is this paper’s central hypothesis. If they exist, one would
expect that an instability resulting from an impulsive addition of CO2
would play out over a time scale with the same order of magnitude as
the damping time scale tx. The upshot is that an unstable trajectory
would reach its maximum extent roughly 104 years after the threshold is
breached. But how that process plays out remains unknown.
MATERIALS AND METHODS
The event database
After a brief discussion of possible sources of error, this section sum-
marizes how each event has been identified and quantified. Geological,
paleoenvironmental, and paleobiological context is also given, as are
relevant references. All data, including the identification of abbrevia-
tions, are provided in Table 1.
When multiple observations are used to compute a mean, the sam-
ple SD, the unaveraged observations, or both are noted in the event’s
description. However, these quantities only indicate the variability of
a particular set of observations rather than an estimate of the obser-
vational error. As with any compilation of isotopic data [for example,
Hayes et al. (43)], inadequate sampling may be a more significant
problem. Here, it takes two forms: the sheer paucity of observations,
and biases deriving from not only the original investigators—who
may favor larger isotopic excursions—but also the needs of the present
compilation, which favors carbon isotopic studies associated with strong
geochronological constraints. An additional source of error derives
from incompleteness of the sedimentary record. For example, a hiatus
in sediment deposition, or erosion after deposition, may make the true
carbon isotopic minimum effectively unavailable to the modern observer,
thereby leading to an underestimate of |Dr|.
The ideal geochronological control would provide absolute dates
at the initiation and termination of isotopic excursions. Only one
event in the database—at the Ediacaran-Cambrian boundary—meets
this standard. Among the others, all but 1 of the 18 events extending
from the Toarcian (182 Ma) to the Miocene (16 Ma) are timed by as-
trochronology, including, in some cases, information obtained from
isotope geochronology. The duration of five of the older, pre-Jurassic
events is determined from at least two dates provided by U-Pb geo-
chronology that bracket less than twice the accumulation of sediment
deriving from the event itself; a sixth event (end-Permian) comes close
to this standard.
The remaining pre-Quaternary events are timed by biostratigraphy
with reference to the Geologic Time Scale (44); their time scales are there-
fore the most poorly constrained. Two of them, the end-Ordovician and
Frasnian-Famennian events, are among the Big Five mass extinctions.
A third mass extinction—the end-Triassic—is constrained by both
isotope geochronology and astrochronology, but these constraints
are not directly associated with the carbonate–carbon isotope excursion.
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
5 of 12
S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
If these three mass extinction events and the other four events with poor
temporal constraints (Cambrian Spice, Tournaisian, mid-Capitanian,
and Albian-Cenomanian) are removed from the database, the trends
in Figs. 2 and 3 remain intact and no new trends appear.
Nevertheless, the unknown errors in the observations render quan-
titative error analysis infeasible. The representative error bars for Dr, M,
and tenv in Figs. 2 and 3 are instead guides for interpretation.
The following list proceeds from the oldest to the youngest event.
Ediacaran-Cambrian boundary, ~541 Ma
The magnitude and duration of this negative isotopic excursion derive
from the Oman carbon isotope data of Amthor et al. (45) and the U-Pb
geochronology of Bowring et al . (46). Other data from Morocco, Siberia,
Mongolia, and China (47, 48) suggest global consistency. The absence of
the Ediacaran biota in the Cambrian has led to the suggestion that the
Ediacaran biota vanished by mass extinction at the end of the Ediacaran,
possibly related to the environmental event signaled by this excursion
(45, 49–52). However, gradual extinction coupled with permanent
environmental changes unfavorable to preservation of the Ediacaran
biota may also be possible (53).
Nemakit-Daldynian-Tommotian boundary, ~525 Ma
Both the carbon isotopic excursion and the U-Pb dates establishing
the duration of this Early Cambrian event are from the Moroccan data
of Maloof et al . (47, 54); data from Siberia, Mongolia, and China con-
firm its global nature (47, 48). No mass extinction has been associated
with this interval. However, diversity (number of genera) and disparity
(number of classes) sharply increase in the Tommotian, primarily be-
cause of the rise of small shelly fossils (55). The Nemakit-Daldynian-
Tommotian boundary is also associated with a transition from seawater
chemistry that favored aragonite precipitation to seawater favoring cal-
cite precipitation (47, 56).
Cambrian Spice, ~497 Ma
The duration and magnitude of this late Cambrian event, also known
as the Steptoean positive carbon isotope excursion, are taken from the
Australian data of Saltzman et al . (57). The event has also been found in
North America, China, Kazakhstan, and Sweden (58, 59). The initiation
of this positive isotopic excursion coincides with an extinction of trilo-
bites (58).
End-Ordovician event, ~446 Ma
The latest stage of the Ordovician—the Hirnantian—simultaneously
exhibits a positive isotopic excursion and one of the Big Five mass
extinctions. I estimated the magnitude of the negative limb of the
excursion from data obtained at Vinni Creek and Monitor Range,
Nevada (60), Copenhagen Canyon, Nevada (61), Blackstone Range,
Yukon (60), Wangjiawan, south China (62), the Kardla drill core, Estonia
(63), and Anticosti Island, Quebec (64). The seven data sets yield a mean
excursion of 5.6 ± 0.9‰. Each data set except for Copenhagen Canyon
has been correlated to graptolite zones by Gorjan et al . (62). By estimating
the fraction of the persculptus and extraordinarius zones in which the
excursion occurs at each site and estimating those time intervals from
Gradstein et al . (44), I found that the negative limb of the excursion lasted
approximately 230 ky.
Silurian Mulde event, ~428 Ma
A sequence of two positive excursions, known as the “Mulde event,”
occurs in the Homerian stage of the Silurian, during a period of grap-
tolite extinctions known as the “Big Crisis” (65). The event is apparently
global (65). Cramer et al. (66) have recently bracketed the earlier of the
two excursions between two U-Pb dates obtained in Gotland, Sweden.
I focused on the negative downswing, which occupies somewhat more
than half of the dated interval, and estimated its duration to be about
260 ky. As shown by Cramer et al . (66), the Gotland excursion correlates
well with observations in the West Midlands, England, with the West
Midlands excursion having a magnitude of about 4‰ and the Gotland
excursion about 2‰, a range that is consistent with other observations
(65, 66). I therefore estimated the magnitude to be 3 ± 1‰.
Silurian Lau event, ~423 Ma
A significant global positive carbon isotopic excursion known as the
“Lau Event” occurs in the Ludlow Epoch of the Silurian (65, 67–72).
The excursion typically peaks at 6 to 8‰ [for example, Lehnert et al.
(71)], though values as small as 3 to 4‰ and as high as 10 to 11‰
have been observed in North America (68) and Sweden (70), respective-
ly. Given that uncertainty, I adopted a conservative mean peak value of
7‰. The negative downswing after the peak typically bottoms out at
about 1‰, resulting in a typical total negative shift of about 6‰.
Cramer et al. (72) have recently bracketed the duration of the excursion
by U-Pb geochronology, finding that it must be less than about 1.17 My.
Following their analysis, I took the upswing and downswing to be
roughly symmetric, each occupying about one-half of an approximately
million-year interval, implying that the duration of the downswing is
about 500 ky.
Frasnian-Famennian boundary, ~372 Ma
The boundary between the Frasnian and Famennian stages in the Late
Devonian is associated with one of the Big Five mass extinctions (5, 6).
But the drop in diversity at the boundary is thought to be a consequence
of a lack of originations rather than an elevated extinction rate (19, 20).
The boundary is associated with a global positive carbon isotopic excur-
sion, known as the Upper Kellwasser Event (73). The magnitude and
duration of the excursion are taken from the compilation of European
Devonian data of Buggisch and Joachimski (74).
Tournaisian event, ~352 Ma
Buggisch et al. (75) presented carbon isotopic data from Belgium,
Ireland, France, western Canada, and the western United States that
collectively exhibit a pronounced positive isotopic excursion in the
Tournaisian stage of the Carboniferous. Using their biostratigraphic
time scale, I estimated the duration of the negative limb of the excursion
to be 1.51 My; the magnitude of the excursion in their compilation is
about 3.5‰.
Mid-Capitanian event, ~262 Ma
A mass extinction in the Capitanian Stage of the Middle Permian is as-
sociated with a negative isotopic excursion of carbonate carbon of about
5‰ over about 500 ky (76, 77). Observations of the event in south China
suggest that the extinction event coincides with the initiation of Emeishan
volcanism (78) and the onset of the excursion (76).
End-Permian extinction, ~252 Ma
The end-Permian extinction is considered the most severe of the Big
Five (8). A significant negative isotopic excursion occurs at the time
of the extinction, just below the Permian-Triassic boundary. Korte
and Kozur (16) reviewed observations of the excursion in 40 localities
worldwide and concluded that the excursion ranged from 4 to 7‰; I
therefore took the magnitude of the excursion to be at the center of that
distribution (5.5‰). The 60-ky duration of the excursion derives from
the well-studied section at Meishan, south China, where Cao et al . (79)
have provided high-resolution carbon isotopic data and U-Pb geo-
chronology provides strong constraints on the time scale (9).
Early Triassic, ~251 Ma
Several significant isotopic excursions of unknown origin follow the
end-Permian extinction (80–83). Galfetti et al. (83) provided carbon
isotopic data tied to U-Pb geochronology from the Loulou Formation,
northwest Guanxi, south China. The only temporally constrained
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
6 of 12
S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
negative isotopic shift occurs in a negative excursion that begins ap-
proximately 251.22 Ma (83) in the mid-Smithian and declines approx-
imately 3.5‰ over approximately 250 ky.
Triassic-Jurassic boundary, ~201 Ma
One of the Big Five extinctions, the end-Triassic event is temporally
associated with the emplacement of the Central Atlantic magmatic
province (40). Geochemically, the most well-studied section is at St.
Audries Bay in southwest England, where the extinction interval coin-
cides with a rapid fall and subsequent rise of the isotopic composition of
organic matter (84) that lasts 20 to 40 ky according to astrochronology
(85). This so-called initial excursion in organic carbon is widely observed
(86). However, its significance for understanding the evolution of marine
dissolved inorganic carbon is unclear (86, 87), partly because it may be
associated with observed changes in flora (87)—and therefore variations
in the organic matter itself—and partly because inorganic isotopic com-
positions need not track organic isotopic compositions (43). There are
unfortunately few well-resolved isotopic studies of carbonate carbon as-
sociated with the end-Triassic event (88–93). These carbonate studies
also reveal a negative excursion associated with the extinction, but it
is often unclear if they represent the “initial excursion” or the later
120-ky-long “main excursion” seen in organic carbon (84, 85). A recent
review (94) cautions that diagenetic alteration may have corrupted car-
bonate data, yet a careful study of lithological effects at an Italian section
suggests that the carbonate excursion may indeed be primary (87). An
important additional problem is the need to estimate the time scale of
the carbonate excursion. The data of Clémence et al. (92) provide a so-
lution to these problems. In their analysis of inorganic and organic car-
bonate at the Tiefengraben section in the Austrian Alps, they found an
initial 2‰ negative excursion in carbonate that extended to the second
significant minimum in their organic data (92 ), which in turn correlates
reasonably well with the astrochronologically calibrated organic data at
St. Audries Bay (85 ). Their initial carbonate excursion was then
found to last approximately 50 ky [which included one precessional
cycle of the “pre-recovery” interval identified by Ruhl et al . (85)]. This
2.0‰ shift is also consistent with isotopic analyses of carbonate in pris-
tine oysters located near St. Audries Bay (91).
Toarcian oceanic anoxic event, ~183 Ma
The early Toarcian oceanic anoxic event (95) is widely observed in
Europe as a negative excursion in the isotopic composition of carbonate
(96–99). To estimate its magnitude and duration, I used the high-resolution
carbonate record of Hesselbo et al . (96 ) obtained at Peniche (Portugal).
The astrochronology of Suan et al . (97) found a duration of 150 ky (their
interval C1) for the negative isotopic excursion, which is consistent with
the subsequent analysis of Huang and Hesselbo (100). Using this chro-
nology, recent U-Pb dating by Burgess et al. (101) found that the intru-
sive magmatism associated with the Ferrar large igneous province is
synchronous with the negative excursion.
Aptian oceanic anoxic event, ~120 Ma
The negative isotopic excursion associated with the Aptian oceanic
anoxic event has been observed worldwide in open ocean records of
carbonate carbon (102–105). My estimate of the magnitude and duration
of the excursion derived from the astrochronology of Malinverno et al.
(106) performed on the carbonate record of Erba et al. (103). The negative
excursion spans the interval C3 [originally identified by Menegatti et al.
(102)] and lasts approximately 47 ky.
Albian-Cenomanian boundary, ~100 Ma
This complex carbon isotope event occurred near the Albian-Cenomanian
boundary (mid-Cretaceous) in several European sections (107). Here,
I focused on the carbon isotopic data obtained at Ocean Drilling
Program (ODP) Site 1050 at Blake Nose, western North Atlantic, as
presented by Ando et al . (108). When plotted using the time scale of
Huber et al . (109), the ODP data show an unambiguous 0.7‰ negative
excursion across the boundary over approximately 110 ky.
Mid-Cenomanian event, ~95.9 Ma
The mid-Cenomanian positive carbon isotopic excursion has been
clearly observed in European sections (107), the North Atlantic (110),
western North America (111), and elsewhere. My analysis follows from
the correlation of the “English Chalk” reference carbonate curve (107)
to carbon isotopic data from the Cretaceous Western Interior Seaway of
North America (111, 112). I used the time scale of Eldrett et al. (112),
which is based on a combination of astrochronology and U-Pb geo-
chronology. As indicated in figure 11 of Eldrett et al. (112), the negative
swing of the excursion begins at about 96.5 Ma, followed by a decline of
about 0.85‰ over roughly 138 ky.
Cenomanian-Turonian boundary event, ~94.2 Ma
Otherwise known as oceanic anoxic event 2, the positive isotopic
excursion at the Cenomanian-Turonian boundary is probably the
most widely observed of the Cretaceous isotopic events (107, 111, 112).
As for the mid-Cenomanian event, my estimate of the time scale and
magnitude of the Cenomanian-Turonian event follows from the cor-
relation of the English Chalk reference carbonate curve (107) to car-
bon isotopic data from the Cretaceous Western Interior Seaway of
North America (111, 112). I took the initiation of the event to be
at the positive peak labeled B in figure 11 of Eldrett et al. (112)
and found that the ensuing negative excursion declines by 1.60‰
in roughly 553 ky.
Cretaceous-Paleogene extinction, ~65.5 Ma
The end-Cretaceous negative isotopic excursion is associated with
one of the Big Five extinctions (5), widely known for the extinction
of dinosaurs and the Alvarez impact hypothesis (39). Its temporal
association with the eruption of the massive Deccan volcanic province
in India is less widely known (38). I obtained the magnitude and time
scale of the carbon isotopic event from the deep-sea bulk carbon iso-
topic data obtained at ODP Sites 1210 (Northwest Pacific) and 1262
(South Atlantic), as presented by Alegret et al . (113) using an orbitally
tuned time scale (114). I took the initiation and termination of the
approximately 1.15 ± 0.03‰ (the mean of 1.16 and 1.13‰) negative
isotopic excursion at each site to be where values of d1 begin to change
by at least 0.1‰, resulting in a mean event duration of roughly 26 ky
(the average of 8.5 and 44 ky).
Early late Paleocene event, ~58.9 Ma
Also known as the mid-Paleocene biotic event, this negative excursion is
synchronous with dissolution of carbonate and changes in the orga-
nization of benthic and planktonic ecosystems (115, 116). The mag-
nitude and time scale of this negative excursion are taken from the
high-resolution bulk isotopic record of Littler et al. (115).
Paleocene-Eocene Thermal Maximum, ~55.5 Ma
The negative isotopic excursion of the Paleocene-Eocene Thermal Max-
imum is perhaps the most studied carbon isotopic event in Earth history
(11), in large part because of its association with significant climate
warming and clear evidence of ocean acidification (10). The event is also
associated with a significant extinction of benthic foraminifera (12);
however, other groups of benthic and planktic microfossils show little
or no extinction (11). The time scales of the isotopic excursion of bulk
carbonate in two deep-sea cores, from ODP Site 1266 on the Walvis
Ridge in the South Atlantic and ODP Site 690 in the Weddell Sea,
Southern Ocean, have each been estimated independently by identifica-
tion of orbital cycles (117) and the estimation of sedimentation rates
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
7 of 12
S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
from the concentration of 3He (118, 119). The mean of the resulting
four estimates [corresponding to the cumulative time between the onset
of the excursion and the termination of its “core” as summarized in
table 1 of Murphy et al . (119)] is 83 ± 23 ky. I obtained the magnitude
of the isotopic excursion, 2.7 ± 1.1‰, from the mean of 33 published
analyses of bulk Paleocene-Eocene Thermal Maximum carbonates re-
viewed by McInerney and Wing (11). The initial value of the excursion
was obtained similarly (11). Although these averages lack prejudice,
they may nevertheless underestimate the excursion’s size (120) and
overestimate its time scale (121).
Eocene Thermal Maximum 2, ~53.7 Ma
The Eocene is punctuated by several “hyperthermal events,” each repre-
sented by a negative excursion in the isotopic composition of carbonate
carbon and dissolution of deep-sea marine carbonates. Eocene Thermal
Maximum 2, one such event, follows the Paleocene-Eocene Thermal
Maximum (also known as Eocene Thermal Maximum 1) by roughly
2 My. My estimate of the magnitude of Eocene Thermal Maximum 2
derived from averaging estimates from four high-resolution bulk isotopic
records presented by Stap et al. (122). The records derive from ODP Sites
1262, 1263, 1265, and 1267, corresponding to paleowater depths ranging
from about 1500 to 3600 m (122); the shallowest site yields an excursion
of about 1.5‰, whereas the others are each about 1.0‰ (122), yielding
a mean of 1.13 ± 0.25‰. The time scale is about 45 ky.
Eocene Hyperthermal H2, ~53.6 Ma
The magnitude and time scale of this Eocene hyperthermal event come
from ODP Sites 1263, 1265, and 1267. Following Stap et al. (122), I took
the excursion of bulk carbon isotopes to be about 0.6‰ and the time
scale to be about 33 ky.
Eocene Hyperthermal I1, ~53.2 Ma
I obtained the magnitude (0.72‰) and time scale (41 ky) of this Eocene
hyperthermal event from the high-resolution bulk carbon isotopic
record of Littler et al. (115), which derives from ODP Site 1262.
Eocene Hyperthermal I2, ~53.1 Ma
This Eocene hyperthermal follows Eocene Hyperthermal I1 by about
100 ky. Using the same source (115) as for Eocene Hyperthermal I1,
I found an excursion of about 0.61‰ and a time scale of 40 ky.
Eocene Thermal Maximum 3, ~52.5 Ma
I estimated the magnitude and time scale of this hyperthermal event from
the high-resolution benthic carbon isotopic record of Lauretano et al.
(123), obtained at ODP Sites 1262 and 1263. The two records are similar,
yielding an excursion of about 0.8‰ and a time scale of about 37 ky.
Early Oligocene Event, ~33.5 Ma
The positive carbon isotopic excursion just above the Eocene-Oligocene
boundary is associated with the initiation of permanent Cenozoic
ice sheets on Antarctica (124–126). I focused on the negative down-
swing to lighter isotopic compositions that followed. The time scale
and isotopic change are derived from the high-resolution data col-
lected at ODP Site 1218, using the astrochronological time scale
presented by Coxall et al. (125). The negative shift extends over rough-
ly 0.6 ± 0.1‰ during a period of about 430 ± 100 ky. The three ODP
records presented by Zachos and Kump (126) are consistent with
these estimates.
Miocene Climatic Optimum 1, ~16.9 Ma
Holbourn et al . (127) provided high-resolution carbon isotopic records
spanning most of the Miocene Climatic Optimum, a period of warm
climates ranging from 17.0 to 14.7 Ma that interrupted the longer-term
trend of Cenozoic cooling. The records are obtained from Integrated
ODP Site U1337 in the eastern equatorial Pacific Ocean, at a paleowater
depth of 3500 to 4000 m. The Miocene Climatic Optimum began with a
negative isotopic excursion, which I designated Miocene Climatic Op-
timum 1, that has also been identified in lower-resolution records from
three other sites, elsewhere in the Pacific and in the Southern Ocean
(127). I obtained the magnitude (0.5‰) and time scale (28 ky) from
the bulk carbonate record. The event is synchronous with shoaling of
the carbonate compensation depth (CCD).
Miocene Climatic Optimum 2, ~16.4 Ma
Holbourn et al. (127) identified three other Miocene Climatic Opti-
mum negative excursions, which I designated Miocene Climatic Op-
timum 2 to Miocene Climatic Optimum 4, during which d1 decreases
sharply and the CCD appears to shoal. Using the same high-resolution
record as for Miocene Climatic Optimum 1, I found that Miocene
Climatic Optimum 2 has a time scale of 22 ky and a magnitude of
0.76‰.
Miocene Climatic Optimum 3, ~16.0 Ma
See Miocene Climatic Optimum 2. For Miocene Climatic Optimum 3, I
found a time scale of 32 ky and a magnitude of 0.66‰. This event also
appears in the records of ODP Site 1146 in the western Pacific Ocean
(128) and ODP Site U1338 in the eastern equatorial Pacific (129).
Miocene Climatic Optimum 4, ~15.6 Ma
See Miocene Climatic Optimum 2. For Miocene Climatic Optimum 4,
the excursion has a magnitude of 0.6‰ and a time scale of 29 ky.
Last Glacial Maximum–to–Holocene deglacial, ~0.021 Ma
Using a compilation of 480 records of benthic foraminiferal carbon iso-
tope analyses, Peterson et al . (130) estimated a whole-ocean increase in
the isotopic composition of dissolved inorganic carbon of 0.34 ± 0.19‰
from the Last Glacial Maximum to the Holocene. I took the time scale to
be 18 ky, given their Last Glacial Maximum averages over 19 to 23 ky
and Holocene averages over 0 to 6 ky. This is the only positive isotopic
shift in the database in which the magnitude and time scale derive from
the upswing; for the others, such as the Ordovician, the calculations cor-
respond to the ensuing downshift.
The boundary value problem
The boundary value problem derives from a model of the global carbon
isotopic system presented by Rothman et al. (131, 132). This model
partitions carbon into two reservoirs. The inorganic reservoir contains
a mass m of carbon with isotopic composition d1; the organic reservoir
contains carbon of isotopic composition d2. In the Phanerozoic, the
isotopic composition of the organic reservoir typically equilibrates
with respect to inputs and outputs at time scales much faster than
the inorganic reservoir (131). The evolution of d1 is then described
by equation 8 of Rothman et al . (132):
m dd1
dt
¼ jiðdi d1Þ þ b2e
ð8Þ
where ji is the flux of carbon of isotopic compositiond i into the surface
carbon cycle (deriving typically from volcanism and continental
weathering), b2 is the burial flux of organic carbon out of the surface
cycle, and e = d1 − d2 is the isotopic fractionation between inorganic
and organic carbon.
Suppose now that a new source of carbon with isotopic composition
d′
i flows into the carbon cycle with flux j′
i , and that the burial flux b2
increases by an amount b′
2 while e remains constant. Equation 8 then
becomes
m dd1
dt
¼ j*
i
ðd*
i
d1Þ þ j′
i
ðd′
i
d1Þ þ ðb*
2
þ b′
2
Þe
ð9Þ
l
D
o
w
n
o
a
d
e
d
f
r
o
m
h
t
t
p
:
/
/
a
d
v
a
n
c
e
s
.
s
c
e
n
c
e
m
a
g
i
.
o
r
g
/
o
n
J
a
n
u
a
r
y
3
,
2
0
1
8
Rothman, Sci. Adv. 2017; 3 : e1700906
20 September 2017
8 of 12