Excitability in ramped systems: the compost-bomb instability

S. Wieczorek, P. Ashwin, C. M. Luke, P. M. Cox

Abstract

The paper studies a novel excitability type where a large excitable response appears when a system’s parameter is varied gradually, or ramped, above some critical rate. This occurs even though there is a (unique) stable quiescent state for any fixed setting of the ramped parameter. We give a necessary and a sufficient condition for the existence of a critical ramping rate in a general class of slow–fast systems with folded slow (critical) manifold. Additionally, we derive an analytical condition for the critical rate by relating the excitability threshold to a canard trajectory through a folded saddle singularity. The general framework is used to explain a potential climate tipping point termed the ‘compost-bomb instability’—an explosive release of soil carbon from peatlands into the atmosphere occurs above some critical rate of global warming even though there is a unique asymptotically stable soil carbon equilibrium for any fixed atmospheric temperature.

1. Introduction

An excitable system remains in a quiescent state, if undisturbed, produces a small response to a ‘small’ stimulus but fires a large transient response when the small stimulus exceeds a certain threshold. The notion of excitability was first introduced in biology and physiology in an attempt to understand the spiking behaviour of neurons (Hodgkin & Huxley 1952; FitzHugh 1955, 1961) and their electronic simulators (Nagumo et al. 1962). Soon after this work, excitability was found in chemical reactions (Zaikin & Zhabotinsky 1970; Ruoff 1982). More recently, there have been a number of theoretical and experimental demonstrations of excitability in liquid crystals (Coullet et al. 1994), optical systems including lasers (Dubbeldam et al. 1999; Yacomotti et al. 1999; Tredicce 2000; Wieczorek et al. 2002; Wünsche et al. 2002; Krauskopf et al. 2003; Goulding et al. 2007) and photonic crystals (Yacomotti et al. 2006). A hallmark of excitability is a genuine or apparent discontinuity in the system’s response versus the stimulus strength (FitzHugh 1955; Hoppensteadt & Izhikevich 1997; Doi et al. 1999; Izhikevich 2006). It has become clear that this strongly nonlinear phenomenon is mathematically intriguing and relevant to different fields of science. In this paper, we demonstrate a new form of excitability that is relevant to the stability of peatlands under global warming.

From the original notion of excitability, which is purely phenomenological, one can conclude that an excitable system has to have the following properties:

  • (P1) A quiescent state.

  • (P2) An excitability threshold above which the system initially evolves away from the quiescent state giving rise to an excitable response.

  • (P3) A return mechanism that specifies the type of excitable response and, when the stimulus is off, brings the system back to the quiescent state so it can be excited again.

A typical candidate for a stimulus that perturbs the system from the quiescent state to above the excitability threshold is a fast and large enough change in one of the system parameters, a short impulse, for example (Wünsche et al. 2002). Other possibilities include stochastic fluctuations in the system variable(s) (Lindner et al. 2004). This work studies excitability owing to parameter ramping—a steady, slow and monotonic change in one of the system parameters referred to as the ramped parameter. In the problem considered here, a (unique) quiescent state exists for any fixed setting of the ramped parameter. However, a very large excitable response appears when the parameter is ramped sufficiently fast from one setting to another.

(a) Summary of main results

In §2, we define the excitability phenomena in rigorous terms of dynamical systems, give a brief survey on classification of excitability, and describe two different types of excitability, namely type A and type B. This work focuses on an in-depth analysis of type B excitable models, and §3 describes excitability properties in these models with state jumps.

In §4, we study a general class of type B excitable models with parameter ramping, of which the compost-bomb problem is a representative. We show that if a suitable parameter is ramped in a slow–fast system with an asymptotically stable equilibrium and locally folded critical (slow) manifold, then there may be a critical value of the ramping rate above which an excitable response appears. This result is obtained in two steps using concepts from singular perturbation theory. In the first step, we show that a necessary and sufficient condition for the existence of a critical ramping rate is a folded saddle singularity in the corresponding desingularized system. In the second step, we derive an analytical condition for calculating the critical ramping rate.

The general analysis in §4 is motivated by a need to understand the response of peatlands to global warming (or atmospheric temperature ramping), which represents a potential tipping point in the response of the climate system to anthropogenic forcing (Lenton et al. 2008). It is estimated that peatland soils contain 400–1000 billion tonnes of carbon, which is of the same order of magnitude as the carbon content of the atmosphere. Peat carbon is increased by plant litter and reduced by microbial decomposition in the soil. Peat decomposition is expected to accelerate under global warming, leading to concerns that carbon could be released to the atmosphere, accelerating the rate of carbon dioxide increase and providing a positive feedback to global warming (Cox et al. 2000; Khvorostyanov et al. 2008b). There is also strong empirical evidence of peat instability in response to warm and dry climate anomalies, such as the peatland fires in Russia in summer 2010.

Although many global climate-carbon cycle models predict loss of soil carbon under global warming (Friedlingstein et al. 2006), few properly deal with organic soils such as peats, and all ignore the effects of biochemical heat release associated with microbial decomposition (Khvorostyanov et al. 2008a). In their recent work, Luke & Cox (in press) define the peatland soil system in terms of its vertically integrated soil carbon content, C (kg m−2) and soil temperature, T (°C). Soil carbon is increased by litter fall from plants, Π=1.055 (kg m−2 yr−1), and reduced by microbial decomposition, which depends on the store of carbon, C, and also the temperature sensitivity of the decomposition rate per unit carbon, r(T), Embedded Image1.1Empirical studies suggest that r(T) is approximately exponential in temperature, such that it increases by a factor of 2 to 5 for each 10°C of warming (Kirschbaum 1995). We choose Embedded Imagewith r0=0.01 (yr−1) and Embedded Image/10 (°C−1). Soil temperature relaxes back to the prescribed atmospheric temperature, Ta, with a time scale dependent on the thermal inertia, μ=2.5×106 (J m−2°C−1), and the soil-to-atmosphere heat transfer coefficient, λ=5.049×106 (J yr−1 m−2°C−1). Soil temperature is increased by microbial respiration, Cr(T), with a constant of proportionality, A=3.9×107 (J kg−1), derived from the enthalpy of the respiration reaction (Thornley 1971). This gives the soil temperature equation Embedded Image1.2with a small parameter ϵ=μ/A≈0.064 (kg °C−1 m−2). Global warming is approximated by an atmospheric temperature ramp Embedded Image1.3with a constant rate v in units of °C yr−1. Equations (1.1)–(1.3) are the climate-carbon cycle model with global warming. The system (1.1)–(1.2) has a unique asymptotically stable equilibrium for any fixed atmospheric temperature Ta. However, numerical simulations of equations (1.1)–(1.3) suggest that biochemical heat release could destabilize peatland above some critical rate of global warming, vc. The result is an explosive increase in the soil temperature accompanied by a catastrophic release of peatland soil carbon into the atmosphere, which has been termed the ‘compost-bomb instability’ (Luke & Cox in press). Until now, the reason for the rate dependence of the compost-bomb instability has not been rigorously understood, but the present paper shows that this arises from a novel form of excitability.

