Presence of a Chaotic Region at the Sleep-Wake Transition in a Simplified Thalamocortical Circuit Model

Sleep and wakefulness are characterized by distinct states of thalamocortical network oscillations. The complex interplay of ionic conductances within the thalamo-reticular-cortical network give rise to these multiple modes of activity and a rapid transition exists between these modes. To better understand this transition, we constructed a simplified computational model based on physiological recordings and physiologically realistic parameters of a three-neuron network containing a thalamocortical cell, a thalamic reticular neuron, and a corticothalamic cell. The network can assume multiple states of oscillatory activity, resembling sleep, wakefulness, and the transition between these two. We found that during the transition period, but not during other states, thalamic and cortical neurons displayed chaotic dynamics, based on the presence of strange attractors, estimation of positive Lyapunov exponents and the presence of a fractal dimension in the spike trains. These dynamics were quantitatively dependent on certain features of the network, such as the presence of corticothalamic feedback and the strength of inhibition between the thalamic reticular nucleus and thalamocortical neurons. These data suggest that chaotic dynamics facilitate a rapid transition between sleep and wakefulness and produce a series of experimentally testable predictions to further investigate the events occurring during the sleep-wake transition period.


INTRODUCTION
The specific anatomy of thalamocortical circuitry along with the interactions of ionic conductances in thalamic and cortical cells give rise to multiple modes of activity that characterize behavioral states such as the sleep-wake cycle and generalized epilepsy. During drowsiness and sleep, large groups of neurons fire synchronously giving rise to oscillatory activity such as spindle (6-14 Hz), delta (1-4 Hz), and slow (<1 Hz) oscillations, all of which may be observed in the constituent cells of the thalamus and cortex; thalamocortical (TC) cells, neurons in the thalamic reticular nucleus (RE) and deep-layer cortical neurons (CX; Steriade, 2005;Huguenard and McCormick, 2007). TC neurons project primarily to layer 4 of cortex, and also send projections to other layers, including infragranular layers Meyer et al., 2010). Corticothalamic projections, which arise from layers 5 to 6, terminate on TC neurons as well as on neurons of the RE (Liu and Jones, 1999;Zikopoulos and Barbas, 2006;Llano and Sherman, 2008). TRN neurons, which receive excitatory input from TC and CX neurons, send GABAergic projections to TC neurons (Pinault, 2004). Interplay of intrinsic ionic properties of single neurons from these structures, modulated and orchestrated by the synaptic interactions of large groups of neurons, likely play a large role in generating coherent oscillatory activity in the brain.
One cellular mechanism which has been postulated to play an important role for the generation of large-scale oscillatory activity is post-inhibitory rebound (PIR) present in TC and RE neurons, as well as a subset of CX neurons (Stafstrom et al., 1984;Steriade et al., 1993b). This mechanism is associated with low-threshold Ca 2+ (I T ) as well as hyperpolarization-activated cation currents (I H ) through hyperpolarization-activated cyclic nucleotide-gated (HCN) channels (Lüthi and McCormick, 1998;Sherman, 2001;Timofeev et al., 2002;Kramer et al., 2008). For the calcium component, after TC neurons have been hyperpolarized for a protracted period of time (∼50-100 ms), I T becomes de-inactivated, releasing the low-threshold calcium conductance for activation (Jahnsen and Llinas, 1984). Elimination of I T greatly disrupts spindle and absence seizure activity (Kim et al., 2001;Astori et al., 2011), suggesting that this current is critical for producing synchronized oscillations in the thalamus. For I H , after being activated by hyperpolarization, the return to baseline voltage is associated with delayed closure of HCN channels, leading to rebound depolarization (Pape, 1996;Slater et al., 2013). The combined effect of these two currents is typically a low-threshold slow rebound depolarization with a superimposed burst of traditional sodium action potentials (Lüthi and McCormick, 1998). In the context of a network containing TC, RE, and CX cells, this post-inhibitory rebound can be triggered after TRN-based GABAergic inhibition in TC cells, producing a TC burst, which would produce excitation in a RE cell. Such a cycle may be initiated by cortical activation of the thalamic reticular nucleus, and can generate oscillations associated with sleep spindles (Steriade, 2005) as well as paroxysmal activity (Meeren et al., 2002;Polack et al., 2007).
What is still unclear is the nature of the state transitions between sleep and wakefulness that are characterized by marked difference in thalamocortical oscillatory behavior. Slow-wave sleep is characterized rhythmic, delta-range bursting in TC cells (Amzica and Steriade, 1998), while the waking state is dominated by irregular isolated action potentials and occasional bursts (Ramcharan et al., 2000;Bezdudnaya et al., 2006). Previous work in relatively simple models has suggested that TC cells near the onset of rhythmic bursting show chaotic behavior (Wang, 1994;Paul et al., 1998); a type of dynamical interaction which may facilitate the rapid transition between discrete states in neural systems (van Vreeswijk and Sompolinsky, 1996). Therefore, the objective of the present study was to characterize state transitions using tools of nonlinear dynamics in a reduced but biologically realistic model of the thalamo-reticulo-cortical circuit that maintains the essential intrinsic conductances in these nuclei and their synaptic connectivity. We assume here that transitions from slow oscillatory behavior seen in model TC neurons to rapid, irregular firing corresponds to a transition between sleep and wakefulness.
The traditional linear spike train analyses treat variability in a biological system as noise and average out this factor. This approach does not take into account that this variability may be an inherent property of an interconnected dynamical system such as the thalamocortical system. Chaos refers to deterministic or non-random behavior which exhibits very rapid growth of errors leading to different system states and thus renders impossible any long-term prediction. Others have argued that the irregularity of spike patterns in large networks is an emergent property of the network that enables it to respond very rapidly to external input, much faster than the time constant of individual neurons, and therefore benefits neural processing (van Vreeswijk and Sompolinsky, 1996;Bertschinger and Natschläger, 2004;Canavier and Shepard, 2009;Sussillo and Abbott, 2009). Previous studies have shown that artificial neuronal networks and real networks of neurons can demonstrate the presence of chaotic regimes (Canavier et al., 1990;Wang, 1994;van Vreeswijk and Sompolinsky, 1996;Siegel and Read, 2001;Bertschinger and Natschläger, 2004;Battaglia et al., 2007;Sussillo and Abbott, 2009;Rajan et al., 2010;Jia et al., 2012). Other studies have demonstrated fractalness in spike trains (Teich, 1989;Teich et al., 1997;Darbin et al., 2006;Gebber et al., 2006). These studies suggest that chaotic dynamics exist in a large number of simulated and real neural systems.
In the present study, we have constructed a reduced 3 neuron thalamo-reticulo-cortical model consisting of one CX corticothalamic neuron, one TC neuron and one RE neuron with intrinsic conductances, based in part on our intracellular recordings, and biologically realistic synaptic connectivity between them. The model is a direct implementation of the circuitry that has been hypothesized to be essential for thalamocortical oscillatory behavior (Steriade et al., 1993a,b;Steriade, 2005;Contreras, 2014; see Figure 1). The implementation of post-inhibitory rebound and associated low-threshold Ca 2+ conductance was based on in vitro data from cortical slices and then fit to the model to obtain realistic rebound excitation in an infragranular CX neuron. In the interconnected model, we demonstrated periodic oscillatory and irregular aperiodic spike burst patterns depending on the strength of synaptic connectivity.
Finally, we demonstrate that within a range of synaptic conductance values representing putative physiological transition states, chaotic behavior in the network model is obtained.

