Model-based analysis and forecast of sleep–wake regulatory dynamics: Tools and applications to data

F. Bahari; J. Kimbugwe; K. D. Alloway; B. J. Gluckman

doi:10.1063/5.0024024

Chaos. 2021 Jan; 31(1): 013139.

Published online 2021 Jan 25. doi: 10.1063/5.0024024

PMCID: PMC7837756

PMID: 33754773

Model-based analysis and forecast of sleep–wake regulatory dynamics: Tools and applications to data

F. Bahari,^1,2,^a) J. Kimbugwe,^1,2 K. D. Alloway,^2,3 and B. J. Gluckman^1,2,4,^b)

Author information Article notes Copyright and License information PMC Disclaimer

Associated Data

Data Availability Statement: Computational examples for the algorithms developed and/or integrated together, an example of animal sleep data on which these algorithms were applied, and a code for generating the figures from Penn State's ScholarSphere at https://doi.org/doi:10.26207/fcc5-7h57.

Abstract

Extensive clinical and experimental evidence links sleep–wake regulation and state of vigilance (SOV) to neurological disorders including schizophrenia and epilepsy. To understand the bidirectional coupling between disease severity and sleep disturbances, we need to investigate the underlying neurophysiological interactions of the sleep–wake regulatory system (SWRS) in normal and pathological brains. We utilized unscented Kalman filter based data assimilation (DA) and physiologically based mathematical models of a sleep–wake regulatory network synchronized with experimental measurements to reconstruct and predict the state of SWRS in chronically implanted animals. Critical to applying this technique to real biological systems is the need to estimate the underlying model parameters. We have developed an estimation method capable of simultaneously fitting and tracking multiple model parameters to optimize the reconstructed system state. We add to this fixed-lag smoothing to improve reconstruction of random input to the system and those that have a delayed effect on the observed dynamics. To demonstrate application of our DA framework, we have experimentally recorded brain activity from freely behaving rodents and classified discrete SOV continuously for many-day long recordings. These discretized observations were then used as the “noisy observables” in the implemented framework to estimate time-dependent model parameters and then to forecast future state and state transitions from out-of-sample recordings.

Sleep and neurological diseases appear to mutually affect each other. Here, we extend tools for using model-based smoothing, observation, and forecasting to observe and forecast state of vigilance (SOV) (sleep/wake state), including innovations in multi-parameter fitting, and for the first time demonstrate the application of these tools to forecast sleep/wake state from experimentally measured time series of animal data. Because the methods utilize mechanistic models of the sleep–wake regulatory system, these tools have the potential to provide insight into the mechanisms linking neurological disease dynamics and the sleep–wake regulatory system.

I. INTRODUCTION

Dysregulated sleep is one of the many comorbidities of a wide variety of neurological disorders. Not only is sleep regulation affected in diseases such as epilepsy, Alzheimer’s, schizophrenia, and dementia, it seems to modulate the disease severity and/or symptoms.^1–3 Improved understanding of the underlying mechanisms that regulate state of vigilance (SOV) transitions will allow us to better understand the role of sleep in neurological disorders.

The sleep–wake regulatory system (SWRS) is built on a variety of neurochemical interactions.^4–6 These interactions could be relatively fast, on the scale of minutes, to bring about state transitions or relatively slow, on the scale of hours, to account for sleep pressure as well as environmental inputs such as light cycle and other time of day cues or loud noises.⁷ Measurements spanning all the neuronal and chemical elements involved are often both technically challenging and highly invasive. This severely limits our knowledge of the underlying mechanisms of sleep–wake regulation and, therefore, affects our understanding of the role of sleep in neurological disorders.

Increasing amounts of our experimental knowledge of sleep–wake physiology is being utilized to design complex mathematical models of sleep regulation. These models can offer a minimally invasive and mechanistic view of the underlying neurophysiological phenomena. Computational models of SWRS are often species specific, high-dimensional with many parameters, and include both non-linear and stochastic components.^8–11 However, validation of these models and then utilization of them to gain further insight into dynamics of normal brain, and what is different in diseased brain, is technically challenging. These challenges include that the regions of interest are deep subthalamic and brainstem regions and must be monitored in unconstrained subjects over long periods.

We aim to contribute computational tools to observe dynamics of the SWRS with details currently unavailable from experimental data. We have been developing a data assimilation framework that utilizes mathematical models of the SWRS along with sparse experimental observations from the SWRS to reconstruct the remaining unmeasured components.

Data assimilation (DA) is an iterative prediction–correction scheme that couples mathematical equations of a computational model to observed measurements from a system. The purpose is to estimate both measured variables and those that are unobserved, or hidden, as well as to predict the future system state. Ensemble Kalman filters are often used in the DA algorithm for nonlinear systems.¹² One such ensemble filter is the unscented Kalman filter (UKF).¹³ Within the iterative prediction–correction scheme, UKF uses the governing equations of the sleep models and experimentally measured observations to fully reconstruct model state and correct model-generated forecasts to agree with the experimental measurements.¹⁴

The objectives of the work presented here are to demonstrate applicability of DA in the context of a relatively high-dimensional nonlinear biological model of the SWRS, to investigate the model capability in forecasting future states, and to identify interactions of the mechanics embodied in the models with neurological disorders.

Our group has previously illustrated the DA framework and validation of the mathematical algorithms in Sedigh-Sarvestani et al.³ Briefly, we used artificially generated data from the Diniz Behn and Booth⁹ model or its extension by Fleshner et al.⁸ We selected a subset of the generated variables as the measured data set and added measurement noise to them. We then utilized the UKF to reconstruct the unobserved variables (model states) and forecast future system states. We validated the model dynamics by calculating the difference between the reconstructed estimates and parameters and the known original data set.

The work presented in this article adds to the previous work on multiple fronts. First, we extended the previously described multiple-shooting method in Sedigh-Sarvestani et al.³ to simultaneous multivariate parameter fitting and state reconstruction. Second, we employed smoothing algorithms to reconstruct an unobservable model variable that represents random excitatory internal and external inputs to promote wakefulness. Third, we improved the previously described observation function to utilize the hypnogram to derive model observables. Finally, Sedigh-Sarvestani et al. described a computational framework to demonstrate the feasibility of using DA to reconstruct physiological data. Here, we utilize our improved computational toolbox to apply DA to experimental measurements from normal rodents to then reconstruct the animal state of vigilance and form state forecasts.

In Sec. II, we introduce the physiology behind the mathematical model. The background mathematics of the model, UKF, and the smoothing algorithm are described in Appendixes A–D. In Sec. III, we present a computational toolbox for simultaneously estimating multiple parameters and reconstructing states. The section includes application of the computational toolbox to model-generated data, discussions of UKF-based state estimation, smoothing, global, local, and time-dependent multiparameter estimation, and extensions to observation from hypnogram. We then apply these tools to reconstruct the model state and forecast behavioral state from experimental measurements.

II. METHODS

The ground truth here is either model-generated data or experimental discretized state of vigilance (SOV) acquired from the animal. A 4th order Runge–Kutta estimator with an integration time of 0.1 s is used for all simulations.

A. Represented sleep physiology in the models of sleep–wake regulatory system

Neurophysiological experiments have confirmed that sleep–wake states are driven by brainstem and hypothalamic circuits.^4–6 The states are predominantly characterized from electroencephalogram (EEG) and electromyogram (EMG) and have been described as discrete states with fast-switching mechanisms. Models of underlying mechanisms of sleep–wake regulation attempt to reproduce these quantized transitions. To do so, these models often invoke the experimentally hypothesized co-inhibitory dynamics analogous to electrical flip-flops. Non-rapid eye movement (NREM) sleep is promoted by gamma-aminobutyric acid (GABA)-ergic ventrolateral preoptic nucleus (VLPO) neurons in the hypothalamus. Wake is promoted by monoaminergic cell groups in the brainstem, including the noradrenergic locus coeruleus (LC) and the serotonergic dorsal raphe (DR) neurons. Mutual inhibition between these two groups creates a winner takes all flip-flop dynamic. Transitions between the monoaminergic cell-groups in LC and DR and the rapid eye movement (REM) promoting cholinergic brainstem cell-groups in laterodorsal tegmentum (LDT) and pedunculopontine tegmentum (PPT) are also analogous to discrete predator–prey like interactions.¹⁵

Transitions between different sleep–wake states occur on faster timescales (over seconds). There exists slower modulatory dynamics including the homeostatic and circadian sleep drives. Orexin producing neurons in the lateral hypothalamus send descending projections to all monoaminergic and cholinergic cell groups described above and are thought to slowly reinforce arousal during sleep.¹⁶ Previous reports have found that extracellular adenosine increases during prolonged wakefulness in several cortical and subcortical regions.¹⁷ Further research has proposed adenosine as a homeostatic accumulator of the need to sleep.¹⁸

The circadian drive is regulated by the suprachiasmatic nucleus (SCN) in the hypothalamus, which sets a roughly 24-h cycle affecting sleep and many other physiological functions. The SCN sends indirect projections to VLPO in the hypothalamus that results in inhibition of sleep during the day.¹⁹ “Day” is subjectively defined by species’ dependent diurnal behavior and refers roughly to 12-h periods consisting mostly of active-wake behavior. The SCN clock can be modulated by a variety of cortical inputs in response to external cues.

We use two mathematical models of sleep–wake regulation: Diniz Behn and Booth (DBB) model of sleep⁹ and the Fleshner, Booth, Forger, and Diniz Behn (FBFD) model of sleep.⁸

The DBB model of the sleep–wake regulatory system consists of five distinct neuronal populations: wake promoting groups of LC and DR, the NREM promoting group at VLPO, and two groups in the LDT/PPT populations, one that is REM promoting (R) and one that is active during both REM and wake (WR). The cell groups’ interactions are modeled via neurotransmitter release dynamics: LC neurons release norepinephrine (NE), DR neurons release serotonin (SR), neurons in the two groups in the LDT/PPT release acetylcholine (ACh), and VLPO neurons release GABA.

Within the DBB model, h represents the homeostatic sleep drive and the random excitatory thalamic input is modeled by the variable $δ$ . The thalamic input represents a random sensory input such as loud noises that could wake the subject up. State of vigilance is then determined by evaluation of the firing rates of LC, DR, VLPO, and LDT/PPT cell groups with respect to a threshold.

The activity of each cell group is described by its firing rate $(F)$ and the accumulation and dissipation of the neurotransmitter $(C)$ that it releases to the post-synaptic populations.

The FBFD model extends the DBB model to include the brain’s circadian drive from SCN by adding a GABAergic cell group representing SCN. The SCN has an enforced 24-h circadian cycle (CIRC), with higher activity during the light period and lower activity during the dark period. In addition to its own cyclic dynamics enforced by CIRC, the variable SYN represents the synaptic inputs to SCN from release of SR and ACh. This combination is modeled by the sum of CIRC and SYN in I $_{S C N}$ .

The DBB and FBFD models are further described, including equations, in Appendixes A and B, respectively.

B. Parameter estimation

