Self-sustained oscillations and global climate changes

The periodic changes of atmospheric CO2 and temperature over the last 5 Myr reveal three features that challenge current climate research, namely: (i) the mid-Pleistocene transition of dominant 41-kyr cycles to dominant 100-kyr cycles, (ii) the absence of a strong precession signal of approximately 20 kyr, and (iii) the cooling through the middle and late Holocene. These features are not directly addressable by Earth’s orbital changes described by Milankovitch. Here we show that a closed photochemical system exposed to a constant illumination source can sustain oscillations. In this simple conceptual model, the oscillations are intrinsic to the system and occur even in the absence of periodic radiative forcing. With proper adaptations to the Earth system, this oscillator explains the main features of past climate dynamics. Our model places photosynthesis and the carbon cycle as key drivers of climate change. We use this model to predict the relaxation of a 1,000 PgC pulse of CO2. The removal of 50% of this CO2 will require one century, and will lead to a warmer and wetter future. However, more pronounced glaciation cycles emerge on the millennial timescale.

did climate cycles evolve from 40-to 100-kyr periodicities ca. 1 Myr ago? What is the origin of the dominant 100-kyr frequency of the last 5 glacial cycles? Abandoning the hypothesis of orbital forcing raises an even more fundamental question: how can a closed system such as Earth trigger and sustain oscillations on the millennia time scale? The answer to this question should point to a physical mechanism capable of sustained oscillations in a closed system such as Earth.
According to the second law of thermodynamics, in an isolated system, exchanging neither energy nor matter with the outside world, entropy increases monotonically to its maximum at equilibrium. An open system can exhibit persistent oscillations of some of its components because there is a constant influx of another component. The first efforts to find open systems with undamped oscillations appeared in chemical sciences and were authored by Lotka 26 , but it was the explanation of Volterra for predator-prey population relations 27 that popularized undamped oscillations. Prigogine demonstrated the possibility of having chemical oscillations of intermediates of an overall reaction A + B → C + D 28 . The mechanism proposed, named Brusselator 29 , is the prototypical model of an autocatalytic chemical reaction providing oscillations and even chaos when two Brusselators are coupled 30 . The first experimental observation of oscillations was the Belousov-Zhabotinsky (BZ) reaction 31 , that undergoes approximately 100 oscillations before achieving the final equilibrium state. Field and Noyes proposed a general kinetic scheme, named Oregonator, to explain such oscillations and showed that sustained oscillations required some chemical components to remain constant 32 , i.e., sustained oscillations required influx of material. This led to the believe that a closed system at constant temperature and pressure should attain equilibrium 33 . Sustained oscillations were only anticipated for open systems described by a set of non-linear chemical reactions with autocatalysis and feedback 34 .
Earth is a (nearly) closed system receiving an influx of energy from the Sun. The oscillators mentioned above cannot represent an internal Earth system mechanism with glaciation cycles. However, a closed chemical reaction system where a net product of other steps is photochemically transformed back in the reactant may have sustained oscillations 35 . Photosynthesis can maintain the state of disequilibrium on Earth, but storage and release of free energy alone are insufficient to drive oscillations. Persistent oscillations also require autocatalysis, where one component X is a catalyst of its own production, D + X → Z + 2X.
We present a simple conceptual model for intrinsic oscillations in a closed system under constant illumination, such as the Earth, that accounts for the major trends of temperature, CO 2 and biomass cycles over the last 5 million years. We design the simplest oscillator for a closed system that can sustain periodic changes of its composition under constant illumination, and then add the minimum complexity required to represent the Earth system. The model only employs two adjustable parameters, which are fitted to historic CO 2 cycles, but incorporates a mechanism that triggers oscillations. This trigger is associated with a change in the global efficiency of photosynthesis in the Pliocene. The oscillations are enabled by an autocatalytic step. The dominant periodicity of 100 kyr arises naturally from the timescales of CO 2 land uptake, ocean invasion, reaction with calcium carbonate Changes in solar energy received at 25°N in the month of March relative to average (ca. 1 GJ/cm 2 , centered in zero) 18,19 . Corresponding frequency distributions calculated for the time interval from -5 Myr to present (C,E), or for the last 500 kyr (D,F).

