Neuronal Coding of pacemaker neurons - A random dynamical systems approach

The behaviour of neurons under the influence of periodic external input has been modelled very successfully by circle maps. The aim of this note is to extend certain aspects of this analysis to a much more general class of forcing processes. We apply results on the fibred rotation number of randomly forced circle maps to show the uniqueness of the asymptotic firing frequency of ergodically forced pacemaker neurons. The details of the analysis are carried out for the forced leaky integrate-and-fire model, but the results should also remain valid for a large class of further models.


Introduction
Already in 1907, long before the molecular mechanisms of neural signal transduction had been clarified, Louis Lapicque proposed a simple model for the firing behaviour of a neuron [1,2,3]. A crucial feature of this so-called integrate-and-fire model (IFM) is the separation of time-scales: the stereotypical and extremely fast generation of an action potential is thought of as being concentrated in a single moment of time, whereas the much slower evolution of the membrane potential in the interspike intervals is modelled as a continuous process. For many questions concerning the behaviour of neural systems this level of abstraction turned out to be exactly the adequate one, such that even nowadays, more than a hundred years after Lapicque's original paper, the different variations of the IFM still play a central role in theoretical neuroscience [4]. One of their great achievements was the explanation of so-called 'paradoxical segments' that were discovered in the experimental investigation of pacemaker neurons in the nervous system of crayfish (Procambarus clarkii), sea slugs (Aplysia californica) and horseshoe crabs (Limulus polyphemus) [5,6]. The counter-intuitive observation that was made in these experiments was that an increase in the frequency of periodic inhibitory presynaptic input can lead to an increase of the post-synaptic firing frequency. This paradoxon was explained by relating the respective theoretical models to monotone circle maps whose rotation number equals the ratio between the input and the output frequency. When the circle map has a stable periodic orbit, then input and output frequency remain directly proportional on a small neighbourhood (the paradoxical segment), disrespective of whether the input is inhibitory or excitatory [5,7,8,9] (see also Section 2). Similar ideas have also been pioneered before by V. Arnold in the study of cardiac cells [10,11].  [5,6]. From left to right: Procambarus clarkii, Aplysia californica and Limulus polyphemus [12].
It must be said, however, that despite the great success of IFMs and the long history of their investigation their rigorous mathematical description is still restricted to a few very special situations. In particular, the only types of external forcing that can be treated analytically so far are either periodic [13,14] or stationary stochastic input [15,16,17]. Even the superposition of the two -noisy periodic input -is mostly accessible only by numerical methods [18]. Our goal here is take up the ideas used in the analysis of the periodically forced IFM an to extend these to more general forcing processes. Thereby, we restrict ourselves to deterministic and/or random forcing, although it should be possible to adapt the approach to models generated by stochastic differential equations as well. In order to state the main results, we first recall the construction of the IFM.
The membrane potential V (t) of a neuron N1 remains between a lower threshold V l and an upper threshold Vu. V (t) can never drop below V l due to physiological constraints, whereas when it reaches Vu the neuron 'fires', meaning that an action potential is triggered and V (t) drops back to a rest potential Vr ∈ [V l , Vu). Between the two thresholds, the potential evolves according to an infinitesimal law The dependence of F on t corresponds to the influence of external time-dependent factors. The reset procedure when V (t) reaches Vu is usually expressed as where t + denotes the right-hand limit. Identifying the interval [V l , Vu) with the circle T 1 = R/Z, this gives rise to a non-autonomous circle flow (see Figure 1.1). For fixed initial values V (t0) = x0, we denote by tn the time of the n-th firing of the neuron N1. Then a very basic and fundamental question is that of the existence and uniqueness of the asymptotic firing frequency: under what assumptions does the limit exist and when is it independent of the initial values t0 and x0?
In the simplest case the external input is periodic in time with period p ∈ R + . As mentioned, this situation is quite well-understood and has been studied for a number of different versions of the IFM [13,14]. The analysis depends on the choice of a suitable Poincaré section for the flow, by which one obtains a circle map g : T 1 → T 1 whose lift G : R → R generates the rescaled sequence tn, that is tn+1/p = G n (tn/p). (We briefly recall the construction in Section 2.) Under suitable assumptions on the function F this map g has good monotonicity properties that ensure the existence and uniqueness of the rotation number (1.4) ρ(G) = lim n→∞ G n (t0/p)/n = ν −1 N 1 /p and thus of the asymptotic firing frequency. As indicated above, the existence of a stable periodic orbit for the map g yields an explanation for the 'paradoxical segments' in situations with inhibitory presynaptic input.
A good framework to study more general forcing processes is provided by the theory of random dynamical systems, as exposed in [19]. In order to model the external input we will assume that there is an underlying forcing process, given by a measure-preserving or metric dynamical system (metric DS), that is, a quadruple (Θ, B, µ, ω) where (Θ, B, µ) is a probability space and ω : R × Θ → Θ, (t, θ) → ωt(θ) and ω is a flow on Θ that leaves the probability measure µ invariant (compare [19]). We usually write ωt(θ) = θ · t. In this setting (1.1) is replaced by where F : Θ × R → R and θ0 is some initial value in Θ. In order to ensure that (1.5) generates a RDS, we have to impose some standard technical conditions on F . We say F is uniformly Lipschitzcontinuous in V if there exists a constant L > 0 such that For the moment, we will just assume that F is bounded and uniformly Lipschitz-continuous in V . These assumptions seem quite reasonable from the physiological point of view. It is also possible to weaken them further to some extent and we will discuss this issue in detail at the beginning of Section 4. Due to the lack of a suitable structure on Θ, it is not feasible anymore to take a Poincaré section in this situation. However, it turns out that instead the flow generated by (1.2) and (1.5) can be analysed directly by applying a result on the fibred rotation number of randomly forced circle flows. (Basically, what we need is a slight modification of statements from [20] that we present in Section 3.) This leads to the following result. A well-known example to which this statement applies is the leaky integrate-and-fire model (LIFM), in which the function F takes the form Here a : Θ → R + corresponds to the membrane permeability 1 whereas b : Θ → R reflects the external input. In order to apply Theorem 1.1 we have to assume that a and b are bounded. Slightly weaker conditions are again discussed at the beginning of Section 3. The LIFM was first introduced by Stein [7] and then further investigated both theoretically and experimentally by Knight [8,6]. In contrast to the so-called 'perfect integrator' or 'perfect IFM' ( ∂ ∂V F ≡ 0) with 'infinite memory' used by Lapicque, it takes into account the exponential decay of the membrane depolarisation in time after excitation.
The most subtle issue in the proof of Theorem 1.1 is the fact that when F takes negative values, then this might lead to discontinuities in the flow V (t) and to a lack of monotonicity with respect to the initial condition x0 (see Figure 2.1). At this point, the monotonicity of F in V is needed in order to recover a certain 'almost-monotonicity' property. The discontinuity does not present a problem in the setting of RDS, but impedes drawing further conclusions in the situation where the base flow ω is a uniquely ergodic 2 system on a compact metric space Θ. All these complications can be avoided by assuming that F is non-negative. In this case, we have the following. (a) The model described by (1.2) and (1.5) has a unique asymptotic firing frequency: There exists a real number ν such that for µ-almost every θ0 ∈ Θ and all x0 ∈ [V l , Vu) there holds (1.9) lim n→∞ tn/n = ν . 1 Usually the membrane permeability is chosen to be fixed, but as we will see it presents no additional cost to assume that it is time-dependent as well 2 A map or flow on a compact metric space is called uniquely ergodic, if it has exactly one invariant probability measure.
This statement can, for instance, be applied to the so-called quadratic IFM, which is given by (1.2) and (1.5) with under the additional assumption that inf θ∈Θ [14] for the analysis of this model in the periodically forced case and further references). The standard example of a uniquely ergodic base flow would be the Kronecker flow θ · t = θ + t(ω1, . . . , ω d ) with d rationally independent frequencies ω1, . . . , ω d , corresponding to the excitation of N1 by d independent pacemaker neurons. Although we will not discuss the topic in detail, we want to mention that mode-locking phenomena may appear in this setting whenever the asymptotic firing frequency is rationally related to the driving frequencies (ω1, . . . , ω d )