The mathematical models we use are based on approximations of the average cell group activities regulating sleep–wake behavior. Therefore, the parameters of these models are inherently estimates of actual physical quantities in brain. Brains from different species and even brains from individuals in the same species vary. Therefore, we expect them to have different representative parameters just as they express different detailed sleep–wake dynamics. These variations introduce the gap between computational models developed using (mostly) simulated conditions and our real and noisy systems (rodents). Therefore, before we can utilize data assimilation to reconstruct and forecast system dynamics using computational models, we have to first estimate the true model parameter values from measured data.

In this section, we describe our multivariate parameter estimation algorithm. This algorithm includes a parameter selection routine where depending on the location of model trajectory in the state space, a subset of parameters are selected for fitting. We then describe our dynamical parameter estimation routine where we consider SCN dynamics as a parameter and use the DBB model along with data generated from the FBFD model to estimate the slowly changing SCN dynamics.

One parameter estimation approach that utilizes the UKF is to track a combination of parameters and states together via the augmented Kalman filter.²⁰ The assumption in this approach is that parameters vary on the same timescales as the state variables. This becomes challenging in high-dimensional models such as the DBB or FBFD due to the increases in the degrees of freedom when neither parameters nor variables are fixed. Furthermore, in these models depending on where in the state space the dynamics are, their sensitivity to parameters can change. A parameter might only modulate state dynamics at certain locations in the state space. Therefore, estimating the parameter with similar fast timescale of the state variables is both an unnecessary computational burden and prone to error.

To address the challenges with parameter fitting, we iteratively estimate parameters over windows of length $T_{w i n}$ . These windows are longer than a typical sleep–wake cycle of the dynamics. Therefore, we can span all the state-dependent parameters and variable dynamics in one estimation cycle. Our method is based on a global expectation-maximization method of least squares introduced by Levenberg and Marquardt.^21,22 We couple this algorithm with the UKF to simultaneously estimate multiple parameters of the model while fully reconstructing the dynamics of the real system.

We first reconstruct hidden states with the UKF using initial given parameters. We will then fully reconstruct model states over $T_{w i n}$ and use the reconstruction in a parameter estimation step to calculate updated parameters. The updated parameters are then used in the UKF for the next iteration of state reconstruction. We will repeat this process until the parameter estimates have converged.

The Levenberg–Marquardt non-linear least squares algorithm combines two minimization methods through a regularization term that is tuned with variable $λ$ : the gradient-descent method (large $λ$ ) and the Gauss–Newton method (small $λ$ ). In the gradient-descent method, the sum of the squared errors is reduced by updating the parameters such that the system moves in the steepest-descent direction to minimize the squared error. In the Gauss–Newton method, the cost function—the sum of the squared errors—is reduced by assuming that the least squares function is locally quadratic and finding the minimum of the quadratic. The gradient-descent algorithm is in effect when the parameters are far from their optimal value. It is replaced by the Gauss–Newton method when the parameters are close to their optimal value. The parameter $λ$ is adjusted on an iteration by iteration basis based on the total cost function value. When the parameters are too far away from their correct value, the algorithm applies small update values (h) to the parameters to slowly improve the approximation such that the system stays in a stable neighborhood. Once the cost function has decreased to the point that there is no danger of divergence, the algorithm switches to the Gauss–Newton method and applies larger update values to the parameters to fully minimize the now locally quadratic cost function. These are achieved through the equation

h = \frac{J^{T} w (y - \hat{y})}{(J^{T} M J + λ I)},

(1)

where $J$ is the Jacobian of the dynamics with respect to parameters, $M$ is a scaling matrix, and $y$ and $\hat{y}$ are the observed and computed dynamics.

Our parameter estimation step is adapted from the “multiple shooting” method.¹⁴ Within each estimation window, an average cost function ( $C F$ ) is calculated from the divergence between short model-generated trajectories ${\hat{x}}_{t r a j}$ and the UKF-reconstructed trajectories for the measured variables $\hat{x}$ . To prevent divergence of the short model-generated trajectories from the reconstructed ones, we reinitialize ${\hat{x}}_{t r a j}$ on the reconstructed trajectories at regular intervals $d T$ ,

{\hat{x}}_{t r a j} (k) = \tilde{x} (k), where k = 1 * d T, 2 * d T, \dots, T_{w i n} .

(2)

We then calculate a cost function averaged over the window $T_{w i n}$ ,

C F = \int_{0}^{T_{w i n}} {({\hat{x}}_{t r a j} - \hat{x})}^{T} M ({\hat{x}}_{t r a j} - \hat{x}) d t,

(3)

where $M$ denotes a matrix with non-zero elements on diagonal positions corresponding to measured states. In order to properly weight the errors for each variable, the non-zero elements of $M$ are set to the inverse of the standard deviation of the associated variable.

1. Parameter selection

The DBB or FBFD model has, respectively, 12 and 14 variables and 40 parameters that govern neuronal interactions, state-transition timescales, and neurotransmitter releases. Only a subset of the variables in these models can be feasibly observed. Within our experiments, these include neuronal firing rates from DR, PPT, LDT, and VLPO. There are parameters that do not directly affect these observable variables. This means that these parameters are either completely independent of the observables or they influence the observed state variables in combinations such that using the current observables, we might not be able to estimate the parameters separately. This is similar to the parameter identifiability concept.¹⁴ In these cases, often the main concern is whether a unique combination of parameter values exist that minimizes the estimation cost function.

To identify parameter combinations and reduce the number of dependent parameters in our parameter estimation algorithm, we calculated the correlations between parameters. Within each fitting window $T_{w i n}$ , we calculate the sensitivity matrix $g$ ,

g = |\begin{matrix} \frac{\partial F_{1}}{\partial θ_{1}} & \dots & \frac{\partial F_{1}}{\partial θ_{p}} \\ \dots & \dots & \dots \\ \dots & \dots & \dots \\ \frac{\partial F_{D}}{\partial θ_{1}} & \dots & \frac{\partial F_{D}}{\partial θ_{p}} \end{matrix}|,

(4)

where $F$ is the model dynamics and $θ_{1}$ to $θ_{p}$ are the parameters. The correlation matrix is then calculated by inverting the matrix $g^{T} g$ and then dividing the $i$ , $j$ element of ${(g^{T} g)}^{- 1}$ by the square root of the product of the $i$ th and $j$ th diagonal elements. Note that in this model, all parameters are at least locally identifiable (i.e., there exists a finite number of solutions to minimize the cost function of the parameters). Therefore, $g^{T} g$ is non-singular and invertible.

The parameter correlation allows us to group a subset of parameters together and only estimate one parameter in each fitting cycle, while the rest are assigned the same value as the representative parameter.

Because the model state is high dimensional, trajectories can span far and visit areas of the state space that are less frequently visited. The REM state is such an area in that it occurs less frequently and is shorter than NREM or Wake. The subset of parameters that govern REM dynamics affect the model state only for a short period of time when the dynamics transition into REM. Estimation of these parameters at all times (even when they are relatively idle) in the presence of a measurement error could cause undesired divergences in the cost function $C F$ . To identify these parameters and limit their estimation to when it is needed, we calculated the Jacobian of the cost function with respect to each parameter, $\frac{\partial (C F)}{\partial (θ)}$ . Here, the Jacobian is calculated in shorter duration than within the entire fitting cycle ( $T_{w i n}$ ). This is because some parameters might only modulate the model dynamics during a state transition,

J = |\begin{matrix} \frac{\partial C F_{1}}{\partial θ_{1}} & \dots & \frac{\partial C F_{1}}{\partial θ_{p}} \\ \dots & \dots & \dots \\ \dots & \dots & \dots \\ \frac{\partial C F_{D}}{\partial θ_{1}} & \dots & \frac{\partial C F_{D}}{\partial θ_{p}} \end{matrix}| .

(5)

Large deviations in the otherwise stable Jacobian indicate sudden changes in the trajectory path such that certain parameters become effective. We estimate these parameters only during this time frame.

The model parameters include the $g_{i, x}$ as coupling constants between presynaptic and postsynaptic populations, $τ_{X}$ representing response dynamics of the postsynaptic population as well as its neurotransmitter release, and $X_{max}$ as the maximum firing rate a cell group is allowed. The model assumes sigmoidal functions’ firing rate activity and neurotransmitter concentration. For firing rate dynamics, $β_{X}$ and $α_{X}$ represent the midpoint and slope of the sigmoid, respectively. For neurotransmitter concentration dynamics, $γ_{i}$ represents the midpoint of the sigmoid, and $C_{i \infty}$ represents the maximum allowed value.

The parameters $X_{max}$ and $C_{i \infty}$ are scaling parameters that can easily be matched to experimentally derived firing rates without disturbing the correct system outcome. We kept these parameters constant at all times. The model also includes parameters from cell groups (i.e., DR and LC) that have similar effects on the system dynamics. We grouped these parameters together and only fitted one parameter from each group. We then used the estimated values for all the parameters within each corresponding group.

C. Experimental data utilized

The experimental data we used for the purposes of this work were previously described in Billard et al.²³ Briefly, we collected hippocampal local field potentials, electrocorticogram (ECoG), and head acceleration from 12 normal, freely behaving rats to assess the behavioral state.

The data were semi-automatically scored for state of vigilance using linear discriminant analysis (LDA) according to methods previously described in Billard et al. and Sunderam et al.^23,24 Here, we utilize the hypnogram with behavioral states as REM, NREM, and Wake.

III. RESULTS

The main objective of the work presented here is to show the applicability of data assimilation in the context of high-dimensional mathematical models of sleep–wake regulation. We extended previous results from Sedigh-Sarvestani et al.³ on two fronts. In Subsections III A A–III D, we present a computational toolbox to simultaneously estimate multiple parameters and reconstruct state using the UKF. In subsection E, we then apply the toolbox to experimental measurements acquired from freely behaving rodents. We demonstrate that within the data assimilation framework, we can forecast the future state of the system with quite high accuracy.

A. Data assimilation-UKF

As in the previously published work of Sedigh-Sarvestani et al.,³ we accurately reconstructed measured and unmeasured variables of the DBB model with the UKF framework. The details of the DA-UKF framework are described in Appendix C. Briefly, we selected a subset of variables from noisy model-generated data as observations representing experimental measurements and then reconstructed the unobserved variables with the UKF and validated the reconstruction by comparing it to the original model-generated data set.

An example result of the reconstruction and validation process is shown in Fig. 1. The time series of sleep–wake data is generated using the DBB model. One variable, the Wake-active LC region, $F_{L C}$ , is chosen as observable. Explicitly, we added to the model-generated $F_{L C}$ values, random, normally distributed, zero mean noise with a variance of 4% of the variance of $F_{L C}$ . The focus here is to demonstrate that the UKF can be used to reconstruct unobserved model variables with sparse noisy measurements. Therefore, in this step, we kept the parameter values constant and the same as what was used to generate the original data. Default values of the covariance inflation (CI)³ were chosen as $10^{- 5}$ times the typical variance of each variable, with the exception of $δ$ and $F_{R}$ for which we chose 4000 and $10^{- 4}$ , respectively. The CI matrix is nonzero only on the diagonal and is added to compensate underestimates of the forecast uncertainty due to process noise and inadequacies in the model (for more details, see Appendix A). The initial values for the model state, $\hat{x}$ , were arbitrarily chosen.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g001.jpg

