Modeling of EEG power spectrum over frontal and occipital regions during propofol sedation 

Get Complete Project Material File(s) Now! »

Linear delay differential equations

Delays or lags occur in many physical and biological systems, due to the limited processing speeds between non-local elements. In neuroscience, the finite propagation speed of action potentials along axonal branches, or in more general terms, the finite-time interaction between two elements in a spatially extended system is associated with a time delay. In this context, it is important to study the dynamics of differential equations with delay because most of the models describing the activity of a large population of neurons are governed by a set of coupled delay differential equations. In the current study, the time delays are treated as a constant in the system equations.
The dynamic of systems with time delay can be described by Delay Differential Equations (DDEs), also known as difference differential equations, which are a special class of differential equations that the rate of changes in the system depends on its past history. Time-delay systems and their applications can be found in many different scientific areas such as control theory, electronic circuits, nonlinear optics and lasers, neuroscience, networks and communication and so on [126, 127]. It has been extensively recognized that delays play a crucial role in mathematical modeling of real-world systems and explaining their dynamical properties. For instance, it has been shown that the delay values can determine the existence, stability and periodicity of the solutions [128, 129, 130, 131, 132, 133, 134, 135], the bifurcations of equilibria [136, 137, 138, 139], and the synchronization in such systems [140, 141, 142, 143, 144, 145]. In general, the mathematical formulation of a linear DDE is given by 8>>< >>: dY (t) dt = AY (t) + BY (t 􀀀 ); t 0; Y (t) = (t); t 2 [􀀀; 0]; (2.6).
where Y is an N1 activity variable vector, A and B are constant quadratic matrices of order N, > 0 denotes the time-delay, and (t) is a continuous function that defines the initial condition [146]. Although a large number of studies have been devoted to obtain the exact stability conditions for these equations, finding analytical solutions for a wide class of DDEs remains as an open problem [147, 148, 149, 150, 151]. The DDEs are often solved using numerical methods such as the Euler or Runge-Kutta method [146, 152], by employing a Lambert function [153, 154], the discretization of the solution operator as implemented in the software package DDE-BIFTOOL [155, 156, 157], or by discretizating the infinitesimal generator as implemented in the software package TRACE-DDE [158, 159].

Stability of delay differential equations

The solution of a DDE is asymptotically stable if and only if all the roots of the characteristic equation referred to as the characteristic roots have negative real parts. The system is called unstable if there exists a root with positive real part. For an ordinary differential equation the characteristic function is a polynomial and the well-known Routh-Hurwitz criterion can be used to determine the negativity of the real parts of characteristic roots and thus the stability of the solution [160, 161, 162]. However, the characteristic functions for DDEs are not ordinary polynomial, instead, they are named as transcendental equation or exponential polynomials or quasi-polynomials, which the number of their roots can be countably infinite [163]. Assuming a solution of Eq. (2.6) takes the form of Y (t) = Cet, where 2 C, yields the following characteristic equation det I 􀀀 A 􀀀 Be􀀀 = 0; (2.7).
where I is an N N identity matrix. This equation can be written as P() + e􀀀Q() = 0; (2.8). where and P() and Q() are polynomial in . This transcendental equation is known to have an infinite number of complex roots which, in general, cannot be expressed explicitly in terms of elementary functions. It is well-known that the system is asymptotically stable if and only if all the characteristic roots have negative real parts, i.e., Re( < 0). Many approaches have been taken to determine the stability of the transcendental equations, however, in general no closed-form solution for Eq. (2.8) is available [164]. An analytic approach to obtain the solution of linear systems of DDEs is based on the concept of the Lambert function. Every function W(x) that satisfies W(x)eW(x) = x is called a Lambert function [165]. The lamber function, W(x) has an infinite number of branches Wk(x), where k 2 􀀀1; :::;􀀀1; 0; 1; :::1. It has been shown that the solution of Eq. (2.6) in terms of the Lambert function is given by [153, 166] Y (t) = k=X􀀀1 k=􀀀1 Cke( 1 Wk(Be􀀀A)+A) ; (2.9).

Thalamo-cortical circuits