Results
Oscillations in a closed system under illumination. Dynamical systems are often described by a set of ordinary differential equations (ODEs) that contain parameters in addition to variables. In some systems, similar sets of parameter values lead to dynamic behaviors of qualitatively different nature. For example, one set of parameters may lead to a stable equilibrium point-an attractor-whereas the other set leads to an asymptotically stable periodic solution-a limit cycle. Such models are said to contain bifurcation points. The case given above, where an attractor becomes a limit cycle as a parameter is varied, is named a supercritical Hopf bifurcation. When such a bifurcation is present, the system may exhibit exponentially damped oscillations leading to equilibrium or to limit cycle oscillations 36 . These limit cycles are periodic orbits when represented as trajectories in the phase space of the variables. Chemical oscillations arise from Hopf bifurcations in the Oregonator 32 , in the oscillations of glycolysis 37 and in the Brusselator 29 .
Following the tradition of naming oscillators based on their geographical origin, we name Coimbrator the set of ODEs represented in Fig. 2. They enable "self-oscillations" of composition variables as a parameter crosses a critical value. Figure 2 emphasizes that the Coimbrator is a closed system with two photochemical steps. The equations are at most bimolecular, contain three composition variables (X, Y, Z) and include an autocatalytic reaction, which are the minimal sufficient conditions for realistic oscillatory behavior 38 .
This closed system can be simplified when the concentrations of dyes A and B are high compared with the photon flux, and can be considered constant. Moreover, if the photon fluxes with energies hν 1 and hν 2 are also approximately constant, the concentrations of C and D will remain approximately constant (steady-state). Hence, for large A and B and constant illumination, A, B, C and D can be treated as constants and the system can be decoupled from the variations of X, Y and Z. "Methods" section shows how to express these species in dimensionless variables (x, y, z) to obtain the following ODEs where the dot notation denotes the derivative with respect to dimensionless time; α = (k 3 D)/(k 1 C) and β = k 4 / (k 1 C) are parameters. We define 1/(k 1 C) as the characteristic time.
The equilibrium (fixed) points of this three-dimensional system are obtained solving ẋ = 0, ẏ = 0 and ż = 0 simultaneously. They are O = (0,0,0) and P = (1, α, α/β). When β < [α-2 + √(α + 4)]/2 the point P becomes a repelling equilibrium point and a periodic orbit appears. "Methods" section presents the treatment of the bifurcation. The Hopf bifurcation is supercritical in the whole region depicted in Fig. 3. Attracting limit cycles emerge when crossing the bifurcation curve with β decreasing or α increasing. Figure 3 illustrates the bifurcation selecting a segment in the parameter space β = 0.6180 and α ∈ [0.80, 1.25]. Figure 4 presents phase space trajectories to illustrate the crossing of the supercritical Hopf bifurcation. The numerical simulations in the upper panel illustrate a trajectory for certain values of the parameters (e.g., α = 1 and β = 0.7) that is attracted to an equilibrium point. This is true independently of the initial values of the variables (x, y, z). The trajectories in the lower panel calculated with another set of parameters (e.g., α = 1 and β = 0.6) go to an  www.nature.com/scientificreports/ asymptotically stable periodic orbit. A Hopf bifurcation occurs for β = 0.618 when α = 1. Figure 4 also represents the values of the variables as a function of time. For α = 1 and β = 0.7 we see exponentially damped oscillations of the variables leading to a final stable value. This closed system sustains oscillations of the variables for α = 1 and β < 0.618. The Coimbrator leads to oscillations when β crosses a critical value.
The Earth system as an oscillator. The Coimbrator shows that the internal Earth system mechanism that drives climate changes can be an oscillator with parameters that in the Pleistocene crossed those of a supercritical Hopf bifurcation and triggered glacial cycles. We explore this hypothesis incorporating the most essential elements of the Earth system in the oscillator. Photosynthesis in plants involves two photons of different energies to produce a high-energy intermediate, This overall reaction of oxygenic photosynthesis represents the buildup of biomass. Hence, we identify CH 2 O with Z. It represents both live biomass and dead organic carbon. Its consumption is the oxidation reaction The amount of oxygen in the atmosphere decreased by less than 1% in the last 5 million years and remained at least one order of magnitude higher than that of H 2 O or CO 2 . Assuming that O 2 is approximately constant and can be assimilated in k 4 , the autocatalysis reaction becomes where D comprises elements from the photosynthetic systems and carbon fixation cycle. A, B and C are also part of this carbon fixation cycle. This leads to the assignment of X as water vapor and, consequently, Y = CO 2 . This approach emphasizes the analogy with the Coimbrator and helps to identify both the photochemical reaction that drives the system and the autocatalysis required for oscillations. These aspects and the presence of a quadratic term, make the representation of the oscillator as a compartmental model that exchanges substances between physical spaces (e.g., atmosphere, sea, biosphere) less insightful than the representation of Fig. 2.
The main control of water vapor in the atmosphere is the vapor-liquid (X ↔ X 0 ) equilibrium. We add the (X ↔ X 0 ) equilibrium to the oscillator to account for the saturation of water vapor pressure in moist air. Similarly, the CO 2 solubility (Y ↔ Y 0 ) equilibrium is included to the oscillator, together with the reactions of CO 2 dissolved in water, CO 2(aq) . "Methods" section show how choices of physically motivated parameters constrain changes in X and Y. Adding these equilibria to the oscillator changes its nature. We name the new oscillator Glaciator. Figure 5 illustrates how the Glaciator incorporates the basic features of the Earth system. "Methods" section shows how to formulate the Glaciator in terms of dimensionless variables. Assuming that X o and Y 0 are constant, the new set of ODEs involves four parameters (ε, σ, ω, ρ) www.nature.com/scientificreports/ The rates of the processes represented by k 3 D and k 4 depend on the biomass, which is not constant, and the corresponding parameters are better treated as time dependent, α(t) and β(t). We recall that in the Coimbrator we assumed a constant irradiation source, and consequently constant α and β, but this is no longer valid when we account for the insolation changes of Fig. 1 in the Glaciator. The rate constants k 5 and k 6 for the (X ↔ X 0 ) equilibrium depend on the temperature and, consequently, ε and σ should also be time dependent. For simplicity, that is not considered in this study. Figure 5 represents the dissolution of carbonates explicitly for CaCO 3 and MgCO 3 . Identically, CO 2 reactions with calcium-or magnesium-containing silicate minerals (CaSiO 3 or MgSiO 3 ) are equally present. The overall reactions timescales of climate changes. We show in "Methods" section that meaningful values of the parameters in Eq. (5) require ε ≈ 33ω and σ ≈ 3.3ρ. The independent parameters α(t), β(t), ω and σ are defined as ratios of rate constants (or time constants). This dependence on the rates places the timescales of the processes in Fig. 5 under careful scrutiny. The major processes for CO 2 removal from the atmosphere are those of land uptake (timescale 1-10 2 yr), ocean invasion (10-10 3 yr), reaction with calcium carbonate (10 3 -10 4 yr) and silicate weathering (10 4 -10 6 yr) 40 . In the Glaciator, CO 2 removal is represented by We associate the rate of land uptake with Considering that globally water vapor pressure is much larger than CO 2 pressure, p H2O ≫ p CO2 , it is reasonable to incorporate the water vapor pressure in the rate constant and relate 1/(k 2 p H2O ) to the timescale of CO 2 removal by land uptake, 130 yr 41 . This, together with present day p H2O ≈ 11 matm, discussed in "Methods" section, allow us to calculate k 2 ≈ 0.7 atm -1 yr -1 . We associate the timescale of silicate weathering with 1/k 7 and make 1/k 7 ≈ 100 kyr. This is also the residence time of inorganic carbon in the ocean 39 . The value of 1/(k 1 C), defined as the characteristic time, must be related to the timescales of ocean invasion and reaction with calcium carbonate. Admittedly, this step of the Glaciator could be divided in two steps, because at least these two physical processes are present. This would increase the complexity and the number of parameters in the calculations. We favored simplicity and use 1/(k 1 C) to describe the two processes and represent the mean atmospheric lifetime of CO 2 . The mean atmospheric lifetime of CO 2 after a large emission pulse is 12-14 kyr 42 . Hence, we make 1/(k 1 C) ≈ 12 kyr, which is in the upper limit of the geometric mean of ocean invasion and calcium carbonate reaction timescales.
These timescales determine the characteristic concentration (k 1 C)/k 2 = 0.12 matm, and allow for the conversion of the dimensionless variables (x, y, z) to absolute values of H 2 O, CO 2 and biomass. The value of the parameter ρ is given by its definition, ρ = k 7 /(k 1 C), and using the timescales discussed above we obtain ρ = 0.11. This restricts the number of adjustable parameters to three: α(t), β(t) and ω. The parameter ω relates the timescale of (7) −→ X 0 , X 0 Scientific RepoRtS | (2020) 10:11200 | https://doi.org/10.1038/s41598-020-68052-9 www.nature.com/scientificreports/ CO 2 release from the surface ocean to the atmosphere with its mean atmospheric lifetime. We employed ω = 5 but the calculations are not very sensitive to this value. The values of the timescales employed in this work are presented in Table 1. They lead to a peak in atmospheric CO 2 in the Holocene before the onset of anthropogenic perturbations 43 . "Methods" section present a study of the Hopf bifurcation in the Glaciator. For the fixed values ω = 5, ρ = 0.11, ε = 33ω and σ = 3.3ρ, the Hopf bifurcation is supercritical for the values of α and β shown in Fig. 6. This figure shows the crossing of the bifurcation curve when β = 0.46 and α increases in the segment α ∈ [0.8,0.9]. The Glaciator may trigger sustained oscillations when β decreases or α increases from a set of parameters corresponding to a stable equilibrium.
Photosynthetic forcing. Given the definition of α and the processes in Fig. 5, α can be regarded as the ratio between the rate of formation of the biomass and the rate of incorporation of CO 2 . Hence, α is a measure of the efficiency of photosynthesis and should not exceed unity.
There are two major types of photosynthesis, C3 and C4 photosynthesis, named after the number of carbon atoms of the first compound in which CO 2 is incorporated. C3 plants (e.g., trees, wheat, rice, soybean) were spread in ancestral atmospheres, characterized by elevated CO 2 and low O 2 levels, whereas C4 plants (e.g., corn, sugarcane, many grasses) became dominant at mid-to-low latitudes in the last 5 Myr 44 , and now represent ¼ of the primary productivity on the planet 45 . C4 plants are able to concentrate CO 2 and minimize photorespiration. This contributes to a specific activity (i.e., mol CO 2 fixed per mass of enzyme per unit time) of C4 plants that can be twice as large as that of C3 plants 46 , and yields 50% greater efficiencies in radiation use 47 . Given the limit of α, we make α C3 ≈ 0.65 for C3 plants and, because the efficiency of C4 plants is 50% higher than that of  www.nature.com/scientificreports/ C3 plants, we make α C4 ≈ 1. "Methods" section describe how to obtain α(t) from the changes in populations C3 and C4 plants. Using β = 0.46 and an increase in α from 0.84 to 0.925 over the last 10 Myr ago, we calculate that the Hopf bifurcation was crossed ca. 5 Myr ago and triggered a trajectory towards the limit cycle oscillation. The value of β depends on production of (CH 2 O). Higher numbers of photons, larger areas covered by biomass and hydrological cycles driven by latitudinal insolation gradients should lead to temporary increases in β. Inter-or intra-hemispheric insolation gradients show the same features as the 65°N insolation curve, namely the dominant 41 kyr period that is the signature of obliquity [48][49][50] . We incorporate photosynthetic-orbital forcing in the Glaciator making β(t) = β 0 [1 + 0.09F(t)], where the zero of F(t) is the average summer energy of the Northern Hemisphere over the last 5 Myr (ca. 5 GJ/m 2 ) 17 , and F(t) is normalized. The scaling factor 0.09 accounts for relative energy changes. A larger scaling factor could be used to account for the other mechanisms contributing to changes in β. Reasonably higher values did not change the nature of our results. This also means that choice of the insolation function is not critical for the results. The parameter β 0 = 0.46 was selected to reproduce climate changes over the last 5 Myr. Its value implies that respiration, burning or other processes associated with the consumption of the biomass are in the same timescale as the incorporation of CO 2 in the biomass.
The Glaciator can be regarded as a case of slow passage through a Hopf bifurcation 51,52 . Parameters α and β change slowly with respect to time and, as exhibited in the simulations, solutions stay near the unstable stationary state after the point (α, β) has crossed the curve of instantaneous Hopf bifurcation; oscillations emerge only after a delay. As already mentioned, Hopf bifurcation is crossed ca. 5 Myr ago, but large oscillations appear ca. 2 Myr later (Fig. 7). These slow passages through bifurcations have been observed in other climate models 6,53 . periodicities of past climate changes. The Glaciator has two adjustable parameters (β 0 and ω) and one parameter with constrained values (0.84 < α(t) < 0.925) over the last 5 Myr. The values of three parameters are imposed by their assigned physical meanings: ε ≈ 33ω, σ ≈ 3.3ρ and ρ = k 7 /(k 1 C) = 0.11. The characteristic concentration allows for the conversion between the adimensional variables x and y to H 2 O and CO 2 pressures. Using the parameters in Table 1, the Glaciator gives x = 95.5 at t = 0 in the absence of a CO 2 pulse (i.e. p H2O = 11.6 matm), and y = 2.19 (i.e., p CO2 = 262 ppm) at t = 0. These absolute values are in reasonable agreement with the pre-industrial CO 2 pressure of 280 ppm and with the water vapor pressure expected for the temperature of pre-industrial times discussed in "Methods" section. The calculated atmospheric CO 2 in the Last Glacial Maximum (LGM), p CO2 = 171 ppm, is in good agreement with the Antarctic ice-core record of 185 ppm 54 . These results come from the assignment of the timescales and not from an arbitrary scaling.  Table 1. The absolute values of CO 2 pressure were obtained using the characteristic concentration. The temperature was estimated adding changes in CO 2 radiative forcing to the pre-industrial global surface land temperature T s = 8°C 55 .
The average global surface temperature change (∆T) is often related to radiative forcing changes (∆R) externally imposed at the top of the atmosphere by the long-wavelength absorption of greenhouse gases. The most important of them is CO 2 , and its concentration changes have been associated with global temperature changes over the last 420 million years 56 . Temperature and radiative forcing are related by the climate sensitivity parameter S The change in radiative forcing when the concentration of CO 2 changes from that of a reference state (C 0 ) to the state studied (C) is given by 40 where the constant is in units of W/m 2 . Studies of past temperature changes showed that changing CO 2 concentrations by a factor of 2 is consistent with ∆T = 2.8 °C, which is within the range 2.3-3.0 °C suggested by various climate models 56 . Taking ∆T = 2.8 °C for ∆R = 3.7 W/m 2 , gives S = 0.76 °C/(W/m 2 ). An alternative approach for the long timescales considered here is to include slow feedbacks in the CO 2 sensitivity parameter and employ a higher value, S = 1.35 °C/(W/m 2 ) 1 . We adopted this approach for temperature calculations of the past.
The Mg/Ca ratios of foraminifera are one of the most reliable proxies of the paleo-temperatures of seawater 9,57 . A simple relation to obtain deep-sea temperature (in °C) was proposed 57 Figure 7 uses this relation and the Mg/Ca ratios available from a 1.5-million-year record 9 to obtain the corresponding Antarctic deep-sea temperatures. The benthic δ 18 O records of the LR stack 7 were aligned with T w from the Mg/Ca ratios to extend the estimated paleo-temperatures to 5 Myr ago.
Calculated global surface temperatures (T s ) are 3 to 5 °C higher than historical Antarctic deep-sea temperatures (T w ), as expected. Remarkably, 5 Myr ago both exhibited a dominant period of 42 kyr, which changed to a dominant 100 kyr period ca. 1 Myr ago. This reflects the change in the main driver of climate from orbital forcing due to obliquity to a limit-cycle oscillation. Figure 8 presents a detail of the MPT. Ca. 2 Myr ago a stronger 100 kyr period started imposing on the older 42 kyr period to become clearly dominant in the last 1 Myr. According to the Glaciator, this change results from the increase of α(t) in the Pleistocene with the expansion of C4 plants, which increased the efficiency of photosynthesis and led to the crossing of a Hopf bifurcation.
The relative changes in terrestrial carbon stock can also be compared with observations. The extreme values of z in the last cycle of the Glaciator are 216 and 278, i.e., CH 2 O changes by ca. 24%. A recent estimate of the total carbon stocks (soils and vegetation) is 2,807 PgC 58 . Data-based estimates of the difference between the LGM and pre-industrial land carbon storage range from 330 PgC 59 to 821 PgC 60 less carbon in the LGM, which correspond to deficits between 10 and 30%. Hence, the change in terrestrial carbon stock given by the oscillator is consistent with the current estimates. climate response to a co 2 pulse. Industrialization released approximately 300 PgC and business-asusual predictions indicate a total release of 1,000 PgC (471 ppm of CO 2 ) by the end of the century 42 . Climate models are often asked to predict the relaxation time of such a pulse on the current p CO2 level. In the Glaciator, the CO 2 pulse comes with a H 2 O pulse. The scenario in Fig. 9 is a simulation of an instantaneous pulse of 1,000 PgC at present time, with the corresponding increase in water vapor, followed by relaxation assuming that all other parameters remain constant. Figure 8 shows that the relaxation time is approximately 100 years but the www.nature.com/scientificreports/ decay is non-exponential. Half of the CO 2 pulse to the atmosphere is removed in 85 years. The IPCC estimate is that, for an emission pulse of about 1,000 PgC, about half is removed within a few decades 40 . Using the climate sensitivity parameter more appropriate for short timescales, S = 0.76 °C/(W/m 2 ), this CO 2 pulse corresponds to ∆T = 4.4 °C. Figure 9 also illustrates a major difference between the Glaciator and simulations that do not take into account the dominant 100-kyr period of climate change. The CO 2 pulse will eventually be absorbed and a new glacial cycle will begin. The next glacial maximum is predicted to occur in 50 kyr and be cooler than the preceding ones. This is due to the approach to the limit cycle shown in Fig. 4. The amplitude of the cycles depends on the final value of α, α(∞). Historical records show that the mid Pleistocene transition from 40-kyr to 100-kyr dominant cycles took more than half a million years. We estimate that the approach to the limit cycle started 5 Myr ago when α(t) became larger than 0.84 for β 0 = 0.46, and the Hopf bifurcation was first crossed. However, as discussed in "Methods" section, the attraction to the limit cycle just after the bifurcation is very weak.

