Extreme phase sensitivity in systems with fractal isochrons

# Extreme phase sensitivity in systems with fractal isochrons

A. Mauroy Department of Electrical Engineering and Computer Science, University of Liège, 4000 Liège, Belgium    I. Mezić Department of Mechanical Engineering, University of California Santa Barbara, Santa Barbara, CA 93106, USA
###### Abstract

Sensitivity to initial conditions is usually associated with chaotic dynamics and strange attractors. However, even systems with (quasi)periodic dynamics can exhibit it. In this context we report on the fractal properties of the isochrons of some continuous-time asymptotically periodic systems. We define a global measure of phase sensitivity that we call the phase sensitivity coefficient and show that it is an invariant of the system related to the capacity dimension of the isochrons. Similar results are also obtained with discrete-time systems. As an illustration of the framework, we compute the phase sensitivity coefficient for popular models of bursting neurons, suggesting that some elliptic bursting neurons are characterized by isochrons of high fractal dimensions and exhibit a very sensitive (unreliable) phase response.

isochrons, fractals, transient chaos, bursting neurons

## I Introduction

Isochrons and asymptotic phase play a central role for the study of asymptotically periodic systems. The isochrons have been introduced in Winfree2 () as the sets of initial states that converge to the same trajectory on the limit cycle. Equivalently, they are the sets of states that share the same asymptotic phase Winfree (). The notions of isochrons and asymptotic phase are of paramount importance to capture the system sensitivity to external perturbations Taylor (); Sacre3 (). Indeed, an external perturbation has a significant permanent effect on the system only if the perturbed and unperturbed trajectories lie on different isochrons (associated with a noticeable asymptotic phase difference). Beyond their applications to sensitivity analysis, isochrons and asymptotic phase also lead to a powerful phase reduction of the dynamics that is widely used for studying synchronization properties of coupled limit-cycle systems (see e.g. Brown (); Ermentrout_book (); Hoppensteadt (); Kuramoto_book ()).

Thanks to a numerical method based on the so-called Koopman operator framework and introduced in Mauroy_Mezic (), it has been observed recently that isochrons may exhibit a fractal geometry MRMM_bursting (). (Note that a similar phenomenon was discussed in MezicBana () for the asymptotic phase of a discrete-time map.) The fractal property of the isochrons is naturally explained in presence of transient chaos Yorke_chaotic_transient (), a regime characterized by sensitivity to initial conditions caused by a nonattracting (chaotic) set chaotic_dyn_book (). Also, this property is related to an extremely high phase sensitivity of the system, which has a dramatic effect on the response to external inputs and perturbations.

The main goal of this paper is to propose a theoretical framework—and associated numerical methods—that complements the brief and empirical observations given in MRMM_bursting () on the fractal properties of the isochrons. To that end, we define the phase sensitivity coefficient that quantifies the overall uncertainty of the asymptotic phase under small perturbations. Using a result of McDonald (), we prove that this coefficient is closely related to the capacity dimension of the isochrons Farmer_fractal_dim (). A consequence of our results is that, when the isochrons are fractal, a significant decrease of the intensity of a noise perturbation can only slightly reduce the average uncertainty on the phase.

In contrast to local measures of sensitivity (e.g. finite-time Lyapunov exponents), the phase sensitivity coefficient defined in the present paper is an invariant of the system that captures a global property. An important motivation of our approach is therefore to provide a framework for comparing the overall phase sensitivity of different asymptotically periodic systems. As an illustration, the framework is applied to several types of bursting neurons. The results show that elliptic bursting neurons tested are highly sensitive (i.e. with fractal isochrons of high dimension) and therefore characterized by unreliable responses to external inputs. Their (finite) phase response curve is also shown to be fractal.

The paper is organized as follows. In Section II, we rigorously define the notions of asymptotic phase and isochrons. Section III presents basic observations and descriptions of the fractal properties of the isochrons. The relationship between phase sensitivity and fractal dimension of the isochrons is discussed in Section IV. We apply the results to bursting neuron models in Section V. Finally, concluding remarks are given in Section VI.

## Ii Asymptotic phase and isochrons

In this section, we introduce the concepts of phase function and isochrons, for both continuous-time and discrete-time systems. A numerical method is also presented for the computation of the phase function.

### ii.1 Continuous time

