## Abstract

Irregular intrinsic behavior of neurons seems ubiquitous in the nervous system. Even in circuits specialized to provide periodic and reliable patterns to control the repetitive activity of muscles, such as the pyloric central pattern generator (CPG) of the crustacean stomatogastric ganglion (STG), many bursting motor neurons present irregular activity when deprived from synaptic inputs. Moreover, many authors attribute to these irregularities the role of providing flexibility and adaptation capabilities to oscillatory neural networks such as CPGs. These irregular behaviors, related to nonlinear and chaotic properties of the cells, pose serious challenges to developing deterministic Hodgkin-Huxley-type (HH-type) conductance models. Only a few deterministic HH-type models based on experimental conductance values were able to show such nonlinear properties, but most of these models are based on slow oscillatory dynamics of the cytosolic calcium concentration that were never found experimentally in STG neurons. Based on an up-to-date single-compartment deterministic HH-type model of a STG neuron, we developed a stochastic HH-type model based on the microscopic Markovian states that an ion channel can achieve. We used tools from nonlinear analysis to show that the stochastic model is able to express the same kind of irregularities, sensitivity to initial conditions, and low dimensional dynamics found in the neurons isolated from the STG. Without including any nonrealistic dynamics in our whole cell stochastic model, we show that the nontrivial dynamics of the membrane potential naturally emerge from the interplay between the microscopic probabilistic character of the ion channels and the nonlinear interactions among these elements. Moreover, the experimental irregular behavior is reproduced by the stochastic model for the same parameters for which the membrane potential of the original deterministic model exhibits periodic oscillations.

## INTRODUCTION

The study of realistic models for the electrical behavior of neurons, where one can freely manipulate the dynamics of individual conductances and test several different hypothesis about the behavior of single cells, is an invaluable tool to address several questions that neurophysiologists have to deal with, such as: how individual conductances shape the cell behavior or what is the role of these individual cellular properties in determining the behavior of a complex network (Marder 1998). Even in the smallest networks these questions are far from being answered because the neurons are themselves highly nonlinear complex systems (Holden 1997; Izhikevich 2000; Rinzel and Ermentrout 1998).

Model neurons, built on electrophysiological data from real neurons, not only constitute an adequate test-bed for hypotheses about the processes that generate a particular observed experimental behavior but can also be used to make predictions and orientate the experimenter to reveal new dynamics (Dayan and Abbott 2001; Marder 1998).

This closed loop of information exchange between experiments and models can reach really exciting new perspectives when we consider the possibilities of interfacing, in real-time, computer, or analog neuron models and living neural tissue in hybrid networks (Pinto et al. 2000; Prinz 2004; Szücs et al. 2000) to better understand on how the living network operates (Selverston et al. 2000), to restore the properties of damaged neural circuits (Szücs et al. 2000), to build circuits with new properties (Ayers 2004), or even to feed a living network with sensory data coming from artificial devices or to control such devices in brain machine interfaces (Nicolelis 2003).

Since the pioneer work introducing the voltage-clamp technique and mathematical modeling of electrical activity of the squid giant axon (Hodgkin and Huxley 1952), deterministic models of the behavior of ionic conductances have been adopted by neurophysiologists as the main paradigm of the macroscopic electrical behavior of nerve membranes.

The main reason for the success of such models is that they were able to put most of the electrical behavior of nerve cells, including the generation of action potentials, in a mathematically rigorous framework. Because most of the parameters of the model could be inferred almost directly from voltage-clamp experiments, where the contributions of many ionic populations to the current passing through the membrane can be separated and studied in detail, deterministic models became very popular among experimentalists (Dayan and Abbott 2001; Kandel et al. 1991).

These models are called deterministic in the sense that once specific initial conditions are chosen, the time integration of a set of differential equations (that model the specific membrane currents and other intracellular processes such as calcium buffering) determines a unique solution that represent the evolution of the membrane potential.

However, it is well known that the intrinsic electrical behavior of many living neurons is very far from deterministic. Experiments where cortical neurons were isolated and repeatedly presented with the same artificially generated current time series (Mainen and Sejnowski 1995) have shown that the spike timing reliability of the neurons depends on statistical properties of the input signal in disagreement with deterministic predictions.

Even in circuits specialized to provide periodic and reliable patterns to control repetitive activity of muscles, such as the well studied pyloric central pattern generator (CPG) of the crustacean stomatogastric ganglion (STG) (Mulloney and Selverston 1974; Selverston and Moulins 1986), most of the motor neurons of the circuit present variable and irregular bursting behavior when synaptically isolated (Elson et al. 1999; Rabinovich et al. 1997; Selverston et al. 2000).