Discussion
Our work addresses the internal Earth system mechanism that drives climate changes. We propose that this mechanism is an oscillator with the following steps: (i) photosynthesis uses CO 2 and H 2 O to generate highenergy intermediates that store free energy for an extended period of time; (ii) the stored energy is released at the same time as CO 2 and H 2 O that re-initialize the cycle; (iii) the minimal sufficient conditions for oscillatory behavior are met due to the presence of an autocatalytic step. According to this mechanism, the expansion of C4 plants, with the concomitant increase in radiation use efficiency of photosynthesis, triggered glaciation cycles.
The Glaciator with meaningful timescales leads to intrinsic oscillations consistent with the periods observed in records of historical climate changes, and to a good agreement with the absolute values of CO 2 and H 2 O in the atmosphere over the last 5 Myr. This simple conceptual model explains major puzzles of past climate dynamics: the absence of a strong 20 kyr precession signal 8 (Fig. 7F), the MPT of 41 kyr cycles to 100 kyr cycles 61 (Figs. 7C, F), the onset of glaciation cycles one million years ago with a dominant 100-kyr frequency 22 (Fig. 8), and the Holocene temperature conundrum 62 (Fig. 7D).
An unexpected prediction of the Glaciator is that further expansion of C4 plants will drive more pronounced glaciations in the future. Climate is a complex dynamic system and the simple conceptual model presented here only captures the most fundamental mechanisms underlying its dynamics. Nevertheless, it seems that extreme climate changes may happen in the next 50 kyr. Feedback mechanisms not included in the Glaciator may change this view, but our present understanding suggests that dampening the amplitudes of climate changes may be achieved recovering land from C4 plants to C3 plants.
The Glaciator is open to improvements as our knowledge of the various processes responsible for atmospheric CO 2 release and removal evolves. Regardless of possible refinements, in their present forms our oscillators reveal that closed photochemical systems are capable of producing sustained oscillations of their chemical compositions when absorbing light from a constant irradiation source.
Methods numerical simulations. Simulations were carried with Matlab R2018a. We used the ode45 solver to integrate each of the systems of differential equations: the Coimbrator and the Glaciator. This solver is based in a Runge-Kutta-Fehlberg method. In both cases we have fixed a relative error tolerance of 10 -12 .
The Coimbrator is a closed system of autonomous differential equations with a Hopf bifurcation. The Hopf bifurcation is supercritical and oscillations arise for parameter values below the bifurcation curve. The Glaciator is described by non-autonomous differential equations because the parameters α and β are time-dependent: α(t) = P C3 (t)α PC3 + P C4 (t)α PC4 and β(t) = β 0 [1 + 0.09F(t)]. Note that the expression for β(t) involves an unknown function F(t) related with the summer insolation in the Northern Hemisphere available in the literature 17 . Values of F(t) are available in steps of 1 kyr along the last 5 Myr. However, a sample of values of F is not enough for numerical simulations. We approach the values of F(t) at arbitrary t with two different methods. On one www.nature.com/scientificreports/ hand, when running the Glaciator to get simulations of the last 5 Myr we use a linear interpolation in each time interval corresponding to consecutive sample points. Alternatively, we have used the same sample points to get a Fourier series approaching F(t) in the whole period of 5 Myr. The good fitting obtained is illustrated in Fig. 1A. Both methods are used to obtain the simulations for the past, but no remarkable sensitivity is observed.
To get extrapolations of the insolation values in the future, we have used the approach of the Fourier series. In the case of the Glaciator, when running the simulation in the past, initial time is fixed in the value − 5 Myr and the numerical integration extends up to present time. Initial conditions for the dependent variables are chosen close to the equilibrium point of the system when α(t) = P C3 (t)α PC3 + P C4 (t)α PC4 and β = β 0 , namely, x = 133.8721, y = 1.710171, z = 245.4826. It is convenient to express these equations in terms of dimensionless variables to obtain a general solution. We follow the procedure of Field and Noyes 32 to cast the concentrations of X, Y and Z in the dimensionless variables x, y, z. We define g = k 1 C as a pseudo-first order rate under continuous irradiation if A is not depleted from the system. Thus, the quotient g/k 2 has the dimensions of a concentration and is defined as the characteristic concentration. This is used to make the following dimensionless measures of the species:

Reduction of variables in the
The characteristic concentration g/k 2 can later be used to obtain the values of X, Y and Z in units of concentration. The same procedure is applied to the conversion of time t in dimensionless variable τ. The derivatives become Making the appropriate replacements we obtain a simplified system of differential equations, in dimensionless units  (25) is satisfied along the curve defined by (24) for all β > 0. These eigenvalues are ± = ±i √ 2β − α . The third eigenvalue λ 3 = -(β + 2) is negative for all β > 0. The distinction between subcritical and supercritical Hopf bifurcation depends on the value of the first Lyapunov coefficient 36 . Although such coefficient can be obtained for the given system, calculations are quite laborious and we prefer to use MATCONT 63 , a MATLAB numerical continuation package for bifurcation analysis, to check that the Hopf bifurcation is supercritical (see Fig. 3). [X]-a step in the biosynthesis of CH 2 O (e.g., condensation reactions in the formation of a glycosidic bond to make cellulose or in the formation of a peptide bond to yield a protein), k 5 [X]-condensation of water, k 6 [X 0 ]-evaporation of water, k 4 [Z]-combustion/respiration, k 7 [Y]-loss of CO 2 into the ocean leading to silicate weathering, k 7 [Y]-release of CO 2 from the ocean in all its forms. As mentioned above and further discussed below, 1/(k 1 C) is the characteristic time and is associated with more than one physical process in this simplified description of the Earth system.