In §5, we use the results of §4 to show that for an initial condition at the unique stable equilibrium of equations (1.1)–(1.2) and ϵ small, the climate-carbon cycle model (1.1)–(1.3) undergoes a very large excitable response (the compost-bomb instability) if the global warming rate exceeds the critical value Embedded Image1.4where T0a is the initial atmospheric temperature and Embedded Imageis a dimensionless coefficient.

2. Excitable dynamical systems

To define the excitability phenomenon in the language of dynamical systems, we consider an n-dimensional autonomous dynamical system referred to as a stimulus-free problem Embedded Image2.1where Embedded Image is the state vector and Embedded Image the (constant) parameter vector. A stimulus is represented by some prescribed time dependence of the parameter vector, Embedded Image, so the corresponding non-autonomous stimulus problem reads Embedded Image2.2To be able to express the phenomenological properties (P1)–(P3) of equations (2.1) and (2.2) in more rigorous terms, let Embedded Image denote a ball of radius r about a(P(t)) such that a(P) is the only attractor for system (2.1) within Br[a(P)], and define:

Definition 2.1

A quiescent state Embedded Image is an asymptotically stable state for the stimulus-free problem (2.1).

Definition 2.2

An excitable response (with respect to suitably chosen δ and σ such that 0<δσ) is a trajectory Embedded Image that starts within Bδ[a(P(0))], leaves Bσ[a(P(t))] for some time, before entering into Bδ[a(P(t))]. As Embedded Image, X(t) may remain within Bδ[a(P(t))] or leave and return to Bδ[a(P(t))] repeatedly.

Definition 2.3

Given P(t), an excitability threshold is a boundary in the phase space Embedded Image separating initial conditions X(0) that give excitable responses from those that give no excitable responses.

Definition 2.4

The stimulus-free problem (2.1) is called excitable if it has properties (P1)–(P3).

A typical candidate for a quiescent state is a stable equilibrium but other possibilities include small-amplitude periodic orbits or even ‘small’ chaotic attractors (Marino et al. 2004, 2007; Al-Naimee et al. 2009). An excitable response in definition 2.2 is a trajectory that (repeatedly) moves far enough from a(P(t)) before coming close to a(P(t)), where ‘far enough’ is quantified by the choice of σ. In definition 2.3, if the time-dependent component of P(t), namely Embedded Image, is a solution to an autonomous differential equation, pr can be treated as an additional state variable. Then, the excitability threshold is a boundary in the phase space Embedded Image separating initial conditions (X(0),pr(0)) that give excitable responses from those that give no excitable responses. Unlike a quiescent state, an excitability threshold can be defined in both the stimulus free (2.1) and the stimulus (2.2) problems. However, the threshold may change in the presence of a stimulus and it may depend on the form of P(t).

(a) Classification of excitability

Here, in the spirit of the original work by FitzHugh (1955), we classify excitability by the type of threshold and the resulting increase in the system’s response versus the stimulus strength. Specifically, we distinguish two excitability types. In type A excitability, there is a unique threshold given by the stable manifold of an unstable state, and an increase in the excitable response is discontinuous (related to the singular-point threshold phenomenon of FitzHugh (1955)). In type B excitability, the threshold is not unique, meaning that it weakly depends on the choice of σ in definition 2.2, and an increase in the response is abrupt but continuous (related to the quasi-threshold phenomenon of FitzHugh (1955)). Figure 1 shows examples of type A and B systems that produce excitable responses following a state jump from below to above the excitability threshold. Both examples have the time-scale separation resulting in slow–fast dynamics near a. However, a different shape of the critical manifold (green in figure 1) approximating the slow dynamics, and the fact that the type A system has an unstable (saddle) state near a while the type B system does not, give rise to different types of excitability threshold.

Figure 1.

Examples of phase portraits illustrating two different excitability types in system (2.1) with a stable equilibrium a: (a) type A excitability near a saddle point s and (b) type B excitability near a fold L of critical manifold S=SaLSr. Shown are trajectories starting at initial conditions below and above the excitability threshold. The excitable response depends on the return mechanism (dashed) specified by the form of f in (a) and the form of g in (b).

In an example of the type A excitability from figure 1a, the unique excitability threshold for the stimulus-free problem is the stable invariant manifold Ws(s) of the saddle s. For state jumps of magnitude Δ, the plot of the maximum (Euclidean) distance between X(t) and a(P(t)) following a jump versus Δ has a discontinuity defining a unique threshold value for Δ. The return mechanism specified by f(x,p) is less obvious. Previous studies focused on systems near to a homoclinic bifurcation (Kuznetsov 1995). After a single and suitably chosen jump from a to above Ws(s), the system follows the upper branch of Wu(s) and produces a single-pulse response near to a 1-homoclinic bifurcation (Plaza et al. 1997; Ventura et al. 2002; Krauskopf et al. 2003), a n-pulse response near to a n-homoclinic bifurcation (Wieczorek et al. 2002), or an irregular and unpredictable response near to the so-called chaotic Shil’nikov case. An example of the type B excitable system from figure 1b does not have a unique excitability threshold, except in the singular limit Embedded Image, as explained in §3. This system is the main focus of this work and its in-depth analysis is presented in §3 for state jumps and in §4 for parameter ramping.