Every time one of these neurons is isolated from the pyloric CPG, it reveals a similar repertoire of irregular behaviors, but it is known that there are, at least, small changes in specific conductances from neuron to neuron and from animal to animal. In a mathematical sense, this means that if it is possible to write down equations to model the experimental behavior, these equations can never be specified exactly, yet similar behavior is expected from these equations when the parameters are allowed to vary. This is exactly the idea of structural stability from dynamical systems theory (Abraham and Shaw 1992).

So the irregular behavior, found to be related to nonlinear and chaotic properties of the cells (Falcke et al. 2000; Rabinovich et al. 1997), is structurally stable in a high-dimensional experimental parameter space, where the neurons are subject to small fluctuations not only in intrinsic properties but also in several other external parameters, such as temperature, ionic concentrations of the saline, etc.

The structural stability of the irregular behavior is a serious obstacle to developing deterministic Hodgkin-Huxley-type (HH-type) conductance models. In simple HH deterministic models of the axonal dynamics, the chaotic behavior was found to be structurally unstable because it was confined to a tiny volume of the parameter space (Guckenheimer and Oliva 2002). In more realistic models, where a dozen nonlinear equations are coupled, it could be intuitively expected that the higher complexity would enlarge the regions in the parameter space that correspond to irregular behaviors, providing structural stability to them. Surprisingly, this is not the case.

Recently, a very popular deterministic HH-type model for cultured stomatogastric neurons (Turrigiano et al. 1995), which is able to well reproduce tonic spiking and several regimes of bursting activity, was extensively studied (Prinz et al. 2003a) and a database of ∼1.7 million different neurons (each 1 corresponding to a particular set of conductance values in a 8-dimensional conductance parameter space) was built to characterize the possible types of dynamical behavior. In the entire database, only 0.1% of the neurons corresponded to truly irregular behaviors (not related to spike adding transitions in periodic bursting behavior). Moreover, this small percentage of irregular neurons are not condensed in a specific region of the parameter space, instead, they are spread over the parameter space (Prinz et al. 2004), suggesting the irregular behavior is not structurally stable in this model.

Structurally stable complex or chaotic behaviors were only reported in a few deterministic models (Falcke et al. 2000; Komendantov and Kononenko 1996) that included a nonlinear model of the dynamics of calcium exchange between endoplasmic reticulum stores and the cytosol. The main appeal of such models is that the modulation of [Ca^{+2}] is a major way to regulate many cellular processes (Ikeda 2004). These models predict that the cell must present slow macroscopic oscillations of the cytosolic levels of calcium. Such [Ca^{+2}] oscillations were observed experimentally in some vertebrates neurons (Parri and Crunelli 2001; Zhang et al. 2003) but not in STG cells (Levi et al. 2003).

The reasons why realistic deterministic models fail to reproduce the variability of the living behavior are the simplifications introduced by the voltage-clamp method and by the modeling. When a neuron is voltage-clamped, the experimenter has no access to the real spatial distribution of its membrane potential but only to a single point measurement that is taken as a good approximation for the membrane potential all over the neuron. Also, the discrete probabilistic character of the opening/closing of the membrane ion channels (Hille 2001; Levitan and Kaczmarek 1997) is replaced by average conductance activations/inactivations. These two simplifications make it possible to write much simpler differential equations, more amenable to analysis or simulations than complex multicompartment models with probabilistic channels. Moreover, because most of the experimental efforts made to characterize the currents present in STG neurons were made in cultured cells submitted to whole cell voltage clamp, there is not enough data to decide, without any speculative assumptions, neither how many compartments are needed, nor how to segregate or how to split the many different conductances among compartments.

However, some of the lost properties can be recovered even in a still limited single-compartment model by using a stochastic model based on Markovian states of ion channels (Desthexhe et al. 1994) where the probabilistic nature of channel operation is an intrinsic property (Hille 2001; Skaugen and Walloe 1979; White et al. 1998, 2000) and no additive noise needs to be included. If we consider a simple Markovian model of a population of channels of a specific type with only two possible and independent states (closed or open), the average number of channels (*N*) that are expected to suffer a transition from closed to open state in a given time interval can be easily obtained from the HH opening rate. The fluctuation (SD) of the number of channels that eventually undergo the transition is of the order of . Considering *N* to be ∼10^{6} (for axonal sodium channels) the fluctuation represents only 0.1% of *N*, and the behavior of the stochastic sodium current is practically the same as the one presented by the deterministic model. However, if one considers that action potentials are generated in the spike initiation zone, a much smaller region than the whole axon, with just ∼10^{3} to 10^{4} sodium channels and that the spikes are then propagated to other regions of the cell membrane, the fluctuation of the number of channels suffering the transition from closed to open can be 2% of *N*, enough to produce macroscopic effects.

