Abstract

Phase resetting curves (PRCs) are phenomenological and quantitative tools that tabulate the transient changes in the firing period of endogenous neural oscillators as a result of external stimuli, for example, presynaptic inputs. A brief current perturbation can produce either a delay (positive phase resetting) or an advance (negative phase resetting) of the subsequent spike, depending on the timing of the stimulus. We showed that any planar neural oscillator has two remarkable points, which we called neutral points, where brief current perturbations produce no phase resetting and where the PRC flips its sign. Since there are only two neutral points, all PRCs of planar neural oscillators are bimodal. The degree of bimodality of a PRC, that is, the ratio between the amplitudes of the delay and advance lobes of a PRC, can be smoothly adjusted when the bifurcation scenario leading to stable oscillatory behavior combines a saddle node of invariant circle (SNIC) and an Andronov-Hopf bifurcation (HB).

1. Introduction

Neural oscillators are excitable cells which means that as soon as a parameter, such as an external bias current or an ionic conductance, crosses a threshold, then the system switches from a stable rest state characterized by small dampened excursions of the membrane potential to a high amplitude excursion of the membrane potential, called action potential (AP). Repetitive firing excitable cells were categorized by Hodgkin [1] as Class 1, that is, the firing frequency is continuously tunable by adjusting a bias current and such neurons can fire with arbitrarily low frequencies, or Class 2, that is, discontinuous curve with a nonzero firing frequency threshold.

Excitable cells generate APs when a certain control parameter crosses a threshold; that is, a bifurcation occurs in their state space [2]. It became clear that threshold phenomena, such as the AP, are related to fundamental properties of the mathematical model describing the neuron [3] and its corresponding bifurcation diagram [4]. Class 1 excitability was ascribed to a saddle node of invariant circle (SNIC) bifurcation and Class 2 excitability to an Andronov-Hopf bifurcation (HB) (see [5, 6] and references therein). Since neurons are excitable cells operating near a bifurcation point, it is mathematically possible to significantly reduce the number of variables describing the state of a neuron to a few essential ones by using the so-called canonical forms (see [2] for an extensive review of existing canonical models in neuroscience). For example, the four-dimensional conductance based Hodgkin-Huxley (HH) model [1], which was originally introduced as a model of the giant squid axon, includes mathematical descriptions of voltage-gated inward (depolarizing) sodium and an outward (hyperpolarizing) potassium currents. However, using additional assumptions, for example, “instantaneous” activation/inactivation of some ionic channels, HH model can be reduced to two-dimensional models, such as FitzHugh [7] or Morris and Lecar [8]. In particular, Morris-Lecar model used in this paper relies on an instantaneously activating and voltage-gated inward current and a voltage-gated outward delayed rectifier channel (see the appendix for detailed equations). While dimension reduction results in a low-dimensional and sometimes even analytically tractable model, the qualitative aspects of neural dynamics are preserved in canonical models since near the bifurcation points some dimensions become trivial (see [6, 9] for a review of methods and applications of dimension reduction). The most significant benefit of dimension reduction is a better understanding of the general principles that govern AP generation through phase space analysis (see [4, 9] and references therein).

The characteristics of repetitive firing of neural oscillators, such as the intrinsic period , the firing threshold, or the duration of the AP, are determined by their biophysical details, such as the morphology and the ionic channels expressed. However, there are circumstances, for example, in the study of phase locked modes and synchrony of neural populations, when a phenomenological understanding of neural activity suffices. A phase resetting curve (PRC) is not concerned with the biophysical details of neurons but rather tabulates the transient changes in the firing period due to an external perturbation (usually a brief current or synaptic conductance pulse). The most significant effect of a perturbation occurs during the cycle that contains the perturbation and is quantified by the first order PRC; that is, , where is the intrinsic period of oscillation of isolated neuron and is the transiently modified duration of the current cycle due to a perturbation applied at the stimulus time or phase (see Figure 1(a)). Higher order PRCs are defined in a similar way in order to account for long-lasting effects of presynaptic inputs on the second and subsequent cycles (see [10] and the references therein).