Open in a separate window

FIG. 1.

Example of DBB model dynamics and UKF reconstruction. Shown are ODE-integrated state dynamics (black traces) and reconstructed state (red traces). State variables shown are firing rates for the wake-active (LC: $F_{L C}$ ), NREM-active ( $F_{V L P O}$ ), and REM-active (LDT/PPT: $F_{R}$ ) cell groups, neurotransmitter GABA concentration ( $C_{G}$ ), thalamic noise ( $δ$ ), as well as the hypnogram. Noise-contaminated measurements of $F_{L C}$ and $F_{R}$ (blue) were passed to the unscented Kalman filter (UKF) framework to reconstruct all other state variables. Here, the framework was given, and the parameter values were used to generate data and random initial conditions. The reconstructed firing rate dynamics (red) quickly approach and track the true (black) dynamics. The stochastically driven $δ$ dynamics are only somewhat reconstructed.

The reconstructed dynamics of NREM-active firing rate variable $F_{V L P O}$ , the REM-active firing rate variable $F_{R}$ , and the stochastic thalamic noise variable $δ$ as well as the true (model generated) dynamics are shown in red and black of Fig. 1. The reconstruction accuracy of the observed variable $F_{L C}$ can be discerned from how close the reconstructed traces are to both the observation, shown in blue, and the true, model-generated values in black. Similarly, reconstructed values of $F_{V L P O}$ are also close to the true state.

The accuracy of state reconstruction is calculated from the mean square difference between the reconstructed ( ${\hat{x}}_{i}$ ) and true ( $x_{i}$ ) values for each variable. This error is then normalized by the variance of each variable’s dynamics to form a normalized mean squared error for each variable,

ϵ_{i}^{2} = \frac{∣ {(x_{i} - {\hat{x}}_{i})}^{2} ∣}{v a r (x_{i})} .

(6)

The $ϵ_{i}^{2}$ approaches zero as the reconstruction improves. The maximum value of $ϵ_{i}^{2}$ depends on the ratio of the full range of a variable to its variance. For variables such as $δ$ that have a large range and a small variance, $ϵ_{i}^{2}$ is much larger than typical model variables.

For the reconstructions presented here, the initial values of the model states were far from the true system state; therefore, a transient period occurs before the estimated state is close to the actual state. Once the model state converges to the true values, the DA framework keeps the error low and model states stable.

B. Smoothing effect on reconstruction of unpredictable variables

The model thalamic noise state $δ$ is generated from rare (few times per hour) randomly generated pulse inputs. The effect of these rare pulses appears in other states with a delay and, therefore, is difficult to reconstruct.

We first used covariance inflation (CI)³ to improve $δ$ reconstruction. The CI improves reconstruction of $δ$ because it expands the ball of sigma points (see Appendix C) and allows for a larger set of available values. However, as shown in Fig. 1, reconstruction of $δ$ using the CI (red trace) is still not optimal. Furthermore, because $δ$ receives no input from other model variables while it only randomly affects them, its reconstruction remains poor regardless of the choice of observable.

The UKF step only updates the state during current estimation time. To accommodate the delay between delta and its effect being felt in observation, we implemented a smoothing step that refines state in the past based on current observations. We utilized fixed-lag unscented Rauch–Tung–Striebel (FL-URTS) smoother to reduce reconstruction error of $δ$ and consequently of the short awakenings. The details of the FL-URTS smoother are described in Appendix D. The fixed-lag smoothing algorithm permits us to initiate smoothing within each filter iteration with minimal computational effort. Fixed-lag smoothers applied at discrete time $k$ utilize information from steps $[k - n, k]$ to correct state estimate at time $t = k - n$ and propagate it forward. This comes at a computational cost that is multiplicative in the lag $n$ with respect to a single UKF step. We optimized the lag $n$ based on the error between true model-generated data $(x)$ and estimated state $\hat{(x)}$ at each step. We found optimized overall performance balanced against computational time with $n = 30$ .

Shown in Fig. 2 is the UKF reconstruction of $δ$ with noisy observations of $F_{L C}$ and $F_{R}$ with and without smoothing. The dynamics of $δ$ , shown in the upper panels, are often lost without smoothing (left panel). Not only does the smoothing improve reconstruction of $δ$ , it also captures its shape and timing more accurately as highlighted in the insets at expanded time base starting at time 0.5 h. Therefore, the average reconstruction error (lower panels) is lower with smoothing.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g002.jpg

Open in a separate window

FIG. 2.

Effect of smoothing on reconstruction accuracy of $δ$ . (a) Default values of the model with no smoother result in poor reconstruction of $δ$ even with optimized values for CI. The black trace represents the original model-generated data as ground truth. The reconstructed $δ$ with optimized CI is shown in red. (b) Using a smoother will reduce the reconstruction error of $δ$ to less than half its original reconstruction metric (blue trace).The reconstructed, smoothed $δ$ with optimized CI is shown in red. The insets in A and B show the dynamics of reconstructed $δ$ using smoothing expanded within the 0.5-0.6 time-range.

C. Parameter fitting

1. Global parameter estimation

The DA framework we present here requires both the mathematical model for the dynamics and the model’s parameters. To optimize the parameters, we have implemented a version of a multiple-shooting method¹⁴ as illustrated in Fig. 3.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g003.jpg

Open in a separate window

FIG. 3.

Parameter estimation with the LM method. We chose three different parameters: $g_{S, R}$ , $g_{A, L C}$ , and $g_{G, W R}$ for fitting while observing F $_{L C}$ and F $_{R}$ . Parameter estimation is performed using a multiple-shooting method:¹⁴ the divergence of short model-generated trajectories that originate on the reconstructed trajectories from the UKF-reconstructed dynamics is minimized. This divergence minimization is done over time windows longer than the cycle time of the dynamics to better sample the state space. Here, we use 20-min long windows, with a 6-min overlap. The trajectories for the short model generated (cyan), reconstructed (red), and true (black) F $_{L C}$ and F $_{R}$ dynamics for different periods of the convergence of the three parameters are shown in the first two panels. We start with random initial values for each parameter and estimate its correct value at each step: $g_{S, R}$ (green), $g_{A, L C}$ (red), and $g_{G, W R}$ (magenta). Yellow and cyan traces indicate the correct value of the parameters. For parameter values significantly different than the true value, the short trajectories diverge quickly from the reconstructed values, which in turn results in a poorly reconstructed state. As parameters approach their true values, both short model-generated and reconstructed trajectories approach their true values as well. Reconstruction metric $ϵ_{i}^{2}$ is computed for each data assimilation window for all variables. Here, in the last two panels, the reconstruction metric for the original noisy observation of F $_{L C}$ as well as an unobserved variable h is shown. The reconstruction metric for both F $_{L C}$ and h decreases as the parameters converge to their correct values. Note that parameter estimation has the effect of optimizing unmeasured state reconstruction.

For illustration purposes, we generated data with fixed parameters and assimilated noisy measurements of $F_{L C}$ and $F_{R}$ to reconstruct the dynamics. Initially, in the UKF step, all model parameters were set to the same values that were used to generate the true data set—except the parameters $g_{A, L C}$ , $g_{S, R}$ , and $g_{G, W R}$ that couple ACh to LC, SR to the REM-active populations, and GABA to the wake-REM active populations, respectively. These parameters were set to arbitrary initial values.

Parameter estimation is performed by minimizing the distance between the UKF-reconstructed traces and short model-generated trajectories that originate on the reconstructed traces. For these computations, we set the length of the short model-generated trajectories at 2 min. This is long enough that differences in parameters result in measurable divergence between the short computed trajectories and the reconstructed dynamics. Here, we define measurable much larger than the measurement noise, but not so large that the distance between the computed and reconstructed trajectories becomes comparable to the range of the state space.

We average the divergence over time windows that are longer than one sleep–wake cycle in order to sample the full state space for each step of the estimation process. This is critical because the sensitivity of the dynamics to any particular parameter is often localized in the state space, and repeated fitting elsewhere results in estimated parameter drift.³

Our global least-square algorithm permits simultaneous multiple parameter estimation. As shown in Fig. 3, our estimation of three parameters converges to the true value. In the two upper panels of Fig. 3, we demonstrate trajectories for the short model-generated (cyan), reconstructed (red), and true (black) $F_{L C}$ and $F_{R}$ dynamics for different periods of the convergence of $g_{A, L C}$ , $g_{S, R}$ , and $g_{G, W R}$ . Note that initially, for parameter values that are significantly different than the true values, the short trajectories diverge quickly from the reconstructed values, and the reconstructed values of $F_{L C}$ and $F_{R}$ are different from the true ones. When the parameters approach the true value, both the short model-generated and reconstructed trajectories approach the true dynamics.

The parameter estimation routine optimizes short model-generated forecasts. To investigate the effect on reconstruction accuracy, we computed the normalized mean square reconstruction error $ϵ_{i}^{2}$ for each variable, averaged over each parameter estimation window. The reconstruction error for variables $F_{L C}$ and $h$ is shown in Fig. 3. As shown for initial values of the parameters, even reconstruction of the measured variable $F_{L C}$ is quite poor. As a reference point, the normalized mean squared error of the measured data—a noisy version of $F_{L C}$ —is $ϵ_{L C, max}^{2} = 0.01$ , shown as a horizontal blue line. As the estimated parameter converges to its optimized value, $ϵ_{L C}^{2}$ falls well below $ϵ_{L C, m a x}^{2}$ , and the reconstruction improves for all variables.

2. Optimization of the parameter selection procedure

The estimation algorithm described here is different from the augmented UKF routine in that it estimates parameters on a much slower timescale compared to the variables. However, we still estimate all selected parameters over all state space regardless of where the dynamics are in the space. The REM-active parameters such as $g_{N, R}, g_{S, R}, g_{G, R}$ , and $τ_{R}$ are examples of parameters that affect the model state only for a short period of time when the dynamics transition into REM. At all other times, these parameters remain relatively idle and estimating them will increase jitter in their value and make them more susceptible to diverge. To identify these parameters and limit their estimation to when it is needed, we calculated the Jacobian of the cost function with respect to each parameter, $\frac{\partial (C F)}{\partial (θ)}$ , in each fitting cycle.

The cost function is calculated from the difference between short model-generated trajectories and the UKF-reconstructed trajectories for the measured variables. The derivative of the cost function with respect to model parameters, $\frac{\partial (C F)}{\partial (θ)}$ , is a feature of the mathematical model alone and not the data. Therefore, it can be computed and applied at any relative point in the state space. Parameters such as $g_{N, R}, g_{S, R}, g_{G, R}$ , and $τ_{R}$ only affect the model state for a short period of time ahead of and during transitions into REM. The Jacobian for these parameters is very small or zero at times other than transitions into REM. Estimation of these parameters during long periods where the Jacobian values are small will still lead, through noise, to random walk of the parameter values. We prevent this random walk by only estimating the parameters at areas in the state space where the Jacobian is not very small or zero.