We consider a nonlinear system , with and analytic, which generates a flow and admits a periodic orbit of period , i.e. with . Each point of the periodic orbit is associated with a phase according to the mapping , where is an arbitrarily chosen point of Winfree ().

If is an asymptotically stable limit cycle with a basin of attraction , the (asymptotic) phase function assigns the same phase to the initial states converging to the same trajectory on the limit cycle, i.e.

 Θ(x)=θ⇔limt→∞∥φ(t,x)−φ(t,xγ(θ))∥=0. (1)

The level sets of —i.e. the sets of states that share the same asymptotic behavior—are the so-called isochrons Winfree2 ()

 Iθ={x∈B|Θ(x)=θ}.

If the limit cycle is normally hyperbolic, the isochrons are co-dimension- manifolds that invariantly foliate the basin of attraction Hirsch (). Figure 1(a) shows the asymptotic phase and isochrons for the Van der Pol model.

Since the phase function is defined only for initial conditions of trajectories converging to the limit cycle, it is not defined outside the basin of attraction . In particular, the boundary of is called the phaseless set and satisfies the following property Guckenheimer_iso (): the values of the phase function evaluated on any neighborhood of span the entire circle (equivalently, the isochrons come arbitrarily close to ). The phaseless set usually corresponds to a (unstable) fixed point and its stable manifold. In Figure 1(a), the phaseless set is the unstable fixed point at the origin.

Primarily developed in the context of computational neuroscience, an important tool using the notion of phase is the well-known (finite) phase response curve

 Ze(θ)=Θ(xγ(θ)+e)−θ, (2)

with and . It represents the phase shift of a trajectory on the limit cycle that is subjected to an impulsive perturbation , where is the Dirac function.

### ii.2 Discrete time

The phase function can also be defined in the case of discrete-time systems. Consider the map , , , which generates the flow . We assume that the system admits an invariant set , a topological circle, on which the dynamics has an irrational rotation number . This invariant set can be characterized as a closure of a dense orbit, i.e. for some . In this case, there exists a sequence such that , with , and for all .

###### Remark 1 (Proof of the existence of the sequence tk).

We have by definition

 ω0=limt→∞Ft(s(xγ))−s(xγ)t,

where defines a coordinate on the circle and is the lifting of the map on the circle to real line. Let the sequence converge to and define

 ωk=Ftk(s(xγ))−s(xγ)tk, (3)

where or (we know such exist due to the fact that every trajectory in the circle is dense). Therefore, it follows from (3) that or and thus .

The phase is defined as follows. Each point is associated with a phase according to the mapping , where is an arbitrarily chosen point of and where the sequence satisfies . If the invariant set is an attractor with a basin of attraction (i.e. is the smallest set such that for all ), then the phase function is defined on by (1). In addition, the isochrons can be defined as the level sets of the phase function, but they might not be connected manifolds. The phase function and 2 isochrons of a simple discrete-time map are shown in Figure 1(b).

### ii.3 Numerical computation

There exist several methods for computing the isochrons of limit cycles (see e.g. Guillamon (); Huguet (); Izi_book (); Langfield (); Osinga (); Sherwood ()). For our purpose, we will use the forward-integration method proposed recently in Mauroy_Mezic (), which is based on the fact that the phase function is related to an eigenfunction of the so-called Koopman operator MezicBana (). (Note that the method presented in Guillamon (); Huguet () precisely solves a partial differential equation which is closely related to the eigenvalue equation for the Koopman operator.) The method is appropriate to deal with the complex geometry of the isochrons that we investigate in this paper. It is also well-suited to the computation of isochrons in non-planar models. Through this framework, the phase function of continuous-time systems is directly given by the argument of the Fourier average evaluated along the trajectories, i.e.

 Θ(x)=∠(limT→∞1T∫T0g∘φ(t,x)e−iω0tdt), (4)

where is an arbitrary function (observable) such that the first Fourier coefficient (first harmonic) of the periodic function (with ) is nonzero. Note that the state associated with the phase is determined by the specific choice of . Similarly, the phase function of a discrete-time map is given by

 Θ(x)=∠(limT→∞1TT∑t=0g∘φ(t,x)e−iω0t). (5)

The Fourier averages (4) and (5) can be easily computed through the numerical integration of trajectories, with the initial conditions on a uniform grid that spans a region of interest in the state space. The isochrons are obtained by plotting the level sets of these Fourier averages. Further numerical details can be found in Appendix C