Here, we consider an excitability problem where the stable state a(P) for system (2.1) exists for any fixed value of the ramped parameter P. More recently, various classifications of excitability have emerged in the context of a different but related problem of excitable bursting where the stable state a(P) for system (2.1) loses stability or disappears at a certain setting of the ramped parameter P (Rinzel & Ermentrout 1989). For example, Hoppensteadt & Izhikevich (1997) and Izhikevich (2006) distinguish between class 1 and class 2 neuronal excitability depending on the bifurcation of the quiescent state and the associated continuous or discontinuous increase in the frequency of bursting. Similarly, Gerstner & Kistler (2002) distinguish between type I and type II excitabilities that are the same as Hoppensteadt & Izhikevich’s class 1 and 2 excitabilities, respectively. More generally, Golubitsky et al. (2001) classify excitable bursters by the codimension of the bifurcation in whose unfolding they first appear. Clearly, class 1 and 2, type I and II as well as Golubitsky’s excitabilities cannot be directly compared to our type A and B excitabilities because they refer to a different dynamical problem. Nonetheless, different excitability problems can be related in the following sense. For example, type A excitability occurs near a saddle-node-homoclinic bifurcation (Kuznetsov 1995) that gives rise to class 1/type I excitability, and type B excitability occurs near a singular Hopf bifurcation followed by a canard explosion (Benoît et al. 1981) that gives rise to class 2/type II excitability (Gerstner & Kistler 2002; Izhikevich 2006).

3. Type B excitability with state jumps

In this section, we fix P and discuss excitability in the (autonomous) stimulus-free problem (2.1) with instantaneous jumps of magnitude Δ in the state vector X. This is a convenient framework for discussing excitability thresholds and the novel return mechanism found in the climate-carbon cycle model (1.1)–(1.2). It is also a good starting point for the analysis of the more complicated problem (2.2) with a time-dependent parameter in §4.

We focus on an example of the type B excitable behaviour shown in figure 1b in a singularly perturbed system of the form Embedded Image3.1with slow variable Embedded Image, fast variable Embedded Image, and parameter vector Embedded Image. There is a singular perturbation parameter 0<ϵ≪1, meaning that the system evolves on fast (t/ϵ) and slow (t) time scales. It is useful to introduce a critical manifold Embedded Imagethat approximates the slow dynamics and is the dz/dt=0 isocline when g does not depend on ϵ. Assume that system (3.1) satisfies the following conditions:

  • (A1) There is an asymptotically stable equilibrium Embedded Image

  • (A2) ∂g(x,z,P,0)/∂x≠0 for any (x,z)∈S(P) so that the critical manifold can be written as x=h(z,P). Furthermore, near to a(P,ϵ), the critical manifold has a fold transverse to the (slow) x direction, L(P), where ∂h/∂z=0 and ∂2h/∂z2≠0. Differentiating the critical manifold condition, g(h(z,P),z,P,0)=0, with respect to z shows that ∂h/∂z=0 is the same as (∂g(x,z,P,0)/∂z)|x=h(z,P)=0, and ∂2h/∂z2≠0 requires ∂2g(x,z,P,0)/∂z2|x=h(z,P)≠0. Therefore, Embedded Image

  • (A3) Initial conditions from within Embedded Image in the (x,z) phase space converge to a(P,ϵ) as Embedded Image.

System (3.1) satisfying assumptions (A1)–(A3) has properties (P1)–(P3) so we call this system excitable. In particular, fold L(P) provides an excitability threshold as shown in figure 1b, and g(x,z,P,ϵ) specifies the return mechanism, that is the shape, magnitude and duration of excitable responses. We will now describe the excitability threshold owing to a locally folded critical manifold, and the new return mechanism found in the climate-carbon cycle model (1.1)–(1.2).

(a) Excitability threshold

In contrast to the example from figure 1a, the excitability threshold in the singularly perturbed system (3.1) from figure 1b is not unique. This property is best understood by looking at the limiting problem Embedded Image on different time scales (Arnol’d 1994; Szmolyan & Wechselberger 2001). On the slow time scale t, one obtains the reduced system: Embedded Image3.2that describes the evolution of the slow variable x on the critical manifold S(P). In figure 1b, S(P) is partitioned into an attracting part Sa(P), a repelling part Sr(P), and the fold point L(P). On the fast time scale τ=t/ϵ, one obtains the layer system: Embedded Image3.3that describes the evolution of the fast variable z for fixed x, in fast time τ=t/ϵ. Sa(P) and Sr(P) correspond to the stable and unstable equilibrium, respectively, of the layer system. A fold L(P) of S(P) transverse to the slow x direction corresponds to a saddle-node bifurcation for the layer system. Therefore, in the limit Embedded Image there is a unique excitability threshold, Embedded Image (figure 2a), where Embedded Image is the left branch of the centre manifold of the saddle-node equilibrium L(P) for the one-dimensional layer system (3.3).

Figure 2.

Sketches of phase portraits near locally folded critical manifold S=SaLSr (green) illustrating excitability threshold in (a) the reduced system (3.2) and (b) the singularly perturbed system (3.1) with canard trajectories (blue). Embedded Image denotes the centre manifold of the saddle-node equilibrium L for the one-dimensional layer system (3.3), a denotes the stable equilibrium for system (3.1), Sa,ϵ and Sr,ϵ (red) denote the attracting and repelling part of the slow manifold for system (3.1), respectively.