Electrophysiology
All procedures were performed under the guidelines of the National Institutes of Health and approved by the Institutional Animal Care and Use Committee. Adult male rats (100-250) g were anesthetized with sodium pentobarbital and decapitated and the brain quickly removed and put in ice-cold artificial cerebrospinal fluid (ACSF) bubbled with 95% O 2 -5% CO 2 . The composition of ACSF was in mM NaCl (126), KCl (3), MgSO 4 (1), NaH 2 PO 4 (1.25), NaHCO 2 (26), Glucose (10), and CaCl 2 (2). Coronal slices (400 µm) were prepared FIGURE 1 | (A) Basic thalamo-cortical circuit consists of reciprocal excitatory connections between TC and CX neurons with collaterals to RE. RE neurons have inhibitory projections to TC (reprinted with permission, Steriade et al., 1993a). (B) Thalamo-Reticulo-Cortical network model developed for simulations in the current study. This computational model with anatomically realistic connections is an implementation of the proposed circuitry for the neural basis of thalamocortical oscillations. through the primary somatosensory area and placed in an interface type slice chamber with a continuous flow of ACSF maintained at a temperature of (37 ± 1) • C. Somatosensory cortex was used since oscillatory phenomena that accompany states of sleep and drowsiness involve large regions of the brain and are found throughout the brain and because spindle oscillations that accompany drowsiness and sleep are found in the somatosensory cortex (Contreras and Steriade, 1995;Khazipov et al., 2004;Rosanova and Timofeev, 2005;Halassa et al., 2011). Blind intracellular recordings employed sharp micropipettes (Brown Flaming Puller; 100-160 M ) filled with 1M K-acetate (pH 7.4). Membrane potentials were recorded with an Axoprobe-1A amplifier. The output was displayed on a digital storage oscilloscope and stored on a videocassette recorder via a Neuro-corder (Model DR-890, Neurodata Instruments) for offline analysis. Analysis was done using custom written software (Dataman PC, Cauller). Neurons in infragranular layers V/VI were classified as either non-adapting or adapting based upon their response to sustained somatic current injection.

Computational Model
A reduced cortical infragranular pyramidal cell model was developed in the GENESIS simulation environment (cortical infragranular cell, CX, in Figure 1). It consisted of four compartments: tuft, neck, soma and base. The tuft (325 µm length and 1 µm diameter) and neck (175 × 1 µm) together represent the apical dendrite. The neck was attached to the apical side of the soma (15 × 15 µm) whereas the base (225 × 2 µm) represented the basal dendrites. The membrane resistivity (R M ), axial resistivity (R A ) and membrane capacitivity (C M ) were set at 0.5 /m 2 , 1.0 /m, and 0.017 F/m 2 , respectively. The model consisted of the following currents: fast sodium (I Na ), potassium delayed rectifier (I K_DR ), cortical low-threshold/rebound (I CX_REB ), hyperpolarizationactivated cation (I H ), potassium after-hyperpolarization (I AHP ), calcium dependent potassium (I K [Ca] ), and the potassium Acurrent (I A ). The channel description for the I CX_REB current was implemented based on data directly obtained from the activation and inactivation curves obtained from in vitro experiments described above. All other channel descriptions were taken from existing GENESIS user libraries. The maximal conductance values used were (in S/m 2 ): g K_DR = 200, g CX_REB = 7, g H = 12.55, g AHP = 2.1, g K[Ca] = 0.28, and g A = 12. To determine rebound parameters in isolation, the fast Na + maximal conductance was set to zero (g Na = 0) to simulate the effect of the intracellular Na + blocker QX-314, as was done in physiological experiments.
The differential equation governing the membrane voltage of the model was given by C dV CX dt = −(I Na + I K DR + I CX REB + I H + I AHP + I A + I K[Ca 2+ ] ) Here, C is the capacitance of the membrane, V CX is the membrane voltage and I j is the current passing through each specific channel. g j is the maximal conductance value for each channel (Bower and Beeman, 1998, Chapters 4 and 6). m and h represent the instantaneous activation and inactivation variables. They are raised to the power M and N, respectively. E j is the reversal potential of each channel j. The steady state activation and inactivation variables m ∞ and h ∞ were either obtained from experimental data or from equations fitted to experimental data as described above. The instantaneous values of these gates change with respect to time as Where, τ m and τ h are the activation and inactivation time constants, α and β are the forward and backward rate constants, respectively (Bower and Beeman, 1998). From first-order kinetics, α and β are of the form Where, A-F are constants and V CX is the membrane potential of the cortical neuron. For Ca 2+ dependent processes, α and β are multiplied by a factor representing the intracellular Ca 2+ concentration. Change in the intracellular Ca 2+ concentration is described by Here, τ Ca is a factor representing the rate of decay of [Ca 2+ ]. B is a factor representing the flux of ions in a thin shell near the membrane surface produced by I Ca (Bower and Beeman, 1998, Chapter 19). The synaptic conductance change is modeled by the "alpha function" (Rall, 1967).
Where, g syn is the instantanteous synaptic conductance value of either of g GLU or g GABA . At t = t p , the function rises to a maximal conductance value of g max for each of the synaptic conductances and then decreases slowly to zero. The synaptic current is then given as Where V m is the membrane potential and E syn is the synaptic reversal potential of the synaptic conductance (i.e., E GLU or E GABA ). For the network model, the simplified 4 compartment model of an infragranular CX cell was synaptically connected to the thalamo-reticular (TC-RE) model (Paul et al., 1998). The TC cell model included I T , I H , I Na , I K , and I GABA and the RE model consisted of I TS , I K[Ca2+] , I CAN , I Na , I K , and I AMPA to implement the three neuron network TC-RE-CX model (Figure 1). This network model is an implementation of the Steriade (Steriade et al., 1993b,c;Steriade, 2005) hypothesis for the neural basis of oscillations in the thalamus and the cortex. The TC-RE network has reciprocal excitatory and inhibitory connections. The TC cell has excitatory AMPA/Kainate projections to the RE cell which, in turn, has GABAergic inhibitory projections to the TC. The CX cell is connected to this network with excitatory glutamate projections to both TC and to RE. The TC cell has reciprocal excitatory projections to the neck compartment of CX. This approximates the thalamocortical excitatory projections at the proximal apical dendrites (layer IV synapses). The maximal conductance values for the CX cell was slightly adjusted to obtain oscillatory behavior in the CX cell in concert with those in the TC and RE cells. These values are given in Table 1. The network model was constructed such that RE acted as the initiator or pacemaker cell, based on findings from lesion experiments in the cat (Steriade et al., 1987).