## Iii Fractal properties of the phase function

We consider particular (continuous-time and discrete-time) dynamical systems and show that their associated phase function exhibits fractal patterns, suggesting that the isochrons are fractal. These systems are characterized by a very high phase sensitivity.

### iii.1 Continuous-time model

The Lorenz system

 ˙x=σ(y−x)˙y=x(r−z)−y˙z=xy−bz (6)

admits a stable limit cycle for the parameters , , and . It also exhibits transient chaos Yorke_chaotic_transient () and has a chaotic saddle (see e.g. chaotic_dyn_book (), also called fractal repeller Gaspard_fractal_repeller ()), i.e. a fractal invariant set containing the union of all unstable periodic orbits. The chaotic saddle and its stable manifold correspond to the phaseless set , since trajectories starting on this set do not converge toward the limit cycle. The chaotic saddle is fractal, and therefore is also characterized by fractal properties. In addition, includes the stable manifold of the saddle node at the origin, which was shown in McDonald () to have a fractal Cantor-like geometry.

For the Lorenz system (and in general, for asymptotically periodic systems that exhibit transient chaos), the particular fractal properties of the phaseless set have a significant effect on the phase function, which exhibits unusually complex patterns (Figure 2(a)-(c)). While it is well-known that the phase function and the isochrons may be complicated near the phaseless set (e.g. near an unstable fixed point, see Langfield (); Osinga ()), the remarkable fact relies here in their fractal properties, which is induced by the fractal geometry of the phaseless set itself (see the close-up in Figure 2(d)). Note that, for the sake of clarity, we show only the phase function in Figure 2, the fractal isochrons being mainly concentrated near the phaseless set .

In presence of transient chaos, the neighborhood of is characterized by a (extremely) high phase sensitivity. Figure 3 shows that two trajectories starting from this region, with very close initial conditions, can have different behaviors reflected in their asymptotic phase. During a time of the order of several limit cycle periods, the two trajectories remain close to the fractal phaseless set , then “escape” and diverge near the axis, subsequently reaching different regions of the limit cycle. Note that it is not sensitivity to initial conditions in classical sense, where exponential divergence is forever (i.e. positive Lyapunov exponent). Instead, it is the popular notion of phase sensitivity where small changes in initial conditions can separate the trajectories so that they are characterized by different asymptotic behaviors on the limit cycle (i.e. different phases).

As shown in Figure 4, this phenomenon can be captured through the computation of the largest finite-time (or local) Lyapunov exponent local_lyap_expo2 (); local_lyap_expo1 (). For given initial condition and time horizon , the finite-time Lyapunov exponents are given by the logarithm of the eigenvalues of the matrix

 Λ=(MT(T)M(T))1/(2T)

where denote the transpose of . For continuous-time systems, is the fundamental matrix solution of

 dMdt=J(φ(t,x))M

with and is the Jacobian matrix of the vector field . For discrete-time maps, we have

 M(T)=T−1∏t=0J(φ(t,x)).

Regions of high finite-time Lyapunov exponent (black regions) are associated with a high sensitivity to initial conditions. By comparing with Figure 2(c), we verify that these regions lie close to the fractal phaseless set.

### iii.2 Discrete-time model

We consider the following map (taken from MezicBana ()):

 x(t+1)=(1−γ)x(t)+asin2(2πy(t))y(t+1)=x(t)+y(t)+asin(2πy(t))mod1 (7)

with and . The map has an attractor (a topological circle on which the dynamics has a rotation number ) near . The corresponding asymptotic phase function is characterized by a complex geometry (see Figure 5(a) or Figure 12(a) in MezicBana ()) with self-similar patterns (Figure 5(b)). In contrast to the continuous-time Lorenz system, these patterns are not observed at very small scales, but have a minimal size that depends on the distance to (Figure 6). (Note that the patterns have an arbitrarily small size as they are observed far from , i.e. for .) The region of self-similar patterns is therefore an almost phaseless set111The term “almost phaseless set” has been coined in Osinga () rather than a phaseless set, i.e. a region of very large—but not arbitrarily large—phase variation. It is characterized by a high sensitivity to initial conditions, as shown by the values of the largest finite-time Lyapunov exponent (Figure 7). It is important to note that there is no precise definition of the almost phaseless set, which depends on the threshold chosen to assess whether the phase variation is large or not. Moreover, the almost phaseless set does not correspond to a specific invariant set known a priori. In contrast to the phaseless set in the continuous-time situation, it cannot be interpreted in terms of a fractal chaotic saddle, since all trajectories—including those with the initial conditions in —converge to the attractor .