Reduction of variables in the Glaciator. The Glaciator is defined by the reactions
Making, as for the Coimbrator, g = k 1 C, the quotient g/k 2 is defined as the characteristic concentration. Additionally, 1/g is defined as the characteristic time. The characteristic time and concentration are used to obtain the dimensionless variables In terms of the adimensional variables, we obtain where (26) www.nature.com/scientificreports/ Considering that k 6 x 0 and k 8 y 0 are proportional to the water evaporation and to the CO 2 desorption fluxes, respectively, the following relation must be obeyed The desorption flux of CO 2 from water at 313 K is j CO2 ≈ 7 × 10 -4 mol/(m 2 s) 64 . At this temperature, the water evaporation flux from seawater ranges from 1.5 × 10 -2 to 4.6 mol/(m 2 s), depending on the air velocity 65 . Using j H2O ≈ 2.3 × 10 -2 mol/(m 2 s) we estimate ε ≈ 33ω.
constrains of the parameters in the Glaciator. The saturation vapor pressure can be calculated from an approximation to the Clausius-Clapeyron equation where the constants where chosen to yield the vapor pressure in mbar, T is in K and a = 0.064 K -1 66 . Using the average Earth surface land temperature T s = 8 °C to represent pre-industrial global temperature (ca. 1900) 55 , we obtain e 0 * ≈ 10.9 mbar, which corresponds to p H2O ≈ 11 matm.
Water vapor equilibrates relatively rapidly with liquid water (k 5 x ≈ k 6 x 0 ). Assuming CO 2 in the atmosphere also equilibrates relatively rapidly with CO 2 dissolved in the sea (k 7 y ≈ k 8 y 0 ), the following relations can be established Given that we reduced the number of variables needed to apply the Glaciator making x ≈ 10y and using ε ≈ 33ω, to set σ ≈ 3.3ρ.
Hopf bifurcation in the Glaciator. The Glaciator, as the Coimbrator, has two equilibrium points. We have to study the following system of nonlinear equations From the third equation we obtain After substituting z in the second equation we get and finally, substituting y in the first equation, we obtain the following equation for x Since [(α -ρ)(ρ + 1) -ω + ε) 2 + 4ε (ρ + 1) (α + σ)] > 0, we always have two equilibria. In this case the analytical characterization of the Hopf bifurcation is plausible but more involved. Using again, as we did for the Coimbrator, the package MATCONT, we obtain the results depicted in Fig. 6, for which we have fixed the values  www.nature.com/scientificreports/ As shown in Fig. 10, the first Lyapunov coefficient is negative in the whole range of parameters considered in Fig. 6 and hence the Hopf bifurcation is supercritical. A periodic orbit in a three-dimensional phase space has Floquet multipliers 1, m 1 and m 2 (m 1 and m 2 are the eigenvalues of the linearization of the first return map defined around the periodic orbit). If |m i |< 1, for i = 1,2, then the periodic orbit is attracting. The right panel in Fig. 10 shows the value of Floquet multipliers for the periodic orbits born at α = 0.84 when β = 0.46. Because the first Lyapunov coefficient is very small, the attraction of the limit cycle just after the bifurcation is very slow. For the computation of the first Lyapunov coefficient, instead of using the approximations obtained with MATCONT, we wrote our own code, designed for this specific model, to obtain more accurate numerical results.
Populations of C3 and C4 plants. The  where N max is the population maximum at t = ∞, N 0 is the population at -t 0 , which is selected to be the origin of the expansion, and r is the rate of population increase. The sudden C4 plants expansion occurred ca. 10 Myr ago 44 , which leads to t 0 = 10 Myr. N max = 0.25 was selected in view of the fact that C4 plants account for 25% of terrestrial photosynthesis. With these values of t 0 and N max , the known profile of C4 plants expansion 44 can be reproduced with the population N 0 = 0.15 and the rate r = 0.3 Myr -1 . The expansion of C4 plants was achieved, in part, at the expense of C3 plants. We modelled the reduction of C3 plants population as where the factor 0.15 was selected to reflect evolution from nearly exclusive C3 vegetation 20 Myr ago to the composition of biomass today: 389.3 PgC of C3 vegetation and 18.9 PgC of C4 vegetation 67 . This corresponds to a fraction of C3 biomass of 0.954, whereas Eq. (43) gives P C3 (0) = 0.96. Equation (43) assumes that only 15% (40) ω = 5, ρ = 0.11, ε = 33ω, σ = 3.3ρ