Despite its limitations, the PRC is an invaluable tool in predicting the phase locked and synchronous modes of neural networks. The PRC-based approaches to treat Parkinson’s disease and epilepsy are at the forefront of practical applications of the extensive theory of phase resetting. It is believe that resting tremor is caused by a cluster of neurons located in the thalamus and the basal ganglia which fire synchronously [1113]. By using the PRC method, it is possible to predict the shape and timing of the necessary stimuli that destabilize a synchronous state of a large population of neurons with a single pulse [14], double pulses [15], or bipolar pulses [16]. Deep brain stimulation uses permanently implanted electrodes into a target area of the brain with the purpose of delivering precise current pulses aiming at correcting abnormal synchronization among neurons through careful phase resetting [1518]. Epileptic seizures can also be suppressed by applying brief current stimuli in order to reset the phases of some neurons and break the synchrony of the entire network [1921].

The shape of the PRC plays a crucial role in the way neurons synchronize or desynchronize when they are part of a neural network. It has been shown that neurons characterized by a Type II PRC, which has an almost sinusoidal shape and comparable sizes of the positive and negative lobes, synchronize more readily due to their ability to both slow down (if the oscillator is ahead) and speed up (if the oscillator is behind) [2224]. In contrast, the neurons endowed with a Type I PRC, which have a disproportionate ratio of the two PRC lobes, can predominantly speed up (or only slow down) their rhythm when coupled with other neurons. Figure 1(b) shows an example of a Type I PRC with a disproportionate ratio between the amplitudes of the positive and negative lobes. Traditionally, if the ratio of the amplitudes of the two lobes was less than [25] or other arbitrary number, then the PRC was classified as “unimodal” (see Figure 1(b)). The distinction between Type I and Type II PRCs is more profound than the simple visual appearance and was linked to the bifurcation mechanisms leading to stable oscillations.

It was long advocated that the PRCs are either of Type I (see Figure 1(b)) when the bifurcation mechanism leading to oscillatory behavior is a saddle node of invariant circle (SNIC) or of Type II when the bifurcation is of Andronov-Hopf (HB) type [2224, 26, 27]. Recent theoretical studies suggested that a SNIC bifurcation can also lead to bimodal PRCs [2830].

In this paper we shown that all PRCs are bimodal and the degree of bimodality of a PRC can be smoothly adjusted if the bifurcation scenario leading to stable limit cycle oscillations combines a SNIC and a HB. Furthermore, our current study extends previous results [28, 30] by introducing the concept of neutral points, that is, phases at which external perturbations produce no phase resetting, and showing that a SNIC that ends with a HB leads to a PRC that is a linear combination of both Type I and Type II PRCs.

2. Theory

2.1. Natural Coordinates for Limit Cycle Oscillators

Biological oscillators with a relatively stable periodic activity are mathematically described as limit cycle (LC) oscillators (see [28] and references therein). In order to gain insight into the oscillatory behavior of nonlinear oscillators and their corresponding PRCs, we used a natural reference frame, that is, the tangent and normal to the LC frame, instead of the usual fixed reference frame [29]. In such a natural reference frame, any external perturbation of the LC oscillator determines either a tangent displacement, that is, a sudden phase jump, or a normal (to the LC) perturbation. While the contribution of the tangent perturbations to the PRC is straightforward to account for, the phase resetting induced by normal perturbations is notoriously harder to account for [30, 31] and, according to [32], “the (A/N: canonical) model above may be of limited value since it fails to reflect the relaxation nature of the dynamics.”

Without reducing the generality of our results, here we only focused on planar neural oscillators whose vector flows are , where are the fixed reference frame variables, which in the case of Morris-Lecar model [27, 33] are the membrane potential, , and the potassium conductance, , and the superscript means transposed (see the appendix for model’s equations).

The analysis of the response of a dynamical system to external perturbations is significantly clearer in the moving reference frame of the figurative point traveling along the LC, , than in the fixed frame of the original system of differential equations (see Figure 2(a)). The choice of the unitary matrix that switches between the reference frames, that is, , is not unique [3436]. However, the unit vector of the tangent direction is uniquely defined by and then we chose one of the two possible orientations of the normal vector to be .

In the fixed reference frame, , the temporal evolution of small perturbations, , from is described by the linearized equations of the dynamical system; that is, , where is the Jacobian matrix evaluated along . Using the unitary transform U, the same perturbation with respect to the moving reference frame, , is given by [29, 34] where , + , , and (we also corrected a typo in of [34]). The coefficient measures the effect along the direction for a perturbation applied along the direction with respect to . For example, the coefficient quantifies the effect of a pure normal perturbation measured along the tangent direction; that is, the amount of phase resetting induced along tangent direction to by relaxation processes. In all planar systems , which means that a perturbation tangent to never induces a normal displacement of the figurative point.