The geometric singular perturbation theory describes how to paste together the dynamical behaviour of the fast and slow limiting problems to obtain dynamics for the singularly perturbed problem with sufficiently small but non-zero ϵ (Fenichel 1979; Jones 1995). Away from L(P), the normally hyperbolic manifolds Sa(P) and Sr(P) of hyperbolic equilibria for the layer system persist for sufficiently small ϵ as invariant (and still normally hyperbolic) attracting Sa,ϵ(P,ϵ) and repelling Sr,ϵ(P,ϵ) slow manifolds (shown in red in figure 2b) for the singularly perturbed system (3.1). The smaller the ϵ, the closer to L(P) is this persistence valid. In figure 2, the non-hyperbolic equilibrium L(P) for the layer system is the only place where normal hyperbolicity fails on S(P). Near to L(P), Sa,ϵ(P,ϵ) and Sr,ϵ(P,ϵ) typically split, allowing for canard trajectories or ducks that follow the repelling slow manifold for a considerable amount of time before rapidly moving away from Sr,ϵ(P,ϵ) in the fast z direction (figure 2b; Benoît et al. 1981; Benoît 1983; Arnol’d 1994; Krupa & Szmolyan 2001; Szmolyan & Wechselberger 2001). Owing to the ϵ-dependent splitting between Sa,ϵ and Sr,ϵ and the associated continuum of canard trajectories, the unique excitability threshold of the limiting problem does not persist in the singularly perturbed system (3.1). Rather, for different choices of σ in definition 2.2, the excitability threshold is given by a different canard trajectory. One may want to choose σ such that the excitability threshold is the canard trajectory that follows Sr,ϵ(P,ϵ) for the longest time. In contrast to type A excitability, the plot of max[z(t)/za] versus the jump magnitude Δ is continuous giving no unique threshold value for Δ, where za denotes the z component of a(P,ϵ).

We illustrate these effects for the excitable system Embedded Image3.4with Embedded ImageFigure 3 shows the response of system (3.4) to jumps of magnitude Δ in x, with dependence on the singular perturbation parameter ϵ and the order N of the specific polynomial function Embedded Image. While the increase in the response remains distinct as an apparent discontinuity for ϵ≤10−2 and N≥2, its position on the Δ-axis varies noticeably with ϵ≥10−3 and N≤4. The unusually large response max[z(t)/za] in figure 3a is shown to scale as ϵ−1 owing to the novel return mechanism described in the following section.

Figure 3.

Response of the excitable system (3.4) to jumps in x of magnitude Δ. We used p1=0.5, p2=0.1, p3=−1 and Embedded Image. (a) Responses for N=3 and different ϵ. (b) Responses for ϵ=10−2 and different N=2,3,4, and Embedded Image (from right to left). Note the logarithmic vertical scale in both panels.

(b) The return mechanism

The most commonly studied return mechanism is found in variants of the excitable van der Pol and FitzHugh–Nagumo equations (van der Pol 1920; FitzHugh 1961; Nagumo et al. 1962; Doi et al. 1999; Izhikevich 2006) and can be modelled using excitable systems of the form: Embedded Image3.5Specifically, for Q(z,P)=p7z3+p6z2+p5z+p4 and Embedded Image, the return mechanism is given by a cubic or an N-shaped critical manifold x=−Q(z,P). This is illustrated in figure 4a for Q(z,P)=−(z−1.2)3+z−1.2. Upon perturbing such a system from the stable equilibrium a to above the excitability threshold where Embedded Image, the trajectory rapidly moves away from a towards larger z at a speed ∼O−1) until it slows down by the attracting part Sa2 of critical manifold S (figure 4a). Then, the trajectory crosses Sa2 towards smaller x and moves along Sa2 at a speed ∼O(1) until it reaches fold L2. Next, the trajectory rapidly moves towards smaller z, slows down by Sa1, crosses Sa1 towards larger x and approaches a along Sa1. Such an excursion in the phase space results in a square-like excitable response (figure 4b) whose magnitude is given approximately by the horizontal distance between L1 and L2, Embedded Imageand whose duration is ∼O(1).

Figure 4.

Different return mechanisms ((a) and (c)) give rise to different types of excitable response ((b) and (d)). Shown are trajectories starting at initial conditions below (blue) and above (red) the excitability threshold, and the two nullclines (green). (a), (b) are obtained using equation (3.5) with p1=0.5, p2=p3=−1, ϵ= 10−3 and Q(z,P)=−(z−1.2)3+z−1.2. (c)–(e) are obtained using equation (3.4) with p1=0.5, p2=0.1, p3=−1, ϵ=10−3 and Q(z,P)=z3/3!+z2/2!+z+1. (e) expanded view around the right-most turning point of the (red) excitable trajectory from (c). Note the logarithmic vertical scales in (c)–(e).

The same type of excitability threshold with a distinctively different return mechanism is found in the climate-carbon cycle model (1.1)–(1.2) and system (3.4). In system (3.4), the critical manifold x=−(p2+p3z)/Q(z,P) grows monotonically with z, then folds and decays monotonically to zero if p4,p6>0, p5≥0 and pj≥0 for j=7,8,…. This is illustrated in figure 4c for Q(z,P)=z3/3!+z2/2!+z+1. Upon perturbing such a system from a to above the excitability threshold where Embedded Image, the trajectory rapidly moves away from a towards larger z but it does not encounter any slow manifold that would abruptly stop its predominantly horizontal motion. Rather, for small enough ϵ, the trajectory converges rapidly to a line x=Kϵz, where K>0, and moves along this line until it arrives at approximately Embedded Image, where Embedded Image. To see this, assume that the trajectory starts at (z0,x0) at time t0, and P is such that Q(z,P)>Kzn for some K>0 and n≥3 to the right of z0. Then, for fixed x≠0, the term (−1,ϵ−1)xQ(z,P) in the vector field for (dx/dt,dz/dt) will dominate for large enough z, to the extent that away from x=0, trajectories will run parallel to a line x=Kϵz. As the trajectory approaches x∼0, it sharply turns around by crossing the two nullclines dz/dt=0 and dx/dt=0 (figure 4e). Next, the trajectory moves towards smaller z, slows down by Sa, crosses Sa towards larger x, and approaches a at a speed ∼O(1) along Sa. Such an excursion in the phase space results in an excitable response that is much shorter and has a much larger amplitude when compared with the response in systems (3.5). As Embedded Image, the response amplitude becomes Embedded Imagewhere K depends linearly on x0. This result explains the ϵ−1 scaling of the amplitude of excitable responses plotted in figure 3a. The time at which maximum is reached can be estimated using the time to blow-up for the singular system and will be Embedded Image.

