Dynamical systems theory in physiology.

Readers of the current issue will find something unfamiliar, perhaps tantalizing, perhaps unsettling. I am referring to the articles by Cha et al. (see “Ionic mechanisms and Ca and “Time-dependent changes in membrane excitability…”). The first of these articles is less exotic; it presents computer simulations of a model for bursting electrical activity in pancreatic β cells. The second uses bifurcation diagrams to analyze the behavior of the model. I will argue that this is relevant far beyond β cells—the leading edge of a wedge driving the methods of dynamical systems theory into the heart of biology.

Mathematical modeling of cell electrical activity has a long history in physiology, going back to the work of Hodgkin and Huxley (Hodgkin and Huxley, 1952; Chay and Keizer, 1983) for action potential generation and propagation in squid giant axon. The model of Cha et al. (2011a,b) is based on the Hodgkin and Huxley formalism, but augmented with mechanisms for maintaining ionic balance (pumps and exchangers for Ca2+, Na+, and K+, and the endoplasmic reticulum). In addition, a nod is given in the direction of metabolism, as β cells are first and foremost metabolic sensors and use ATP-dependent K+ (K(ATP)) channels, to transduce the rate of glucose metabolism into intensity of electrical activity.

Cha et al. (2011a,b) follow the path blazed by Chay and Keizer (1983). Their model was based on the simple idea, first proposed by Atwater et al. (1980), that bursting results from slow modulation of spiking by calcium. That is, during the active spiking phase of the burst, calcium builds up and turns on calcium-activated K+ (K(Ca)) channels until membrane potential falls below the threshold level and spiking terminates. During the ensuing silent phase, calcium would be pumped out of the cell, lowering the spike threshold and allowing the next active phase to begin. Rinzel (1985) formalized this mathematically, recognizing that the key element of Chay–Keizer was a fast spiking system modulated by slow negative feedback. This led to a profusion of models, different biophysically but essentially equivalent mathematically, with alternate proposals for the source of the negative feedback. These included inactivation of the L-type Ca2+ channel (Chay, 1990), indirect activation of K(ATP) channels by Ca2+ via its effects on ATP consumption or production (Keizer and Magnus, 1989), and activation of electrogenic Na–K exchange (Fridlyand et al., 2003). The model of Cha et al. (2011a,b) is most similar to that of Fridlyand and colleagues, including all of the above mechanisms as well as some others, paying close attention to balancing the ion fluxes and including a prominent role for the Na+/K+ exchanger.

Time-based simulations of models are well recognized as useful. They provide a quantitative test of the hypothesis posed by the model, showing whether the hypothesized mechanisms can in fact work together quantitatively to produce the target behavior. They also allow perturbations to be tested that may not be possible in the experimental system. One limitation, however, is that each simulation tests only one set of parameter values. With an analytic solution (a formula for the solution as a function of time), one would know everything there is to know as parameters vary, but this is not possible for complex, nonlinear models like those for β cells.

One solution to this problem is to construct a bifurcation diagram, which is a plot of the solutions versus parameters. This may seem no better than time integrations, as such solutions also have to be obtained by a numerical computation. However, the diagram, in part, is a plot just of the steady states, which is a much easier math problem, requiring only the solution of the system of nonlinear algebraic equations obtained by setting the right-hand sides of the differential equations to 0 and not the time-dependent solution of the differential equations themselves.

An example of such a bifurcation diagram is Fig. 1 in Cha et al. (2011b), computed using the popular public domain program XPPAUT (Ermentrout, 2002). The figure shows the behaviors as the parameter for glucose concentration are varied, namely that as glucose is raised, the system goes from silent to bursting to continuous spiking, with membrane potential and cytosolic calcium, and, implicitly, insulin secretion increasing. Even more important than convenience and efficiency, the diagram shows solutions that are unstable, which cannot be achieved by solving the differential equations in time because the system moves away from those solutions, not toward them. The standard example for illustrating this is a marble on a surface with hills and valleys. The marble will roll away from the hills (the unstable solutions) into the valleys (the stable solutions).

Figure 1.

Linear oscillations. Sinusoidal solutions of (1) for a pendulum with friction (b < 0), negative friction (b > 0), and no friction (b = 0). In the last case, amplitude is constant but depends on the initial velocity of the pendulum. The dashed blue curve is started at the same position but with twice the velocity as the solid black curve.