While the Jacobian matrix describes the stability of a solution to small perturbations with respect to the fixed reference frame , the matrix (see (1)) determines the stability with respect to normal and tangent perturbations at different phases along . A negative value of any coefficient means that the perturbation decays exponentially along the respective direction, that is, the LC is stable to perturbations applied at that phase.

In the PRC theory, the actual perturbation only affects the membrane potential equation; that is, with respect to the fixed reference frame, (see thick horizontal arrow in Figure 2(a)). The corresponding tangent and normal projections of the perturbation can be found using the unitary matrix : Equation (2) gives the initial conditions with which we integrated the perturbation equations (1). Since for any planar dynamical systems, , the differential equation for normal perturbations, that is, , has the analytic solution: where the initial condition is determined by the injected current of amplitude ; that is, . The differential equation for tangent perturbation, which actually determines the total amount of phase resetting due to the current perturbation , contains both pure tangent perturbations (through coefficient) and normal perturbations (through ); that is, which has the analytic solution: where .

Those fundamental theoretical results highlighted by (3) and (5) relate to the local stability coefficients of the normal and tangent geometrical displacements from the unperturbed limit cycle . In the subsequent section, we will confer these geometric (phase space) displacements into temporal advances or delays of the perturbed figurative point, that is, phase resetting.

3. Results

3.1. Why Are All PRCs Bimodal?

The fact that all PRCs are bimodal results from the existence of two remarkable points along —the neutral points or zero phase resetting points. A perturbation applied at a neutral point induces no phase resetting. For planar systems, the existence of the two neutral points can be explain in terms of the stability coefficients of the LC. Indeed, let us assume for a moment that the crossover coefficient vanishes. This means that the only contribution to the PRC is due to tangent perturbations through . When the applied perturbation is perpendicular to , its corresponding tangent projection vanishes; therefore, its contribution to the PRC also vanishes because the initial tangent displacement in (5) is zero. Since the normal perturbation does not contribute to the tangent displacement (because ), the points where the brief (infinitesimal) voltage perturbation is normal to are the neutral points, that is, the points “” and “” in Figure 2(b). In other words, if , then the maximum and the minimum of the membrane potential are the neutral points of the limit cycle.

If , the exact positions (or phases) of the two neutral points and in Figure 2(b) can be determined from the requirement that the contribution of the pure tangential perturbations to (5), that is, , is exactly balanced by the contribution of the normal relaxation processes through . Furthermore, for planar systems, there are two and only two neutral points along ; that is, the PRC is only bimodal. Indeed, assuming again for a moment that there is no normal-to-tangent coupling between perturbations (), then the only contribution to the PRC is due to . On the other hand, contribution vanishes at local minima or maxima of the membrane potential because the tangent projection of the perturbation, , vanishes. For a planar system, there is only one maximum and only one minimum of the membrane potential; otherwise at least one more independent variable (in addition to the membrane potential and one slow variable ) is required in order to further modulate the membrane potential and produce more local minima and maxima, such as bursting activity. Furthermore, at neutral points the sign of the PRC flips. Indeed, assuming for a moment no normal-to-tangent perturbation coupling (), then the tangent projection of a membrane potential perturbation (see Figure 2(b)) always delays the next spike if applied during the repolarization of the AP around the two neutral points of the limit cycle (see point in Figure 2(b)). On the other hand, outside the repolarization region of the AP, a brief depolarization always advances the next spike. Therefore, the PRC must flip its sign at neutral points.

3.2. Are PRCs Unimodal Near a SNIC Bifurcation?

The theoretical arguments of our conjecture regarding the bimodality of all PRCs are supported by numerical simulations. We found that a brief depolarization applied during the repolarization phase of the AP always delays the next spike of ML model neuron (lengthening the duration of the current cycle), whereas a brief hyperpolarization during the repolarization phase always advances the next spike (see Figure 3). All PRCs in Figure 3 were numerically obtained by integrating the differential equations of the ML model (see the appendix) with a current perturbation in the shape of a rectangular pulse of amplitude (in absolute units, which is 2.45% of the bias current at the SNIC bifurcation point) and a duration of .