4. Type B excitability with parameter ramping

This section discusses excitability in the (non-autonomous) stimulus problem (2.2), where one component of the parameter vector P(t), namely Embedded Image, is ramped at a rate v, Embedded Image4.1and all the remaining components of P(t), including Embedded Image and a small parameter Embedded Image, do not vary in time. Henceforth, we refer to pr(t) as the ramped parameter, to the stimulus-free problem (2.1) as the unramped system, and to the stimulus problem (2.2) as the ramped system. In addition to assumptions (A1)–(A3), we make two extra assumptions:

  • (A4) v is constant and, without any loss of generality, v>0.

  • (A5) In the unramped system (2.1), given by setting v=0, a stable equilibrium a(P) and excitability threshold, e.g. due to a fold L(P) in type B systems from figure 1b, exist for any fixed setting of the ramped parameter pr(t).

Whether the ramped system (2.2) produces an excitable response depends on the initial condition X(0) and on pr(t). This brings us to an interesting question: what is the excitability threshold in ramped systems that preserve a quiescent state and, given X(0) and pr(0), what is the critical rate, vc, required to trigger an excitable response? In the context of the climate-carbon cycle problem with global warming (1.1)–(1.3), this question has relevance to the notion of dangerous rates of climate change.

We focus on type B excitable systems from figure 1b with parameter ramping. If v is small compared with the fast time scale of system (3.1), that is v≪ϵ−1, then pr(t) becomes an additional slow variable and one has to consider a singularly perturbed problem Embedded Image4.2with two slow variables Embedded Image and one fast variable Embedded Image. System (4.2) has no equilibrium points when v≠0. Rather, a one-dimensional manifold Embedded Imageindicates the position of the stable equilibrium of the unramped system (3.1) for different fixed settings of pr(t) in the (x,z,pr) phase space.

In §3, we showed that the excitability threshold in the unramped system (3.1) is unique only in the limiting case Embedded Image. Guided by this result we seek answers to the questions about the excitability threshold and critical rate in the ramped system (4.2) by analysing the ‘slow’ reduced system given by Embedded Image4.3With q denoting f or g and k denoting x, z or pr, we simplify notation by introducing Embedded ImageIn the (x,z,pr) phase space, the reduced system evolves on a folded two-dimensional critical manifold Embedded ImageNear to a(p,ϵ), S(p) is partitioned into the attracting part Sa(p), repelling part Sr(p), and a one-dimensional fold transverse to the (slow) x direction Embedded Imagewhich follows from assumptions (A2) and (A5). Furthermore, S(p) is a graph over pr and z, Embedded Image4.4which follows from Embedded Image in assumption (A2) and the implicit function theorem.

Since excitable responses appear in the fast variable z (figure 4cd), we need to study the evolution of z on S(p) in slow time t, and this is obtained by differentiating the algebraic constraint in equations (4.3) with respect to t, Embedded Image4.5Then, the reduced system (4.3) becomes Embedded Image4.6Furthermore, it is useful to project (4.6) onto the (z,pr) plane using equation (4.4) to get the projected reduced system Embedded Image4.7Note that system (4.7) is singular at a fold, L(p), where Embedded Image. As a result, z(t) blows up (diverges off to infinity in finite time t) when trajectories reach typical points on L(p). However, trajectories that approach special points on L(p) where Embedded Image in a (hyperbolic) direction so that dz/dt remains finite may cross a fold with finite speed (Benoît 1983). Such special points Embedded Imageare called folded singularities (Takens 1976; Arnol’d 1994), pseudo singularities (Argemi 1978) or canard points (Szmolyan & Wechselberger 2001). In graphical terms, folded singularities are intersection points between one-dimensional manifolds L(p) and Embedded Imagein the (x,z,pr) phase space. To eliminate certain degenerate cases, we assume that

  • (A6) Embedded Image is a graph over pr and can be expressed as z=u(pr).

The associated trajectories that cross from Sa(p) to Sr(p) via F(p,v) are called singular canards. Folded singularities are absolutely essential to the understanding of the excitability threshold in slowly ramped systems that preserve stable equilibrium and, to make progress, we need to analyse the flow near folded singularities of the projected reduced system (4.7). The analysis is greatly facilitated by the following scaling also known as desingularization (Dumortier & Roussarie 1996; Krupa & Szmolyan 2001), Embedded Image4.8which preserves the direction of time on Sa(p) but reverses it on Sr(p). The scaling gives the desingularized system, Embedded Image4.9Clearly, a folded singularity of the projected reduced system (4.7) is a regular equilibrium of the desingularized system (4.9). This means that we can study the phase portrait of the desingularized system (4.9) in a neighbourhood of F(p,v) using standard techniques. Then, we obtain the phase portrait of the projected reduced system (4.7) simply by changing the direction of time on Sr(p) in the phase portrait of (4.9). In general, folded singularities are classified as folded nodes, folded foci, folded saddles and folded saddle-nodes, based on their classification as equilibria of the desingularized system. The relevant problem of possible types of folded singularities and singular canards in Embedded Image with a folded two-dimensional critical manifold was addressed by Szmolyan & Wechselberger (2001).

In the phase portrait of the projected reduced system (4.7), one should be able to identify up to three different types of initial condition within Sa(p). The first type consists of initial conditions for which trajectories approach L(p) and then blow-up along a fast direction. They represent initial states of the ramped system that lead to excitable responses. The second type consists of initial conditions for which trajectories remain within Sa(p) for all time. If these trajectories additionally approach and remain near a(p,ϵ) as Embedded Image, we say that the ramped system adiabatically follows the stable but changing equilibrium of the unramped system. The third type is initial conditions that give rise to canard trajectories crossing from Sa(p) to Sr(p) via the folded singularity F(p,v). According to definition 2.3, a boundary separating regions of the first and second type is the excitability threshold within Sa(p).