Nonlinear Dynamical Analysis
To identify chaotic regimes within the dynamic system four methods of analyses were used. These are elaborated below: Phase space analysis: Two time-dependent variables were plotted against each other over time. Such a plot provides a graphical picture of the dynamical system being investigated and differentiates between a fixed point attractor, a limit cycle and strange attractors. Presence of a strange attractor suggests the presence of a chaotic regime within the dynamical system. Bifurcation plots: A bifurcation plot illustrates the transition from periodic to aperiodic or chaotic states thus plotting the transition from order to chaos in a dynamical system. These are usually plotted with a system parameter on the x-axis and a representation of an attractor on the y-axis. At a bifurcation point, the attractor in the plot splits into two if the attractor changes from a period of one to a period of two. Fractal dimension: Two measures were used to compute the fractal dimension of neuronal spike train pattern. These are the Fano Factor and the Allen Factor. Fano Factor refers to the ratio of the variance in the number of spikes within a time window T to the mean number of spikes within that same time window. This is given as Experimentally, the Fano Factor is calculated by computing the above metric repeatedly for different time windows (T) over the range of interest. For a true Poission process, F(T) is equal to 1. F(T) > 1 imply that the variance in the number of events grows faster than the mean, suggesting the presence of long-term correlations in the data and fractal properties (Teich, 1989;Lowen and Teich, 1996;Koch, 1999). Due to the nature of computation of the Fano Factor, it cannot increase faster than ∼T 1 (Lowen and Teich, 1996). Therefore, a second measure, the Allan Factor, which has been used to detect selfsimilarity with fractal exponents greater than unity (Allan, 1966;Barnes and Allan, 1966), was used here. The Allan Factor is derived from the Allan variance which is given as the average variation in the difference of adjacent counts. The Allan Factor is then the ratio of the Allan variance of the even count to twice the mean. It is related to the Fano Factor by A(T) = 2F(T) -F(2T) where, F(T) is the Fano Factor for counting time T (Scharf et al., 1995). The Allan Factor can rise as fast as ∼T 3 and therefore, can be used to estimate fractal exponents in the range 0-3 (Lowen and Teich, 1996). Lyapunov exponents: The method for computing Lyapunov Exponents involved a direct approach which was possible within a simulation environment. This algorithm is graphically shown in Figure 2. Two identical networks of the TC-RE-CX reciprocally connected model were created and were run for a period of 15 s simulation time to eliminate transients. An iteration time step (dt) of 10 µs was used.
After 15 s of simulation run time, a perturbation of 0.005 mV introduced into the copy of the original system. The two parallel simulations were now run for 5 time steps to allow the initial perturbation introduced into the membrane voltage to filter through to all other parameters. The error E0 was measured at this time and the simulations were run for an evolve time of 0.1 s during which the introduced "error" was allowed to amplify (Ek). The parameter of interest was the somatic membrane potential which was recorded from both the original system and the perturbed copy and the log of the absolute difference ("error") was computed. All parameters of the perturbed system were reset to the values of the original system and this procedure was repeated for 1000 iterations. The average rate of divergence or convergence (i.e., Lyapunov exponent) was then computed using the formula: This procedure of perturb, run, collect data, reset and perturb again prevented the errors from becoming unbounded, therefore enabling the algorithm to detect any local "stretches" while avoiding the global "fold" (Paul et al., 1998).
For all methods (except Lyapunov Exponents), the analysis was done over 20 s of data after removal of the first 5 s in order to avoid the initial transients. The characterization of network activity as periodic or aperiodic was based on the spike pattern of the CX cell. The two parameters that were systematically varied while computing the Lyapunov exponent were the strength of TC→CX synaptic excitation and the strength of RE→TC synaptic inhibition. Two representative data sets are used for detailed analysis for each of the methods to distinguish between periodic state (as observed during sleep) and chaotic state (transition state). Additionally, the bifurcation plot and Lyapunov exponents were computed for a range of data points that demonstrate a transition from one mode to the other. Physiologically, transition from sleep to wakefulness was simulated by decreasing a potassium leak conductance (g LEAK ) to the TC and RE cells from its initial setting. This approach has been used previously in thalamic Willis et al., 2015) and cortical (Hill and Tononi, 2005) models, and simulates excitatory brainstem input (e.g., cholinergic) that accompanies the transition between these behavioral states and modulates the resting membrane potential of these cells. Lyapunov Exponents were computed for a range of g LEAK values.
FIGURE 2 | Algorithm for computing Lyapunov Exponents. Two identical networks of the TC-RE-CX reciprocally connected model were created and an error was introduced in the membrane potential of the copy system and allowed to amplify for a period of time (0.1 s) after which the difference between the original system and the copy system was measured and all parameters of the copy system was reset to the values of the original system. The average rate of divergence or convergence was then computed using the formula shown in the figure.