Unstable solutions are also physically unrealizable by the cells, so why are they of interest? The reason is that oscillations, the behavior of β cells that we really want to understand, are born out of steady states that become unstable. The change of stability is called a bifurcation, which is why such diagrams are called bifurcation diagrams and not just state diagrams. In other words, the diagram is not merely a list of the states the system can be in but shows how they arise from an underlying dynamic process. Indeed, Fig. 1 in Cha et al. (2011b) contains more than the steady states; it also displays the maximum and minimum values attained during some of the oscillations undergone by the system and the bifurcations that lead to them.

The idea of oscillations arising as a parameter changes can be understood with the aid of another simple mechanical example, a pendulum, which can be represented in a simplified form by the linear equationwhere θ represents the angle of the pendulum from a vertical line through the pivot, and b represents friction (valid only for small deflections from the vertical). This equation has only one steady state, θ = 0, and solutions take the formor a similar equation with cosines, or both. (If one substitutes expression 2 into Eq. 1, one can solve for λ and ω, but we don’t need to do this because we are only interested in how the behavior changes qualitatively as b varies.) If b < 0, as would be normal for friction, then λ < 0, and the solutions are decaying oscillations with an exponential envelope (Fig. 1). These go to 0 as time goes to infinity, so 0 is a stable solution; if the system starts at 0 it stays there, and if it is moved off 0 it returns, like the marble rolling into the valley. If, on the other hand, b > 0, which would correspond to some kind of “negative friction” that injects energy into the system in proportion to the speed of the pendulum, then the solutions oscillate with ever increasing amplitude. In this case, we would say that the solution 0 is unstable. (One may wonder where negative friction comes from; we will return to that below.) At the border between these stable and unstable parameter regimes, b = 0, λ would also be 0, and the solution would be a pure sinusoidal oscillation. In a nonlinear system, something similar can happen, and it is called a Hopf bifurcation (HB). Programs like XPPAUT detect this by linearizing the system and calculating constants of the linearized system, called eigenvalues, from which λ and ω can be obtained. When λ changes sign, an oscillation of the full system is born that is initially sinusoidal but distorts because of nonlinear effects as the parameter is moved away from the bifurcation and can become highly asymmetrical, like a neuronal or β-cell action potential.

There is one critical difference between the oscillations of the nonlinear system and those of the pendulum. A pendulum has no fixed amplitude; the amplitude depends on the initial position and velocity of the pendulum. The solid and dashed traces in Fig. 1 for b = 0 have smaller and larger initial velocity, respectively. The amplitude can also be changed by briefly tapping the bob (not depicted). In the nonlinear system, if the oscillation is perturbed, say, by a brief injected current in the case of the β cell, then the oscillation will either return to its original amplitude or run away to some other behavior. Thus, nonlinear oscillations can also be characterized as stable or unstable and are usually referred to as “limit cycles” because the system is either attracted to or repelled by them in the limit as time goes to infinity. The linear oscillation, in contrast, is neutrally stable; like a marble on a flat surface, it does not revert to its old position if it is moved.

Going back to Fig. 1 in Cha et al. (2011b), we find stable oscillations for high glucose concentrations that become unstable through a more complex bifurcation, called a torus bifurcation, as glucose is reduced <18.84 mM. (Note that XPPAUT has calculated the unstable oscillations as well as the stable ones, which again cannot be done by integrating the differential equations in time.) A complex sequence of other bifurcations that are not shown ensues that culminates in the bursting oscillations shown in Fig. 2 of the simulation paper (Cha et al., 2011a). That figure shows that the scenario outlined by the bifurcation diagram is realized in simulations. At 20 mM glucose the system undergoes continuous spiking; at 16, 12, and 8 mM glucose it bursts; and at 6 mM glucose it is quiescent.

Figure 2.

Bifurcation diagrams and time simulations of bursting. Each panel is made with a simplified Chay–Keizer model (Tsaneva-Atanasova et al., 2010). Bifurcation diagrams with overlaid burst trajectories (blue) in upper subpanel. Solid lines, stable solutions; dashed lines, unstable equilibria; dotted line, unstable oscillations. (A) Pure square wave without spikes. (B) β cell–like square-wave bursting. (C) Pituitary-like bursting. (D) Bursting without bistability requires a second slow variable, as in this caricature of bursting in the R15 neuron of Aplysia.

