COMPUTATION OF THE TOPOLOGICAL ENTROPY IN CHAOTIC BIOPHYSICAL BURSTING MODELS FOR EXCITABLE

One of the interesting complex behaviors in many cell
membranes is bursting, in which a rapid oscillatory state
alternates with phases of relative quiescence. Although there is
an elegant interpretation of many experimental results in terms of
nonlinear dynamical systems, the dynamics of bursting
models is not completely described. In the present paper, we study
the dynamical behavior of two specific three-variable models from
the literature that replicate chaotic bursting. With results from
symbolic dynamics, we characterize the topological entropy of
one-dimensional maps that describe the salient dynamics on the
attractors. The analysis of the variation of this important
numerical invariant with the parameters of the systems allows us
to quantify the complexity of the phenomenon and to distinguish
different chaotic scenarios. This work provides an example of how
our understanding of physiological models can be enhanced by the
theory of dynamical systems.


Introduction and preliminaries
The study of bursting is crucial in physical and biological systems.This complex behavior is particularly interesting in neural systems where it plays an important role in information processing [10][11][12]15].Bursting oscillations are ubiquitous in the experimental and modeling studies of excitable cells.Biological significance and dynamical complexity of bursting oscillations stimulated mathematical investigations of the mechanisms of bursting behaviors (see [6,16,20], and references therein).
It is of biological interest to ask how the rhythm of bursting can come about the action potentials, starting from the "building blocks" of the neural electrical activity.This question leads to challenging mathematical problems.The bursting activity of individual excitable cells is the result of high-dimensional dynamics of nonlinear events responsible for variations in the ionic currents across the membrane.The numerical studies of such neural activity are usually based on either realistic ionic-based models or phenomenological models.The ionic-based models proposed for a single neuron are designed to represent the biophysical mechanisms of the membrane, with the parameters and functions derived from experimental data.Some of these models contain many nonlinear differential equations.The strong nonlinearity and high dimensionality of the phase space is a significant obstacle in understanding the qualitative behavior of such dynamical systems [9].An alternative approach is to construct a polynomial model that retains the important qualitative features but is simpler to analyse and understand [8].Such a model has the same relationship to the above models as the FitzHugh-Nagumo model does to the Hodgkin-Huxley model: it is phenomenological in nature and based on the fact that the underlying behavior of biophysical models can be distilled into a simpler model containing only polynomials.
In order to isolate essential aspects of bursting dynamics, there have been proposed systems that contain just three dynamical variables.Indeed, the study of dynamical principles and mechanisms underlying the bursting behavior is hardly possible without numerical studies of low-dimensional models.In this work, we focus our attention on two classical three-variable systems: the ionic model of Chay for the bursting of the pancreatic β-cell [2], and the phenomenological model of Hindmarsh-Rose [8].These models have the form of a three-variable autonomous system: with fast variables x and y, and slow variable z.In order to understand the dynamics of models, it is necessary to explain the observed behavior from basic features of the differential systems.For example, we can gain some important qualitative insights by studying representative return maps, considered as one-dimensional interval maps (e.g., see [2,22]).The aim of this paper is to provide a contribution for the detailed analysis of the chaotic behavior of the bursting models.More precisely, using results of symbolic dynamics theory, we compute the topological entropy of iterated maps (with the shape of a logistic map) which incorporate the important dynamical properties of the attractors.The variation of this measure of the amount of chaos in a dynamical system with the parameters gives us a finer distinction between different states of complexity.