## Iv Phase sensitivity and fractal dimension

Phaseless sets are always characterized by a high sensitivity of the asymptotic phase. But when they have a fractal geometry, they occupy an important portion of the state space, so that the overall phase sensitivity of the system is accentuated. The largest Lyapunov exponent cannot capture this overall phase sensitivity, since it is always equal to zero (the finite-time Lyapunov exponent tends to zero as the time horizon increases). Also, although the computation of finite-time Lyapunov exponents can be used to highlight sensitive regions of the state space (see Figures 4 and 7), it only provides a local measure of the system sensitivity which depends on the chosen time horizon and initial condition. Therefore, a new notion is required to measure the overall phase sensitivity of the system. To this end, we define a phase sensitivity coefficient and we show that this coefficient is closely related to the fractal dimension of the isochrons. This implies that the coefficient is an invariant of the system and can be used to compare the phase sensitivity of different systems.

### iv.1 Phase sensitivity

Consider the geodesic distance on the circle

 d(θ,θ′)=mink∈Z|θ−θ′+k2π|. (8)

We define the phase sensitivity function by

 f(x,ϵ)=maxx′∈B(x,ϵ)∩Bd(Θ(x),Θ(x′)), (9)

where is the basin of attraction of the limit cycle and is a ball with center and radius . If there is an uncertainty (e.g. induced by external perturbations or noise) on the initial condition of a trajectory, then the asymptotic behavior of the trajectory will be associated with an uncertainty on the phase. Note that is the worst-case uncertainty.

For a given compact set , the average phase sensitivity function is computed as

 ⟨f(x,ϵ)⟩A∩B=1μ[A∩B]∫A∩Bf(x,ϵ)dx,

where is the Lebesgue measure on . If (unless ), we define

 α=limϵ→0ln⟨f(x,ϵ)⟩A∩Blnϵβ=1−limϵ→0ln⟨f(x,ϵ)⟩A∩Blnϵ, (10)

or equivalently, for . We refer to as the phase sensitivity coefficient. If , reducing the uncertainty on the initial condition by a certain amount (e.g. reducing the noise intensity) reduces the average phase uncertainty by the same amount (at least when is small). This is the usual situation observed with globally asymptotically stable periodic systems (e.g. Van der Pol model, see Figure 1). But if the phase function has fractal properties, we observe that . In this case, a reduction of the uncertainty produces only a slight reduction of the average phase uncertainty. For instance, in the Lorenz model considered in Section III.1, a logarithmic plot of the average phase sensitivity function with respect to shows that the phase sensitivity coefficient is equal to (see Figure 8). In this case, a reduction of the uncertainty by a factor only reduces the average phase uncertainty by a factor . The phase sensitivity coefficient is therefore directly related to the overall phase sensitivity of the system.

### iv.2 Fractal dimension

We now show that the phase sensitivity coefficient (10) is an invariant of the system that is closely related to the fractal co-dimension of the phaseless set and the isochrons. For a given , consider the set

 MS(ϵ)={x∈A∩B|B(x,ϵ)∩S≠∅},

i.e. the sets of points lying within a distance of the phaseless set . Since the values of the phase function span on any neighborhood of , we have

 f(x,ϵ)=π∀x∈MS(ϵ). (11)

If the limit cycle is normally hyperbolic, the isochrons are as smooth as the vector field, and so is the phase function in , since smoothly increases from one isochron to another. It follows that, for small enough and , we have

 f(x,ϵ)=max∥e∥=ϵ|∇Θ(x)⋅e|+O(ϵ2)=ϵ∥∇Θ(x)∥+O(ϵ2)∀x∈A∩B∖MS(ϵ+δϵ), (12)