Periodic forcing revisited
In this section we briefly recall the analysis of the periodically forced IFM and in particular the construction of the circle map g mentioned in the introduction. This will also allow to discuss the modifications needed to extend the approach to randomly forced models. In order to simplify notation, we will from now on always assume that V l = 0 and Vu = 1.
Suppose that the function F in (1.1) is p-periodic in the second variable, that is Equivalently, we may assume that Θ = T 1 and ωt(θ) = θ + t/p mod 1 in (1.5). For simplicity, we also suppose that In this case, we may assume without loss of generality that Vr = V l = 0, since the interval (V l , Vr) is not accessible from outside and does not play a role in the description of the dynamics. The ODE (1.1) generates a flow in the sense that x(t) = Φt−t 0 (x0) is the solution of (1.1) with initial conditions x(t0) = x0. The first firing time t1 can then be expressed as t1 = t0 + inf{t ≥ 0 | Φt−t 0 (x0) = 1}. Fixing the initial value x0 = Vr, this allows to define a mapG : R → R, t0 → t1 that generates the sequence of spiking times tn, that is, tn+1 =G(tn) ∀n ∈ N0. RescalingG, we let G : The periodicity assumption (2.1) implies thatG(t + p) =G(t) + p and therefore G(t + 1) = G(t) + 1.
Consequently G is the lift of a circle map g : T 1 → T 1 . Furthermore, using (2.2) it is possible to show the map G is continuous and strictly monotonically increasing, such that g is an orientationpreserving circle homeomorphism. Such maps have a well-defined rotation number defined by the limit which always exists and is independent of t (see, for example, [22]). Of course, this immediately implies the existence and uniqueness of the asymptotic firing frequency νN 1 = 1 p · ρ(G) −1 . There is also an alternative way of constructing the map g that yields some additional insight and which we want to discuss on an informal level. Suppose we let ThenF defines a vector field on R 2 that is invariant under integer translations and hence projects to a vector field on the two-torus T 2 . Consequently, the flowΨ : projects to a flow Ψ on T 2 . It is now easy to see that the map g defined above is simply the return map of Ψ to the Poincaré section T 1 × {0}. (Note that (2.2) ensures that this section is transversal to the flow direction.) However, from this point of view we also see that taking this particular Poincaré section has something arbitrary. For example, one might just as well have chosen to define another circle homeomorphismĝ with liftĜ : R → R by taking the Poincaré map of the vertical section {0} × T 1 . This simply corresponds to stopping the flow at timestn = np. The asymptotic firing frequency could thus be obtained as νN 1 = ρ(Ĝ)/p. One can even arrive at the same conclusion without resort to Poincaré sections at all. AsF is constant in the first coordinate, the flow Ψ has skew product structure over the simple base flow ωt(t0) = t0 + t/p mod 1. In other words, Ψ is a periodically forced circle flow with liftΨ. As circle homeomorphisms, such flows have a well-defined vertical rotation number and one obtains Now, in the case of periodic forcing all this is surely a mere tautology. However, things become quite different as soon as one wants to consider more general types of forcing. When the driving space Θ is an arbitrary measurable space, then taking a vertical Poincaré section does not make sense anymore, whereas the Poincaré return map to Θ × {0} just yields a self-map of Θ that might be very difficult to analyse. In this context, the fact that we can also directly consider the forced circle flow Ψ generated by (1.2) and (1.5) turns out to be a great advantage. In the next section, we will present a result on randomly forced circle flows and their lifts that ensures the existence and uniqueness of the vertical or fibred rotation number of such skew product flows. The remaining sections are then dedicated to the construction of suitable lifts for the potential flow in the situation of Theorems 1.1 and 1.2, which immediately leads to proofs of the respective statements.
Finally, before turning to the rigorous analysis in the next sections, we briefly want to indicate why further technical problems appear when condition (2.2) does not hold anymore. First of all, this evidently leads to discontinuities of the flow when the potential 'jumps' from Vu = 1 to Vr ∈ (0, 1) or, when considering lifts, from Vu + n to Vr + n + 1 for some n ∈ Z. This would in itself not present a serious problem, since most results on the existence and uniqueness of rotation numbers only require monotonicity, whereas continuity is less important. However, it turns out that dropping condition (2.2) may at the same time lead to a lack of monotonicity of the circle flow. An example of how this can happen is indicated in Figure . This remains true even when F is non-increasing in V , as assumed in Theorem 1.1. However, in this situation it is nevertheless possible to make the argument work. In the case of periodic forcing, one can show that lift G of the circle map g constructed above is monotonically increasing on the image of one of its iterates (as done in [13,14]). On the level of forced flows, one may similarly show that while two different orbits may reverse their order, they will always remain within a bounded distance of each other. (This is the point of view we will adopt in Section 5.) In both cases, this is still sufficient to guarantee the existence and uniqueness of the rotation number.