(a) Necessary and sufficient condition for critical rate

Given assumptions (A1)–(A6), a phase portrait of the projected reduced system (4.7) has the following properties:

  1. a(p,ϵ)∩L(p)=∅ or a(p,ϵ)⊂Sa(p) (an intersection between a(p,ϵ) and L(p) would correspond to a Hopf bifurcation (Benoît et al. 1981) of the stable equilibrium for the unramped system (3.1) and violate (A5)).

  2. Embedded Image contains no folds transverse to the x direction nor self-intersections by assumption (A6).

  3. A folded singularity F(p,v) corresponds to an intersection between L(p) and Embedded Image; this follows from equations (4.7) and definitions of L(p) and Embedded Image.

  4. The vector field points in the direction of increasing pr for all (z,pr)∈S(p)\(L(p)\F(p,v)); this follows from equations (4.7) and (A7).

  5. If v≠0, given Embedded Image in (A2) and definitions of a(p,ϵ), S(p) and Embedded Image, it follows that Embedded Image intersects a(p,ϵ) only at special points (z,pr), such that Embedded Image. Also, Embedded Image in the limit Embedded Image (figure 5a).

One can check that properties (i)–(v) allow for two qualitatively different phase portraits of the projected reduced system (4.7) without a folded singularity, shown in figure 5b,c, one of the eight different portraits with a single-folded singularity shown in figure 5dk, or combinations of these.

Figure 5.

Phase portraits of the projected reduced system system (4.7) satisfying assumptions (A1)–(A6) corresponding to the projections of critical manifold S onto the (z,pr) plane in the reduced system (4.6). Sa,L and Sr denote the attracting part, fold and the repelling part of S, respectively. Initial conditions from Sa that remain in Sa for all time are shaded in grey. a denotes the stable equilibrium of the unramped system (dashed curve), Embedded Image denotes the Embedded Image isocline of desingularized system (4.9) (green curve) and F denotes folded singularities (black dots). Canard trajectories that cross from Sa to Sr through F are shaded in light blue and unique singular canards tangent to the strong stable manifold of the equilibrium for the desingularized system are in dark blue.

An excitable portrait requires an excitability threshold defining a critical value of v. In the two simple cases, either Embedded Image remains within Sa(p) and can intersect a(p,ϵ) if Embedded Image (figure 5b), or Embedded Image lies entirely within Sr(p) (figure 5c). Neither case involves an excitability threshold within Sa(p) because all initial conditions from Sa(p) remain in Sa(p) for all time (figure 5b) or reach L(p) and blow-up (figure 5c). More interesting are transition cases between figure 5b and c, because these involve a folded singularity F(p,v) (figure 5dk). A necessary condition for a folded singularity Embedded Image4.10requires the ramped parameter pr(t) in the fast dz/dt component of the vector field; see appendix A(a) for details. Guided by the grey-shaded regions indicating subsets of Sa(p) that remain in Sa(p) for all time, we identify only two cases that involve an excitability threshold within Sa(p). These are a folded saddle in phase portrait (figure 5d) and a folded saddle-node in figure 5e. Since a folded saddle-node is not structurally stable (Kuznetsov 1995), it will unfold under typical changes in the parameter vector p into either phase portrait in figure 5b or a pair of folded singularities including folded saddle from figure 5d and folded node from figure 5j. This leaves only one phase portrait involving an excitability threshold that is structurally stable.

Hence, under assumptions (A1)–(A6), the ramped system (4.2) has an excitability threshold defining a critical ramping rate, vc, if the desingularized system (4.9) has an equilibrium F(p,v), and Embedded Image4.11meaning that F(p,v) is a saddle; see appendix A(b) for a derivation of equation (4.11). Furthermore, in the limit Embedded Image the excitability threshold is unique. Within Sa(p), this threshold is typically given by the singular canard trajectory through a folded saddle (shown in dark blue in figure 5d). In the three-dimensional (x,z,pr) phase space, the excitability threshold consists of three components: (I) the singular canard trajectory within Sa(p) and all trajectories that converge to it, (II) the part of the fold curve L(p) past the folded saddle F(p) and all trajectories that converge to it, and (III) the subset of Sr(p) dividing initial conditions that end up on the grey-shaded part of Sa(p) from those that blow-up (figure 6a). For sufficiently small but non-zero ϵ and away from L(p), the (normally hyperbolic) attracting Sa(p) and repelling Sr(p) parts of the critical manifold S(p) perturb to nearby invariant (and still normally hyperbolic) slow manifolds Sa,ϵ(p,ϵ) and Sr,ϵ(p,ϵ), respectively (Fenichel 1979; Jones 1995). At the fold curve, however, Sa,ϵ(p,ϵ) and Sr,ϵ(p,ϵ) typically split except for where they intersect along a special trajectory that is called a maximal canard (figure 6b). Specifically, Szmolyan & Wechselberger (2001) prove that a singular canard through a folded saddle for the limiting problem perturbs to a maximal canard for sufficiently small ϵ. As the excitability threshold is no longer unique in analogy to the two-dimensional problem in figure 2, one may want to choose σ in definition 2.2 such that the threshold component (I) in the singularly perturbed system is the maximal canard trajectory and all trajectories that converge to it.

Figure 6.

A sketch of the phase portrait of the ramped system (4.2) near a folded saddle F. (a) Three components of the excitability threshold (I)–(III) in the limit Embedded Image. Perturbations that take trajectories from the shaded region across (I), (II) or (III) will cause an excitable response. The trajectory through F is the singular canard (dark blue). Sa, L and Sr denote the attracting part, fold and repelling part of critical manifold S, respectively (cf. figure 5d). (b) Phase portrait for 0<ϵ≪1. The slow attracting manifold, Sa,ϵ, and the slow repelling manifold, Sr,ϵ, intersect along the maximal canard trajectory (dark blue).