where is the gradient of the phase function and is the Euclidean norm. Then (11) and (12) imply that

 ⟨f(x,ϵ)⟩A∩B=1μ[A∩B](∫MS(ϵ)f(x,ϵ)dx+∫A∩B∖MS(ϵ+δϵ)f(x,ϵ)dx+∫MS(ϵ+δϵ)∖MS(ϵ)f(x,ϵ)dx)=μ[MS(ϵ)]μ[A∩B]π+(1−μ[MS(ϵ+δϵ)]μ[A∩B])ϵ⟨∥∇Θ∥⟩A∩B∖MS(ϵ+δϵ)+1μ[A∩B]∫MS(ϵ+δϵ)∖MS(ϵ)f(x,ϵ)dx+O(ϵ2).

Since by definition , one has

 0≤∫MS(ϵ+δϵ)∖MS(ϵ)f(x,ϵ)dx≤μ[MS(ϵ+δϵ)∖MS(ϵ)]π

and it follows that

 μ[MS(ϵ)]μ[A∩B]π+(1−μ[MS(ϵ+δϵ)]μ[A∩B])ϵ⟨∥∇Θ∥⟩A∩B∖MS(ϵ+δϵ)+O(ϵ2)≤⟨f(x,ϵ)⟩A∩B≤μ[MS(ϵ+δϵ)]μ[A∩B]π+(1−μ[MS(ϵ+δϵ)]μ[A∩B])ϵ⟨∥∇Θ∥⟩A∩B∖MS(ϵ+δϵ)+O(ϵ2).

Considering and taking the limit yield

 limϵ→0ln(μ[MS(ϵ)]+(1−μ[MS(2ϵ)]/μ[A∩B])ϵ⟨∥∇Θ∥⟩A∩B∖MS(2ϵ))lnϵ≥α≥limϵ→0ln(μ[MS(2ϵ)]+(1−μ[MS(2ϵ)]/μ[A∩B])ϵ⟨∥∇Θ∥⟩A∩B∖MS(2ϵ))lnϵ, (13)

where we used and the fact that for . If (i.e. ), we have simply

 α=limϵ→0ln(ϵ⟨∥∇Θ∥⟩A∩B)lnϵ=1β=0.

Otherwise, according to McDonald (), scales as

 μ[MS(ϵ)]∼ϵN−D, (14)

with the dimension of the state space and the capacity dimension of the fractal set , i.e. (see Farmer_fractal_dim ())

 D=limd→0lnN(d)ln(1/d) (15)

where is the number of boxes of size required to cover . In addition, if is a normally hyperbolic repeller of co-dimension (at most) one, the gradient scales as when it is evaluated at a small distance of (see Appendix A). It follows that its integral scales as when the boundary of the domain of integration is at a distance of , so that for 222If is of co-dimension greater than one, intuitive arguments suggest that . However, a rigorous result is beyond the scope of this paper.. Finally, (13) leads to

 limϵ→0ln(ϵN−D+(1−(2ϵ)N−D)ϵlnϵ)lnϵ≥α≥limϵ→0ln((2ϵ)N−D+(1−(2ϵ)N−D)ϵlnϵ)lnϵ

or equivalently

 α=N−Dβ=1−(N−D) (16)

(for of co-dimension ). This important relationship implies that the phase sensitivity coefficient is directly related to the fractal dimension of . More precisely, the coefficient will be strictly greater than only if is fractal. (Note that a phaseless set of co-dimension greater than is associated with a phase sensitivity coefficient equal to , even if it is fractal.) This result also implies that the phase sensitivity coefficient is an intrinsic property of the system. It has a unique value that does not depend on the choice of the set .

The fractal properties of and the phase sensitivity can also be characterized by considering the set

 Mδθ(ϵ)={x∈A∩B|∃x′∈B(x,ϵ) s.t. d(Θ(x),Θ(x′))>δθ}

for given and . Since it is clear that as , (14) yields the additional relationship

 N−D=limϵ→0lnμ[Mδθ(ϵ)]lnϵ (17)

for any value , provided that . If is fractal with , then the fraction of trajectories for which an uncertainty smaller than on the initial condition induces a phase uncertainty greater than on the asymptotic phase will only be slightly reduced by a significant decreasing of . We remark that this property is independent of the value .

The definitions and equalities (10)-(16)-(17) summarize the relationships between the overall phase sensitivity of the system and the fractal dimension of the isochrons and phaseless set. Figure 9 shows that (10) and (16)-(17) yield similar results.

#### The isochrons are fractal.

When the phaseless set is fractal, the isochrons themselves are also fractal. Since any neighborhood of intersects an isochron Guckenheimer_iso (), every box of size used to cover intersects . Hence, we have

 NI(d)≥NS(d),