Uniqueness of the fibred rotation number
As mentioned before, we say (Θ, B, µ, ω) is a metric DS if (Θ, B, µ) is a probability space and ω : R × Θ → Θ, (t, θ) → ωt(θ) is a flow that preserves µ, meaning that µ • ω −1 t = µ ∀t ∈ R. As before, we write ωt(θ) = θ · t. A random dynamical system (RDS) over ω is a measurable mapping that satisfies the cocycle property We say Ψ is a continuous RDS if for all θ ∈ Θ the mapping (t, x) → ψt(θ, x) is continuous. Evidently, what we are interested in is the existence of uniqueness of the fibred rotation number when Ψ is obtained as the lift of a randomly forced circle flow. For this, we will use the following assumptions: • There holds In particular, this implies that Ψ projects to a random dynamical system on the circle T 1 . with V1(t0) = x1 < x2 = V2(t0) have changed order at time s.
• There exists a constant C ≥ 0 such that In other words, the mappings ψt(θ, .) are 'almost monotone' in x, up to a uniform constant C.
• There exists a constant η ≥ 0 such that Remark 3.2. (a) When Θ is a compact metric space and ω a uniquely ergodic flow on Θ, then this result is well-known and due to Herman [23]. (See also [24] for a precursor in the context of quasiperiodic Schrödinger operators.) For the more general case of ergodically forced monotone circle maps, it was proved more recently by Li and Lu [20]. The only new aspect here is that we work with continuous time replace monotonicity by the slightly weaker property (3.5). It is not surprising that this can be proved along the same lines as the previous results, but since the proof is rather short anyway we include the details for the convenience of the reader.
(b) When Ψ is a continuous RDS on the circle, then assumption (3.5) can be replaced by the much more general one that In this case Herman's original proof, which uses the existence of a Ψ-invariant probability measure that projects down to µ, remains valid with hardly any modifications. However, since we do not assume Ψ to be continuous (and this is crucial for the application to the forced LIFM), such an invariant measure does not necessarily exist.
(3.6) and (3.9) together now now imply that for all (θ, It is easy to see that the functionρ is invariant, that isρ(θ · t) =ρ(θ) ∀t ∈ T, and the ergodicity of the base flow ω therefore implies thatρ(θ) = R Θρ dµ =: ρ on a set of full measure Θ ′′ . If we now let Θ0 = Θ ′ ∩ Θ ′′ then all the assertions of the theorem are satisfied.

Construction of the potential flow
The aim of this section is to formalise the model described by (1.2) and (1.5) in the introduction and to prove Theorem 1.2. (The proof of Theorem 1.1 will then be given in Section 5.) More precisely, we will construct a lift b V for the circle flow V that describes the evolution of the membrane potential (recall that we identify the interval [V l , Vu) with the circle. We assume without loss of generality that V l = 0 and Vu = 1. Unfortunately, the construction has a somewhat technical flavour, which basically comes from the need to treat the discontinuities produced by (1.2) in a formally precise way. However, in order to treat the problem in a rigorous way a certain amount of detail seems unavoidable, in particular since the 'almost-monotonicity property' needed for the proof of Theorem 1.1 is a rather subtle issue.
Suppose that (Θ, B, µ, ω) is a metric DS and consider the non-autonomous differential equation with F : R × Θ → R, which is equivalent to (1.5). As mentioned in the introduction, we first discuss the precise technical conditions on F needed for the construction. Given any f : R → R, we let 1] |f (x)| and (The terminology follows [19,Appendix B].) We then require that there exists a constant C > 0, such that for all θ ∈ Θ there holds The restriction to the interval [0, 1] in the definition of the norms is explained by the fact that due to the reset procedure described by (1.2) this is the only part of the phase space we are interested in. We could equally assume that the conditions are satisfied on all of R and just modify the function F if they are violated outside of [0, 1]. Assumption (4.2) is needed to ensure that the hypothesis of Theorem 3.1 in the appendix are met -roughly spoken, it just ensures that action potentials or spikes cannot be generated at an arbitrary rate (there will be at most C/(1 − Vr) + 1 spikes in any interval of length 1). (4.3) is just a standard technical condition which ensures that (4.1) generates a RDS over the base flow ω (see [19,Theorem 2.2.2]). 3 From now on, we will always assume that Φ is the skew product flow generated by (4.1).
In the following, we will give a precise definition of the lift of the potential semiflow corresponding to the model described by (1.2) and (1.5). The membrane potential at time t with initial values t = 0, θ and x0 is then obtained as V (t) = b Vt(θ, x0) mod 1. The only assumption on F needed for this construction is (4.3). (4.2) and the additional assumptions made in Theorems 1.1 and 1.2 will be required only later in order to ensure that b V meets the assumptions of Theorem 3.1. Given any x ∈ R, we denote by [x] ∈ Z its integer part and by {x} = x − [x] ∈ [0, 1) its fractional part. For any (θ, x) ∈ Θ × R we define with τ (θ, x) = ∞ if the set on the right is empty. This is the time when the first spike is generated. For any s ∈ (R + ) N we denote by S = S(s) the sequence given by Sn = P n i=0 si. Then, given (θ, x) ∈ Θ × R, we recursively define a sequence sn = sn(θ, x) and the corresponding sequence Sn = Sn(θ) by (4.7) s0 = 0 , s1 = τ (θ, x) and sn+1 = τ (θ · Sn, Vr) .
(Recall that Vr ∈ [0, 1) = [V l , Vu) is the rest potential introduced in (1.2).) Sn is the time when the n-th action potential is triggered, whereas sn is the length of the time interval between the n − 1-th and the n-th spike. Given (t, θ, x) ∈ R + × Θ × R we define n(t, θ, x) as the unique integer n such that Sn(θ, x) ≤ t < Sn+1(θ, x), that is In other words, n(t, θ, x) is just the number of spikes generated until time t. The lifted potential flow is then defined by (4.5) with   2) and (1.5) is typically not continuous and therefore a priori has many different lifts that do not only differ by an integer constant. However, the particular lift b V defined above is the unique one that 'counts' the number n(t, θ, x) of spikes (action potentials) generated up to time t, in the sense that n(t, θ, x) equals the number of integers in the interval [x, b Vt(θ, x)]. In this way, the asymptotic firing frequency is obtained as the inverse of the rotation number provided this limit exists. To show that this is the case under the assumptions made in Theorem 1.2 is the aim of the remainder of this section.
(b) We note that Further, when (4.2) holds then Finally, we note that by definition we have Consequently the mapping t → b Vt(θ, x) is continuous from the right. The only discontinuities occur exactly at times t = Sn, except in the case Vr = 0 in which t → b Vt(θ, x) is continuous.
This shows that in the situation of Theorem 1.2 the semi-flow b V satisfies all assumptions of Theorem 3.1, which proves part (a) of the theorem. Part (b) for uniquely ergodic base flows then follows from Herman's classical result on the fibred rotation number [23,Theorem 5.4].
Continuity and monotonicity. Suppose that F is strictly positive and Vr = 0. As mentioned in Remark 4.1, the mapping (t, x) → b Vt(θ, x) is continuous in t in this case. In order to show the monotonicity, fix θ ∈ Θ and x1 < x2 ∈ R and let t0 = inf{t ∈ R + | b Vt(θ, x1) ≥ b Vt(θ, x2)}. Suppose for a contradiction that t0 < ∞. Then continuity in t yields V0 := b Vt 0 (θ, x1) = b Vt 0 (θ, x2). We distinguish two cases. First, assume that V0 / ∈ Z + Vr. In this case the orbits of (θ, x1) and (θ, x2) coincide with integer translates of orbits of the flow Φ on a small interval I0 around t0. They are therefore either equal or distinct on all of I0, but cannot merge exactly at time t0. Secondly, assume that V0 = k + Vr for some k ∈ Z. Since F is strictly positive, this would imply that limt→t 0 b Vt(θ, x1) = limt→t 0 b Vt(θ, x2) = k, which would again mean that two distinct orbits of the flow Φ would have to merge at time t0. Hence, in both cases we arrive at a contradiction.
Now assume in addition that ω is a continuous flow on a topological space Θ, F is continuous and Then F is bounded on Θ × [0, 1], which together with the uniform Lipschitz continuity of F in V implies that for any compact interval I ⊆ R the flow Φ generated by (4.1) is uniformly continuous on I ×Θ×[0, 1]. Now, the flow b V is obtained by concatenating integer translates of finite trajectories of Φ in Θ × [0, 1]. b V will therefore inherit the uniform continuity of Φ, provided that no discontinuities are created by this concatenation in (4.9).
In order to see this, observe that (4.18) together with the continuity of F implies that τ defined in (4.6) is continuous in (θ, x), except when x ∈ Z. By induction, this yields that the spiking times Sn(θ, x) depend continuously on (θ, x) as well unless x ∈ Z. Consequently n(t, θ, x) is continuous in (t, θ, x) unless x ∈ Z or t = Sn(θ, x) for some n ∈ Z. Furthermore [x] is obviously locally constant when x / ∈ Z. Thus, it only remains to check that b Vt(θ, x) defined by (4.9) is continuous in (t, θ, x) when x ∈ Z or when t = Sn(θ, x). However, this can be seen quite easily by having a careful look at (4.9). When x passes an integer, then [x] will jump up by one, but at the same time n(t, θ, x) will drop down by one, such that the two discontinuities cancel each other. Similarly, when t and Sn(θ, x) change order then n(t, θ, x) has a discontinuity of size 1, but at the same time the quantity ϕt−S n(t,θ,x) (θ · S n(t,θ,x) (θ, x), Vr) jumps by one in the opposite direction.  (4.5) and (4.9) satisfies
This shows that (5.2) holds on the interval [S ′ n , Sn], and for the remaining interval (Sn, S ′ n+1 ] we can now proceed in exactly the same way as in the case n = 0. This proves (5.3) and (5.4) for all n ∈ N and hence (5.2).
In order to treat the case where x2 ∈ (1, 1 + x1) we have to obtain some information on m. More precisely, we claim that (5.8) m ≤ Vr 1 − Vr + 1 .