Based on this mechanism of enhancement of the effects of the fluctuations of smaller numbers of ion channels, a stochastic model of an axon membrane patch (Schneidman et al. 1998) was compared with its deterministic version, and it was shown that the stochastic model is able to reproduce the same spike time reliability properties found in living neurons (Mainen and Sejnowski 1995) where the deterministic model fails.

Inspired by the success of the stochastic axon patch and by the fact that in stomatogastric neurons the slow ionic conductances, responsible for the starting and stopping of a burst of action potentials, are several orders of magnitude smaller than the conductances for axonal sodium and potassium (Prinz et al. 2003b), we present a single compartment stochastic HH-type model for a whole stomatogastric neuron.

The stochastic model presented here was obtained by translating to stochastic an up-to-date deterministic model for a stomatogastric neuron (Prinz et al. 2003b; Turrigiano et al. 1995). Our stochastic algorithm was adapted from a method developed for simulation of chemical reactions (Gillespie 1977) and later implemented for nonlinear membrane conductances (Skaugen and Walløe 1979), and although computationally intensive, it is considered an exact algorithm (Mino et al. 2002).

The translation of the deterministic model was carefully done in order not to change or include any new dynamics or ad hoc noise sources that were not present in the original model, so rather than comparing our results to the behavior of real living stomatogastric neurons, we can directly compare them to the results of the deterministic model when using the same parameters (maximum conductances of the different ion channel types), and the nonperiodic dynamics can be fully attributed to the role of the stochastic nature of the channels' behavior in generating macroscopic fluctuations of the membrane potential.

The present paper deals with data analysis using several tools from nonlinear dynamics to characterize and compare the results between the models and the ones reported in experiments. Bifurcation diagrams, return maps of interspike intervals (ISIs), reconstruction of dynamic attractors in phase-space, and Lyapunov exponents are used for identification and analysis of several types of neuronal activity (Abarbanel 1996; Hegger et al. 1999).

The beauty of the model presented here is that it allows the recovery of the irregularities found in real neurons by building up chaotic like global events from local stochastic dynamics that were directly translated from a realistic single-compartment deterministic model that does not itself exhibit any structurally stable chaotic behavior. Preliminary reports on the progress of our findings have been previously published in abstract form (Carelli and Pinto 2004; Carelli et al. 2004).

## METHODS

### Deterministic model

The deterministic HH-type model used in this paper is a single-compartment version based on voltage-clamp experiments of cultured STG neurons (Prinz et al. 2003b; Turrigiano et al. 1995). With appropriate tuning of the ion conductances this model is able to mimic very well both periodic bursting and tonic firing patterns. The behavior of different cells of the pyloric CPG as well as general properties of other bursting neurons can be, in principle, mimicked by fine tuning the maximal conductances of the model.

We used a sixth-order Runge-Kutta algorithm, with a fixed time step of 10 μs, implemented in a program written in C++ language on an Athlon XP-2400 PC-compatible computer to integrate the set of nonlinear differential equations of the model.

The membrane potential dynamics is given by where *V* is the membrane potential, *C*_{m} is the membrane capacitance, *I*_{ext} is an external current that can be either used to tune the model to a desired behavior or to simulate synaptic input, *V*_{i} is the reversal potential of the ion of the conductance of type *i* calculated from the Nernst equation. *g _{i}* is the instantaneous value of the conductance of type

*i*, expressed by where ḡ

*is the maximum conductance of type*

_{i}*i*,

*m*and

_{i}*h*are the HH activation and inactivation variables that assume values in the interval [0,1], γ

_{i}_{i}and δ

_{i}are integers chosen to fit the experimental data.

The model includes two types of spike generating conductances: sodium (*i =* Na) and rectifier potassium (*i = K*_{d}), as well as five types of slower conductances: *Ca*_{T} (transient calcium), *Ca*_{S} (slow calcium), *K*_{[Ca]} (calcium-dependent potassium), *A* (adaptation current), and *H* (hyperpolarization-activated current).

The dynamics of *m _{i}* and

*h*(for all conductances but

_{i}*K*

_{[Ca]}) is determined by where

*m*

_{i,}_{∞}, τ

_{mi},

*h*

_{i,}_{∞}, τ

_{hi}are the sigmoidal functions (shown in Table 1) inferred from voltage-clamp experimental data. The activation variable of the calcium-dependent potassium conductance, different from the others that depend only on the membrane potential, also depends on the calcium concentration.