Shown in Fig. 4 are the time series of the REM-active cell group firing rate, $F_{R}$ , the parameters $g_{S, R}$ and $g_{A, L C}$ , and their respective Jacobians for three fitting cycles. As the firing rate of the REM-active cell group increases and the dynamics transition into REM, $\frac{\partial (C F)}{\partial (g_{S, R})}$ peaks because the REM-related parameters are now in effect. As the parameter estimation routine initiates, we use $\frac{\partial (C F)}{\partial (θ)}$ to update the parameter value. We then use the updated parameter value for the next filter iteration and reduce the reconstruction error, $ϵ_{i}^{2}$ . The selective fitting of the parameters lowers the computational cost of simultaneous parameter estimation and state reconstruction.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g004.jpg

Open in a separate window

FIG. 4.

Changes in the Jacobian of the cost function with respect to $g_{S, R}$ during REM transitions. The true firing rate of the REM-active cell group (F $_{R}$ ) is indicated by the black trace. The cyan and red traces indicate the model-generated trajectories and UKF reconstruction, respectively. The parameters $g_{S, R}$ and $g_{A, L C}$ are chosen for fitting. As F $_{R}$ peaks, $\frac{\partial (C F)}{\partial (g_{S, R})}$ (green trace), which is a directly REM-coupled parameter, increases. However, the $\frac{\partial (C F)}{\partial (g_{A, L C})}$ trace (blue) directly linked to the wake cell group LC does not change. The change in Jacobian marks the parameters that are going to change given where in the dynamical state the trajectories are. This will reduce unnecessary computational burden to fit all parameters at all times.

3. Dynamical parameter tracking

With the Levenberg–Marquardt parameter estimation algorithm, we can also estimate variables that change much slower than other model states. Because of their much slower timescale, these variables can be regarded as parameters. Here, we used our parameter estimation algorithm and a slightly modified DBB model to reconstruct dynamics from the FBFD model. The DBB model lacks any circadian dynamics, while the FBFD model is an expansion of the DBB model that specifically includes the circadian oscillations driven by SCN.

We used the DBB model, in the UKF process, to assimilate noisy measurements of $F_{L C}$ and $F_{R}$ generated from the full FBFD model and then use them within the parameter estimation routine to track the value of $C_{G (S C N)}$ , the concentration of GABA released from SCN. An example of the output of this procedure is shown in Fig. 5 for a three-day period. We have omitted the initial 12 h of reconstruction because it is the transient period of convergence for both the UKF and parameter estimation.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g005.jpg

Open in a separate window

FIG. 5.

Parameter estimation for tracking of circadian dynamics. Noisy measurements of F $_{L C}$ from the full FBFD model were assimilated with a version of the DBB model that represented input from the SCN as a quasi-static parameter $C_{G, S C N}$ whose value was estimated and tracked in overlapping 30-min long windows (6 minute overlap). Two upper panel: Short excerpts of reconstructed dynamics for various phases of the circadian cycle are shown for F $_{L C}$ and F $_{V L P O}$ where black, blue, and red traces indicate true model values, the observation, and the UKF reconstruction, respectively. Third panel: Estimated (red) and true (black) value of the tracked parameter $C_{G, S C N}$ are shown. Note that the tracked parameter value is estimated with inherent smoothing on the timescale of a 30 min and does not reconstruct the feedback of the sleep–wake regulatory cell groups to the SCN. Normalized reconstruction error for various variables is shown in the lower three panels. The reconstruction of unobserved variables F $_{W R}$ and homeostatic sleep drive h is quite good as indicated by small $ϵ^{2}$ values.

SCN modulates the overall sleep cycles such that the dynamics have more frequent sleep periods in the light cycle and dominant longer wake bouts in the dark period. Short-time example time series for $F_{L C}$ is shown in the upper panel for different phases of the circadian cycle. The original DBB model lacks the SCN associated variables and their effect on the sleep network. However, here, we modified the DBB model by adding a single quasi-static parameter $C_{G, S C N}$ to represent the input contribution of SCN’s GABAergic projections to the sleep network. This parameter is added to other neurotransmitter variables in Eq. (A2). We then use observations of $F_{L C}$ from the full FBFD model to estimate the parameter.

The estimated value for $C_{G, S C N}$ (red) and the true model-generated time series for SCN (black) are shown in the third panel of Fig. 5. The reconstructed $C_{G, S C N}$ is an estimate with inherent averaging over window length including at least one sleep–wake cycle (half-hour periods). Therefore, it does not reproduce the fast dynamics of the model-generated SCN input to the sleep network. However, our estimate tracks the mean value of the true $C_{G, S C N}$ quite well. In addition, all other model variables are reconstructed accurately. Examples of the normalized reconstruction error, averaged over the fitting windows, are shown in the three bottom panels of Fig. 5 for sample variables. Note that even the homeostatic sleep drive $h$ , a variable with no direct coupling to the observed variables, is reconstructed well over most of the day. This is because of the simultaneous parameter estimation and state reconstruction. As we estimate $C_{G, S C N}$ , the Kalman gain weighs all the variables and improves estimates of all the variables including the unobserved ones.

4. Simultaneous estimation of fast and slow parameters

Experimental measurements from normal or diseased rodents are acquired in normal conditions with 12-h light–dark cycles. The true sleep–wake regulation in the animals includes both the faster variables governing state transitions, neuronal firing rates, and neurotransmitter concentration of the cell groups, as well as the slower circadian drive from SCN. To utilize the model with experimental measurements, we need to fit the parameters for both of these mechanisms. Therefore, we combined the Levenberg–Marquardt parameter estimation algorithm (demonstrated in Fig. 3) with the dynamical parameter tracking routine (demonstrated in Fig. 5).

We first generated 24 days of data from the FBFD model with the SCN circadian drive and used $F_{L C}$ and $F_{R}$ as observables. From this dataset, we selected a fixed hour from each day when SCN is relatively constant (12–1 p.m.), accumulating a 24-h long ensemble. We used the 24-h long dataset with constant SCN to optimize the fast timescale parameters that govern state transitions, as described in Sec. III C 1.

Based on the parameter correlation results, we optimized parameters listed in Fig. 6 with the simultaneous state and the multiple parameter estimation algorithm.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g006.jpg

Open in a separate window

FIG. 6.

Optimized fast timescale parameters. We used a subset of data with relatively constant SCN dynamics to estimate fast timescale parameters. Based on parameter correlation calculations, described in Sec. III C 2, a subset of these parameters were optimized. Parameters included in the steady state firing rate function, $V L P O_{max}$ , $L C_{max}$ , $D R_{max}$ , $R_{max}$ , $W R_{max}$ , $α$ , and $β$ values were fixed using the original model values. The parameters determine relative maximum firing rates of each cell group and are not significantly modulated by sleep dynamics.

Next, we used 3 days of the original 24 days of data when circadian drive and effects of SCN are present with the optimized fast timescale parameters to estimate the SCN dynamics, as described in Sec. III C 3.

This procedure, similar to fast–slow decomposition,²⁵ provides us with both fast and slow timescale parameters as a function of time of day, which can then be used for DA state reconstruction and forecasting.

In Sec. III E, we combine this algorithm with an observation function to estimate parameters and forecast animal SOV from experimental recordings.

D. Reconstruction from hypnogram

In the DBB or FBFD model of sleep–wake regulation, the variables are firing rates from small and deeply embedded cell groups such as VLPO and LC. To use the DA framework introduced here with real experimental data, we need to acquire measurements of neuronal activity from the cell groups represented in the model. However, instrumenting freely behaving animals to simultaneously measure neuronal activity in such brain regions is technically challenging.²³ To remedy these circumstances, we can still utilize DA but with measurements that only remotely resemble variables of the model: measurements of state of vigilance (SOV). The feasibility of this framework is shown using model-generated SOV. In Sec. III E, we apply it to reconstruct and forecast system dynamics with real data. The procedural steps of the simultaneous state and the multiple parameter estimation algorithm are depicted in Fig. 7.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g007.jpg

Open in a separate window

FIG. 7.

Procedural steps of the DA-UKF algorithm. The data assimilation framework receives the observed data, whether noisy model-generated data or experimental SOV, makes forecasts, corrects its predicted SOV trajectory and model parameters according to observed data, and then repeats the process.

We first sleep-score the model-generated data by assigning an SOV to each time point. We determine the SOV based on relative values of $F_{L C}$ , $F_{V L P O}$ , and $F_{R}$ . From the SOV data, we compute the probability distribution functions (PDFs) of the variables $F_{L C}$ , $F_{V L P O}$ , and $F_{R}$ conditioned on each SOV and time spent in it, $P (F_{i} | (S O V, t - t_{S O V_{o n s e t}}))$ .

In the DA-UKF routine, we generate a new dataset that we also sleep-score. The UKF reconstruction step requires that the observation function derived from the measured values—here SOV as a function of time—to include values and error estimates for each variable. To translate observed SOV to inputs for model variables in the UKF, we use the mean of conditional PDFs $P (F_{i} | (S O V, t - t_{S O V_{o n s e t}}))$ for each variable, $i$ , and each SOV and then use the standard deviation of the distribution as uncertainties. Therefore, we utilize observations of SOV to infer observations of the model variables. The SOV-derived observations for the DBB model dynamics are shown in Fig. 8.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g008.jpg

Open in a separate window

FIG. 8.

Hypnogram derived observations for the FBFD model dynamics. Mean (blue traces) and variance (red traces) for $F_{L C}$ , $F_{V L P O}$ , and $F_{R}$ from model-generated data as a function of time from state transitions into identified SOV. Because state and time since last transition can be causally measured from hypnograms, we use these means and variances as inferred observation functions to translate hypnogram into state- and state-transition conditioned observations.

1. Parameter estimation from hypnogram

We first demonstrate the feasibility of SOV-derived observation function by simultaneous estimation of model state, multiple fast timescale parameters, and SCN dynamics using discrete SOV measurements as observations. The procedure is as described in Sec. III C 4: a 24-hour long dataset with relatively constant SCN dynamics is utilized to estimate fast timescale parameters. Then, a multiple day long dataset and the optimized fast timescale parameters are used to estimate SCN dynamics.

Shown in Fig. 9 are the results of the multi-parameter estimation procedure with the same initial conditions for unknown parameters $g_{A, L C}$ , $g_{G, W R}$ , and $g_{S, R}$ as in Fig. 3; only here, we used SOV-derived observations of the firing rate. After a short convergence period, the estimated parameter values approach the parameters used to generate the data. Likewise, the reconstruction error in variables directly linked to estimated parameters, i.e., $F_{L C}$ , decreases as the parameters approach their correct value.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g009.jpg

Open in a separate window

FIG. 9.