Close to the SNIC bifurcation point (see Figure 3(A1)) the positive portion of the PRC is hardly visible in print unless it is magnified 1000 times (see the solid dots in Figure 3(A2)). When the magnified PRC is plotted together with the AP (Figure 3(A2)), we notice that the delay is only induced during a small portion of the repolarization phase that starts shortly after the membrane potential maximum (the first neutral point in the absence of normal-to-tangent coupling) and ends before the membrane potential minimum (the second neutral point in the absence of normal-to-tangent coupling). We also notice from Figure 3 that, very close to the SNIC bifurcation point, small perturbations induce significant phase resettings. Indeed, stimulus amplitude has a very large effect (over maximum phase resetting) close to the SNIC bifurcation point (Figure 3(A1)), whereas for a stimulus amplitude that is only about away from the SNIC bifurcation the amount of phase resetting is reduced 100 times (Figure 3(A2)). The first column of Figure 3 clearly supports the conjecture that Type I PRC, which is characteristic to stable oscillations emerging through a SNIC bifurcation, is always bimodal and its bimodality becomes more apparent while the system moves away from the SNIC bifurcation point. At the same time, the positive portion of the PRC, that is, the delay of the subsequent spike in response to a brief depolarization, always occurs during the repolarization phase (second column in Figure 3), as we conjectured.

The curve shows the characteristic Class I excitability behavior (Figure 4(a)); that is, the firing frequency can be continuously adjusted arbitrarily close to zero [1], when the bias current approaches the SNIC bifurcation, up to about 120 Hz (Figure 4(a)). We used the freely available software package AUTO [37] to obtain the bifurcation diagram (Figure 4(b)) for a fixed value of calcium conductance () and variable bias current. The bifurcation diagram has a branch of hyperpolarized stable steady states for low values of bias current (continuous line), which becomes unstable at the SNIC bifurcation point (see point in Figure 4(b)). The emerging LC at SNIC bifurcation is stable with the peak-to-peak amplitude represented by the pairs of solid circles for each value of the bias current up to the maximum that corresponds to the limit point in Figure 4(b). To further support our conjecture that all PRCs are bimodal, we numerically computed the ratio of the PRC lobes (see Figure 1(b)). Very close to the SNIC bifurcation point the amplitude of the smaller lobe is a few orders of magnitude smaller than the magnitude of the other lobe (Figure 4(c)), which justifies the name “unimodal” for Type I PRC, associated with SNIC bifurcations.

3.3. Why the Bimodality of the PRCs Can Be Smoothly Tuned?

While it is clear that every PRC is bimodal, one question is still open: why does the bimodality of the PRCs change with the control parameters? For example, in Figure 4(c) we notice that the ratio of the amplitudes of the two lobes smoothly increases with the bias current until they are comparable.

How is it possible for the PRC of a system for which the oscillatory behavior emerges through a SNIC bifurcation (and whose PRC was supposed to be “unimodal”) to be so bimodal that we cannot distinguish it from a Type II PRC (which is characteristic for HB)? We believe that the answer must take into account the entire structure of the bifurcation diagram and not only the fact that at some point the stable oscillatory behavior was generated by a SNIC. For example, the complete bifurcation diagram for ML model shown in Figure 4(b) has two bifurcation points, both of which generate stable limit cycle oscillations: on the hyperpolarized stable steady state branch there is indeed a SNIC bifurcation point for low values of the bias current, and on the depolarized branch there is an Andronov-Hopf (HB) point (Figure 4(b)). The oscillatory behavior emerges through a SNIC bifurcation by slowly increasing the bias current along the hyperpolarized branch ( point in Figure 4(b)), which produces a PRC that can be approximated with the “unimodal” Type I expression in the close proximity of the SNIC bifurcation [2, 22, 26]. On the other hand, stable oscillations emerge as soon as the HB point is reached (HB point in Figure 4(b)) by slowly decreasing the bias current along the depolarized branch. As such, the close proximity of the HB produces a PRC that can be approximated with the “bimodal” Type II expression with being a constant phase [2, 22, 26]. We conjectured that no PRC is purely of Type I or Type II but a (linear) combination of the two types. To check this conjecture, we fitted the numerically generated PRCs with a linear combination of the and ; that is, , where and are fitting coefficients. Ideally, we expect that the PRC is entirely captured by a Type I iPRC very close to a SNIC bifurcation; that is, .