There is a simple calcium buffer based on the competition between the mechanism of Ca^{2+} extrusion and the Ca^{2+} influx through Ca_{T} and Ca_{S} channels, expressed by where *f*, [Ca^{2+}]_{0}, and τ_{Ca} are constants. More details of the deterministic model can be found in Prinz et al. (2003b) and Turrigiano et al. (1995).

The parameters used in our simulations were *1*) *C*_{m} = 1.0 μF/cm^{2}; *2*) maximal conductances (ḡ* _{i}*) in mS/cm

^{2}: 200 for Na

^{+}, 100 for

*K*

_{d}, 2.5 for

*Ca*

_{T}, 4.0 for

*Ca*

_{S}, 5.0 for

*K*

_{[Ca]}, 50 for

*A,*0.01 for

*H,*0.01 for the leak current;

*3*) reversal potentials in mV: +50 for Na

^{+}, −80 for

*K*

_{d},

*K*

_{[Ca]}, and

*A,*−50 for the leak current, and −20 for

*H.*Since the intracellular concentration of calcium is not constant,

*V*

_{Ca}is computed as a function of [Ca

^{2+}] considering an external calcium concentration of 3,000 μM, according to the Nernst equation And

*5*) calcium buffer: τ

_{Ca}= 200 ms,

*f*= 14.96 μM/nA, [Ca]

_{0}

*=*0.5 μM.

The initial conditions used in the simulations were: *V*_{m} = −55.0 mV, [Ca^{2+}] = 0.5 μM*, m*_{Na} = 0.9. All other activation and inactivation variables were initialized to 0.1, and an initial transient of 10 s was disregarded in all simulations to ensure stability of the time series.

### Stochastic version

The stochastic version of the model was also implemented in a C++ written program and ran on the same computer used for the simulations of the deterministic model. An Euler algorithm with a time step of 10 μs was used to integrate the differential equations for *V* and [Ca^{2+}]. Although these equations were kept unchanged from the deterministic model, *g _{i}* now is expressed as where

*n*(

_{i}*A*

_{m}*, V, t*) is the total number of channels of type i in the open and activated state at a given time for a neuronal membrane with surface area

*A*, and φ is the conductance of a single ion channel, here for simplicity considered to be 20 pS for all channel types. The total number of channels of type

_{m}*i*,

*N*is thus obtained from the maximal conductances of the deterministic model as the closest integer to The average opening and closing rates, and are obtained directly from the sigmoidal expressions given in Table 1, according to The same method is used to obtain and , the rates of activation and inactivation.

_{i}A microscopic Markovian kinetic description is considered for representing the states that each ion channel can assume. Based on the current interpretation given to the exponents γ_{i} and δ_{i}, empirically inferred by HH, each channel of type *i* is modeled as a microscopic protein formed by γ_{i} gating subunits and δ_{i} inactivating subunits. All subunits in a channel act independently, according to the probabilities of transition given by their α's and β's. To allow the flow of current, an ion channel must have all its γ_{i} gating subunits in the open state and all its δ_{i} inactivating subunits in the activated state. *K*_{[Ca]} channels have more complex microscopic structures and more than one open state because of the [Ca^{2+}] dependence of the channels, but here we are not considering the effects of the stochastic binding of ions Ca^{2+} to the *K*_{[Ca]} channel binding sites. To keep the fidelity to the deterministic version the [Ca^{2+}] dependence is only included in the evaluation of .

According to this formalism, the Markovian kinetic scheme of states for *H* channels is the simplest one [*m _{i}*] represents the number of

*H*channels in the state

*m*and the open state is

_{i}*m*

_{1}. The total conductance provided by all open

*H*channels is given by

*g*(

_{H}*V, t*) = [

*m*](

_{i}*V, t*)φ.

The kinetic scheme for *K*_{d} and *K*_{[Ca]} channels is The open state is *m*_{4} and *g _{i}*(

*V, t*) = [

*m*

_{4}](

*V, t*)φ.

Because the Na*, Ca*_{T}, *Ca*_{S}, and *A* channels also have inactivation dynamics, a more complex scheme was used Here the open state is *m*_{3}*h*_{1} and *g _{i}*(

*V, t*) = [

*m*

_{3}

*h*

_{1}](

*V, t*)φ.

The key feature of the translation of the model consists in recovering the stochastic character of ion channel operation by giving a probabilistic interpretation to the original HH rates of transition between open and closed states (Schneidman et al. 1998) that is reflected in the number of open and activated channels *n _{i}*(

*V, t*). If we know the population

*N*of channels in state

_{j}*j*and the transition rate

*r*from state

_{jk}*j*to state