Hypnogram derived simultaneous state tracking and multiple parameter fitting. We chose three different parameters: g $_{S, R}$ , g $_{A, L C}$ , and g $_{G, W R}$ for fitting while observing F $_{L C}$ , F $_{V L P O}$ , and F $_{R}$ derived from the hypnogram. Parameter estimation is performed with the multiple-shooting method described in Sec. III C 3 over 20-min long windows, with a 6-min overlap. The trajectories for the short model generated (cyan), reconstructed (red), and true (black) F $_{L C}$ and F $_{V L P O}$ dynamics for different periods of the convergence of the three parameters are shown in the first two panels. We start with random initial values for each parameter and estimate its correct value at each step: g $_{S, R}$ (cyan), g $_{A, L C}$ (magenta), and g $_{G, W R}$ (green). Dashed lines indicate the correct value of the parameters. Note that initially, for parameter values significantly different from the true value, the short trajectories diverge quickly from the reconstructed values, and the reconstructed values of F $_{L C}$ and F $_{V L P O}$ are different from the true values. When the parameters approach their true values, both short model-generated and reconstructed trajectories approach their true values as well and the cost function (CF) approaches zero (black trace). Reconstruction metric $ϵ_{i}^{2}$ is computed for each data assimilation window and is shown in the last two panels for F $_{L C}$ (green) and F $_{V L P O}$ (red). As a reference, the variance of the noisy observations of F $_{L C}$ and F $_{V L P O}$ is also shown in dashed green and red lines. $ϵ_{i}^{2}$ decreases as the parameters and states approach their correct values.

Once the DBB-relevant parameters were estimated correctly, we utilized them and the SOV-based observations to estimate the SCN dynamics. As demonstrated in Fig. 10, the estimated SCN (magenta) approached the true values (black). However, the convergence for SCN dynamics—using the SOV-based observations—was not as accurate as UKF estimation using direct observations of model variables. The normalized reconstruction error for variables directly linked to the estimated parameters was similarly affected. This error might be originated from the UKF algorithm attempting to constrain the observed variables to the PDF values mapped from the SOV. Similarly, the parameter estimation algorithm attempts to minimize the error between model forecasts and the UKF reconstructions of the observed variables. Therefore, any difference between model reconstruction and observed variables will certainly affect the reconstruction and forecasting.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g010.jpg

Open in a separate window

FIG. 10.

Hypnogram derived simultaneous state tracking and multiple parameter fitting for SCN. The observations are the hypnogram derived F $_{L C}$ and F $_{V L P O}$ . After estimating the fast parameters with a constant SCN, we utilized the estimated parameter values to estimate SCN over longer time-periods. The estimated $C_{G, S C N}$ (magenta) follows the true value relatively well (black). This results in good reconstruction of the other model variables such as the firing rate during REM (indicated by REM) and the homeostatic sleep drive (H). Observation of F $_{L C}$ (Wake) and F $_{V L P O}$ (NREM) are indicated by green traces.

E. Application of SOV reconstruction and forecasting to experimental data

We developed a multi-variate parameter estimation algorithm and utilized a smoothing routine to simultaneously estimate multiple model parameters and reconstruct model state including the random thalamic input $δ$ . We have shown the capability of our data assimilation framework to forecast model state using model-generated data. In this section, we demonstrate the performance of our framework when applied to experimental data. This framework once validated can be utilized in a variety of experimental settings to forecast SOV, detect hidden dynamics, for example, in a diseased brain, and provide educated guesses at suitable interventions.

The hypnogram used in this section as observation is derived from experimental physiological and behavioral measurements from chronically implanted rodents.²³ Here, we use the model derived observation function described in Sec. III D that translates hypnogram data to model variables to observe and forecast SOV from hypnograms measured from animals in long-term recording experiments. As in Sec. III D 4, we used a subset of experimental SOV data where SCN was relatively constant to estimate parameter values for the representative fast timescale parameters (listed in Fig. 6). Then, we used multiple day long subset of the experimental SOV to reconstruct SCN dynamics. The cost function (CF) of the multiple fast timescale parameter estimation algorithm over the first 48 h of simultaneous state reconstruction and parameter fitting is demonstrated in Fig. 11. As with model-generated data, the CF decreases and settles at a low value as the fitting converges.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g011.jpg

Open in a separate window

FIG. 11.

Reconstruction results using DA-UKF and experimental SOV from rodents. The model predictions are corrected at every measurement point. The cost function, similar to Fig. 9, is calculated as the difference between original model values and the ones derived from animal SOV. The model predictions converge and seem to stabilize after about 12 h.

We tracked model parameters—including the SCN dynamics—for multiple days and then averaged across days to get time-of-day averages for these tracked values. We then used these time-of-day dependent parameter sets to assimilate and forecast SOV for out-of-sample hypnogram data from the same animal. Shown in Fig. 12(a) is the time-of-day dependence of the extracted SCN value.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g012.jpg

Open in a separate window

FIG. 12.

Reconstruction and forecast of SCN and SOV from animal data. (a) SCN dynamics were predicted within the DA-UKF framework with experimentally measured animal SOV. The SCN was reconstructed and averaged over 5 days to form a 24-h long time-series. The yellow and gray backgrounds indicate the animal day and night, respectively. (b) Model forecast for out-of-sample data is shown in red for 1 h in the animal’s day when the animal mostly sleeps (left panel) and for 1 h in the animal’s night when the animals are mostly awake (right panel). The original experimental SOV and firing rate traces derived from the SOV are shown in black. The firing rates shown are $F_{L C}$ from wake-active LC, $F_{V L P O}$ from NREM-active VLPO, and $F_{R}$ from REM-active PPT/LDT cell groups. The model is set to forecast every 5 s with a fixed update rate of 5 s. (c) Forecast error (FE) as a function of forecast time. FE was computed as the mean value of the binary test that the forecast SOV equaled the experimentally measured SOV. FE was computed for 24 h of out of sample data. FE increases as the forecast horizon becomes longer.

The model forecasts 5 s ahead of the most recent observation. The observation (red trace) overlaid on the ground truth (black trace, experimental SOV, and observations derived from it) is shown in Fig. 12(b) for the NREM-active, REM-active, and wake-active cell groups as well as the SOV. The two panels represent 1 h long traces during light and dark periods. The model forecast compares favorably to the experimental observation of SOV. The forecast accuracy is good during any given state away from state-transition points. However, the forecast can be off by a few seconds during state-transition times.

After any measurement time, the model state is updated, and forecasts can change. We expect that forecasts would degrade as forecast time increases. To quantify that, we computed the mean forecast error (FE) as the average of the binary test that the forecast SOV differed from the actual SOV, as a function of forecast time, averaged over a full day of out-of-sample data. FE is shown in Fig. 12(c) and is degraded to $\approx 50 %$ for forecasts $\approx 10$ min in the future.

IV. DISCUSSION

We introduced the application of data assimilation in the study of sleep–wake regulation. Sleep–wake regulation is a complex phenomenon that includes a variety of chemical and neuronal interactions, some of which are still unknown. Data assimilation (DA) is a valuable tool in the study of such complex systems when measurements are often incomplete, noisy, or both. It utilizes all available information including experimental observations and short-term model-generated predictions of a system to build more accurate forecasts.

Although DA has been successfully applied to study neuronal dynamics before,^26–32 other than the previous work from our group³ confirming its applicability, DA framework and UKF have not been utilized in sleep modeling research and have not been applied to experimental data within the sleep research field. Our main goal is to use physiologically based mathematical models of sleep–wake regulatory system (SWRS) to gain a better understanding of the underlying mechanisms of sleep–wake regulation as well as form forecasts of the dynamics in normal and diseased brain. Within our framework, we take advantage of both computational and experimental data available to us. These include the model dynamics that represent the interactions within the SWRS and the experimental observations of the sleep–wake regulatory network.

We utilized DA with a computational model of SWRS to (1) use information from unmeasured variables, (2) estimate parameters specifically for a subject even when a parameter is associated with an unobserved variable, (3) allow for uncertainty or noise in a model structure or measurements, and (4) forecast future dynamics.

A. Application of smoothing algorithms for unmeasured variables

In models of sleep–wake regulation, as with many biological models, there exists an input that is not deterministic within the model. In the presence of random thalamic noise, $δ$ , the DA routine fails to synchronize the model with experimental observations. This is because we cannot measure $δ$ and observations only reflect its effect. Previously, Sedigh-Sarvestani et al. attempted to increase reconstruction accuracy of $δ$ by adding a large CI value to the covariance of the sigma points assigned to $δ$ .³ This approach essentially enlarges the ball of sigma points and allows for a larger span of model trajectories. However, if the model dynamics are completely blind to $δ$ , widening the sigma points for the current estimate will only result in small improvements, if any.

To reduce the reconstruction error of random thalamic noise variable, $δ$ , we utilized fixed-lag URTS smoother. The smoothing algorithms are concerned with the problem of finding maximum likelihood estimates of $x (t)$ over time intervals of $t_{0} < t < t_{k}$ .³³ The interpolation with delay can reduce the covariance of the error because more variables have been observed during the delay.³⁴ Because the fixed-lag URTS smoother can use the state dynamics after $δ$ occurrence, we can interpolate the prior $δ$ variable and re-adjust filter and forecast trajectories. By embedding a smoothing algorithm within the UKF routine, we can thus provide a better estimate of a fully unobservable variable, i.e., $δ$ and improve reconstruction fidelity of other observed and unobserved states. Although smoothing algorithms have been used in many control theory problems, we are unaware of any within a neuronal dynamics field.

B. Application of parameter estimation tools for subject-specific estimations

Mathematical models of systems almost always depend on a set of parameters that govern state interactions and dynamics. In most cases, certainly with models of biological systems, parameter values are either unknown or only known imprecisely. Correct parameter estimation is crucial for accurate model reconstruction or forecast. We adopted a multiple-shooting parameter estimation method that computes the difference between the short model-generated forecasts and the UKF-reconstructed trajectories over time windows that are long enough to explore the entire state space. The estimation step involves the minimization of a least-squared error, which also maximizes likelihood. Because we simultaneously update parameter estimates by minimizing the divergence of trajectories that are reconstructed using previous parameter estimates, our simultaneous parameter and state-estimation algorithm become an expectation-maximization method, with all the associated global optimization implications.^35,36

Any parameter estimation method is inherently limited by the identifiability of the model. Identifiability is a structural property of any dynamical system that is defined as the ability to find a unique set of correct parameters, given noise-free observations of the dynamics.³⁷ For structurally unidentifiable parameters, no parameter estimation method will prevail. We detected subsets of unidentifiable parameters by utilizing the model’s sensitivity matrix. We found that if we keep those parameters constant, multiple-shooting methods will converge reasonably well for combinations of identifiable parameters.

C. Application of DA to reconstruct hidden dynamics

Data assimilation permits us to use models that may not represent full dynamics of the system and successfully assimilate experimental data and estimate both observed and unobserved dynamics. In the example of Fig. 11, we used a model without any circadian drive (DBB model) to correctly estimate a 24-h long sleep–wake cycle as well as the underlying interaction with the SCN. Our results indicate that our DA framework can utilize inadequate models to provide insight into unobserved dynamics that are outside the scope of the model’s governing equations. We anticipate this to be of use for identification of patterns that might be present as neurological co-morbidities in a diseased brain.