where and are the number of boxes of size required to cover and , respectively. It follows that, for ,

 NI(d)log(1/d)≥NS(d)log(1/d)

and taking the limit , we obtain from (15) that the capacity dimension of is greater or equal to the capacity dimension of .

#### The phase response curve is fractal.

The (finite) phase response curve (2) is fractal provided that the curve satisfies . The fractal dimension is obtained as follows. We assume that the isochrons—i.e. the level sets of —have a fractal dimensional equal to (see (16)). Then the hypersurface is of dimension and it follows that the curve is of dimension . This implies that the dimension of the phase response curve (2) is also equal to . An example of fractal phase response curve will be given in Section V (Figure 13).

#### Almost phaseless set.

The discrete map (7) has an almost phaseless set characterized by self-similar patterns scaling down to , at best. This implies that the phase sensitivity coefficient, as defined in (10), is equal to zero, so that the dimension of is one. But for large values (i.e. ), the phase sensitivity function behaves as if the fractal dimension were more than one: the average phase sensitivity function remains close to one and slowly decreases as decreases (see Figure 10). In other words, the “infinitesimal” phase sensitivity coefficient

 1−dd(lnϵ)ln⟨f(x,ϵ)⟩A∩B=1−ϵddϵln⟨f(x,ϵ)⟩A∩B (18)

where denotes the derivative with respect to , is large for large . In contrast, for maps that do not have an almost phaseless set with these self-similarity properties, the average phase sensitivity function scales as for all values of (black stars in Figure 10).

### iv.3 Numerical computation

For efficient numerical computations of the phase sensitivity coefficient, it is necessary to reduce the number of evaluations of the phase function, which can be computationally expensive. Toward this aim, we propose the following guidelines.

#### Choice of the set A.

In order to capture the overall phase sensitivity of the system, it seems natural to consider a set that contains a large region of the state space. However, the phase sensitivity coefficient is unique and does not depend on the size and location of the set , provided that this set has a non-empty intersection with the phaseless set , as required by our definition. This implies that accurate results are obtained with small sets that cover a small portion of . From a practical point of view, regions of interest containing can be located through the use of finite-time Lyapunov exponents (see Figure 4) or simply by detecting regions of high variation of the phase function (see Figure 2). In addition, in the generic case where the phaseless set is of co-dimension (or less), it is more efficient to choose a set of dimension . Such a set has an intersection with so that for all . It follows that the results of Section IV.2 are still valid (with the measure being the one-dimensional Lebesgue measure).

#### Approximation of the ball B.

It is efficient to approximate the ball by considering a few sample points on its boundary. In particular, two points on opposite sides of are enough to obtain the exact value of the phase sensitivity coefficient. In this case, the phase sensitivity function is approximated by

 ~f(x,ϵ)=maxx′∈{xk−ϵe,xk+ϵe}d(Θ(xk),Θ(x′))≤f(x,ϵ)

where is a unit direction. For , it is clear that we have instead of (11). However, in the generic case where the phaseless set is of co-dimension (or less), the two points and on the boundary of lie on both sides of , so that

 limϵ→0~f(x,ϵ)≠0∀x∈MS(ϵ). (19)

 ~f(x,ϵ)=ϵ|∇Θ(x)⋅e|+O(ϵ2)∀x∈A∩B∖MS(ϵ+δϵ) (20)

instead of (12). Using (19) and (20) in the computations of Section IV.2, it is easy to see that the result still holds when is replaced by . This has been confirmed by numerical simulations. Note also that the value of the phase sensitivity coefficient does not depend on the location of the two opposite points on the ball.

#### Method based on the sets Mδθ.

As shown in Figure 9, the phase sensitivity coefficient can also be computed with the ratio (17). However, the value in (17) is usually underestimated when the ball is approximated with a few points (see above) so that the numerical computation may yield a zero value for small (but positive) values . The computation of the phase sensitivity coefficient is therefore easier and more accurate when using the ratio (10).

## V Application to bursting neurons

The notions of phase and isochrons play a central role in the study of neuron models, showing the sensitivity of neurons to external inputs. Motivated by preliminary observations presented in MRMM_bursting (), we illustrate the framework based on the phase sensitivity coefficient on popular (periodic) bursting neuron models (see Appendix B). Our comparison of different models shows that some elliptic bursting models exhibit strong fractal properties associated with very high phase sensitivity.