By the same token, we expect that the PRC is entirely captured by a Type II iPRC very close to an Andronov-Hopf bifurcation point; that is, . Since the amplitudes and depend on the characteristics of the injected current perturbation (amplitude and duration), we computed their ratio , which is insensitive to such variations (see Figure 4(d)). Our numerical simulations show that indeed the contribution of is one order of magnitude larger than that of near a SNIC bifurcation point whereas the situation is reversed near a HB point (see Figure 4(d)). Furthermore, by setting an arbitrary ratio for classifying the PRCs into “unimodal” or “bimodal”, for example, 0.175 [25] or , we can actually use Figure 4(d) to estimate the range of the bias currents for which such an arbitrary classification is valid.

3.4. How Robust Is the SNIC-HB Scenario That Produces Tunable PRC Bimodality?

The complete bifurcation diagram shown in Figure 4(b) was obtained for fixed calcium conductance (). Since the pair SNIC-HB of bifurcation points is crucial to a tunable bimodality of the PRC, we also explored the reproducibility of such a bifurcation diagram while changing . For example, for large values of calcium conductance, , the bifurcation scenario is always the same: while increasing the bias current from very low values (below 0.08), stable LC oscillations emerge through a SNIC bifurcation (point in Figure 5(a)), but while decreasing the bias current from large values (above 0.3), stable limit cycle oscillations emerge through an Andronov-Hopf bifurcation (HB point in Figure 5(a)). It is this particular bifurcation scenario, with a SNIC bifurcation at one end and an Andronov-Hopf bifurcation at the other end of the bifurcation diagram, that leads to a PRC which is a linear combination of a Type I iPRC (typically near a SNIC bifurcation) and a Type II iPRC (typically near an Andronov-Hopf bifurcation) as shown above. For a critical value of calcium conductance (close to 0.6), both ends of the stable LC oscillations are marked by Andronov-Hopf bifurcation points (Figure 5(b)). The SNIC-HB bifurcation scenario changes to a HB-HB when the two limit points and of the SNIC fold collide and give way to a HB point (see Figure 5(b)). The bifurcation diagrams (Figure 5) also show that large calcium conductances lead to larger peak-to-peak membrane potential oscillations, especially through a significant increase of the depolarization overshoot. This effect is due to the fact the a larger calcium conductance drives the depolarization of the cell more effectively as soon as the membrane potential exceeds the half activation threshold of calcium channels (the parameter of the model described in the appendix) and leads to a rapid and strong overshoot (Figure 5(a)). Large calcium currents, due to a large value of , require a relatively modest depolarizing bias current compensation to sustain oscillations. This is the reason why the bias current range is narrower for larger in the SNIC-HB bifurcation diagrams. On the other hand, for small calcium conductances (Figure 5(b)) the oscillations have small amplitude since calcium channels do not reach their half activation voltage and the dynamics of the system are dominated by potassium channels. Indeed, small hyperpolarizations induced by potassium immediately and dramatically decrease the number of open voltage-gated calcium channels (since the calcium channels already work below their half activation threshold). At the same time, the bias current required to compensate the dynamics described by a HB-HB scenario, which is dominated by the outward potassium current, is significantly larger than the bias required for a SNIC-HB scenario (compare bias current values in Figure 5(a) with those in Figure 5(b)).

4. Discussions

Network models showed that selectively blocking some channels in critical inhibitory neurons can lead to switches from unsynchronized high frequency to synchronized low frequency network oscillations [40]. General anesthetics are selective blockers that could also induce dramatic change in the PRC shape and switch the brain from firing patterns associated with consciousness to states associated with general anesthesia [4143]. Therefore, the conductances of voltage-gated channels, such as calcium in our ML example, could be pharmacologically manipulated to produce a dramatic change in the bimodality of the PRC. Although this paper focuses on understanding the structure of PRCs and the effect of biophysically relevant control parameters, such as the bias current and calcium conductance, at the single cell level, the results are relevant actually at network level. In the end, the PRC is a tool used for predicting neural population behavior [14, 44].