*k*, then in the time interval Δ

*t*, the mean number of channels expected to suffer this transition is Δ

*N*

_{jk}=

*N*Δ

_{j}*t r*. Because the probability of transition for one channel is

_{jk}*P*

*=*Δ

*t r*, we can find the number Δ

_{jk}*N*of channels that jump from the state

_{jk}*j*to the state

*k*by choosing a random number from the binomial distribution of probabilities

The main drawback of the stochastic model is that the binomial distribution requires intensive computation and the simulation of long time series can become infeasible. To surmount this constraint, we used an approximation for the exact binomial distribution. When the mean of the binomial distribution is bigger than five, the distribution has a symmetric shape, and a Gaussian approximation is possible. However when 0 ≤ mean ≤ 5, the binomial is asymmetric, and the Poisson distribution becomes a good approximation, as can be seen in Fig. 1. Our simulation program automatically detects and uses the best approximation for the current mean value of the distribution.

All values of parameters used in the simulations of the stochastic version are either the same as in the deterministic model or obtained from them.

The initial conditions used in the simulations were: *V _{m}* = −55.0 mV, [Ca

^{2+}] = 0.5 μM. The total number of channels of each channel type was obtained using the conductances of the deterministic model and considering a membrane area of 6.28 × 10

^{−4}cm

^{2}, a typical value for a stomatogastric neuron. The total number of channels of each type (

*N*) and the initial number of channels in each population were chosen as follows: Na

_{i}^{+}:

*N*

_{Na}

^{+}= 6.28 × 10

^{6}, initial number of channels in each one of the eight possible states =

*N*

_{Na}

^{+}/8;

*K*

_{d}:

*N*

_{Kd}= 3.14 × 10

^{6}, initial number of channels in each one of the five possible states =

*N*

_{Kd}/5;

*Ca*

_{T}:

*N*

_{CaT}= 7.85 × 10

^{4}, initial [

*m*

_{0}

*h*

_{0}] = initial [

*m*

_{0}

*h*

_{1}] =

*N*

_{CaT}/2, all other initial [

*m*

_{i}*h*] = 0;

_{j}*Ca*

_{S}:

*N*

_{CaS}= 1.25 × 10

^{5}, initial [

*m*

_{0}

*h*

_{0}] = initial [

*m*

_{0}

*h*

_{1}] =

*N*

_{CaS}/2, all other initial [

*m*

_{i}*h*] = 0;

_{j}*K*

_{[Ca]}:

*N*

_{K[Ca]}= initial [

*m*

_{0}] = 1.57 × 10

^{5}, all other initial [

*m*] = 0;

_{i}*A*:

*N*= 1.57 × 10

_{A}^{6}, initial [

*m*

_{0}

*h*

_{0}] = initial [m

_{0}h

_{1}] =

*N*/2, all other initial [

_{A}*m*

_{i}*h*] = 0;

_{j}*H*:

*N*= [

_{H}*m*

_{0}] = 314, all other initial [

*m*] = 0.

_{i}Using these initial conditions, the simulation of the stochastic model neuron starts with a spike that is followed by a fast convergence of the ionic populations to dynamical “steady-state” values.

When not specified in the text, an initial transient of 10 s was disregarded.

### Mixed (stochastic/deterministic) model

A third program in which it is possible to choose each individual current to be stochastic or deterministic was also written. We used this program to test the role of each ion channel type in producing irregular behavior. When a current is chosen to be stochastic, the simulation method described in *Stochastic version* is used for that current, and when the current is deterministic the method described in *Deterministic model* is used, however, in this case, the program integrates the activation/inactivation variables using an Euler method instead of the Runge-Kutta. The use of different integration methods, in this case, do not impact the results shown. We performed several tests using both Euler and Runge-Kutta methods for the deterministic model, and it behaved in a perfectly periodic fashion in both cases.

### Stomatogastric isolated neuron

To allow a direct comparison between the results of the simulations and the behavior of a real neuron, the membrane potential of an isolated lateral pyloric (LP) neuron from the stomatogastric ganglion of the spiny lobster *Panulirus interruptus* was recorded. The nervous system was dissected from the animal, the ganglion was prepared for intracellular recordings, and the cells were identified as described in Mulloney and Selverston 1974. The synaptic isolation of the LP neuron was done by using standard pharmacological and cell-killing techniques (Bal et al. 1988; Miller and Selverston 1979). In short, glutamatergic synapses were blocked by using lobster saline with 10 μM picrotoxin, and the two pyloric dilators (PD) and the ventral dilator (VD) neurons were photoinactivated. The LP was impaled with two electrodes for separated membrane potential recording and current injection.

