Abstract

The dynamical complexity of a system of ordinary differential equations (ODEs) modeling the dynamics of a neuron that interacts with other neurons through on-off excitatory and inhibitory synapses in a neural system was investigated in detail. The model used Morris-Lecar (ML) equations with an additional autonomous variable representing the input from interaction of excitatory neuronal cells with local interneurons. Numerical simulations yielded a rich repertoire of dynamical behavior associated with this three-dimensional system, which included periodic, chaotic oscillation and rare bursts of episodic periodicity called the transient periodicity.

1. Introduction

In this paper, we studied a three-dimensional extension of the Morris-Lecar (ML) system [1, 2]. The ML system is a two-dimensional representation of the four-variable Hodgkin-Huxley system [3]. It assumes that the conductivity of Na+ ions is very high. The ML system supports only two types of dynamics, resting and sustained spiking activity and does not exhibit chaotic oscillations. In this paper, we explored the dynamics of a model system that had an extra variable added to the ML system. The model described a condition in which excitatory principal cells interacted with local inhibitory inter neurons. These cells are present in a small volume of physiological space. Hence, the effect of the spatial structure of individual cells on the model dynamics was ignored.

The emergent behavior of a population of neurons was studied by considering population of neurons as the basic unit of an excitatory pathway connecting to inhibitory interneurons. The feedback inhibition present in (1) and (3) is an important part of the system architecture. Thus, by defining a population of neurons as a dynamical system [4], we studied the temporal behavior of neurons interacting with each other as well as with the other population. The 𝑉-π‘Š subsystem was governed by a faster time scale compared to that of state variable 𝑍. All interesting behavior of the model system resulted from an interaction between the faster subsystem and a slower one with the population of inhibitory inter-neurons. The differential equations that represent the dynamical system were derived by considering the behavior of prototypical single neurons (ML neurons).

The mean field equations describing the time evolution of neural system are given as follows:𝑑𝑉𝑑𝑑=βˆ’π‘”πΆπ‘Žπ‘šβˆž(π‘‰βˆ’1)βˆ’π‘”πΎπ‘Šξ€·π‘‰βˆ’π‘‰πΎξ€Έβˆ’π‘”πΏξ€·π‘‰βˆ’π‘‰πΏξ€Έ+πΌβˆ’π›Όinh(𝑍)𝑍,(1)π‘‘π‘Š=πœ™ξ€·π‘Šπ‘‘π‘‘βˆžξ€Έβˆ’π‘Šπœπ‘Š,(2)𝑑𝑍𝑑𝑑=𝑏𝑐𝐼+𝛼excξ€Έ,(𝑉)𝑉(3) where 𝑉 and 𝑍 are the mean membrane potentials for the excitatory and inhibitory cells. π‘Š is the fraction of open potassium channels in the population of principal cells. 𝑉𝐾 and 𝑉𝐿 are the Nernst potentials for the potassium ions and leakage channels. πœπ‘Š, a voltage dependent constant, is given by πœπ‘Š=1coshξ€·ξ€·π‘‰βˆ’π‘‰3ξ€Έ/2𝑉4ξ€Έ(4) where 𝑉3 and 𝑉4 are the midpoint potentials at which the potassium and calcium currents are half-activated. The π‘”πΆπ‘Ž and 𝑔𝐾 are the total conductances for the population of calcium and potassium channels, respectively. The constant 𝑔𝐿 denotes β€œleakage conductance,” and it measures the effect of the channels that are always open. The parameter 𝑐 represents the relative strength of current into the inhibitory interneurons to the current into principal cells. The parameter 𝐼, the applied current, is maintained as a constant stimulus. The changing polarity of the mean membrane potential of the principal cells ensured that the mean membrane potential of the inhibitory cells did not grow beyond manageable limit.