(b) A framework for calculating the critical rate

To calculate the critical rate of ramping, vc, consider a phase portrait as in figure 5d with folded saddle at Embedded Image and assume the system to be at Embedded Image at t=0. Then, the critical rate is the value of v at which the (v-dependent) excitability threshold crosses Embedded Image. In most cases, the threshold can be computed only numerically as the stable invariant manifold Ws(F) of the saddle equilibrium F(p,v) for the desingularized system (4.9). This means that vc too has to be obtained numerically. However, if Embedded Image is sufficiently close to F(p,v), the threshold can be approximated by its linearization at F(p,v), that is a straight line through F(p,v) in the direction of the eigenvector w(p,v)=(w1(p,v),w2(p,v))T corresponding to the negative eigenvalue of saddle equilibrium F(p,v). Then, vc can be calculated from the condition for the threshold line to cross Embedded Image: Embedded Image4.12In some cases, this implicit general condition simplifies to an explicit formula for the critical rate vc. For example, later equation (5.5) gives such a formula for the climate-carbon cycle model (1.1)–(1.3).

5. The compost-bomb instability

In this section, we explain the compost-bomb instability reported in Luke & Cox (in press). Specifically, we derive the critical speed of global warming above which the mechanism of type B excitability causes an abrupt increase in peatland soil temperature that is accompanied by a potentially catastrophic release of peatland carbon into the atmosphere.

To use the general framework developed in §4, we recognize that the climate-carbon cycle model with global warming (1.1)–(1.3) is a singularly perturbed system with two slow variables Embedded Image, one fast variable Embedded Image and singular perturbation parameter ϵ. The slow dynamics of equations (1.1)–(1.3) can be approximated by the reduced flow, obtained by setting ϵ=0, evolving on a two-dimensional critical manifold Embedded Image5.1that has a unique fold transverse to the (slow) C direction Embedded ImageThe ramped system (1.1)–(1.3) has no equilibrium points but the unramped system, obtained by setting v=0 in equations (1.1)–(1.3), has a unique equilibrium at Embedded ImageWe find from linearizing (1.1)–(1.2) that (Ceq,Teq) is asymptotically stable if Embedded Imagethat is for all Embedded Image. In the phase space of the ramped system (1.1)–(1.3), Embedded Imageindicates the position of the unique asymptotically stable equilibrium of the unramped system for different but fixed settings of Ta. It follows that, for the given parameter values, the climate-cycle carbon model (1.1)–(1.3) satisfies requirements (A1)–(A5) for all Embedded Image.

To describe the evolution of the fast variable T on S in slow time t, we set ϵ=0 in equations (1.1)–(1.3) to get the reduced system, differentiate the resulting algebraic equation with respect to t, and rewrite the reduced system as Embedded Image5.2Then, we use the condition for S in equation (5.1) to obtain the projection of the reduced system onto the (T,Ta) plane Embedded Image5.3that is singular at the fold L where α(TTa)−1=0. This singularity is removed by the following scaling Embedded Imagewhich gives the desingularized system Embedded Image5.4In the desingularized system (5.4), the Embedded Image isocline Embedded Imageintersects the Embedded Image isocline L at a point Embedded ImageFor the given parameter values, condition (4.11) is satisfied so F is a saddle equilibrium of the desingularized system (5.4) (and, hence, a folded saddle of the reduced system (5.3)) with eigenvalues Embedded Imageand the eigenvector w=(w1,w2)T corresponding to the negative eigenvalue s is given by Embedded Imagewhere Embedded Image

We now have all the necessary ingredients to extract the critical rate of global warming, vc, from the general condition (4.12). Given an initial condition (C0,T0,T0a)∈Sa sufficiently close to F, using equation (4.12), and as far as Embedded Image, the climate-carbon cycle model (1.1)–(1.3) will exhibit a very large excitable response (the compost-bomb instability) if the rate of global warming exceeds the critical value Embedded Image5.5where Elin accounts for the deviation of the actual stable manifold of F from its linearization at the initial condition (C0,T0,T0a) used in the derivation of equation (4.12), and Eϵ is a correction for non-zero ϵ. In the special but realistic case with the initial condition at the unique stable equilibrium of the unramped system, that is Embedded Image, the critical rate formula (5.5) simplifies to Embedded Image5.6where ‘≈’ arises from neglecting Elin and Eϵ. Furthermore, if 2B2≫1, equation (5.6) simplifies to Embedded Image5.7that is independent of the initial conditions C0 and T0. A comparison in figure 7a between numerically calculated vc (solid curve) and its approximation given by equations (5.6) (dashed line) and (5.7) (dotted line) reveals a very good agreement for sufficiently small ϵ. A small discrepancy between the solid curve and the dashed line is expected even as Embedded Image because of Elin. Given a value of vc, it is interesting to calculate the critical time, tc, defined as the time to reach the temperature of 100 (°C). Results in figure 7b suggest that about 20 years of global warming at a constant rate vc≈0.09 (°C yr−1) may already cause spontaneous combustion of peatlands.

Figure 7.

(a) Critical rate of global warming, vc, as a function of ϵ=μ/A with fixed A=3.9×107 (J kg−1) and varied μ, obtained using (solid curve) numerical solutions to equations (1.1)–(1.3) as well as analytical formulae (dashed line) (5.6) and (dotted line) (5.7) derived for Embedded Image. (b) Critical time to reach 100°C, tc, calculated for v=vc+0.001 (°C yr−1). The initial condition is the stable equilibrium of the unramped system with T0a=0 (°C). The dots indicate ϵ≈0.064 (kg °C−1 m−2) for μ=2.5×106 (J m−2°C−1).