We have established a framework to mechanistically study sleep–wake regulation by incorporating experimental measurements into physiologically based mathematical models. To effectively utilize this DA-based method, several points must be considered. First, initial values for filter parameters such as process noise and covariance matrices should be estimated off-line via iterative state reconstruction and parameter estimation. Second, for validation purposes, often, we do not have access to all model variables. At minimum, we can utilize the hypnogram to validate the UKF’s predictions of state transitions. Previously, we have developed a system to semi-automatically score the behavioral state of a freely moving animal in real time and form the hypnogram with a resolution of a few seconds.²⁴

Understanding the role of sleep–wake regulation in a normal or diseased brain can drastically improve and even change our approaches in proposing effective treatments. We have established a framework to mechanistically study these dynamics by incorporating sparse experimental measurements into a relatively high-dimensional physiologically based mathematical model. We demonstrated that once the error between the tracked and its corresponding true state is minimized, we can reliably reconstruct the unobserved system variables (Fig. 3). Furthermore, we confirmed that smoothing algorithms can improve reconstruction of dynamics in the presence of unobservable and uncontrollable inputs such as $δ$ . We have demonstrated a parameter estimation method that allows us to track non-stationary model parameters and accommodate slow dynamics that are not included in the UKF model such as circadian drive input from SCN (Fig. 5). In addition, we have introduced a Jacobian based routine embedded in the parameter estimation method to allow for more efficient estimation of parameters that are only active in specific areas of the state space. Finally, we have demonstrated that we can even use discretized SOV, which is not one of the model variables, as the observable to successfully reconstruct model state.

The simultaneous parameter estimation and the state tracking framework we have presented here can be applied to any system identification and validation problem. Both our parameter estimation and reconstruction methods are based on minimization of prediction errors. We can use a similar method to compute the divergence of the predictions from UKF-based tracking of system dynamics with different models.

Our efforts will introduce new avenues in model validation, as well as improved understanding of network interactions involved in sleep–wake transitions. We anticipate that when coupled with mechanistic models that embody physiological dynamics of brain state, data assimilation methods will allow for more robust neural prosthetic controllers. Although initial validation of these models will require costly measurements (i.e., single unit activity), once validated, they can be inspected to identify the least costly observables that would yield sufficiently accurate estimations and hence control algorithms. The elimination of highly invasive and damaging measurements together with the capability of data assimilation methods to provide a mechanistic explanation of the underlying dynamics can then offer a platform to further identify unstable phenomena in different pathologies and thus propose effective intervention.

ACKNOWLEDGMENTS

We thank Myles W. Billard for assistance with the design and development of the acquisition hardware and experiments with rodents. This work was supported by the National Institutes of Health (NIH) grant (No. R01EB019804) and a doctoral Academic Computing Fellowship from the Pennsylvania State University to F.B.

APPENDIX A: DINIZ BEHN AND BOOTH (DBB) MODEL OF SLEEP

The DBB model of the sleep–wake regulatory system consists of three sets of neural masses representing activity in five distinct neuronal populations: wake promoting groups of LC and DR, the NREM promoting group at VLPO, and two groups in the LDT/PPT populations, one that is REM promoting (R) and one that is active during both REM and wake (WR). The cell groups’ interactions are modeled via neurotransimitter release dynamics: LC transmits norepinephrine (NE), DR transmits serotonin (SR), the two groups in the LDT/PPT transmit acetylcholine (ACh), and VLPO transmits GABA.

Information processing of any sort in brain, i.e., sleep–wake regulation or memory consolidation, is carried out by large groups of interconnected neurons. The activity of these groups can be modeled either by ensemble models that comprise thousands of formal neuron models, such as the Hodgkins and Huxley model of a neuron or by neural mass models. Neural mass models explain the activity of populations of neurons by reducing their dynamics to a probabilistic distribution function. This function captures the likely distribution of neuronal states at a given time. In turn, this can be further reduced to a single variable describing the mean firing rate. For further detail on neural mass models, see Ref. ³⁸.

Within the DBB model, h represents the homeostatic sleep drive, and the random excitatory thalamic input is modeled by the variable $δ$ . The thalamic input represents a random sensory input such as loud noises that could wake the subject up. State of vigilance is then determined by evaluation of the firing rates of LC, DR, VLPO, and LDT/PPT cell groups with respect to a threshold.

The activity of each cell group is described by its firing rate $(F)$ and the accumulation and dissipation of the neurotransmitter $(C)$ that it releases to the post-synaptic populations and are modeled as first-order chemical kinetics with time constant $τ_{X}$ ,

{\dot{F}}_{X} = \frac{F_{X \infty} (I_{X}) - F_{X}}{τ_{X}},

(A1)

where the first-order time derivative is indicated by a dot. Here, $I_{X}$ is a weighted sum of neurotransmitter release $i$ that is affecting cell group firing rate $F_{X}$ , with coupling constants $g_{i, X}$ ,

\begin{aligned} I_{L C} = g_{A, L C} C_{A} - g_{N, L C} C_{N} - g_{G, L C} C_{G} + δ, \\ I_{D R} = g_{A, D R} C_{A} - g_{S, D R} C_{S} - g_{G, D R} C_{G} + δ, \\ I_{V L P O} = - g_{N, V L P O} C_{N} - g_{S, V L P O} C_{S} - g_{G, V L P O} C_{G}, \\ I_{R} = g_{A, R} C_{A} - g_{N, R} C_{N} - g_{S, R} C_{S} - g_{G, R} C_{G}, \\ I_{W R} = g_{A, W R} C_{A} - g_{G, W R} C_{G}, \end{aligned}

(A2)

where $A$ represents Ach, $N$ represents NE, $G$ represents GABA, and $S$ represents SR.

The steady state firing rate, $F_{X \infty} (I_{X})$ , is given by maximum firing rate parameter $X_{max}$ times a sigmoidal function with midpoint $β_{X}$ and slope $α_{X}$ ,

F_{X \infty} (I_{X}) = X_{max} (0.5 \{1 + \tanh [\frac{I_{X} - β_{X}}{α_{X}}]\}) .

(A3)

The sigmoidal halfpoint for VLPO $β_{V L P O}$ is proportional to the homeostatic sleep drive $h$ . For all other cell groups, $β_{X}$ is constant. The concentration of the neurotransmitter in the post-synaptic space is also modeled to evolve according to a first-order kinetics with time constant $τ_{i}$ and steady state neurotransmitter concentration $C_{i \infty}$ given by

\begin{aligned} {\dot{C}}_{i} = \frac{C_{i \infty} (F_{X}) - C_{i}}{τ_{i}}, \end{aligned}

(A4)

\begin{aligned} C_{i \infty} = \tanh (\frac{F_{X}}{γ_{i}}), \end{aligned}

(A5)

where $γ_{i}$ is an adjustable scale parameter and $F_{i}$ denotes the sum of the firing rates of cell groups that produce species $i$ .

Projections from thalamocortical circuits are thought to contribute noisy excitatory input to the wake-active populations LC and DR. This is modeled through state variable $δ$ and implemented as Poisson distributed impulses $P (t)$ with a rate of 0.003 Hz and leaky integrator dynamics with decay constant $τ_{δ} = 10$ s,

\dot{δ} = P (t) - δ / τ_{δ} .

(A6)

The homeostatic drive variable $h$ regulates the duration of sleep and wake bouts by changing the sigmoidal midpoint $β_{V L P O}$ for firing of the NREM-active VLPO cell group and is thought to derive from cortical metabolite and accumulation and dissipation⁴ and is modeled as

\begin{aligned} \dot{h} & = H [(F_{L C} + F_{D R}) - θ_{w}] \frac{1 - h}{τ_{h w}} \\ - H [θ_{w} - (F_{L C} + F_{D R})] \frac{h}{τ_{h s}}, \end{aligned}

(A7)

where $H$ is the Heaviside function, $τ_{h s}$ and $τ_{h w}$ determine accumulation and dissipation rates, and switching from accumulation to dissipation occurs by comparison of the sum of wake-active firing rates $F_{L C}$ and $F_{D} R$ to threshold $θ_{w}$ .

APPENDIX B: FLESHNER, BOOTH, FORGER, AND DINIZ BEHN (FBFD) MODEL OF SLEEP

The DBB model does not include the brain’s circadian drive from SCN. In an extension of the DBB model, Fleshner et al. introduced an additional GABAergic cell-group representing SCN.⁸ The firing rate dynamics of SCN follow the same dynamics as the original DBB model, with input $I_{S C N}$ composed of synaptic input $S Y N (t)$ and 24-h periodic circadian input (CIRC). The only difference is that the SCN has an enforced 24-h circadian cycle (CIRC), with higher activity during the light period and lower activity during the dark period. Direct excitatory projections from LC, DR, and LDT/PPT cell groups to the SCN increase its activity during wake and REM. VLPO indirectly affects SCN by inhibiting wake and REM-active cell groups, decreasing its firing rate. In a dynamical feedback loop, SCN activity inhibits LC, DR, and LDT/PPT while exciting VLPO.

In addition to its own cyclic dynamics enforced by CIRC, the variable $S Y N$ represents the synaptic inputs to SCN from release of SR and ACh. This combination is modeled by the sum of $C I R C$ and $S Y N$ in $I_{S C N}$ . Note that the amplitude of $S Y N$ is smaller than the amplitude of $C I R C$ . However, because it represents a core sleep–wake regulatory system responsible for minute timescale state transitions, its oscillation timescale is faster than that of $C I R C$ ,

\begin{aligned} I_{S C N} = C I R C (t) + S Y N (t), \end{aligned}

(B1)

\begin{aligned} C I R C (t) = \sin (2 π (t - t_{d a y b r e a k}) / τ_{d a y}), \end{aligned}

(B2)

\begin{aligned} S Y N (t) = g_{A, S C N} C_{A} + g_{S, S C N} C_{S} . \end{aligned}

(B3)

Here, $τ_{d a y} = 24$ h. In order to match the nocturnal behavior and imposed light–dark cycle of our animals, we have shifted the phase of $C I R C$ from its original value in Ref. ⁸ by adding $t_{d a y b r e a k} = 6$ . GABAergic mediated projections from the SCN on the sleep–wake network are modeled by the neurotransmitter concentration denoted as $C_{G (S C N)}$ , which adds to the input of the LC, DR, VLPO, and R firing rates and modify Eq. (A2) from the DBB model to

\begin{aligned} I_{L C} & = g_{A, L C} C_{A} - g_{N, L C} C_{N} - g_{G, L C} C_{G} - g_{G (S C N), L C} C_{G (S C N)} + δ, \\ I_{D R} & = g_{A, D R} C_{A} - g_{S, D R} C_{S} - g_{G, D R} C_{G} - g_{G (S C N), D R} C_{G (S C N)} + δ, \\ I_{V L P O} & = - g_{N, V L P O} C_{N} - g_{S, V L P O} C_{S} - g_{G, V L P O} C_{G} \\ - g_{G (S C N), V L P O} C_{G (S C N)}, \\ I_{R} & = g_{A, R} C_{A} - g_{N, R} C_{N} - g_{S, R} C_{S} - g_{G, R} C_{G} - g_{G (S C N), R} C_{G (S C N)}, \\ I_{W R} & = g_{A, W R} C_{A} - g_{G, W R} C_{G} . \end{aligned}