In vitro Determination of Cortical Neuron Post-inhibitory Rebound Parameters
We examined PIR properties mediated by low threshold Ca 2+ conductance (I T ) in real infragranular cortical neurons in brain slices. These neurons are hypothesized to contribute to the dynamic behavior of cortico-thalamic networks and PIR plays a critical role in network oscillations. The rebound depolarization properties obtained from in vitro experiments provide the activation and inactivation variables for the implementation of I T conductance in infragranular cortical neurons. PIR is mediated by a low threshold calcium conductance or low threshold calcium spike (LTS). The mean resting membrane potential of the rebound cells was −65.56 ± 4.27 mV. In the presence of intracellular blocker QX-314 (50 mM), the mean LTS amplitude 13.57 ± 4.05 mV, and the LTS width was 67.6 ± 21.55 ms. Figures 3Ai,ii demonstrates typical voltage responses to current stimuli in a neuron for each of the protocols used to obtain low-threshold conductance properties, namely, the activation parameter (m), and the inactivation parameter (h). The Ca 2+ dependence of the post-inhibitory rebound is shown by block and recovery with Co 2+ (Figure 3Aiii). The current stimulus protocols used for obtaining the voltage responses for activation and inactivation parameters are shown in Figures 3Bi,Bii, respectively. Figure 4 shows the activation, inactivation, and deinactivation variables based on the in vitro data. The activation curve ( Figure 4A) was obtained by normalizing the low-threshold spike (LTS) amplitude for each depolarizing current stimulus to the maximal LTS amplitude of that cell. The initial activation of LTS occurred at −76 mV and complete activation at −42 mV with the greatest variability between neurons occurring in the linear range between −55 and −45 mV. The inactivation curve ( Figure 4B) was obtained by normalizing the rebound potential at the end of each hyperpolarizing current prepulse to the maximum rebound potential for a neuron generated by the most hyperpolarized prepulse. The time course for deinactivation curve ( Figure 4C) was obtained by eliciting hyperpolarizing pulses of constant amplitude from a holding potential of −40 mV and varying the width of the hyperpolarizing prepulse. The amplitude of the rebound potential was normalized to the maximum amplitude for that cell and plotted as a function of hyperpolarization time. The activation and inactivation curves for the rebound conductance obtained from the experimental data were reduced to 30 data points by a fit to the nearest 5 mV value in the range −90 to −40 mV. These data were used to interpolate a table of 3000 voltage states across that range to obtain the smoothed and fitted activation and inactivation curves ( Figure 4D) to be used in the cortical neuron model.

Matching In vitro Data to Single Neuron Model
Details of the reduced infragranular pyramidal cell model have been described in the Methods Section (also see cortical component CX in Figure 1). The model consisted of the FIGURE 3 | Voltage responses to current stimulus protocols. (A) Typical voltage responses to current stimulus for (i) activation and (ii) inactivation of the low-threshold conductance. Blockade of PIR response by Co 2+ (iii), demonstrates the Ca 2+ dependence of the PIR response. (B) (i) Experimental protocol for determination of the activation parameter consisted of 500 ms step pulses of increasing current amplitude. The holding current was adjusted such that the membrane potential was at −80 mV at the start of the step pulse. (ii) Protocol for determination of the inactivation parameter consisted of hyperpolarizing current steps of increasing amplitude and of 500 ms duration. The holding current at the start of the step pulses is adjusted such that the membrane potential is at −45 mV.
following currents: fast sodium (I Na ), potassium delayed rectifier (I K_DR ), cortical low-threshold/rebound (I CX_REB ), hyperpolarization-activated cation (I H ), potassium afterhyperpolarization (I AHP ), calcium dependent potassium (I K[Ca] ), and the potassium A-current (I A ). The channel description for the I CX_REB current was based on in vitro data directly obtained from the activation and inactivation curves ( Figure 4D). The experimentally determined activation and inactivation curves (Figures 4A,B) were reduced to 30 data points by "eyeball" fit to nearest 5 mV value from −90 to +5 mV. In our experimental dataset, this ranged from −80 to −40 mV and all values outside this range were set to 0 or 1. The GENESIS function tweaktau was then used to extrapolate the data to 3000 data points within the range −80 to −40 mV to obtain a smooth curve. The function scaletabchan was used to left shift the curves by −6 mV to give a more accurate depiction of channel properties and obtain realistic LTS and PIR responses.
An extensive parameter search was performed to obtain physiologically realistic responses to current stimulus protocols to demonstrate post-inhibitory rebound in the model. We found it was necessary to left shift the curves in the model by 6 mV (i.e., in a more hyperpolarized direction) to generate rebound that matched experimentally observed results. Figure 4D shows the final smoothed and shifted activation and inactivation curves for the I CX_REB conductance implemented in the model. The maximal conductance values used were (in S/m 2 ): g K_DR = 200, g CX_REB = 7, g H = 12.55, g AHP = 2.1, g K[Ca] = 0.28, and g A = 12. To determine rebound parameters in isolation, the fast Na + maximal conductance was set to zero (g Na = 0) to simulate the effect of the intracellular Na + blocker QX-314. Figure 5 shows the simulated LTS and post-inhibitory rebound in the reduced 4-compartment cortical neuron model. With the membrane potential set at −80 mV and gNa = 0 (to simulate the action of QX-314), incremental LTS spikes and rebound potentials are observed with increasing depolarizing and hyperpolarizing current pulses, respectively (Figures 5A,B). The near all-or-nothing nature of the LTS evoked by stimuli between 0.7 and 0.8 nA reflects the steep slope of the activation curve. In these simulations, the maximum width of the LTS at its base was 70 ms which is close to the average value observed in the experimental data (67.5 ms). The maximum amplitude of the simulated LTS was 15 mV, near the average experimental value of 13.7 mV. The maximum rebound potential followed a hyperpolarizing current pulse current pulse of −0.7 nA which was well within the physiological range. The rebound width was 110 ms which exceeded the mean but was in the range of experimentally observed values. The amplitude of the rebound spike was close the mean observed value ( Figure 5B). Figure 5C shows the increasing amplitude of rebound response to a constant hyperpolarizing current pulse of −0.8 nA with the maximal low-threshold conductance g CX_REB varied from 0 to 7 S/m 2 . This finding shows that the post-inhibitory rebound response in the cortical cell model arises primarily as a result of the low-threshold conductance implemented from experimentally observed activation and inactivation data.