Discrete dynamics, topological entropy, and chaotic maps
In this paragraph, we describe techniques of symbolic dynamics, in particular some results concerning to Markov partitions associated to unimodal maps (family of continuous maps on the interval with two monotonic subintervals and one turning point).For more details see [7,14,17,18].
A unimodal map f on the interval I = [c 0 ,c 2 ] is piecewise monotone and I is subdivided into two subintervals: Jorge Duarte et al. 3 in such a way that the restriction of f to interval L is strictly increasing and the restriction of f to interval R is decreasing.Every such maximal interval on which the function f is monotone is called a lap of f , and the number = ( f ) of distinct laps is called the lap number of f .Denoting by c the turning point (relative extremum) of f and with the purpose of studying the topological properties, we associate to each orbit O(x) a sequence of symbols, itinerary (i(x)) j = S = S 1 S 2 ...S j ..., where S j ∈ Ꮽ = {L, C,R} and The point c plays an important role as well as its orbit (2. 3) The dynamics of the interval is characterized by the symbolic sequence associated to the orbit of point c, that is, the turning point itinerary.When O(c) is a k-periodic orbit, we obtain a sequence of symbols that can be characterized by a block of length k, the kneading sequence [17].We introduce, in the set of symbols, an order relation L < C < R.
The order of the symbols is extended to the symbolic sequences.Thus, for two of such sequences P and Q in Ꮽ N , let i be such that P i = Q i and P j = Q j for j < i. Considering the R-parity of a sequence, meaning odd or even number of occurrence of a symbol R in the sequence, if the R-parity of the block And if the R-parity of the same block is odd, we say that The ordered sequence of elements x ij of O(c) determines a partition ᏼ of the interval x 1 ] into a finite number of subintervals labeled by I 1 , I 2 ,...,I k−1 .
To this partition we associate a (k − 1) × (k − 1) transition matrix M = [a i j ] with (2.4) Take, for example, the period-5 kneading sequence S = RLRRC.Its successive points of the orbit are obtained by shifting the periodic sequence by one symbol at a time, that is, 4 Computation of the topological entropy in bursting models Partition of the interval corresponding to the kneading sequence RLRRC.
These points are ordered in the following way: The dynamical invariant range is now divided into four subintervals, as shown in Figure 2.1.By inspecting the partition ᏼ of the interval I = [x 2 ,x 1 ], we write down the transition matrix (2.6) Now we consider the topological entropy.This numerical invariant allows us to quantify the complexity of the phenomenon.A possible definition of chaos in the context of one-dimensional dynamical systems states that a dynamical system is called chaotic if its topological entropy is positive.Thus, the topological entropy can be computed to express whether a map has a chaotic behavior.
The topological entropy of a unimodal interval map f , denoted by h top ( f ), is given by where λ max (M( f )) is the spectral radius of the transition matrix M( f ) and s( f ) is the growth rate, of the lap number of f k (kth-iterate of f ) (see [13,17,19]).Now, regarding the previous considerations we discuss dynamical properties of the three-variable bursting models of Chay and Hindmarsh-Rose.In order to facilitate the analysis and make this paper self-contained, we outline briefly some important aspects of these models pointed out in previous papers (for further information see [2,8]).

The Chay β-cell model. One of the first models for bursting was proposed by
Atwater et al. [1].It was based on extensive experimental data, incorporating the important cellular mechanisms that were thought to underlie bursting.Following this experimental work, Chay and Keizer developed a mathematical model for the ionic and electrical behavior of the pancreatic β-cell [3].Chay reduced the model to three variables [2], in the context of the general model (1.1), presented in the first section.The Chay's βcell model was recently studied using qualitative methods for singular perturbed systems of differential equations [16].
Jorge Duarte et al. 5 The ionic events in excitable membranes may be represented by the following system: where (2.10) An excitable membrane is considered to contain voltage-sensitive channels which allow Na + and Ca 2+ ions to enter the cell, voltage-sensitive K + channels which allow K + ions to leave, and voltage-insensitive K + channels known to be activated by intracellular calcium ions [1].
The excitable cell model consists of the following three dynamical variables: (1) the membrane potential, denoted by V , is the difference of external and internal voltages; (2) the variable n is the probability of opening the voltage-sensitive K + channel; (3) the variable C is the intracellular concentration of Ca 2+ ions.The fast variables are V and n.The slow variable is C, with ρ the time scale parameter.For further information concerning the physiological significance of the variables, the reader is referred to the papers [1,2].The values of the parameters used in the numerical computations are as a bifurcation parameter.
In order to illustrate the chaotic nature of the model, it is interesting to exhibit some numerical results about the behavior of the dynamical variables, namely, (1) the dynamics of membrane potentials, which reveal the transition from a regular to a nonperiodic mode; (2) the time variation of C, which closely follows the mode of the potentials; (3) iterated maps on the interval, which confirm that the irregular feature of the action potentials is indeed deterministic in nature.
6 Computation of the topological entropy in bursting models      in Figure 2.12, the data from the chaotic time series appear to fall on a logistic curve.Indeed, there is so much to be gained by treating the graph as a function At this point, we direct our attention to the study of the topological entropy.In order to illustrate the outlined formalism about the computation of this numerical invariant, we discuss the following example (for more details see [4]).