Bursting is the alternation between a relatively quiescent state and a succession of rapid spikes in a system. This phenomenon is observed with slow-fast dynamics and is explained by particular bifurcations in the fast subsystems Rinzel_Lee (). According to the types of bifurcations involved in the bursting mechanism, bursting models can be classified in different categories: elliptic (E) bursting (subcritical Hopf bifurcation/fold limit cycle bifurcation), square-wave (SW) bursting (saddle-node bifurcation/homoclinic bifurcation), parabolic (P) bursting (saddle-node on a limit cycle bifurcation); see Izhikevich_burst_classif () for more details.

As shown in Figure 11, each type of bursting is associated with a particular range of values for the phase sensitivity coefficient, an observation which reflects different values of fractal dimension and overall phase sensitivity. Elliptic (E) bursting models considered in our analysis (in blue) are characterized by a high phase sensitivity coefficient , and therefore by strong fractal properties. This is in agreement with the preliminary observations of MRMM_bursting (). In contrast, parabolic (P) bursting models (in red) have a phase sensitivity coefficient equal to , thereby exhibiting no fractal properties. Square-wave (SW) bursting models (in green) roughly correspond to an intermediate situation with a low phase sensitivity coefficient. They are characterized by no fractal properties (HR model) or weak fractal properties (ML model).

Our analysis based on the phase sensitivity coefficient is confirmed by the following numerical experiment. We consider a network of neurons with random initial conditions on the limit cycle (uniform distribution in phase) and we assume that these neurons receive a common impulsive input , where is a vector in the direction (i.e. membrane voltage). The state of the neurons instantaneously jumps to , which corresponds to an asymptotic phase that we compute. Then we perform the same simulation with identical initial conditions but with a slightly perturbed input for each neuron, where is a small random variable. In this case, the neurons jump to the state and we compute their corresponding phase . Depending on the pulse size , the neurons may (or may not) reach a region of high phase sensitivity associated with a high phase error . In order to cover a large part of the state space, we consider different pulse sizes and obtain statistical results on that are consistent with the values of the phase sensitivity coefficient (Figure 12, see also Appendix D for detailed results obtained with different pulse sizes). It is also noticeable that the phase sensitivity coefficient, computed on a very small subset of the state space (Figure 11), captures well the sensitivity of the network, which is computed globally in a large region of the state space (Figure 12).

Although applied to a small collection of models, the results tend to show that elliptic bursting models are characterized by a very high overall phase sensitivity. For the elliptic bursting neurons tested, a small uncertainty on the input signal (or a small noise perturbation) may induce high variation in phase (see Figure 12), and a reduction of this uncertainty only slightly reduces the uncertainty on the phase. These neuron models are therefore characterized by sensitive and unreliable responses to inputs. In particular, it is clear from the results of Figure 12 that these neurons with identical initial phases would not remain synchronized under the effect of a common (impulsive) input with a (very small) additive noise. Also, the phase response curve of these models has strong fractal properties (Figure 13(a)) and cannot be used in practice since it is impossible to know exactly its values (in the fractal regions, see Figure 13(b)). Note that the sensitivity observed here only results from the properties of the neuron dynamics and does not depend on the type of forcing (or the type of coupling in a network). This is in contrast to high sensitivity to initial conditions usually reported in neuron models with an external input (e.g. shear-induced chaos in the periodically kicked Morris-Lecar neuron Lin_shear_chaos ()).

## Vi Conclusion

We have discussed the fractal properties of asymptotically periodic systems and we have provided a theoretical framework to quantify these properties. Through the notion of phase sensitivity coefficient, the fractal capacity dimension of the isochrons has been related to the overall phase sensitivity of the system.

The main implication of the results is that there exist systems—with isochrons of high fractal dimension—that are characterized by a (extremely) high phase sensitivity. For these systems, reducing the intensity of a noise perturbation only slightly decreases the average uncertainty on the phase. This is for instance the case for some elliptic bursting neuron models, whose response to external inputs is unreliable.

We will finally note that the rich geometric properties of asymptotically periodic systems illustrate the importance of developing efficient methods for computing the isochrons and the phase function. It is only when these methods are well-suited to complex dynamics in high dimensional spaces that they may unveil new properties, such as the fractal isochrons described in this paper.