Network Model
Details of the network model are outlined in the Methods Section. Figure 6A shows the development of oscillations in the interconnected network simulation (see Figure 1). For this simulation, the synaptic parameters were adjusted to keep the CX and TC cells quiescent while RE was bursting intrinsically at 2.5 Hz (Figures 6Ai-iii). In the first second of simulation time, there are two bursts in the RE neuron whereas the TC and CX neurons show no spiking. This initial condition is based on the theory that a few RE or TC cells may be subgrouped as initiator cells that may start oscillating intrinsically and then gradually recruit other RE and TC cells as well as in vivo

FIGURE 5 | Simulations of low-threshold spike (LTS) and post-inhibitory rebound (PIR) in cortical layer V cells. (A) Simulated LTS response for increasing depolarizing 500 ms step pulses (left panel). Close up view of the LTS is seen on the right panel showing realistic responses
compared to experimental data. The resting membrane potential was set at −80 mV and g Na = 0 S/m 2 to simulate action of QX-314 and block action potentials. (B) Simulated PIR response to increasing hyperpolarizing step pulses of 500 ms. Close up view of the PIR is seen on the right panel showing realistic PIR responses compared to experimental data. The resting membrane potential was set at −80 mV and g Na = 0 S/m 2 to simulate action of QX-314 and block action potentials. (C) Relationship between the low-threshold conductance channel (g CX_REB ) implemented from experimental data and the post-inhibitory rebound response observed in the cortical cell model. The maximal value of g CX_REB is varied from 0 to 7 S/m 2 . Close-up view of PIR response and its dependence to the low-threshold channel conductance (g CX_REB ) is seen on the right panel. Figures on the right for (A-C) are drawn to different scale to highlight LTS and rebound response. and computational data that show oscillations to be present in isolated RE nuclei (Steriade et al., 1987;Destexhe et al., 1994;Fuentealba et al., 2004). To simulate the RE initiator condition, the synaptic connectivity was removed between RE and TC cells for the initial 1 s period. Addition of inhibitory synaptic connectivity (RE→TC, g GABA = 1 S/m 2 ), produced hyperpolarization in the TC neuron leading to post-inhibitory rebound potential in the TC neuron along with burst spike activity. Burst spiking in the TC neuron also led to reciprocal excitatory EPSPs in the RE neuron and feed forward EPSPs in the CX neuron and the entire circuit oscillated at 6.5 Hz, which is in the lower range of spindle oscillations. Therefore, the RE initiator neurons entrained the entire TC-RE-CX circuit into the spindle oscillations frequency range. We also found that the low threshold conductance in the CX cell was essential for its rhythmic burst firing since setting g CX_REB = 0 abolishes bursting activity in the CX neuron, though subthreshold oscillations remain ( Figure 6B).
We also examined the effect of cortical feedback and resting thalamic leak conductances on the oscillation frequency. Across the range of excitatory corticothalamic synaptic conductances, the oscillation frequency reduced from 6.4 to 5.2 Hz for feedback ranging from 0 to 2 S/m 2 (Figure 6C). For the simulations in Figures 6A,B, g FEEDBACK was set at 0.4 S/m 2 . Figure 6D shows the effect of varying the TC resting potassium conductance value, g LEAK , on oscillation frequency. At the initial value of 11 nS, the frequency of oscillation was at the lower end of the spindle range of 6 Hz. Increasing the g LEAK caused a decrease in the frequency, whereas decreasing g LEAK caused an increase in the frequency. The range of frequencies was from 6.8 Hz at a g LEAK value of 0 nS −4.8 Hz at a g LEAK value 20 nS. Therefore, there is an inverse relationship between this passive potassium conductance value and the frequency of oscillation. Altering the g LEAK conductance represents changes in the resting membrane potential of TC neurons that may arise from excitatory brainstem inputs during different physiological states Willis et al., 2015).