### Data analysis

Data analysis was performed in the original simulated time series as well as in the interspike intervals (ISIs) sequences obtained from them. ISIs were detected in the time series using a simple threshold method implemented in a home-made C++ program: every time the membrane potential crossed a 0.0-mV threshold level with positive derivative a spike was attributed to the model neuron. Nonlinear analysis methods were applied to the original time series and in the ISI series: we used the method of false nearest neighbors (Abarbanel 1996; Kennel et al. 1992) to calculate the embedding dimension of the dynamical attractor reconstructed from the simulated data and calculated the positive Lyapunov exponent (Sano and Sawada 1985), which is a measure of the presence of sensitivity to initial conditions and chaos in a system. We used algorithms from the package TISEAN (Hegger et al. 1999) in our nonlinear analysis.

## RESULTS AND ANALYSIS

A first comparison of the dynamical behavior of the deterministic and the stochastic model is shown in Fig. 2 where the time series obtained from both models for the same values of parameters are plotted.

As can be seen in Fig. 2, irregular activity is obtained in the stochastic model, where the original deterministic model presented only regular oscillations. However, the activity of the stochastic model still resembles the main characteristics of the deterministic model as shown in Table 2.

There is a tendency of the stochastic model to present shorter bursts and smaller number of spikes/burst than the deterministic model. This fact was found to be related to the behavior of the Ca^{2+} and *K*_{[Ca]} stochastic conductances as will be further discussed later.

As a first attempt to identify how much determinism there is in the time series, we followed the time evolution of the membrane potential when the model neurons start from close initial conditions. After a hyperpolarization, the conductances of the neuron are reset to similar values at every burst, so we detected the first spike of the neurons after the hyperpolarizations and used these first spikes to trigger a plot of superimposed bursts. In Fig. 3 we plotted 50 superimposed bursts for each model. The bursts of the deterministic model are all identical, so the superimposed bursts appear just as a single curve that repeats over and over. However, in the stochastic model, only the first spikes overlap, then the precision of the spiking time decreases with the higher order of the spike in the burst and the last spike is very unreliable, consequently the burst end is highly variable. The picture for the stochastic model clearly reproduces the behavior of isolated biological neurons and suggests the small differences present in the afterhyperpolarization are amplified in the stochastic model as is the case in a nonlinear system with sensitivity to the initial conditions. Further experimental data for comparison can be found in Selverston et al. (2000).

Another way used to analyze the behavior of the models during the bursting activity was to detect their ISIs and to plot their ISI return maps (ISI_{n+1} vs. ISI_{n}) as shown in Fig. 4. Because these maps associate the value of an ISI with the previous one, they are a straight way to detect determinism in the ISI patterns.

The ISI patterns are precise for all ISIs of the deterministic model, but only the first ISIs repeat in the bursts presented by the stochastic model, as we can see in the small clouds of scattered dots produced by the first ISIs. This initial and more precise behavior of the ISIs of the stochastic model closely resembles what happens in the deterministic model and is the signature of determinism that is still present in the stochastic model, even though it was built from stochastic elements. Again, only the stochastic model reproduces the behavior found in ISI return maps of STG neurons activity. More detailed analysis of ISI return maps from STG neurons can be seen in Szücs et al. (2003).

We have also studied and compared the behavior of the models when a depolarizing DC current is injected in the neuron. The results obtained with the models are shown in Fig. 5 can also be directly compared with the behavior of a living LP neuron. The stochastic model clearly is a better approximation of the experimental behavior than the deterministic model. However, there is still some similarity between the models, such as the change of the behavior from bursting to tonic firing as the current is increased.

To characterize the way the transition from bursting to tonic firing takes place in both models, we built bifurcation diagrams (Fig. 6) where several ISIs were plotted for each value of injected current. The only irregularities found in the deterministic model (Fig. 6, *left*) are due to the intermittent switching between bursts with *N* spikes and bursts with *n* + 1 spikes, which is characteristic of a period adding bifurcation. For *I*_{ext} ≈ 0.16 nA, there is a sudden transition between bursting behavior and tonic spiking. In the stochastic model, the period adding bifurcation is not clearly present and the sudden transition does not take place, instead the interburst ISI intervals get shorter and mix with the intraburst ISIs.