It is known that the cortex, the thalamus, and the thalamo-cortical feedback loops contribute in the EEG oscillations. However, the precise role of each structure to the dynamics of EEG is unknown. While the macroscopic cortical activities are directly reflected in EEG data, the thalamus acts as both the gate for sensory input to the cortex and the site for feedback from cortical cells [226]. It has been shown that the specific changes in EEG patterns during the brain  transition states can be determined by neural activities of the cortex, the thalamus and the interactions between them [35, 36, 37]. Therefore, a population-level model of thalamo-cortical system is a plausible candidate for investigating the mechanisms responsible for the various features observed in EEG data [186, 29, 30, 23, 91, 227].
In the last years, a large number of studies using mean-field models have been devoted to understand how EEG recordings are related to underlying brain activity [228, 23, 24, 176, 229, 35]. However, the exact neurobiological mechanisms that generate the EEG oscillatory signatures along with its relation to behavioral states remain to be fully understood. It has been shown that depend on the functional state of the brain, thalamo-cortical circuits are involved in generating different spontaneous oscillatory patterns such as oscillatory tonic and burst firing, across a wide range of frequencies from slow and spindle rhythms to gamma-frequency band [230, 228, 231, 232].
The reciprocal thalamo-cortical loops has been proposed to be responsible for normal EEG oscillations during wake and sleep states [233]. Other researches illustrate that the thalamocortical circuits crucially involve in the formation of sleep, normal physiologic activities [234, 235, 236], consciousness processing and anesthetic-induced unconsciousness [237, 44, 81, 238, 239], processing of sensory information [240, 241, 242], attention regulation [243, 244, 245, 246], and also synchronization of cortical activity [247, 248, 249, 250, 251, 252]. There are further evidences showing that a change in these circuits can impair consciousness and induces spike-and-slow wave discharges (SWD), which are essential features of the absence seizures [253, 254, 230, 228, 255]. It has also been reported that mediating of T-type calcium channels and HCN-channels, which could influence the excitation/inhibition balance in these circuits play an important role in regulating thalamo-cortical rhythmicity and absence seizures [233, 256, 236]. In addition, clinical studies have revealed that enhancing oscillations in the thalamo-cortical loops could lead to epileptic seizures [257].

Robinson model

In the 1990’s David Liley and Jim Wright started developing partial differential equation models for neural population activity [96, 97]. Their works have inspired other teams, e.g., Peter Robinson and colleagues, who developed a similar neural model which has been proven to be very successful in reproducing key features of EEG signals recorded during different behavioral states such as sleep [29, 34], anesthesia [35, 36], and seizures [176, 258].
The model incorporates key anatomic connectivities within the cortex, within the thalamus, and between these structures: cortical, thalamic, cortico-thalamic and thalamo-cortical pathways. In addition, the model is furnished with synaptic and dendritic dynamics, nonlinearity of the firing response, and finite axonal transmission speeds [259].
The body of the model is based on a population-level description of a single thalamo-cortical module [173, 29] comprising four populations of neurons, namely excitatory (E) and inhibitory (I) cortical population, a population built of thalamo-cortical relay neurons (S) and of thalamic reticular neurons (R), as shown in Fig. 2.4. This model is based on the original idea of Lopes da Silva et al., stating that the -rhythm represents the noisy thalamic input signal band-pass filtered by feedback-connected cortical and thalamic neural populations [180]. The details of the model and the nominal parameter values are taken from the previous works [29, 30, 23]. Here we just briefly describe the key concepts of the model. The average soma membrane potential denoted by Va, for a = E; I; S;R is modeled by Va(t) = X b=E;I;R;S hb(t) abb(t 􀀀 ab); (2.31).

A thalamo-cortical model to reproduce the EEG rhythms