(B4)

Because SCN has a much slower oscillation timescale on the order of 24 h, the output of the FBFD model on short timescales is similar to that of the DBB model in terms of typical bout durations as well as cycle times through states. However, as is typical for real rats, these characteristics change within diurnal timescales. As observed with nocturnal rats, REM and NREM are primarily observed during the putative light phase, while long periods of wake are observed during the putative dark phases.

APPENDIX C: DATA ASSIMILATION—UNSCENTED KALMAN FILTER

Data assimilation is the process where observed data are combined with mathematical models to produce an estimate of the system state. The Kalman filter embodies a data assimilation framework. It uses a prediction–correction scheme, where the predictions come from the mathematical equations that govern the model, available to the Kalman filter, and the correction comes from sparse, noisy measurements. Kalman’s initial algorithm was designed for linear systems.³⁹ The unscented Kalman filter (UKF) is an ensemble version developed to address nonlinear systems without linearization.^13,40,41

The UKF algorithm is described in detail in many standard textbooks.²⁰ For the purposes of this article, we present a brief overview. As shown in Fig. 13 based on the model equations, we start from a set of initial conditions ${\hat{x}}_{0}$ at time $t_{0}$ . Then, the UKF forms a series of points surrounding each state value, called sigma points. Given a $D$ dimensional state space for the model dynamics $F$ , we choose $2 D_{x}$ sigma points (using Cholesky decomposition) with cross-covariance ${\hat{P}}_{x x, k}$ to represent state uncertainty. UKF then generates a series of forecast trajectories by iterating the sigma points through the nonlinear model equations. The $i$ th sigma point is denoted $χ_{i, k}$ prior to iteration and ${\tilde{χ}}_{i, k + 1}$ after iteration. The model forecast for the states ${\tilde{x}}_{k + 1}$ is then the mean of the forward iterated sigma points,

{\tilde{x}}_{k + 1} = \frac{1}{2 D_{x}} \sum_{i = 1}^{2 D_{x}} {\tilde{χ}}_{i, k + 1} .

(C1)

Similar to ${\hat{P}}_{x x, k}$ , the forecast uncertainty ${\tilde{P}}_{x x, k + 1}$ is the cross-covariance of sigma points after iteration, ${\tilde{χ}}_{i, k + 1}$ , plus an additive covariance inflater matrix $C I$ ,

{\tilde{P}}_{x x, k + 1} = \frac{1}{2 D_{x}} \sum_{i = 1}^{2 D_{x}} ({\tilde{χ}}_{i, k + 1} - {\tilde{x}}_{k + 1}) {({\tilde{χ}}_{i k + 1} - {\tilde{x}}_{k + 1})}^{T} + C I .

(C2)

FIG. 13.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000031-013139_1-g013.jpg

Open in a separate window

Unscented Kalman Filter (UKF) algorithm and output. UKF uses a prediction–correction scheme, where the predictions come from the mathematical equations that govern the model, given to the UKF, and the correction comes from measurements. UKF starts with a set of initial conditions that it propagates through the model equations and builds a set of forecast trajectories (shown in red). This process will continue until a measurement point. At that point, the UKF compares the forecast trajectories with the measurement value. Based on the error, it either widens or tightens the cloud of points and builds the next set of forecast trajectories (magenta). As time goes on, the error in the trajectories diverges unless they are corrected based on the measurements.

The CI matrix is nonzero only on the diagonal and is added to compensate underestimates of the forecast uncertainty due to process noise and inadequacies in the model.^42,43 The UKF forecast is then corrected with respect to the measurement $y_{k + 1}$ at time $t_{k + 1}$ . $y$ need not contain the same number of variables as $x$ . All variables will be corrected using a correction factor that weighs the observation and forecast according to the Kalman gain $K_{k + 1}$ ,

{\hat{x}}_{k + 1} = {\tilde{x}}_{k + 1} + K_{k + 1} (y_{k + 1} - {\tilde{y}}_{k + 1}),

(C3)

where ${\tilde{y}}_{k + 1}$ denotes the mean of estimated sigma points for the observed variables, ${\tilde{Y}}_{i, k + 1}$ ,

{\tilde{y}}_{k + 1} = \frac{1}{2 D_{x}} \sum_{i = 1}^{2 D_{x}} {\tilde{Y}}_{i, k + 1} .

(C4)

The Kalman gain $K_{k + 1}$ is the ratio,

K_{k + 1} = {\tilde{P}}_{x y, k + 1} {({\tilde{P}}_{y y, k + 1})}^{- 1},

(C5)

where ${\tilde{P}}_{x y, k + 1}$ and ${\tilde{P}}_{y y, k + 1}$ are cross-covariances of the sigma points either in the full dimensional space of $x$ or in the subspace spanned by the measurements,

\begin{aligned} {\tilde{P}}_{x y, k + 1} = \frac{1}{2 D} \sum_{i = 1}^{2 D} ({\tilde{χ}}_{i, k + 1} - {\tilde{x}}_{k + 1}) {({\tilde{Y}}_{i k + 1} - {\tilde{y}}_{k + 1})}^{T}, \end{aligned}

(C6)

\begin{aligned} {\tilde{P}}_{y y, k + 1} = \frac{1}{2 D} \sum_{i = 1}^{2 D} ({\tilde{Y}}_{i, k + 1} - {\tilde{y}}_{k + 1}) {({\tilde{Y}}_{i, k + 1} - {\tilde{y}}_{k + 1})}^{T} + R . \end{aligned}

(C7)

Measurement uncertainty is represented here with $R$ . The forecast uncertainty ${\hat{P}}_{x x, k + 1}$ is computed using the Kalman gain,

{\hat{P}}_{x x, k + 1} = {\tilde{P}}_{x x, k + 1} - K_{k + 1} {\tilde{P}}_{x y, k + 1} .

(C8)

Within this recursive prediction–correction scheme, the UKF synchronizes the model state to the measurements, and the forecast uncertainty collapses. Because the correction factor weighs all the variables, the UKF can improve the estimate of the experimentally inaccessible variables.

APPENDIX D: SMOOTHING ALGORITHMS FOR RECONSTRUCTION

Kalman filters aim to estimate model state at a given time point based on observations available up to and at that time point, as well as the mean squared error of the estimate. These filters can then be used to compute the likelihood of certain observations and forecast future state dynamics. In contrast, smoothing algorithms aim to estimate the past state dynamics given all the available data. Smoothing allows us to interpolate missing observations/events that filtering did not previously predict and to obtain the mean squared error of the interpolated estimates. Smoothing and filtering are essentially similar algorithms one moving backwards in the dynamical time series describing the system, while the other moves forward.

In the DBB and FBFD model, the state variable $δ$ is representative of the random external sensory information that are relayed through thalamic regions that can modulate sleep–wake cycles. In both models, $δ$ is described by a random Poisson process inputted to LC and DR that causes short awakenings whenever it occurs. Because $δ$ does not receive any inputs from any model variable or parameter, it is impossible to reconstruct using our current observations. Sedigh-Sarvestani et al. attempted to reconstruct $δ$ by increasing its associated covariance inflation ( $C I$ ) value so that during reconstruction, the Kalman filter relied heavily on model equations instead of observations.³ Although that work slightly improved the model performance, the reconstruction error for short awakenings still remained large.

To reduce reconstruction error of $δ$ and thus the short awakenings, we utilized fixed-lag unscented Rauch–Tung–Striebel (FL-URTS) smoother. The FL-URTS smoothing algorithm is described in detail in various state-estimation and control theory references.^44–46 Here, we present an overview and equations needed to understand the details presented in this chapter. The smoother is embedded in the UKF loop so that the lag between smoothing and reconstruction is minimized. With this fixed-lag smoother, at time point $k$ , a fixed length history of dynamics of length $n$ of both observations and UKF reconstruction are processed to better reconstruct the history of state. The UKF is then applied to this history to further constrain or estimate the current state estimate and improve forecasts.

We start with the UKF procedure its estimates ${\tilde{x}}_{k + 1}$ and ${\hat{P}}_{x x, k + 1}$ from ${\tilde{x}}_{k}$ and ${\hat{P}}_{x x, k}$ . Then, once sufficient time has passed ( $k \geq$ fixed-lag $n$ ), we initialize the smoother,

\begin{aligned} {\hat{x}}_{k | k - 1}^{s} = {\hat{x}}_{k}, \\ {\hat{P}}_{k | k - 1}^{s} = {\hat{P}}_{k} . \end{aligned}

(D1)

We then utilized the UKF-generated sigma points for ${\tilde{x}}_{k}$ , covariance ${\hat{P}}_{x x, k}$ , and cross-covariance ${\tilde{P}}_{x y, k}$ and propagate them backwards (from time $k$ to time $k - n$ ) through the nonlinear model equations. The smoother gain $D_{k - 1}$ , similar to the Kalman gain, is the ratio between the smoother cross-covariance and covariance. The smoothed estimate is then calculated as

{\hat{x^{s}}}_{k - 1} = {\tilde{x}}_{k - 1} + D_{k - 1} ({\hat{x}}_{k}^{s} - {\hat{x}}_{k}) .

(D2)

This process continues for the entire fixed-lag duration.

Because the smoothing gain affects all the variables, smoothing can improve the estimate of the unobserved variables that are also affected by $δ$ as well.

DATA AVAILABILITY

Computational examples for the algorithms developed and/or integrated together, an example of animal sleep data on which these algorithms were applied, and a code for generating the figures from Penn State's ScholarSphere at https://doi.org/doi:10.26207/fcc5-7h57.

Notes

Note: This paper is part of the Focus Issue on Dynamical Disease: A Translational Perspective.

REFERENCES

1. Wulff K., Gatti S., Wettstein J. G., and Foster R. G., “Sleep and circadian rhythm disruption in psychiatric and neurodegenerative disease,” Nat. Rev. Neurosci. 11, 589–599 (2010). 10.1038/nrn2868 [PubMed] [CrossRef] [Google Scholar]

2. Pritchett D., Wulff K., Oliver P. L., Bannerman D. M., Davies K. E., Harrison P. J., Peirson S. N., and Foster R. G., “Evaluating the links between schizophrenia and sleep and circadian rhythm disruption,” J. Neural Transm. 119, 1061–1075 (2012). 10.1007/s00702-012-0817-8 [PubMed] [CrossRef] [Google Scholar]

3. Sedigh-Sarvestani M., Schiff S. J., and Gluckman B. J., “Reconstructing mammalian sleep dynamics with data assimilation,” PLoS Comput. Biol. 8, 1002788 (2012). 10.1371/journal.pcbi.1002788 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

4. Datta S. and MacLean R. R., “Neurobiological mechanisms for the regulation of mammalian sleep-wake behavior: Reinterpretation of historical evidence and inclusion of contemporary cellular and molecular evidence,” Neurosci. Biobehav. Rev. 31, 775–824 (2007). 10.1016/j.neubiorev.2007.02.004 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