Thus, like a Shakespeare play in which the battle scenes take place offstage, the interesting behaviors of the β-cell model are only hinted at, but those hints are enough for those in the know. In simpler systems the bursting solutions and the bifurcations that lead to them can be calculated with XPPAUT if one has sufficient determination (see Tsaneva-Atanasova et al., 2010). The mathematics of this process, which includes intervals of chaos, is not yet completely understood. Nonetheless, we should pause to admire the achievement; a nontime-dependent calculation of steady states, both equilibria and oscillations, has predicted where to find particular behaviors in the full, time-dependent system. This is the hallmark of a successful theory—it is able to capture the key features of a system with a compressed description, much as Newton’s unitary central, square-law force was able to account for Kepler’s multiple empirical laws.

We have seen so far that bifurcation diagrams can provide useful summaries of a complex system, but Rinzel’s profound innovation (Rinzel, 1985) was to introduce a different kind of diagram that led to insights into how the bursts are generated and modulated, not just where they occur. He formalized the biophysical idea of slow negative feedback by decomposing the system of equations into a set of (at least two) fast equations to represent the dynamics of the voltage-dependent Ca2+ and K+ channels, and a slow equation (cytosolic calcium) to turn the spiking on and off. Because the two subsystems operate on well-separated timescales, they can each be analyzed semi-independently. From the point of view of the slow variables, the fast variables can be viewed as in quasi–steady state, and from the point of view of the fast variables, the slow variable can be viewed as a constant, that is, a parameter. Thus, the program was to construct a bifurcation diagram of the fast subsystem with the slow variable as a parameter. The general idea of subdividing systems by timescales is ubiquitous and highly successful in applied mathematics, but Rinzel was the first to apply it to bursting. The decomposition simplifies the analysis tremendously by reducing the whole burst trajectory to a sequence of lower dimensional problems. The same strategy was implemented by Cha et al. (2011b) in Figs. 2, 4, and 5, but with a much larger number of fast variables and multiple slow variables.

This changes some aspects and introduces new challenges, but the essentials are already present in the simplest case. In fact, we can get most of the way to our goal if we simplify further by eliminating the spikes (see Fig. 2 A). This will allow us to concentrate on how one fast and one slow variable can produce an oscillation via an HB. This mechanism lies at the heart of both the fast spikes and the slow alternation between spiking and silence in the β-cell model. Indeed, the same concept applies to a wide variety of biological oscillations, including neuronal and cardiac action potentials, circadian rhythms, calcium oscillations driven by release from the endoplasmic reticulum, oscillating chemical reactions, the eukaryotic cell cycle, and many more (Keener and Sneyd, 1998; Fall et al., 2002).

Biophysical intuition is adequate to understand action potentials; inward (Na+ or Ca2+) current builds up autocatalytically (i.e., provides positive feedback), and the rise in Vm activates outward current (K+) and/or inactivates inward current (i.e., provides negative feedback). It is much harder to intuit whether one gets one spike or a train, as this depends on a subtle quantitative balance between the depolarizing and hyperpolarizing currents. As already is the case for an action potential, the negative feedback must develop more slowly than the positive feedback or the two would cancel out, but how much slower can make the difference between a single spike and repetitive spiking.

Fig. 2 shows four different oscillation patterns that can be obtained from a simplified version of the Chay-Keizer model (Tsaneva-Atanasova et al., 2010) by varying a few parameters. The model consists of equations for membrane potential (Vm), activation of the delayed rectifier K+ channel, and free cytosolic calcium (Ca). In each panel, the upper subpanel shows the bifurcation diagram of the fast subsystem, and the lower subpanel shows the Vm as a function of time. In addition, the Vm–Ca trajectory is overlaid on the bifurcation diagram to show that the analysis done with cytosolic calcium fixed and viewed as a parameter predicts the time evolution. Moreover, once one is familiar with such diagrams, it is often easy to sketch out the solution using just the geometric information they contain, without needing to solve the equations numerically. In Fig. 2 A, the spikes have been eliminated by speeding up the delayed rectifier K+ channels, leaving only a square wave–like oscillation. Once we understand how this comes about, the other three panels follow suit.

The Z-shaped curve represents the steady-state values of Vm for each value of Ca. Generally, Vm decreases as Ca increases because of the activation of K(Ca) channels. However, for an intermediate range of Ca values, Vm has three steady states, upper and lower states that are stable and an intermediate state that is unstable, indicated by dashing. Bistability, like oscillations, is ubiquitous in biology. An interesting application is the genetic engineering of a toggle switch (Gardner et al., 2000).