To provide an easy way to compare the results of our models with other results found in the literature (such as the Figs. 6 and 9 of Rowat and Elson 2004), we measured the spike time in burst (STB) of the time series of both models and plotted the bifurcation diagram of the STBs as a function of the external current, as shown in Fig. 7. The STBs are measured using as reference the first spike of the burst. When the model neuron is tonic, the STB simply is taken equal to the ISI. In the deterministic model (Fig. 7, *left*) a period adding bifurcation (for 0.10 < *I*_{ext} <0.16 nA) is observed. This means that the burst duration is increased by the addition of an extra spike at the end of each burst. Finally, there is a sudden transition when the burst length tends to infinity (*I*_{ext} ≈ 0.16 nA), the model stops bursting and a tonic firing periodic regime takes place. In the stochastic model (Fig. 7, *right*), however, the transition from bursting to tonic firing is not so clear. For a range of current values in which the deterministic model is already tonic (0.17 < *I*_{ext} <0.28 nA), the stochastic model switches between irregular tonic firing and bursting behavior.

Some works in the literature (Falcke et al. 2000; Selverston et al. 2000) indicate that the LP neuron presents statistical dynamical properties compatible with deterministic chaos. This analysis assumes that the neurons have deterministic dynamical evolution rules. We applied the same analysis to our stochastic model time series, as if they were experimentally obtained, to show that our results are compatible with those from real neurons.

The phase space of the time series was reconstructed by using delayed coordinates (Abarbanel 1996), and the embedding dimension (the minimum number of dynamical variables needed to fully detangle the dynamical attractor) was estimated to be three by the method of false nearest neighbors. This is in agreement with experimental results (Falcke et al. 2000; Selverston et al. 2000).

The Lyapunov spectrum of the stochastic model time series, calculated by using the package TISEAN, has one positive exponent (λ_{1} *=* 0.050 ms^{−1}), one exponent that lies very close to zero (λ_{2} *=* −0.007 ms^{−1}) and a negative exponent (λ_{3} *=* −0.070 ms^{−1}). The sum of all exponents is negative, as expected for a dissipative system where the trajectories are confined in a finite region of the phase space. From the Lyapunov exponents, the Kaplan-Yorke fractal dimension was estimated to be 2.6. Considering that we are using a single-compartment model, with large stereotyped spikes and without deep hyperpolarizations, these results are compatible to the ones found in Falcke et al. (2000).

Because most of the works where the stochastic nature of the ion channels was considered relied on the effects of the fluctuations produced by a small population of channels (Rowat and Elson 2004; Schneidman et al. 1998; Skaugen and Walloe 1979; White et al. 1998), there is another important question about the irregularities reproduced by our stochastic model: How robust are these irregularities to changes in the number of channels considered in the model (surface of the neuronal membrane). To address this question, we ran the stochastic model for many different sizes of the membrane and plotted the ISI bifurcation diagram as a function of the membrane area, as shown in Fig. 8. As the area is increased, the ISIs of the stochastic model become more reliable and the model behavior tends to be more periodic. In the limit of the area going to infinity, the periodic behavior of the deterministic model is recovered, as one can see by comparing the scattered ISIs of the stochastic model with the deterministic model ISIs (which are area independent) shown in the right of the figure as small horizontal segments. The stochastic nature of ion channels is indeed important for the behavior of an STG neuron, as we can see in the irregular behavior pointed to by the arrow at the left of Fig. 8.

To identify which are the conductances that play more important roles in generating the irregularities showed by the stochastic model we implemented mixed models, as explained in *Mixed (stochastic/deterministic) model,* where we could choose each individual current either to be stochastic or deterministic. The results of several mixed models are resumed in Fig. 9 where we plotted statistical properties of models where each one of the currents was made stochastic alone and the other currents were kept deterministic. We also compare the results with those obtained when all currents are stochastic and when all the currents but Na^{+} and K_{d} are stochastic.

The main currents responsible for the variability found in the stochastic model are: *I*_{CaT}, *I*_{CaS}, *I*_{K(Ca)}, known to be the “slow wave” or “plateau” currents in living neurons and *I*_{A}, the adaptation current. This is not surprising because the two calcium currents have a special relation with their competing *K*_{[Ca]} current: they change the cytosolic [Ca^{2+}] which, in turn, regulates the gating of the *K*_{[Ca]} channels. The key feature here is that there are much fewer of these channels than the spike generating Na^{+} and K_{d} channels, so their relative statistical fluctuations are more important than the fluctuations of spiking channels in determining when a burst is about to start, in increasing the activity of the neuron, or in ending the burst. The reason why the *H* current (which has the smallest number of channels) doesn't play an important role in the irregularities present in the stochastic model is that this current is practically not activated in this isolated model except when neuron is strongly hyperpolarized. Our finding is that the fluctuations in the number of spiking channels is not important for the macroscopic irregular behavior: the spiking channels tend to be modulated by the slower irregular oscillations produced by the slow wave channels (present in smaller number in the membrane), as can be seen in Fig. 9 in the very similar behavior presented by the models built either with all currents stochastic or with only the “slow” currents stochastic.