## Acknowledgments

This work was funded by Army Research Office Grant W911NF-11-1-0511, with Program Manager Dr. Sam Stanton. It was completed while A. Mauroy was with the Department of Mechanical Engineering, University of California Santa Barbara. A. Mauroy is currently supported by a BELSPO (Belgian Science Policy) return grant.

## Appendix A Scaling of ∥∇Θ∥

The (infinitesimal) phase difference between two trajectories

 Θ(φ(−t,x+dx))−Θ(φ(−t,x))=∇ΘT(φ(−t,x))(φ(−t,x+dx)−φ(−t,x))

is constant. It follows that we have, for all infinitesimal ,

 ∇ΘT(x)dx=∇ΘT(φ(−t,x))M(t)dx,

where is the fundamental matrix solution of

 dMdt=J(φ(−t,x))M

with and is the Jacobian matrix of the vector field .

If we assume that is normally hyperbolic and of co-dimension (at most) , we can choose in the direction tangent to the fibers of and we obtain

 ∥∇Θ(x)∥∥dx⊥∥cosβ(0)=∥∇Θ(φ(−t,x))∥∥M(t)dx⊥∥cosβ(t) (21)

with the angle between the gradient and the direction tangent to the fiber. Since and for all , it follows that either for all or for all . The first case (i.e. ) is not possible since the isochrons are not the fibers of (otherwise, would not be phaseless). For the same reason, the isochrons are not tangent to the normal bundle of in the neighborhood of , so that one cannot have as . Also, one has

 ∥M(t)dx⊥∥=O(eλ⊥t)t→∞

where is the (negative) Lyapunov exponent of the system in backward time, associated with the normal direction . (Note that the value is the generalized Lyapunov type number along the fiber Fenichel1 (); Wiggins ().) It follows from (21) that

 ∥∇Θ(φ(−t,x))∥=O(e−λ⊥t)t→∞.

Since the distance between the trajectories and satisfies as , we have finally

## Appendix B Bursting neuron models

### b.1 Modified Morris-Lecar (ML) model (square-wave, elliptic, or parabolic bursting)

 C˙V = −gCam∞(V)(V−VCa)−gKn(V−VK)−gL(V−VL) −gKCaz(h)(V−VK)−gCaSs(V−VCa)+I ˙n = ϕ(w∞(V)−n)/τ(V) ˙h = ϵ1(−μgCam∞(V)(V−VCa)−h) ˙s = ϵ2(s∞(V)−s)/τs

with

 m∞(V) = 0.5(1+tanhV−V1V2) w∞(V) = 0.5(1+tanhV−V3V4) s∞(V) = 0.5(1+tanhV−V5V6) z(h) = hCa0+h, τ(V) = cosh−1V−V32V4.

The parameters are given in Table 1.

### b.2 Hindmarsh-Rose (HR) model (square-wave bursting)

 ˙V = n−aV3+bV2−h+I ˙n = c−dV2−n ˙h = r(σ(V−V0)−h)

The parameters are , , , , , , , .

### b.3 FitzHugh-Rinzel (FR) model (elliptic bursting)

 ˙V = V−V3/3−w+y+I ˙w = δ(a+V−bw) ˙y = μ(c−V−dy)

The parameters are , , , , , , and .

### b.4 Plant model (parabolic bursting)

 C˙V = −gNam3∞(V)h(V−VNa)−gCax(V−VCa) −(gKn4+kKCac0.5+c)(V−VK)−gL(V−VL) ˙n = (h∞(V)−h)/τh(V) ˙h = (n∞(V)−n)/τn(V) ˙x = (x∞(V)−x)/τx ˙c = f(k1x(VCa−V)−c)

with

 w∞(V) = α∞(V)α∞(V)+β∞(V)for w=m,h,n τ∞(V) = 12.5α∞(V)+β∞(V)for% w=h,n x∞(V) = 1exp(−0.15(V+50))+1

where

 αm(V) = 0.150−Vsexp((50−Vs)/10)−1 βm(V) = 4exp((25−Vs)/18) αn(V) = 0.0155−Vsexp((55−Vs)/10)−1 βn(V) = 0.125exp((45−Vs)/80) αh(V) = 0.07exp((25−Vs)/20) βh(V) = 1exp((55−Vs)/10)+1

and . The parameters are , , , , , , , , , ,