In the β cell, bistability results from the balance of inward Ca2+ and outward K+ currents. As Ca decreases from the right, Vm rises gently until the voltage threshold of the L-type Ca2+ channels is reached, which increases Ca2+ entry, which raises Vm further, which opens more Ca2+ channels. This positive feedback loop, which creates the bulge to the right, is the bioelectrical analogue of negative friction. As Vm increases further, voltage-dependent K+ channels open, which eventually balances the Ca2+ current and dominates it as the driving force of K+ efflux increases and the driving force for Ca2+ decreases. Vm continues to rise slightly, but Ca decreases, along the upper branch of the Z-curve. Note that the Z-shape and bistability are not guaranteed by the mere presence of Ca2+ and K+ currents but depend on the parameter values. For example, if the activation curve of the Ca2+ channels is shifted too far to the right, then they will activate at voltages where the K+ channels are already active and lead to a monotonic Ca–V relation. Given bistability, however, the stage is set for the oscillation depicted in the lower panel.

The diagram was calculated on the assumption that Ca is fixed, but in reality, it varies slowly, and Vm adjusts quasi-statically according to the Ca–V relation. Suppose the system starts on the lower branch, with high Ca and low Vm. In this condition, outward pumping of Ca2+ will (again assuming properly chosen parameters, as discussed below) exceed influx through channels and Ca will decrease. Once the left knee of the Z-curve, marked by limit point (LP), is reached, however, the lower branch comes to an end. If the system passes even infinitesimally to the left of LP, the fast dynamics of Vm are unleashed, and the only available equilibrium is the upper one. Once the upper branch of the Z-curve is reached, Vm is once again in equilibrium, but at a higher value of Ca. On this branch, Ca2+ influx exceeds efflux because the Ca2+ channels are wide open. This pulls the system to the right, with Ca increasing and Vm gently decreasing quasi-statically. The denouement is clear. When the system passes the right knee (upper LP), the system again finds itself with a single equilibrium, in this case, the lower one. Once on the lower branch, the cycle repeats. The reader should trace out the time course in the lower subpanel to confirm that every feature corresponds to the events described above in the bifurcation diagram in which time is only present implicitly.

Another shorthand way to describe the oscillation cycle is in terms of thresholds of the fast subsystem, which are generated by the physically unrealizable middle states. These points separate the upper and lower states, not only geometrically but dynamically. If Vm lies below one of these states (again thinking of Ca as fixed), it is repelled downward to the lower steady state. If Vm lies above the middle steady state, it is repelled upward. We can thus think of the lower LP as the place where the threshold falls to reach the rising voltage because of slowly decreasing activation of the K(Ca) channels, and once Vm crosses this threshold, it jumps up rapidly. The same happens in reverse at the end of the plateau, when the falling voltage meets the rising threshold at the upper LP and a rapid repolarization ensues. A valuable innovation in Fig. 4 of Cha et al. (2011b) is to overlay the thresholds coming together on the time course, which may be easier to understand than the Vm–Ca plane and is also applicable when there are many slow variables.

The mere existence of the Z-curve is not a guarantee of oscillations. For example, if the Ca2+ pump is too weak to push Ca down to the level of the LP, the system gets stuck on the lower branch, and if the pump is strong enough to overwhelm the entry of Ca2+ in the high voltage state, the system gets stuck on the upper branch. It is also necessary for Ca to be much slower than Vm. See Rinzel and Ermentrout (1998) for a nice demonstration that in a system with fast activation and slow inhibition, oscillations are assured when certain geometric conditions are satisfied and the slow variable is sufficiently slow.

These are also threshold phenomena, but with a difference—it represents a threshold of the full system, with Ca as a dynamic variable, not of the fast subsystem alone. In their original model, Chay and Keizer (1983) in fact proposed that this was the way that increasing glucose brought the β cells above threshold for bursting, namely by increasing the activity of the plasma membrane Ca2+ ATPase. Later, when the K(ATP) channel was discovered, it was suggested that the threshold was instead set by a sufficient degree of closure of K(ATP) by the increased ATP/ADP ratio. Rinzel et al. (1986) showed that both hypotheses could be framed in terms of the geometry of the Z-curve and whether the full system got stuck on the bottom branch or could escape. With either mechanism, bursting is born via an HB of the full system. This shows the power of mathematical abstraction to reveal the underlying unity of disparate mechanisms. The model of Cha et al. is much more complex than Chay–Keizer, but it too crosses the burst threshold via an HB, indicated by the HB point of their Fig. 1 (Cha et al., 2011b). The fast–slow bifurcation diagrams of Cha et al. (2011b) differ in some details from Fig. 2 B, but the main features are the same, showing the power of simplified models, even wrong ones, to get at the heart of the matter.