The important variables in the model are π‘šβˆžξ‚Έξ‚΅(𝑉)=0.51+tanhπ‘‰βˆ’π‘‰1𝑉2π‘Šξ‚Άξ‚Ή,(5)βˆžξ‚Έξ‚΅(𝑉)=0.51+tanhπ‘‰βˆ’π‘‰3𝑉4𝛼,(6)exc(𝑉)=𝛼excξ‚Έξ‚΅1+tanhπ‘‰βˆ’π‘‰5𝑉6𝛼,(7)inh(𝑍)=𝛼inhξ‚Έξ‚΅1+tanhπ‘βˆ’π‘‰7𝑉6ξ‚Άξ‚Ή.(8)π‘šβˆž and π‘Šβˆž are nondimensionalized functions to describe the voltage-regulated calcium and potassium ion channels, respectively. The parameters 𝑏 and πœ™ are temperature scaling factors. The hyperbolic tangent function describes the collective behavior of a large number of channels. On the right-hand sides of (7) and (8), 𝛼exc and 𝛼inh are dimensionless synaptic strengths of excitatory principal cells and inhibitory inter-neurons. On the left of these equations are functions that are parameterized by these dimensionless constants. These functions were assumed to be of tangent hyperbolic form. 𝑉6 is the steepness parameter for both the functions. 𝑉5 and 𝑉7 are the thresholds for the synaptic strength functions of the excitatory principal cells and inhibitory interneurons. Both the functions have the same functional properties. The negative sign in (1) signifies an inhibitory effect from the interneurons with mean membrane potential 𝑍 on the principal cells. Parameter values for 𝑉1, 𝑉2, 𝑉3, and 𝑉4 are chosen in such a way that nondimensionalized functions π‘šβˆž and π‘Šβˆž reach their equilibrium values almost instantaneously. The aim of the paper is to report the rich dynamical repertoire displayed by the system of ODEs.

The key parameters of the model are 𝑏, 𝑐, and 𝐼. The parameter 𝑏 ensures that the inherent time scale of the group of inhibitory cells is different from the rest of the system. The parameter 𝑐 measures the strength of feedforward inhibition. When 𝑉>𝑉5 (a fixed threshold potential), the connection between the principal cell and the interneuron is opened. The parameter 𝐼 represents constant external current. It turned out to be useful for studying the effects of disinhibition factors on model dynamics. Yet another key parameter of this model system is 𝑉6, which controls the steepness of the β€œon-off” switch in 𝛼exc(𝑉) and 𝛼inh(𝑍). For a single presynaptic terminal, excitatory or inhibitory, one would expect β€œall or none” behavior. In that case, one could have a simple representation of the on/off switch via a Heaviside step function. In a population of cells, the connection strengths vary smoothly as distributions of thresholds exist and the values of parameters 𝑉5 and 𝑉7 represent the end points of this distribution.