5. Saper C. B., Scammell T. E., and Lu J., “Hypothalamic regulation of sleep and circadian rhythms,” Nature 437, 1257–1263 (2005). 10.1038/nature04284 [PubMed] [CrossRef] [Google Scholar]

6. Sakai K., “Sleep-waking discharge profiles of dorsal raphe nucleus neurons in mice,” Neuroscience 197, 200–224 (2011). 10.1016/j.neuroscience.2011.09.024 [PubMed] [CrossRef] [Google Scholar]

7. Booth V., Xique I., and Diniz Behn C. G., “One-dimensional map for the circadian modulation of sleep in a sleep-wake regulatory network model for human sleep,” SIAM J. Appl. Dyn. Syst. 16, 1089–1112 (2017). 10.1137/16M1071328 [CrossRef] [Google Scholar]

8. Fleshner M., Booth V., Forger D. B., and Diniz Behn C. G., “Circadian regulation of sleep-wake behaviour in nocturnal rats requires multiple signals from suprachiasmatic nucleus,” Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 369, 3855–3883 (2011). 10.1098/rsta.2011.0085 [PubMed] [CrossRef] [Google Scholar]

9. Diniz Behn C. G. and Booth V., “Simulating microinjection experiments in a novel model of the rat sleep-wake regulatory network,” J. Neurophysiol. 103, 1937–1953 (2010). 10.1152/jn.00795.2009 [PubMed] [CrossRef] [Google Scholar]

10. Tamakawa Y., Karashima A., Koyama Y., Katayama N., and Nakao M., “A quartet neural system model orchestrating sleep and wakefulness mechanisms,” J. Neurophysiol. 95, 2055–2069 (2006). 10.1152/jn.00575.2005 [PubMed] [CrossRef] [Google Scholar]

11. Rempe M. J., Best J., and Terman D., “A mathematical model of the sleep/wake cycle,” J. Math. Biol. 60, 615–644 (2010). 10.1007/s00285-009-0276-5 [PubMed] [CrossRef] [Google Scholar]

12. Simon D., Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches (John Wiley & Sons, Inc., 2006). 10.1002/0470045345 [CrossRef] [Google Scholar]

13. Julier S. J. and Uhlmann J. K., “Unscented filtering and nonlinear estimation,” in Proceedings of the IEEE (IEEE, 2004).

14. Voss H. U., Timmer J., and Kurths J., “Nonlinear dynamical system identification from uncertain and indirect measurements,” Int. J. Bifurcat. Chaos 14, 1905–1933 (2004). 10.1142/S0218127404010345 [CrossRef] [Google Scholar]

15. Mccarley R. W. and Hobson J. A., “Neuronal excitability modulation over the sleep cycle: A structural and mathematical model,” Science 189, 58–60 (1975). 10.1126/science.1135627 [PubMed] [CrossRef] [Google Scholar]

16. Ohno K. and Sakurai T., “Orexin neuronal circuitry: Role in the regulation of sleep and wakefulness,” Front. Neuroendocrinol. 29, 70–87 (2008). 10.1016/j.yfrne.2007.08.001 [PubMed] [CrossRef] [Google Scholar]

17. Porkka-Heiskanen T., Strecker R. E., and McCarley R. W., “Brain site-specificity of extracellular adenosine concentration changes during sleep deprivation and spontaneous sleep: An in vivo microdialysis study,” Neuroscience 99, 507–517 (2000). 10.1016/S0306-4522(00)00220-7 [PubMed] [CrossRef] [Google Scholar]

18. Huang Z.-L., Qu W.-M., Eguchi N., Chen J.-F., Schwarzschild M. A., Fredholm B. B., Urade Y., and Hayaishi O., “Adenosine A

_{2 A}

, but not A

_{1}

, receptors mediate the arousal effect of caffeine,” Nat. Neurosci. 8, 858–859 (2005). 10.1038/nn1491 [PubMed] [CrossRef] [Google Scholar]

19. Deurveilher S. and Semba K., “Indirect projections from the suprachiasmatic nucleus to major arousal-promoting cell groups in rat: Implications for the circadian control of behavioural state,” Neuroscience 130, 165–183 (2005). 10.1016/j.neuroscience.2004.08.030 [PubMed] [CrossRef] [Google Scholar]

20. Schiff S. J., Neural Control Engineering (The MIT Press, 2011). [Google Scholar]

21. Levenberg K., “A method for the solution of certain non-linear problems in least squares,” Q. Appl. Math. 2, 93–101 (1944). 10.1090/qam/10666 [CrossRef] [Google Scholar]

22. Marquardt D. W., “An algorithm for least-squares estimation of nonlinear parameters,” J. Soc. Indus. Appl. Math. 11, 0111030 (1963). 10.1137/0111030 [CrossRef] [Google Scholar]

23. Billard M. W., Bahari F., Kimbugwe J., Alloway K. D., and Gluckman B. J., “The systemdrive: A multisite, multiregion microdrive with independent drive axis angling for chronic multimodal systems neuroscience recordings in freely behaving animals,” eNeuro 5, ENEURO.0261-18.2018 (2018). 10.1523/ENEURO.0261-18.2018 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

24. Sunderam S., Chernyy N., Peixoto N., Mason J. P., Weinstein S. L., Schiff S. J., and Gluckman B. J., “Improved sleep-wake and behavior discrimination using MEMS accelerometers,” J. Neurosci. Methods 163, 373–383 (2007). 10.1016/j.jneumeth.2007.03.007 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

25. Behn C. G. and Booth V., “A fast-slow analysis of the dynamics of REM sleep,” SIAM J. Appl. Dyn. Syst. 11, 212–242 (2012). 10.1137/110832823 [CrossRef] [Google Scholar]

26. Toth B. A., Kostuk M., Meliza C. D., Margoliash D., and Abarbanel H. D. I., “Dynamical estimation of neuron and network properties I: Variational methods,” Biol. Cybern. 105, 217–237 (2011). 10.1007/s00422-011-0459-1 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

27. Hu X., Nenov V., Bergsneider M., Glenn T. C., Vespa P., and Martin N., “Estimation of hidden state variables of the intracranial system using constrained nonlinear Kalman filters,” IEEE Trans. Biomed. Eng. 54, 597–610 (2007). 10.1109/TBME.2006.890130 [PubMed] [CrossRef] [Google Scholar]

28. Quach M., Brunel N., and D’alché-Buc F., “Estimating parameters and hidden variables in non-linear state-space models based on ODEs for biological networks inference,” Bioinformatics 23, 3209–3216 (2007). 10.1093/bioinformatics/btm510 [PubMed] [CrossRef] [Google Scholar]

29. Eberle C. and Ament C., “The unscented Kalman filter estimates the plasma insulin from glucose measurement,” Biosystems 103, 67–72 (2011). 10.1016/j.biosystems.2010.09.012 [PubMed] [CrossRef] [Google Scholar]

30. Ullah G. and Schiff S. J., “Assimilating seizure dynamics,” PLoS Comput. Biol. 6, e1000776 (2010). 10.1371/journal.pcbi.1000776 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

31. Schiff S. J., “Towards model-based control of Parkinson’s disease,” Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 368, 2143–2146 (2010). 10.1098/rsta.2010.0050 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

32. Freestone D., Aram P., Dewar M., Scerri K., Grayden D., and Kadirkamanathan V., “A data-driven framework for neural field modeling,” NeuroImage 56, 1043–1058 (2011). 10.1016/j.neuroimage.2011.02.027 [PubMed] [CrossRef] [Google Scholar]

33. Brayson A. E. and Frazier M., “Smoothing for linear and nonlinear dynamics systems,” in Proceedings of the Optimum System Synthesis Conference (Flight Control Laboratory Aeronautical Systems Division Air Force Systems Command Wright-Patterson Air Force Base, Ohio, 1963), pp. 353–364.

34. Rauch H. E., Striebel C. T., and Tung F., “Maximum likelihood estimates of linear dynamic systems,” AIAA J. 3, 1445–1450 (1965). 10.2514/3.3166 [CrossRef] [Google Scholar]

35. Wu C. F. J., “On the convergence properties of the EM algorithm,” Annal. Stat. 11, 95–103 (1983). 10.1214/aos/1176346060 [CrossRef] [Google Scholar]

36. Vaida F., “Parameter convergence for EM and MM algorithms,” Technical Report, 2005.

37. Jacquez J. A. and Greif P., “Numerical parameter identifiability and estimability: Integrating identifiability, estimability, and optimal sampling design,” Math. Biosci. 77, 201–227 (1985). 10.1016/0025-5564(85)90098-7 [CrossRef] [Google Scholar]

38. Deco G., Jirsa V. K., Robinson P. A., Breakspear M., and Friston K., “The dynamic brain: From spiking neurons to neural masses and cortical fields,” PLoS Comput. Biol. 4, e1000092 (2008). 10.1371/journal.pcbi.1000092 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

39. Kalman R. E., “A new approach to linear filtering and prediction problems,” Trans. ASME J. Basic Eng. 82, 35–45 (1960). 10.1115/1.3662552 [CrossRef] [Google Scholar]

40. Julier S. J., “The scaled unscented transformation,” in Proceedings of the 2002 American Control Conference (IEEE Cat. No. CH37301) (IEEE, 2002), Vol. 6, pp. 4555–4559.

41. Julier S. J. and Uhlmann J. K., “Reduced sigma point filters for the propagation of means and covariances through nonlinear transformations,” in Proceedings of the 2002 American Control Conference (IEEE Cat. No. CH37301) (IEEE, 2002), Vol. 2, pp. 887—-892.

42. Anderson J. L. and Anderson S. L., “A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts,” Mon. Weather Rev. 127, 2741–2758 (1999). 10.1175/1520-0493(1999)127<2741:AMCIOT>2.0.CO;2 [CrossRef] [Google Scholar]

43. Miyoshi T., “The Gaussian approach to adaptive covariance inflation and its implementation with the local ensemble transform Kalman filter,” Mon. Weather Rev. 139, 1519–1535 (2011). 10.1175/2010MWR3570.1 [CrossRef] [Google Scholar]

44. Kohn R. and Ansley C. F., “Filtering and smoothing algorithms for state space models,” Comput. Math. Appl. 18, 515–528 (1989). 10.1016/0898-1221(89)90104-1 [CrossRef] [Google Scholar]

45. Briers M., Doucet A., and Maskell S., “Smoothing algorithms for state-space models,” Ann. Inst. Stat. Math. 62, 61–89 (2010). 10.1007/s10463-009-0236-2 [CrossRef] [Google Scholar]

46. Pei T., Jiachuan L., Li Z., Zhou Y., and Choi Y., “A fixed-lag unscented Rauch-Tung-Striebel smoother for non-linear dynamic state estimation,” Int. J. Digital Content Technol. Appl. 7, 769–777 (2013). [Google Scholar]

Articles from Chaos are provided here courtesy of American Institute of Physics