In this section we consider a thalamo-cortical neuronal population model based on the work of Hutt & Longtin [28], which extends previous cortical models [173, 174, 23] by distinguishing the excitatory and inhibitory synapses. The model consists of a network of four populations of neurons: cortical pyramidal neurons (E), cortical inhibitory neurons (I), thalamo-cortical relay neurons (S) and thalamic reticular nucleus (R), as shown schematically in Fig. 2.8. Note that in this model the connection from population S to I is neglected for simplicity reasons. In addition, the long-range propagation of the signals has been considered by the connection between cortex and thalamus associated with a constant time delay and the model closely follows [229, 91] showing that long connections between cortico-cortical populations through white matter are not necessary to describe experimental EEG dynamics.
The underlying model considers an ensemble of neurons on a mesoscopic scale which includes two types of neurons, namely excitatory and inhibitory cells [28]. The connection elements between neurons are excitatory and inhibitory chemical synapses in which both types of synapses may occur on dendritic branches of both cell types. The present model considers AMPA=NMDA and GABAA synaptic receptors and assumes spatially extended populations of both neuron types. In the following the key concepts of the model are described and we derive a formula for the EEG power spectrum. Let us consider spatially extended populations of excitatory and inhibitory neurons on the order of a few hundred micrometers. Assume that the neural population of type a receives the incoming firing rates Pe(x; t) and Pi(x; t) at spatial location x and time t, which are originated from excitatory and inhibitory cells (or terminate at excitatory and inhibitory synapses), respectively. Then, presynaptic population firing rates arrived at excitatory (e) or inhibitory (i) synapses convert to the mean excitatory and inhibitory postsynaptic potentials V e a (t) and V i a (t), respectively, by the convolution V e;i(x; t) = Z t 􀀀1 he;i(t 􀀀 t0)Pe;i(x; t0)dt0; (2.47).

Table of contents :

Chapter 1 Neuroscience of general anesthesia 
1.1 General anesthesia
1.2 General anesthetics
1.3 Anesthetics effects on single neurons (microscopic level)
1.3.1 Molecular targets of general anesthetics
1.3.2 Neuroanatomical targets of general anesthetics
1.4 Anesthetics effects on EEG signals (macroscopic level)
1.5 Current theories of general anesthesia
1.6 Conclusion
Chapter 2 Studying neural population models of EEG activity 
2.1 Damped harmonic oscillator
2.2 Linear delay differential equations
2.2.1 Stability of delay differential equations
2.2.2 Power spectrum of delay differential equations
2.3 Neural population models
2.3.1 Neural field models
2.3.2 Thalamo-cortical circuits
2.3.3 Robinson model
2.3.4 A thalamo-cortical model to reproduce the EEG rhythms
2.4 Conclusion
Chapter 3 Modeling the anesthetic action on synaptic and extra-synaptic receptors 
3.1 Function of extra-synaptic receptors
3.2 Effect of propofol on cortical and thalamic neural populations
3.3 EEG acquisition and the experimental observations
3.4 Theoretical power spectrum
3.5 Reproducing the experimental observations
3.5.1 The role of synaptic inhibition
3.5.2 The role of extra-synaptic inhibition
3.6 Induction of delta activity in EEG rhythms
3.7 Discussion
3.8 Conclusion
Chapter 4 Modeling of EEG power spectrum over frontal and occipital regions during propofol sedation 
4.1 GABAergic inhibition in thalamic cells
4.2 Theoretical power spectrum
4.3 Experimental power spectrum
4.4 Power spectrum of full thalamo-cortical model
4.5 Reduced thalamo-cortical model
4.5.1 The role of different populations and anatomical loops
4.5.2 Frontal spectrum
4.5.3 Occipital spectrum
4.6 Discussion
4.6.1 Origin of spectral peaks
4.6.2 Multiple resting states
4.6.3 Effective sub-circuits in the cortico-thalamic model
4.6.4 Model limitations
4.7 Conclusion
Chapter 5 Spectral power fitting using stochastic optimization algorithms 
5.1 Parameter estimation (inverse problem)
5.2 Optimization problem
5.2.1 Formulating an optimization problem
5.2.2 Objective function
5.2.3 Types of optimization methods
5.3 Optimization algorithms
5.3.1 Levenberg-Marquardt algorithm (LM)
5.3.2 Particle Swarm Optimization (PSO)
5.3.3 Differential evolution (DE)
5.3.4 Metropolis-Hastings (MH)
5.3.5 Simulated Annealing (SA)
5.4 Precision of the estimates
5.4.1 Confidence regions
5.4.2 Correlation analysis
5.4.3 Sensitivity analysis
5.5 Case Studies
5.5.1 Case Study I: A stochastic damped harmonic oscillator
5.5.2 Case Study II: A stochastic linear delay differential equation
5.5.3 Case Study III: A thalamo-cortical model reproducing the EEG rhythms .
5.6 Results of the parameter estimation in Case Studies I, II, and III
5.7 Discussion
5.8 Conclusion
References

GET THE COMPLETE PROJECT