The end of the burst is determined by the competition between the calcium influx and the outflow of potassium through *K*_{[Ca]} channels. In the beginning of a burst, the calcium conductances are very high, so the stochastic fluctuations do not influence the behavior. However, from the middle to the end of a burst, as shown in Fig. 10, the number of conducting channels in both populations is small (∼10^{3} open channels), making the stochastic fluctuations more important, which causes an increase in the probability of ending the burst before the time deterministically predicted.

### Conclusion

We presented a single-compartment whole cell stochastic model for a stomatogastric neuron based on a deterministic HH model built on realistic measures of conductances from voltage-clamp experiments.

The main question we addressed was: because the behavior of stomatogastric neurons has been studied and known for more than three decades and there is a lot of experimental data from voltage-clamp experiments as well as very good and complex detailed models built on these experimental data, why do the models not present the same kind of irregularities easily found in isolated STG neurons?

Several approaches can be used to build new models: one can include new dynamics in the existing conductances, or a myriad of new complex ion currents, or mechanisms such as the calcium exchange with internal stores. However, we opted for simply using an existing realistic deterministic model and to include in this model the ingredient we believe is essential to reproduce the experimental irregularities: the stochastic nature of the ion channels operation.

The single-compartment deterministic model used as a basis of our work had the traditional spike-generating conductances (Na^{+} and *K*_{d}), several slower conductances (*Ca*_{T}, *Ca*_{S,} *K*_{[Ca]}, *H, A*), a leak current, and a simple Ca^{+2} buffer to account for the competition between the Ca^{+2} influx and the mechanisms of Ca^{+2} extrusion, most inferred from detailed voltage-clamp experiments.

We made the translation of the deterministic model to stochastic in the simplest possible way: the numbers of channels of each type were determined from the deterministic maximum conductances, the channels were modeled as a chain of Markovian states with transition probabilities given directly by the HH voltage-dependent transition rates, which constitute one of the most exact stochastic algorithms for describing the microscopic fluctuations the ion channels are subject to. No external noise was added to any gating variable, and all the stochasticity came directly from the channel transition dynamics. The Ca^{+2} buffer was kept deterministic to ensure that the irregularities obtained from the model are not produced by any other mechanism not present in the original deterministic model.

We found that our stochastic model is able to reproduce the irregular behavior presented by living neurons. Using tools from nonlinear dynamics, we showed that it also presents sensitivity to initial conditions. However, it is well known that the original deterministic model does not present structurally stable deterministic chaos for any set of parameters. So, the chaotic behavior must arise in the stochastic model as the interplay between the stochastic nature of the channels and the nonlinear dynamics already present in the deterministic version. We interpret this emergence of chaotic behavior as the effect of the high dimension of the real state space. The oscillations in the deterministic model lie in a 12-dimensional state space (12 dynamical variables). However, because of the strong attractive periodic orbit, the steady-state behavior can be reduced to a lower dimensional cycle. In the vicinity of this cycle, there are regions where the local dynamics present sensitivity to initial conditions, but these regions are never visited by the periodic orbit. When the stochastic character of ion channels is included, the steady-state solution is perturbed in such a way that the resulting orbit still remains close to the attractive periodic orbit, but now the regions with richer dynamics are accessible, revealing the sensitivity to initial conditions.

Mixed models, where some of the conductances are deterministic and some are stochastic, were used to track which conductances are more important in generating the irregular macroscopic behavior. Using this approach, we found that the main conductances responsible for generating irregularities are the Ca^{+2} and the *K*_{[Ca]} conductances.

Earlier stochastic models used to consider a patch of the membrane where a smaller number of channels could produce relative fluctuations big enough to impact the overall neural behavior. However, here we found that the effective number of conducting channels of a type in a given time is several orders of magnitude smaller than the total number of channels of that type, and the effects of the stochastic nature of channels can be found even in a single-compartment whole cell model, without the need of implementing more compartments or axon patches.

## GRANTS

The authors acknowledge the financial support from the Brazilian agencies Fundação de Amparo à Pesquisa do Estado de São Paulo and Conselho Nacional de Desenvolvimento Científico e Tecnológico.

## Acknowledgments

We thank I. Segev and M.A.L. Nicolelis for encouragement during our poster presentation in the first Neuroscience Symposium of the International Institute for Neurosciences of Natal, Brazil.The authors also acknowledge the hospitality and the help of M. I. Rabinovich and A. I. Selverston, University of California, San Diego, where the experimental data shown in this paper was collected.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 by the American Physiological Society