Inhibitory inter-neurons are found in many brain regions. The ML system does not consider the effect of the inhibitory inter-neurons on the excitatory principal cells. Since it is based on two primary system variables, the M-L system cannot support deterministic chaos. We explored regular and chaotic dynamics, which resulted from an interaction of principal cells with inhibitory inter-neurons. When the values for parameters of the M-L system of equations (as described above) were suitably chosen, we observed Ca++ and K+ oscillations. The action of inhibitory inter-neurons on the excitatory system (equations (1) without the last term and (2)) generated chaotic dynamics as the third variable (Z) was added to the system. The base parameter values were chosen in such a way that the M-L system displayed sustained spiking activity, which has a stable limit cycle as the underlying attractor. The remaining parameters were set in such a way that dynamics of the model system unfolded on a chaotic attractor. The base parameter set for which the model system (equations (1), (2) and (3) displayed oscillations is given in Table 1.

The Morris-Lecar equations are applicable to a spatially isopotential patch of membrane injected with a constant current. Transition from single shot firing to a stable limit cycle (tonic firing) occurred as the excitation current was varied. At very low or very high values of the stimulus, it had a single stable fixed point. The neurons of the inhibitory system provided an additional degree of freedom as it is well known that a minimum of three degrees of freedom are required for a system to support complex dynamics including chaos. We chose the set of values presented in Table 1 for the system parameters as base values. At these values of parameters the system displayed oscillatory behavior. Table 1 gives the chosen parameter values.

2. Simulations

Two-dimensional parameter scans were performed for (𝑏,𝑉𝐾) parameter space with 𝑏 in the range of 0.1 to 0.15 with a step size of 0.01 and π‘‰π‘˜ in the range of βˆ’0.9 to βˆ’0.6 with a step size of 0.05, and the results of these scans are presented in Figure 1. The figure plots points at which the system of differential equations exhibited chaotic solutions and nearby points in the system exhibited other dynamical behavior. We used the characteristic signature of deterministic chaos: sensitive dependence on initial conditions to detect chaos in the system. Chaotic dynamics was detected when two time series (collected after transients were allowed to die out) generated at the same set of parameter values, but, for two different initial conditions which were nearby, did not overlap. Computations of the largest Lyapunov exponents were also performed to confirm the regular or chaotic behavior observed at a chosen set of parameters. Figures 2(a) and 2(b) demonstrate the onset of chaos through period doubling route. We have also examined bifurcation diagrams varying the parameter π‘‰π‘˜ and confirmed that the system exhibits period doubling route to chaos. This portrays sustained spiking activity displayed by the neuronal system for a given combination of temperature scaling factor and the Nernst potential for potassium in the node. The chaotic oscillations of the system are shown in Figures 3(a) and 3(b), and the largest Lyapunov exponent for the attractor is 0.0621.

The dynamical behavior of the system was sensitive to changes in both the parameters. It changed from simple periodic oscillations to complex periodic ones and at times transitioned to chaos. Two-dimensional parameter scans showed how sensitive the dynamics of the model system was when it encountered a change in an intrinsic attribute of the system. At certain points (e.g., 0.12, βˆ’0.7), the system exhibited bistable behavior; a periodic attractor coexisted with a chaotic attractor. A few points represented bistability, especially those at the outer edges of the parameter regions, and chaotic behavior was observed at discrete points in closed intervals.

The other two-dimensional parameter scan was carried out in (𝑏,𝑉6) with 𝑏 range from 0.10 to 0.15 with a step size 0.01 and 𝑉6 range from 0.3 to 0.8 with a step size 0.1. We found that chaotic behavior exists at (0.1, 0.40), (0.11, 0.3), (0.12, 0.3), (0.13, 0.3), (0.14, 0.4), (0.15, 0.35), and (0.15, 0.4). At a few points, neuronal dynamics displayed resting behavior.

Figures 4(a) and 4(b) show stable limit cycle attractor and corresponding time-series at the point where bi-stability was detected.

The third parameter scan was performed in (𝑐,𝑉𝐾) with the following ranges and step sizes. The range for 𝑐 was from 0.165 to 0.25 with step size 0.002, and for 𝑉𝐾 the range was from βˆ’0.85 to βˆ’0.7 with a step size of 0.05. Chaos was observed at the following points: (0.236, βˆ’0.85) (0.238, βˆ’0.7), (0.238, βˆ’0.85), (0.238, βˆ’0.80), and (0.24, βˆ’0.80). Figures 4 and 5 show period-2 cycle and suggest that the system exhibited period doubling route to chaos.

The fourth and the final parameter scan was performed in (𝛼exc,𝛼inh), and it showed chaos at the following points: (1, 0.9), (1.0, 1.1), (1.1, 0.8), (1.1, 0.9), (1.1, 1.0), (1.2, 0.8), (1.2, 0.9), (1.2, 1.0), (1.3, 0.8), (1.3, 0.9), and (1.3, 1.0). Figure 6 shows the presence of β€œedge of chaos” in the neural system. Edge of chaos is a dynamical phenomenon that is characterized by the existence of chaotic behavior in finite but arbitrarily small parameter ranges. This suggests that smooth changes in the strengths of excitatory and inhibitory synapses cause dynamical transitions from chaotic behavior to regular oscillatory dynamics. The model system also show chaotic saddle which is shown in Figure 7.

Results of time series comparisons were verified by computing the largest Lyapunov exponents (πœ†max) using the Wolf et al. algorithm [5]. We present πœ†max with standard error computed from 20 trials of random initial conditions of the system in Table 2. The choice of the delay for reconstruction of the dynamics from the time series was made by choosing the first minimum of mutual information.

We present an interesting finding of transient periodicity displayed by this neural system, which was detected in the (𝑏,𝑉6) parameter space which is shown in Figure 8.

3. Results

Chaotic dynamics is characterized by sensitive dependence on initial conditions. The error committed in fixing initial conditions grows exponentially with respect to time. The largest Lyapunov exponent is a quantity that characterizes the rate of separation of infinitesimally close trajectories. The time series is termed as chaotic if it returns a positive value for the Lyapunov exponent. Computer simulations of our model system suggested that deterministic chaos manifests itself as short-term recurrent phenomenona.

We discovered rare bursts of β€œperiodic” motion at 𝑏=0.12, 𝑉6=0.3 (cf. 8b). This phenomenon is known as β€œtransient periodicity” in the dynamical systems literature. This is caused by the presence of semiperiodic saddles (unstable invariant sets) in the chaotic attractor. Motion in the vicinity of these sets is periodic. In the case of EEG time series obtained by recording averaged action potentials of epochs of epileptic seizure [6], burst pattern often contains nearly periodic dynamics. In chaotic systems, such transient periodicity can reflect the existence of semi-periodic saddles contained in the attractor. Motion in the vicinity of such objects has a prominent periodic component as trajectories are temporarily attracted to these neighborhoods before exiting.

4. Discussion

Model simulations and animal experiments suggest that seizure activity occurs when a critical mass of neurons join and get involved in time-linked high frequency discharging [4]. An asynchronous-to-synchronous transition (longer-time state) is caused either by deterministic changes in system parameters or spontaneous dynamical shifts from chaos to order caused by transient periodicity in chaotic saddles. This situation has a parallel in our model system wherein the dynamics switched to sustained spiking activity or to resting behavior from chaotic behavior when small changes in crucial parameters occurred. This dynamical behavior is known as β€œedge of chaos.” Chaotic saddles are semiperiodic semiattractors that are distinct from the better known attractors in nonlinear dynamical systems which are created during boundary crisis [7]. Presence of these unstable invariant sets generates strong episodic periodic components embedded in chaotic dynamics. This has been found in epidemiological models and also in EEG time series from patients of epilepsy [8, 9].

We envision that chaos is vital for functioning of a healthy brain and synchronization of the neuronal system occurs when all its regions are in transient periodicity represented by chaotic saddles in state space. This is how the intermittent pathology of epileptic seizures is created. It is understood that lack of inhibition [6, 10] caused either by a reduction in the number density of excitatory cells that synapse on the inhibitory cells in the focal region or by the diminished release of neurotransmitter 𝛾-amino butyric acid in a nearby region leads to seizures (transient periodicity). The ictal period is represented by periodic oscillations and the interictal period by chaotic phase. Stable limit cycles and complex periodic dynamics observed in a narrow range of model parameter, 𝑐, correspond to a biophysical state of partial dis-inhibition. Biophysical experiments performed previously [11, 12] have attempted to simulate dynamics observed in this range. The transient periodicity, which was detected for higher values of this parameter, has been observed in electrographic recordings from depth and subdural electrodes performed in epilepsy patients [8].

Acknowledgments

The authors thank Dr. Pravitha Ramanand for discussions and Professor Raima Larter for making available a few chapters of B. Speelman's doctoral thesis.