We analyzed the mechanisms leading to a smooth and continuous change of PRC’s bimodality in a planar model. We used a ML model neuron not only because it represents a biophysically realistic model of an actual neuron [8, 27, 33], but also because, by changing only two parameters, that is, the calcium conductance and the bias current, it displays both SNIC and Andronov-Hopf bifurcations. Some authors used the ratio of the positive and the negative amplitudes of the PRCs to classify them as unimodal or bimodal (see [26, 28] and references therein). In this paper, we showed that all PRCs are bimodal regardless of the bifurcation mechanism that led to oscillatory behavior. Moreover, we showed that there are two remarkable points on any planar limit cycle, which we called the neutral or zero resetting points, at which the PRC flips its sign (see Figure 2). Although the phases of the two neutral points slightly change with the parameters of the external perturbation (amplitude and duration of stimulus), their location can be inferred from the model’s equations using the stability coefficients of the LC. The phase at which a brief depolarization can induce a delay (positive resetting) is situated along the repolarization region of the AP (see Figure 3). Since for fast spiking neurons the repolarization covers only a minute fraction of the [30], it may raise concerns regarding real-world applicability of such weak bimodality. Experimental studies showed that a brief pulse can desynchronize a fully synchronized cluster of neural oscillators if the pulse hits the cluster in a very narrow phase range (5% ) [14, 16, 45]. It is also possible for the two neutral points to actually collide, in which case the PRC becomes truly unimodal. We showed that a ML model can mimic Class I excitability (with a continuous frequency versus current curve—see Figure 4(a)) with a SNIC-HB bifurcation diagram (Figure 4(a)). The tunable bimodality of the PRC is a result of the presence of both SNIC and HB, which allows fine mixing of Type I character (specific to a SNIC bifurcation) with Type II character (specific to a HB) into the final PRC (Figure 4(c)).

Oscillatory activity is a hallmark of neuronal network function in various brain regions, including the olfactory bulb [4648], thalamus [49, 50], hippocampus [51, 52], and neocortex [53]. The existence and stability of synchronous and phase locked modes can be predicted using PRCs (see [10, 31, 5456] and references therein). Neural networks oscillate over more than three orders of magnitude and include delta (0.5–3 Hz), theta (3–8 Hz), gamma (30–90 Hz), and ultrafast (90–200 Hz) oscillations [57, 58]. Stable oscillatory activity at single cell or population levels was associated with higher brain functions [59, 60], such as gamma rhythm that is believed to be the basis of temporal encoding [61, 62], sensory binding of features into a coherent perception [59, 63], and storage and recall of information [6466]. Disruption of synchronous rhythms underlies some psychiatric disorders, such as schizophrenia [67, 68]. The ability to smoothly change the PRC from a Type I to a Type II through experimentally available control parameters, such as bias currents or pharmacologically controlled ionic channels, has potential application in shaping the activity of neural cells and networks.

Appendix

We used a dimensionless Morris-Lecar model [27, 33] described by the following equations: where is the membrane potential, is the slow potassium activation, and all ionic currents are described by , where is the conductance of the voltage-gated channel and is the corresponding reversal potential. In particular, the calcium current is , the potassium current is , and the leak current is . The reversal potentials for calcium, potassium, and leak currents are , , and , respectively. The steady state activation function for calcium channels is , where , , the steady state activation function for potassium channels is , where , , the inverse time constant of potassium channels is , the potassium and leak conductances are , , respectively, and the .

The two control parameters that can switch the bifurcation diagram from a SNIC-HB type to a HB-HB type (see Figure 5) are the calcium conductance and the bias current . For example, if and , (A.1) describe what was classified by Hodgkin as Class I excitable cells. For example, if and , (A.1) describe what was classified by Hodgkin as Class II excitable cells.

The two nullclines of (A.1) are two-dimensional curves and , respectively. A fixed point of (A.1) is any crossing of the two nullclines, that is, phase space points where both the rate of change of membrane potential and the potassium activation vanish. The linear stability of the fixed points is determined by the two eigenvalues of the Jacobian matrix of (A.1). Stable steady states have eigenvalues with negative real part and are represented by a solid line in any bifurcation diagram (see Figures 4(b) and 5). Saddle or hyperbolic points have real eigenvalues with opposed signs, whereas Andronov-Hopf points have pairs of purely imaginary eigenvalues.

Acknowledgment

This research was supported by the National Science Foundation Grant CAREER Award IOS 1054914.