The Hindmarsh-Rose model.
We begin this paragraph presenting a three-variable model for the bursting of neurons developed by Hindmarsh and Rose [8], which is a modified version of the FitzHugh-Nagumo model [5].Much of their paper concerns the xy-subsystem, which describes the action potential.Hindmarsh and Rose added an adaptation current z governed by a linear equation which gradually hyperpolarizes the cell and Jorge Duarte et al. 11 approaches a steady-state value.With this procedure, the model did not fire indefinitely and has the nice property of generating oscillations with a long interspike interval.As pointed out in [8], the resulting autonomous system admits a periodic behavior.We consider this system in the context of the general model (1.1).
The Hindmarsh-Rose model with adaptation has three variables x, y, z and three equilibrium points, satisfying the following polynomial equations: (2.15) For details concerning the physiological significance of the variables, the reader is referred to the original paper [8].The bursting behavior is accomplished in this model by the presence of the variable z that modulates the applied current I on a slower time scale.In this context, the third variable is the slow variable and r is the time scale parameter.In this model (x 1 , y 1 ) are the coordinates of the leftmost equilibrium point of the model without adaptation (i.e., without the variable z).This has the result that (x 1 , y 1 ,0) is an equilibrium point of the model with adaptation.
For numerical investigation we will use a = 1, b = 3, c = 1, and d = 5 as in [8].We obtain where )/2 (the x-coordinate of the resting state in the xy-subsystem) and I is the applied current.The responses of this model to a short depolarizing current pulse depend on the values given to the parameters r and s.
To see what the trajectories do in the long run we use numerical integration.After an initial transient, a structure emerges when the solution (x(t), y(t),z(t)) is visualized as a trajectory in three-dimensional space (Figure 2.15).Some examples of projections of the three-dimensional trajectory onto a two-dimensional plane are presented in Figures 2.16 and 2.17.We show the burst pattern for the current pulse I = 3.25 (Figure 2

.18).
With the purpose of understanding the main features of the three-dimensional flow, we can use the Poincaré map technique to reduce the dimensionality of the phase space and so make the analysis simpler.Now we focus on the construction of a Poincaré map.
Consider an n-dimensional system dx/dt = f (x).Let P be an (n − 1)-dimensional surface, called a Poincaré section.P is required to be transverse to the flow.The Poincaré map   F is a map from P to itself, obtained by following trajectories from one intersection with P to the next.If x n ∈ P denotes the nth intersection, then the Poincaré map is defined by x n+1 = F(x n ).In our particular case, we choose the Poincaré plane defined by x = −0.85.After allowing the initials to decay, we record the successive intersections of the trajectory with the plane, which are specified by two coordinates (y n ,z n ).The iterated map of Figure 2.19 consists of pairs (z n ,z n+1 ), obtained from the successive second coordinates of the points defined by the Poincaré map.
In order to see the long term behavior for different values of the parameters at once, we plot typical bifurcation diagrams (Figures 2. 20-2.22).
With the same procedure used in the case of Chay's model, we compute the topological entropy of the interval maps.Several situations of the variation of the topological entropy with each of the parameters are plotted in Figures 2.23    We note that the Hindmarsh-Rose model exhibits a chaotic behavior in a certain range of the parameters.The topological entropy allows us to quantify the orbit complexity and to extract the order from chaos.

Final considerations
In this paper, we have studied the chaotic dynamics of two specific three-variable bursting models.
A detailed analysis of the iterated maps, that incorporate the salient dynamical properties of the systems, became possible by the study of the variation of the topological entropy with the parameters.
Indeed, the biophysical models exhibit positive topological entropy, which means that in certain conditions the bursting behavior of excitable cells has a chaotic nature.
In the presence of the evidence for the existence of chaos, we are invited to meditate about its role in neurophysiology, in particular, in neural systems.Why evolution has selected chaos as an apparently typical pattern of behavior in neural systems?What might be accomplished by this choice?And perhaps slightly elusively, what could not have been accomplished by the choice of regular, predictable behavior?These delicated questions, of incontestable pertinence, are object of attention of several authors (e.g., see [21]).A central aspect of the discussion consists of understanding how chaos is employed by neural systems to accomplish biologically important goals.These types of questions are quite different than a classical approach associated with the usual inquiry where in the structure of a neural or any other system the origin of dynamical chaos lies.The questions presented go beyond the purely technical and it is plausible that they have no definitive answer.
Chaotic motions explore a broad sector of the system state space, which give a means to explore the opportunities available to the system when the environment changes.Therefore, chaos can act as a precursor to adaptive, reliable, and robust behavior for living systems.This behavior differs from regular, rigidly confined motions which lie on a particular region of the state space.
Jorge Duarte et al. 17 The nonlinear dynamics of systems with chaos facilitates the extraordinary ability of neural systems to adapt, make transitions from one pattern of behavior to another when the environment is altered, and consequently create a rich variety of patterns.
The demonstration that chaotic neurons can synchronize and transmit information tells us that these states of motion, useful for evolution, also allow collective activity directed toward useful functions and goals [21].The investigations suggest that nature uses complex dynamics of neural assemblies in promoting principles of adaptability and reliability as well as in providing rapid response to changing external stimuli for information processing and response.These arguments state the physiological relevance of the qualification of chaos in neural models.

Figure 2 . 20 .
Figure 2.20.Bifurcation diagram for z n as a function of I, with r = 0.005 and s = 4.0.

Figure 2 . 21 .
Figure 2.21.Bifurcation diagram for z n as a function of r, with I = 3.25 and s = 4.0.