In figure 8, we fix the initial condition (red dot) and explain three qualitatively different responses of the climate-carbon cycle model (1.1)–(1.3) for different rates of global warming v. To avoid numerical problems encountered in Luke & Cox (in press) owing to the model stiffness illustrated in figure 4e, we used Embedded Image rather than Embedded Image. For v=0.075 (°C yr−1)<vc (figure 8a), the initial condition is below the excitability threshold approximated by the singular canard trajectory (blue), so the system does not produce any excitable responses. Rather, the (red) trajectory moves towards a straight away, meaning that the ramped system quickly approaches and then adiabatically follows the stable but changing equilibrium of the unramped system. However, as the rate of global warming is increased to v=0.09 (°C yr−1) >vc (figure 8b), the folded singularity F and the excitability threshold shift their position such that the same initial condition is now above the excitability threshold. As a result, the ramped system reaches the fold of the slow manifold and exhibits an explosive increase in the soil temperature, T, associated with a catastrophic release of soil carbon into the atmosphere. After the excitable spike, the system returns to the attracting part of the slow manifold approximated by Sa, finds itself below the excitability threshold, approaches and then adiabatically follows the stable but changing equilibrium a of the unramped system. Clearly, for large enough v, the system can exhibit more than one excitable spike. In figure 8c for v=0.3 (°C yr−1)>vc, when the system returns to the attracting part of the slow manifold following the first excitable spike, it still finds itself above the excitability threshold and produces the second spike only after which the trajectory approaches a. Note that the phenomenon of multi-pulse excitability was first described by Wieczorek et al. (2002) in connection with n-homoclinic orbits to a saddle focus. The multiple-spike excitable response described here manifests in a similar way but has a rather different underlying dynamical mechanism.

Figure 8.

Spiky excitable responses in the climate-carbon cycle model with global warming (1.1)–(1.3) with ϵ≈0.064 (kg °C−1 m−2) and Embedded Image. The solution starting at the equilibrium of the unramped system (C,T,Ta)=(50,8.15,0) (red dot) is shown as (left column) a trajectory in the (C,T,Ta) phase space and (right column) a response of the soil temperature T in time t. Different rows show (a) sub-threshold response for v= 0.075 (°C yr−1)<vc, (b) single-spike excitable response for v= 0.09 (°C yr−1)>vc and (c) double-spike excitable response for v=0.3 (°C yr−1) >vc. For reference, included are: the two-dimensional critical manifold S=SaLSr (grey surface), the unique fold L of S and the unique asymptotically stable equilibrium a of the unramped system (black curves), the folded saddle F (black dot) and the singular canard trajectory (blue) indicating the excitability threshold within Sa (for Embedded Image). Note the logarithmic T-scale.

6. Conclusions

This work identifies and analyses a novel excitability type in systems with ramping—a steady, slow and monotonic change in one of the parameters, pr, called the ramped parameter. The unramped system considered here is of slow–fast nature and its (unique) asymptotically stable equilibrium exists for any fixed setting of pr. When pr is ramped sufficiently slowly from one setting to another, the system can adiabatically follow the stable but changing equilibrium. However, very large excitable responses appear when pr is ramped above some critical rate of ramping.

We showed that such an excitable system forms a singularly perturbed problem with at least two slow variables and focused on the case with locally folded critical (slow) manifold. Using concepts from singular perturbation theory, we studied the system dynamics in terms of folded singularities and associated canard trajectories corresponding to the intersection of an attracting and repelling slow manifold. In this way, we uncovered possible phase portraits of the ramped system, identified those that give rise to excitability, and gave a necessary and sufficient condition for the existence of a critical ramping rate. Furthermore, we identified a novel type of excitability threshold related to a singular canard trajectory through a folded saddle, and explained the very large amplitude of excitable responses by revealing a novel slow–fast return mechanism that allows for single- or multiple-spike excitable responses. Finally, we derived a general condition for calculating the critical rate of ramping.

The general analysis was motivated by the recently reported ‘compost-bomb instability’ (Luke & Cox in press)—a potentially catastrophic explosive release of peatland soil carbon into the atmosphere as the greenhouse gas carbon dioxide, which could significantly accelerate anthropogenic global warming (Khvorostyanov et al. 2008b). Such apparent discontinuities in the response of the climate system to forcing are commonly termed ‘climate tipping points’. Previous studies of climate tipping points have tended to focus on bifurcations or the levels of equilibrium global warming beyond which each tipping point is likely to be excited (Lenton et al. 2008; Thompson & Sieber in press). In contrast, we have shown here that there is a general class of dynamical systems, including the climate-carbon cycle model (1.1)–(1.3), which define a dangerous rate rather than a dangerous level per se. We suspect that such rate-dependent tipping points are much more common in the climate system than is typically assumed, and suggest that deriving the associated critical rates of global warming, as we have done here for the ‘compost-bomb instability’, would provide valuable guidance for climate change policy.

Acknowledgements

S.W. thanks Mark Holland for useful remarks. P.A. thanks Kiyoyuki Tchizawa for helpful discussions about canards in singular systems.

Appendix A

(a) Folded singularity condition

A point (pr,z)∈S is a folded singularity for system (4.7) iff (pr,z) is an equilibrium for system (4.9) or Embedded Image. By definition of L, (pr,z) is an equilibrium for system (4.9) if (pr,z)∈L. Because fS(z,pr,p)=0 iff (pr,z)∈a by definition of a, and aL=∅ by property (i), we have fS(z,pr,p)≠0 if (pr,z)∈L. Also, Embedded Image is required by assumption (A2). Hence, (pr,z) is an equilibrium for system (4.9) iff (pr,z)∈L and Embedded Image.

(b) Folded saddle condition

Linearizing system (4.9) at F gives a 2×2 Jacobian matrix with entries Embedded Image, Embedded Image, Embedded Image, and Embedded Image. If J11J22J12J21<0, the Jacobian matrix has two real eigenvalues of opposite sign so F is a saddle equilibrium of the desingularized system (4.9) and a folded saddle of the projected reduced system (4.7).

  • Received September 17, 2010.
  • Accepted October 29, 2010.

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

References

View Abstract