“Wait a minute,” you say. “I haven’t seen any bursting yet, just a square wave.” This can be profitably viewed as a degenerate burst without any spikes, as we turn to Fig. 2 B, in which spikes have been restored by slowing down the voltage-dependent K+ channels. We again have a Z-curve of steady states, and the overlaid trajectory exhibits spikes that surround the upper steady state. The new element is that the upper steady state is unstable. Looking to the left, we see that this instability is again the result of an HB but, in this case, of the fast subsystem alone. The HB gives rise to an oscillation, which is the fast spiking oscillation.

The circulation of the system follows the same pattern as in Fig. 2 A, but with spikes adorning the plateau. The end of the active phase of the burst is more complicated with spikes and beyond our scope, but Fig. 2 B is suggestive that the active phase ends when the spike minima intersect the rising threshold. In addition to explaining the bursting time course shown below the diagram, it accounts for the fact that at high glucose, the system undergoes continuous spiking, as we already noted in Fig. 2 of Cha et al. (2011a). Now, getting stuck in the plateau means a state of permanent oscillation with Ca high but nearly constant.

Note that in Cha et al. there are S-curves instead of a Z-curve in some of their bifurcation diagrams because the parameter represents inhibition of an inward current rather than activation of an outward current (for example, conductance of the Ca2+ channel in Fig. 2 in Cha et al., 2011b). The dynamics are nonetheless identical, another example of the power of mathematical abstraction.

In Fig. 2 (A and B), the existence of oscillations is tied to bistability between the upper and lower states, whether the upper state is an equilibrium or an oscillation. This implies immediately that it should be possible to shift the system between the active and silent phases with brief electrical perturbations. In fact, this was shown by Cook et al. (1981) for β cells by injecting positive current in the silent phase or negative current in the active phase or by applying pulses of KCl.

This has been cited as evidence that β-cell oscillations are ionic in origin, rather than metabolic. However, we will see from Fig. 2 C that this is not necessarily the case. The same model is used, but the activation curve of the Ca2+ channels has been shifted to the left, resulting in a somewhat different bifurcation diagram. Again, we have a Z-curve, but now the HB of the fast subsystem is shifted to the right, so that it occurs between the knees. This happens because higher Ca and more K(Ca) channel activation are needed to overcome the enhanced inward current. More important, however, is that the oscillations that emerge from the HB are unstable, indicated by dotting. This means that the spikes in the time course below are not stable oscillations but only transients. If Ca were fixed at any value, Vm would come to rest, not oscillate. Nonetheless, if Ca increases during the active phase faster than Vm reaches equilibrium, small but distinct spikes appear. (The increased speed of Ca is revealed by the wide spacing of the spikes compared with that in Fig. 2 B.)

Models like this turn out to represent well the oscillations seen in pituitary somatotrophs, lactotrophs, and corticotrophs, which are endocrine cousins of the β cell (Tabak et al., 2007; Stern et al., 2008) that burst but have smaller spikes and shorter plateaus than β cells. In some cases, the spikes can continue well to the right of the HB, when paradoxically there is neither a stable steady state nor a stable oscillation. This can happen because the repulsion from the unstable steady state is slow compared with the rate of rise of Ca. (Think of how Road Runner hangs in the air until he realizes that he has run off a cliff.) In this condition, it is very difficult to reset the system from the silent phase to the active phase (Stern et al., 2008). Thus, we have a case in which bursting is ionic but not resettable. (It is also possible to have metabolic oscillations that are resettable by electrical perturbations [Goel and Sherman, 2009], but these are the much slower β-cell oscillations, with a period of ∼5 min, not treated by Cha et al., 2011a,b.) Note that the prediction of the nonresettability of pituitary-type bursters would not be at all apparent from the biophysics or even from a time-based simulation. (This prediction remains to be tested.)

As a final example, consider Fig. 2 D. Here, the activation of the delayed rectifier has been slowed down, resulting in large amplitude spikes. As a consequence, the spiking solutions of the fast subsystem fall below threshold at the left knee of the Z-curve. This means that there is no longer bistability—the system can either spike continuously or be silent, but it cannot alternate between the two, at least not by the mechanism of Fig. 2 (A–C).