Phase Plots
In a computational model, phase plots or phase portraits provide a way to plot multiple dynamic variables that vary with time (e.g., activation and inactivation gates of conductances, Ca 2+ concentration, etc.) as a function of each other at the same time during continuously evolving neuronal behavior. These plots provide us with information about the dynamics of these variables that ultimately determine the evolution of the membrane potential. Figure 7 shows two representative data sets for 15 s of simulated data in the putative chaotic condition and in the periodic condition for the CX neuron. The periodic oscillations at 6-7 Hz were obtained by setting TC→CX synaptic excitation to 2.0 S/m 2 and RE→TC, and synaptic inhibition to 1.0 S/m 2 ( Figure 7A). The irregular aperiodic bursts observed in the putative transition region were obtained by setting TC→CX synaptic excitation: 1.4 S/m 2 and RE→TC synaptic inhibition: 0.08 S/m 2 (Figure 7B). To discern the underlying mechanisms between the two states, we constructed four phase plots for each of the periodic and putative chaotic conditions focusing on the Ca 2+ concentration and the calcium channel activation/inactivation dynamics. These include ( Figure 7B. In contrast to Figure 8A, these plots show a rich banded structure and overlapping trajectories which characterize the presence of strange attractors within a chaotic regime.

Bifurcation Plot
Whereas, phase plots reveal information about the underlying dynamics by plotting variables against each other, bifurcation plots show the possible transition from periodic to aperiodic or chaotic states. Figure 9 shows the bifurcation plot of the instantaneous frequency of spikes (1/ISI) in Hz vs. the TC→CX excitatory synaptic conductance with RE→TC synaptic inhibition set at 0.08 S/m 2 . At low TC→CX values a single period is observed. The spread of points in the lower range increases (0.7-1.2 S/m 2 ) until at 1.4 and 1.6 S/m 2 where an apparent period-three limit occurs (indicated by arrow). This implies the presence of a chaotic regime (Li and Yorke, 1975;Canavier et al., 1990). Data points between TC→CX values of 1.4 and 2.2 S/m 2 occupy a spectrum of values in both frequency ranges with no points in the middle region between 4 and 5.5 Hz. We have also shown the abscissa point TC→CX = 2.2 S/m 2 with RE→TC synaptic inhibition increased to 1.0 S/m 2 which is the periodic oscillatory bursting regime (see Figure 7A). At this point the lower range of frequencies abruptly disappears and only the higher frequency points remain with some spread in this range. Thus, the bifurcation plot demonstrates the route to chaos and back to periodicity.

Fano and Allan Factor
Fano factor and Allan Factor were computed for the chaotic and periodic conditions and plotted on a log-log scale (Figure 10). Both measures show power law growth between 10 0 and 10 1 for longer counting times which indicate fractal nature of burst occurrences (Lowen and Teich, 1996). The fractal dimensions estimated from the slope of these plots were 1.16 for the Fano Factor (solid squares) and 2.31 for the Allan Factor (open circles). For the periodic condition, at longer counting times in the same range, the burst pattern shows negative slope.

Lyapunov Exponents
We calculated the Lyapunov exponents for the two data sets shown in Figure 7. With the variable measured being the membrane potential of the CX neuron, the data set representing the transition region the Lyapunov exponent value was computed to be +0.24 for 1000 iterations, which, being a value greater than zero, indicates a chaotic regime (Abarbanel et al., 1993). In contrast, in the periodic condition, a negative Lyapunov exponent value of −0.05 was obtained.
Given the importance of RE→TC inhibition in generating post-inhibitory rebound and oscillatory states, we further characterized LE as a function of RE→TC synaptic inhibition ( Figure 11A). With the CX membrane potential once again as the variable of interest, positive Lyapunov exponent values are obtained for RE→TC inhibition 0.02 S/m 2 to just below 1.0 S/m 2 . At ≥1.0 S/m 2 the Lyapunov exponents abruptly become negative indicating a transition to non-chaotic regime. Therefore, the cortical membrane potential demonstrates sensitivity to initial conditions between 0.02 and 1.0 S/m 2 , which is a hallmark of a chaotic system. Given the importance of the strength of CX→TC feedback on oscillation frequency (Figure 6C), we also examined Lyapunov exponents at varying CX→TC excitatory feedback synaptic strength while keeping other synaptic connectivity strengths constant (RE→TC = 0.08 S/m 2 ; TC→CX = 1.4 S/m 2 ). Varying CX→TC synaptic strength from 0.05 to 1 S/m 2 produces Lyapunov exponents in the range 0.2-0.6 ( Figure 11B). However, removal of the cortical feedback loop (CX→TC = 0.0 S/m 2 ) produces a larger Lyapunov exponent of 1.4 (not shown in semilog scale in Figure 11B). In our previous simulations with only an interconnected TC and RE neurons, positive Lyapunov exponents were obtained in the range of 1-2 ( Figure 5 in Paul et al., 1998). There is a wide spectrum of Lyapunov Exponent magnitudes reported in neural simulations with positive values in the range of 10 −4 to show chaotic behavior in a single neuron model of R15 bursting cell of Aplysia (Canavier et al., 1990). Figure 11C shows the Lyapunov exponents for the CX cell as a function of the potassium leak conductance (K_leak). Decreasing the leak conductance from its initial value simulates the effect of excitatory brainstem neuromodulatory input that accompanies the transition from sleep to wakefulness. The other control parameters are set at TC→CX = 1.4 S/m 2 and RE→TC = 1.0 S/m 2 . With these settings, the LE value at initial g LEAK of 11.25 nS is −0.05-only just negative. For lower g LEAK , Lyapunov exponents were positive (except at g LEAK = 2.25 nS), suggesting FIGURE 8 | (A) Phase Plots for various parameters representing low-threshold conductance dynamics and Ca 2+ concentration in CX and TC cell of network simulation. All phase plots exhibit a 2 period limit cycle. "Control" parameter settings are TC→CX synaptic excitation: 2.0 S/m 2 and RE→TC synaptic inhibition: 1.0 S/m 2 . Computed Lyapunov Exponent is −0.05. (B) With TC→CX synaptic excitation: 1.6 S/m 2 and RE → TC synaptic inhibition: 0.08 S/m 2 , phase plot analysis shows strange attractor dynamics for the same parameters as in (A). Lyapunov Exponent is +0.24 and fractal dimension estimate from the Allan Factor time curve is 2.3. that brainstem modulators, which modulate the state of arousal (Williams et al., 1994;Jones, 2003), may produce a transition from periodic to chaotic dynamics.

DISCUSSION Summary
A minimal network model consisting of one TC cell, one RE cell, and one CX cell was developed based in part on realistic biophysical parameters as well as measured values from physiological experiments. These cells were interconnected with realistic anatomical projections and implemented as previously hypothesized circuit for the neural basis of thalamocortical FIGURE 9 | Bifurcation plot of CX cell in network model in transition region. RE→TC synaptic inhibition is set at 0.08 S/m 2 . Frequency of burst occurrences are plotted vs. increasing TC→CX synaptic excitation. Bursts are defined as groups of spikes with interspike intervals of less than 10 ms. On the x-axis, the point 2.2 represents a TC→CX synaptic excitation of 2.2 S/m 2 and RE →TC synaptic inhibition of 0.08 S/m 2 as well as TC→CX: 2.2 S/m 2 and RE →TC: 1.0 S/m 2 as distinct points on the x-scale.. The latter represents the fully periodic oscillation state. Generation of apparent period-three cycle is indicated in the route to chaos.
oscillations (Steriade et al., 1993a,b;Steriade, 2005;Contreras, 2014). The network model was constructed such that the RE cell acted as an initiator cell and subsequently entrained the TC and CX cells into a bursting rhythm. We observed that chaotic dynamics were present between thalamocortical states that accompany the physiological transition from wakefulness to drowsiness. The presence of chaos was dependent on a number of physiological parameters, such as the resting leak conductance of TC neurons, and the strength of thalamocortical transmission. These data add to the growing body of literature suggesting that chaotic dynamics are present in real and simulated neural systems (Canavier et al., 1990;Wang, 1994;van Vreeswijk and Sompolinsky, 1996;Siegel and Read, 2001;Bertschinger and Natschläger, 2004;Battaglia et al., 2007;Sussillo and Abbott, 2009;Rajan et al., 2010;Jia et al., 2012) and suggest that the transition between sleep and wakefulness may be characterized by chaotic dynamics. The limitations and implications of this study are discussed below.

Nonlinear Dynamical Analysis
The present study examined the region between the full scale oscillation observed during sleep and drowsiness and the transition region characterized by irregular burst occurrences. Simulation runs with phase plots of various parameters representing the low-threshold conductance dynamics revealed the presence of strange attractors in the transition region, in contrast to a two-period limit cycle observed in the periodic condition. Previous studies have shown the presence of strange attractors in the transition region between beating and bursting modes in simple models that incorporate a single isolated neuron with biologically realistic conductances and a variable externally applied current that determines the state of the neuron (Plant and Kim, 1976;Canavier et al., 1990;Wang, 1994;van Vreeswijk and Sompolinsky, 1996;Carelli et al., 2005). Other studies have similarly shown the presence of a chaotic region in a pair of coupled oscillators (Paul et al., 1998;Iqbal et al., 2014). Paul et al. (1998) demonstrated a novel method of calculating Lyapunov Exponents which was possible within the GENESIS environment. We have now extended those findings to a more complex and anatomically realistic circuit with the addition of a cortical component and used multiple lines of analysis that converge on the finding that a chaotic regime occurs during the transition between wakefulness and sleep. In this study, we have taken a proposed hypothesized circuit for thalamocortical oscillations, experimentally determined the activation, and inactivation parameters of the critical postinhibitory rebound conductance in the cortical component of the network model and then demonstrated the presence of a chaotic regime in the cortical neuron whose network dynamics are solely under the control of network synaptic input to CX neuron or output from CX neuron. The major synaptic strength being varied is the RE→TC inhibition. Therefore, with no direct externally applied current, we show that the CX neuron dynamics demonstrate the presence of strange attractors and a chaotic regime. Therefore, the thalamo-reticulo-cortical loop with its closed organization and architecture and multiple time scales, isolated from sensory inputs, may also generate chaotic complexity (Cauller, 2003). The physiological significance of this is elaborated further in the discussion of nonlinear dynamics in sensory processing.
In the bifurcation plot we show a route to chaos and back to period-one regime with an apparent period-three stable regime at TC→CX excitation strength of 1.4 and 1.6 S/m 2. The wide spread of points show multiple limit cycles within the putative chaotic regime and is clearly different from what is observed at very low TC→CX excitatory synaptic strength as well as when RE→TC inhibition is increased to 1.0 S/m 2 . Although the period-three orbit is somewhat arbitrary in our simulations, the wide spread of points within the transition region tend to fall into two distinct clumps and may be described as a bifurcation between dynamical states (Ciszak and Bellesi, 2011).
Presence of a period-three regime is known to be a hallmark of the transition to chaos (Canavier et al., 1990). Both the Fano Factor and the Allan Factor demonstrate a power law growth for longer counting times in the range 0-1 and indicate that the irregular burst pattern observed in the transition region is fractal in nature.
Finally, simulation runs were carried out to estimate the Lyapunov Exponent of the dynamical system. The Lyapunov Exponent measures the sensitivity of the dynamical system to just slightly different initial conditions. Positive Lyapunov exponents is considered to be the quintessential indicator of a chaotic regime (Abarbanel et al., 1993). Simulation runs in the representative transition region produced a positive Lyapunov exponent. In the periodic case, the Lyapunov Exponent was negative. All of these findings strongly indicate a chaotic regime in the transition region of network activity. These findings are consistent with a body of literature suggesting that sleep-wake transitions may be characterized by power law dynamics and as a bifurcation between different states (Lo et al., 2004;Chu-Shore et al., 2010;Ciszak and Bellesi, 2011), though this has not yet been measured using tools that provide resolution at the level of synaptic currents. In addition, it is well known that paroxysmal states such as seizures are more prone to occur during the sleep-wake transition (Da Silva et al., 1984; , 1991). These data suggest that the presence of strange attractors at the transition state allows rapid global transitions from one state to another (e.g., from wakefulness to sleep), but may also lead to rapid transitions to pathological states (e.g., seizures). Lo et al. (2004) noted power-law "fractal-like" behavior in sleep-wake transitions in mice and rats and have suggested that the transitions between sleep and wakefulness were not random fluctuations but related to underlying neural mechanisms of sleep control. The development of neural models for wakefulness-sleep transitions would be helpful for the understanding and regulation of sleep and wakefulness.

Limitations of Current Study
The network model employed in the present study was minimal, in the sense that it consisted of one TC, one RE, and one CX cell, each with simplified morphologies, an approach based on the assumption that the behavior of large networks of neurons can be explained by the behavior of smaller networks of neurons. Clearly other elements could potentially influence the behavior of the model, such as dendritic calcium currents (Crandall et al., 2010), modifications of the calcium dynamics in RE and TC cells to reflect their physiological differences (Talley et al., 1999), or the addition of more cells to introduce potential cross-channel effects, due to either dendro-dendritic synapses, gap junctions, or nonreciprocal interactions between TC and RE cells (Deschênes et al., 1985;Pinault and Deschênes, 1998;Crabtree and Isaac, 2002;Landisman et al., 2002). Furthermore, in our model we do not distinguish between the diffuse matrix thalamocortical pathway projecting to the superficial layers of CX and the spatially selective core pathway which project to the granular layer (Piantoni et al., 2016). It is likely that such modifications would alter the quantitative behavior of the model, though it is not yet known if they would alter it qualitatively (i.e., would alter whether chaotic behavior would be observed). Despite these simplifications, we note that simple models have been used successfully to gain insights into network behavior (Grillner et al., 1988;Destexhe et al., 1993;Ching et al., 2010) though the ability of the current model to scale up to populations of neurons is yet unknown. Thus, a future version of this model could include multiple instances of each type of cell which may increase the validity of the model and provide opportunities to examine the impact(s) of the modifications described above.

Experimental Predictions Made by the Model
The results from the current study predict that spike trains in cortical and thalamic neurons should reflect chaotic dynamics during the transitions from the sleeping to the waking state, and that non-chaotic dynamics should be seen during sleep and during waking state. To determine if chaotic dynamics definitively exist in real spike trains would involve computation of Lyapunov exponents from time series data from an in vivo preparation, which is extremely difficult (but not impossible) to compute on real biological neurons. Therefore, the more feasible approach would be to determine if fractal patterns are seen in spike trains during the sleep-wake transition. However, difficult, it may still be possible to determine Lyapunov exponents from time series data based on reconstructed attractor maps and determining the long term evolution of vectors defined by points on separated by at atleast one mean orbital period (Wolf et al., 1985).
Another prediction made by the model is that corticothalamic feedback has significant effects on thalamocortical function, which could be measured. For example, in the current study, the effect of increasing g FEEDBACK was a lowering of the network oscillation frequency (Figure 6C). In addition, weakening cortical feedback (to zero), elevated the Lyapunov exponent to 1.4, suggesting a greater tendency for chaotic dynamics to be seen. Strengthening or weakening corticothalamic feedback could be done using a range of modern genetic tools (Stroh et al., 2013;Denman and Contreras, 2015), and the effects of this manipulation can be measured on the firing properties of thalamic neurons during sleep-wake transitions.
Similarly, in the current study, there is an inverse relationship between the leak conductance value (g LEAK ) and the frequency of oscillation (Figure 6D). Varying g LEAK on the TC and RE cells simulates the effect of neuromodulatory input (e.g., cholinergics) on the thalamus during transition from sleep to wakefulness. Excitatory brainstem neuromodulatory input causes a block of K LEAK channels thus, causing the membrane potential to become more depolarized (McCormick and Prince, 1987). This effect is simulated by reduction of g LEAK which results in an increase of oscillation frequency. Such an effect should be readily testable using application of cholinergic or other modulators while measuring thalamocortical oscillation frequency.

Nonlinear Dynamics in Sensory Processing
In this study we have explored the chaotic nature of irregular aperiodic bursting that occurs in the transition from sleep to wakefulness. Irregular bursts may also occur during and due to sensory processing of stimuli in the thalamus and cortex (Krahe and Gabbiani, 2004;Desbordes et al., 2008), which must merge with the ongoing thalamocortical activity characterized by spontaneous tonic activity (Cauller, 2003). We hypothesize that the coupling of bottom-up sensory input bursts with top-down feedback may also constitute chaotic regimes with infinite limit cycles which allow an organism to quickly switch attention or respond to external inputs. This hypothesis may be tested in more detailed large scale neuronal models that incorporate a layered cortical structure.
Chaotic behavior has been previously obtained in minimal reciprocal models of neuron pairs with small changes in connection strength (Jackson et al., 1996;Paul et al., 1998). In the current study, we have shown that the simulation of a minimal thalamo-reticulo-cortical loop, isolated from sensory inputs, also generates chaotic complexity. Therefore, the thalamo-reticulo-cortical loop with its closed organization and architecture and multiple time scales, may also generate chaotic complexity with a system of multiple self-organized attractors and attractor sequences where each attractor space or unit refers to a particular spatio-temporal firing pattern (Cauller, 2003). From a physiological perspective, the cortex may contain specific instances self-organized attractor spaces with the rich complexity of chaotic dynamics which contain an internal model based on the sum of experiences of an organism. Top-down influences are continuously probed and tested by bottom up sensory inputs which may modify the attractor sequence. Even in the development phase a newborn infant spends much of its time in REM sleep which may correspond to dynamic self-organization of attractor spaces with chaotic complexity and is primed for further refinement as the infant begins to explore the environment (Cauller, 2003).

Conclusions
The current data extend previous work showing that chaotic dynamics are present in neural spike trains at the transition points between stable states. Chaotic dynamics may permit rapid transitions between states-increasing behavioral and cognitive flexibility, possibly at the cost of transitions to pathological states, such as seizures. The use of a computational model permits quantification of certain values, such as the Lyapunov exponent, which are very difficult to compute in biological data. Calculation of the Lyapunov exponent permitted a systematic exploration of the impact of corticothalamic feedback and thalamic depolarization on chaotic behavior. This exploration will allow new hypotheses to be tested about the transition between sleeping and waking states.

AUTHOR CONTRIBUTIONS
KP and LC designed the study. KP conducted electrophysiological experiments and developed computational models. KP, LC, and DL analyzed data and interpreted results. KP and DL prepared manuscript.

FUNDING
This work was supported by the NIH grant NIDCD-R21DC014765 awarded to DL.