It is possible to convert a diagram like the one in Fig. 2 C to one like Fig. 2 D by blocking large-conductance voltage- and calcium-dependent K+ (BK) channels. This effectively makes the summed K+ current through delayed rectifier and BK channels slower by reducing or eliminating the faster component. Van Goor et al. (2001) showed that this could account for the paradoxical shortening of plateaus by blocking BK channels in pituitary somatrophs and also suggested that pituitary gonadotrophs do not show plateau bursting because they lack the BK channels found in other pituitary cells.

If a second slow variable is added, however, it is possible to get another kind of bursting. In fact, this is the mechanism proposed by Rinzel and Lee (1987) for the bursting in the R15 neuron of the sea snail Aplysia. The fact that there is no bistability means that these bursts would not be resettable, either upward or downward, by brief perturbations, and the fact that there are two slow variables, one providing positive feedback and one providing negative feedback to the fast subsystem, means that the slow variables can oscillate autonomously, even if the spikes are blocked. Thus, the analysis of Rinzel and Lee (1987) accounted for two key experimental observations.

The above examples show some of the successes of bifurcation theory and dynamical systems approaches more generally in solving biological puzzles. They provide insights that are not possible from a biophysical or simulation approach. Beyond that, Fig. 2 hints at a deeper level of theory than the study of particular bursting systems. All of the examples we have considered arise from a common substrate with modest changes in parameters. With the first few particular examples in hand, Rinzel (1987) showed that they could be classified based on the bifurcations that mediated the initiation of the active phase and the termination of the active phase rather than the channels or other mechanisms involved. Izhikevich (2010) showed that many more types of bursters could be anticipated if one considered all the bifurcations that can lead to or abolish oscillations, many of which have now turned up in modeling studies. Indeed, the pattern shown in Fig. 2 C was predicted before it was instantiated by a model. Deeper mathematical analysis shows further that all of these bifurcations can emerge from a single fundamental bifurcation that spins off multiple configurations of subsidiary bifurcations as several parameters are varied (Bertram et al., 1995; Golubitsky et al., 2001). Although it can be very difficult to establish definitively that a particular model or proposed bifurcation scenario is implicated in a given biological phenomenon, we can be confident that as cells, tissues, and organisms randomly explore parameter space, they will encounter these bifurcations and make use of them. With advances like this, mathematical biology leans in the direction of realizing its aspirations to achieve the kind of theoretical unification found in more mature traditions like physics.

We have seen that exploiting multiple timescales is a powerful tool for dissecting the behavior of complex systems by subdividing the variables. Multiple timescales also have biological significance; the semi-independence of the subsystems makes them modular, which is beneficial from an evolutionary point of view because it allows innovation without destroying well-tested subcomponents (Tyson et al., 2003).

Why are fast–slow oscillations ubiquitous in biology? We can think of this as homeostasis of a higher order. Biology is replete with negative feedback. Almost every reaction is inhibited directly or indirectly by its product, which prevents concentrations from running out of narrow bounds. However, if the feedback is slow and some source of positive feedback is available, the system can undergo transients before returning to rest, and under the right conditions, this can result in repeated oscillations. When these oscillations are useful, they can be fixed by evolution.

None of the ideas presented here are new. They are old hat to mathematical biologists although little known to nonmathematical biologists. Physiology, especially electrophysiology, has had a long symbiotic relationship with dynamic modeling because of the early development of techniques for monitoring time-dependent behavior with high time resolution. There was not much of a field of calcium modeling before the invention of imaging techniques, which revealed a wealth of dynamic phenomena such as oscillations and waves. When experimentalists turned to theorists for help in understanding these phenomena, a large repertoire of models was ready-to-hand to help out that was further enriched by new examples and by the challenge of integrating the calcium and electrical subsystems in cells.

As biology forges ahead and live cell–imaging techniques reveal the temporal complexity of more and more cell-signaling mechanisms, dynamical systems theory will be essential. Luckily, there are many accessible sources aimed at bringing this theory within the grasp of experimentalists (Keener and Sneyd, 1998; Rinzel and Ermentrout, 1998; Fall et al., 2002; Tyson et al., 2003; Izhikevich, 2010) and many theorists with the deep grounding in biology needed to lend a hand. Even if your own work has yet to be touched by these developments, be on the lookout: they are coming to a biological system near you.