Delayed closed-loop neurostimulation for the treatment of pathological brain rhythms in mental disorders: a computational study

Mental disorders are among the top most demanding challenges in world-wide health. A large number of mental disorders exhibit pathological rhythms, which serve as the disorders characteristic biomarkers. These rhythms are the targets for neurostimulation techniques. Open-loop neurostimulation employs stimulation protocols, which are rather independent of the patients health and brain state in the moment of treatment. Most alternative closed-loop stimulation protocols consider real-time brain activity observations but appear as adaptive open-loop protocols, where e.g., pre-defined stimulation sets in if observations fulfil pre-defined criteria. The present theoretical work proposes a fully-adaptive closed-loop neurostimulation setup, that tunes the brain activities power spectral density (PSD) according to a user-defined PSD. The utilized brain model is non-parametric and estimated from the observations via magnitude fitting in a pre-stimulus setup phase. Moreover, the algorithm takes into account possible conduction delays in the feedback connection between observation and stimulation electrode. All involved features are illustrated on pathological α- and γ-rhythms known from psychosis. To this end, we simulate numerically a linear neural population brain model and a non-linear cortico-thalamic feedback loop model recently derived to explain brain activity in psychosis.

1. Introduction

Electrical neurostimulation is an old human idea, and has been a well-established therapy for mental disorders for few decades. Caius Plinius during Antiquity and Scribonius Largus, who lived in the first century AD, proposed respectively contacts with the Electric ray (Torpedo Fish) for the treatment of post-partum pain and severe headaches. In the 19th century, electrical stimulation was commonly prescribed by neurologists for nervous disease (Edel and Caroli, 1987). Today, various electrical stimulation techniques exist to modulate neuronal systems and novel techniques for an optimal clinical treatment of a specific pathology gain more and more attention (Sun and Morrell, 2014; Chen et al., 2022). They could be used as an additional therapeutic lever or as an alternative to pharmacological medication, thus representing a hope for pharmaco-resistant forms of disease.

Brain oscillations result from coordinated electrical neuronal tissues activity within and between structures and networks. Implicated in various neural processes, such as perception, attention and cognition, their disruption yields pathological rhythms, which reflect abnormal activity of the implicated brain network, notably at the cellular and molecular level (Basar, 2013). These pathological rhythms serve as good biomarkers for neuropathologies. For instance, neurophysiological studies have revealed that a large number of mental disorders exhibit pathological rhythms, which do not occur in healthy patients (Schulman et al., 2011). Neurostimulation techniques have identified such pathological rhythms as good stimulation targets for the treatment of brain oscillatory disorders. Neurostimulation induces electric currents in neuronal tissue. Depending on the stimulation protocol, i.e. the temporal stimulation current shape, its duration and pause and the number of repetitions, neurostimulation can lead to neural plasticity effects or to pacemaker-like brain stimulation, respectively.

For example, Deep Brain Stimulation (DBS) is an invasive technique and proposed for patients suffering from severe pharmaco-resistant Parkinson's disease (PD) or obsessive-compulsive disorders. In PD patients aberrant hypersynchronicity and hyperactivity in the β-frequency band (12–30 Hz) of the basal ganglia-thalamocortical network can be addressed by the pharmacological medication (e.g. Levodopa) or DBS. The conventional DBS protocols focus on the subthalamic nucleus or globus pallidus stimulation continuously at a temporally constant frequency about 130 Hz. The suppression of the pathological beta oscillations was correlated with improving motor symptoms (Kühn et al., 2008). Recent techniques (Hosain et al., 2014; Fleming et al., 2020) propose to apply an adaptive closed-loop stimulation protocol based on observed intracranial brain activity. In addition to this intracranial neurostimulation technique, transcranial electrical stimulation (TES) and transcranial magnetic stimulation (TMS) are non-invasive neuromodulation approaches in which, respectively, a low electrical current and a magnetic field are applied to the cortical tissues. The TES current modalities include direct currents (tDCS), i.e. constant currents, alternating current (tACS), i.e. typically oscillatory currents, and random noise-shape currents (tRNS), which typically includes frequencies above the β-frequency band. It was shown that tDCS can improve cognitive performance in healthy subjects (Brunelin et al., 2012) and patients (Stagg et al., 2018) and it is applied as a therapeutic means to target brain network dysfunctions, such as Attention-Deficit/Hyperactivity Disorder (Nejati et al., 2020) and major depressive disorder (Bennabi and Haffen, 2018).

Although the neurostimulation techniques mentioned above may permit to alleviate mental disorder patients from symptoms, the success rate of these treatments is still limited (Nasr et al., 2022). This underperformance results from non-optimal choices of the stimulation protocol originating from the lack of understanding of the underlying neural response to stimulations and the non-patient specific stimulation protocol. In other words, typically the stimulation protocol (including size, duration, repetition cycle of the stimulation signal) is open-loop, i.e. pre-defined by the clinician based on heuristic criteria before the stimulation starts (Paulus, 2011). This non-optimal approach is inferior to so-called closed-loop techniques, which automatically adapt to the patients current brain/health state. Such an adaptive, or closed-loop, approach has been introduced for intracranial (Hartshorn and Jobst, 2018; Prosky et al., 2021; Stanslaski et al., 2022) and transcranial stimulation (Tervo et al., 2022) and has been shown to improve neurostimulation in major depression patients (Scangos et al., 2021), epilepsy (Haeusermann et al., 2023) and affective and anxiety disorders (Guerrero Moreno et al., 2021). These proposed closed-loop methods are adaptive in the sense that a pre-defined stimulation signal is applied when observed brain activity fulfills certain criteria, such as passing an amplitude or power threshold. While this adaptive approach improves existing open-loop methods, the pre-defined stimulation signal may still be non-optimally chosen. Recently proposed methods produce better results by using reference signal tracking control schemes such as proportional integral (PI) controller (Westover et al., 2015; Bolus et al., 2018; Su et al., 2019; Zhu et al., 2021), linear quadratic regulator (LQR) (Yang et al., 2018, 2019; Bolus et al., 2021) or model predictive control (MPC) (Fang and Yang, 2022, 2023) which uses an LQR in a MPC framework. However this form of control requires to pre-define a reference signal which is often non-trivial to provide in a patient specific manner. Furthermore, because of the stochastic nature of brain signals, forcing the resting state signal to follow a reference signal with its own independent noise creates an unnecessary constraint for the stimulation signal when we only want to regulate the power in given frequency bands.

We propose to estimate a stimulation signal on the basis of observed brain activity without the need to track a reference signal. The target stimulation signal is computed directly via a linear controller synthesized using a user-defined filter that encodes the desired frequency-domain modifications. We argue that it is the natural choice for a closed-loop optimization in the presence of pathological rhythms: typically the pathology is identified by an abnormal power in a certain frequency band and the closed-loop control aims to modify this power value in such a way that the final brain activity power spectral distribution resembles the distribution of a healthy subject. Examples are pathological too strong β-rhythm magnitudes in Parkinson's disease (Martin et al., 2018) and too weak α-rhythm (Howells et al., 2018) and too strong γ-rhythm in schizophrenia (Y et al., 2015). This approach implies the hypothesis that modifying the observed pathological brain rhythms of a patient to resemble brain rhythms of a healthy subject renders the patients brain state and improves the patients health situation. This assumption was motivated by the impressive improving impact of DBS in psychiatric disorders (Holtzheimer and Mayberg, 2011).

Technically, the proposed method aims to reshape the spectral distribution of observed data, such as electroencephalographic data (EEG). For illustration, we consider pathological brain rhythms observed in psychosis in the α- (Howells et al., 2018) and γ-band (Leicht et al., 2015). Our method relies on the extraction and the filtering in real-time of the brain resting state activity signal, using the EEG and an estimated brain response model. The underlying brain model is fully non-parametric and estimated from observed resting state EEG. Moreover, we consider the fact that the closed-loop feedback exhibits a certain conduction delay between measurement and stimulation. This conduction delay results from the transmission delay in the hardware and the numerical computation time of the stimulation signal. Very first estimates of this delay time are in the range of few tens of milliseconds (Private communication, Isope, 2020), i.e. in the range of EEG signal time scales. Consequently, the present feedback delay in real-world systems may affect the methods performance. Fang and Yang (2022) presents a method to increase the robustness of an adaptive closed-loop controller against delay by reducing the sensitivity of the closed-loop to high frequency disturbances. However, while this decreases generally the error in the closed-loop output, this also prevents to apply the control signal specifically to high frequency ranges like the γ-range. To our best knowledge, the present study is the first providing a method to effectively compensate the frequency-domain errors created by the feedback delay in closed-loop neurostimulation systems without sacrificing the controllability of high frequency ranges.

The remaining article is organized as follows: Section 2 presents the neurostimulation setup and the closed-loop circuit studied in the rest of this paper. Then, we propose a model-based controller design to apply desired modifications to the observed activity signal. Subsequently, we propose a model estimation method to extract the brain input response model needed for the controller design. Later, we address the problem of the closed-loop delay by designing an additional system to approximate the future values of the observations. Finally, we present two brain models, which illustrate and validate the proposed method. Then, Section 3 presents the simulation results of our circuits, including the accuracy of the model estimation step and the delay compensation. Lastly, in Section 4, we discuss the results of the method presented in the paper compared to the state of the art, mention limitations and pinpoint some perspectives and possible experimental tests.

2. Materials and methods

2.1. Simulated neurostimulation

We build a theoretical plant as a circuit containing a stimulation element and an observation element, both connected to the model brain system under study. In real practice, the stimulation element corresponds to the neurostimulation device, such as a TES system or a TMS coil. In contrast, the observation element may represent electro-/magneto-encephalographic electrodes (in the following called EEG) or electrodes observing Local Field Potential. We define the time-dependent functions u:ℝ → ℝ and y:ℝ → ℝ as the input stimulation current and the output EEG signal, respectively.

If no input current is applied, the output is a non-zero stochastic signal y₀ corresponding to the measured resting state EEG activity and a non-zero neurostimulation current alters the output signal as a linear response. This alteration is caused by a change in the brain activity in response to the neurostimulation input and a direct measurement of the input current. The latter is undesirable as it is not correlated with brain dynamics but only with neurostimulation and measurement devices. In the following, we assume that observations include brain dynamics correlated output only while direct current measurements are filtered out. A method to remove the direct current measurement from the EEG signal is discussed in Section 4.

Then, we define the plant $\mathcal{P}$ as the system that takes u as its input and generates an output y which is equal to y₀ when no input is applied. By modeling the dynamics of $\mathcal{P}$, our goal is a neurostimulation signal u that causes predetermined changes in the spectral power amplitude of the output signal y. In our case, the goal is to increase the activity in the alpha band (8 − 12 Hz) and decrease the activity in the gamma band (25 − 55 Hz) motivated by aberrant power spectrum magnitudes in schizophrenia (Howells et al., 2018; Martin et al., 2018). A possible experimental setup involving our method is sketched in Figure 1A.

FIGURE 1

Figure 1. Closed-loop neurostimulation setup. (A) The different steps of our method in a possible experimental setup. First, the clinician specifies the desired frequency domain modifications to apply e.g. increase α-activity and decrease γ-activity. Second, the brain input response model is fitted to the measured EEG signal in order to match as closely as possible the current brain dynamics of the patient. Finally, the clinician's specifications and the fitted model are used to synthesize the closed-loop controller $\mathcal{K}$ that is directly used for closed-loop neurostimulation. (B) The controller directly produces the stimulation signal u in function of the measured total brain output y in real-time.

2.2. Linear time invariant model

We assume that the observed output response to a small neurostimulation input u is linear and time-invariant (LTI). This assumption is supported by multiple results across literature (Popivanov et al., 1996; Liu et al., 2010; Kim and Ching, 2016). Thus, there is an underlying LTI system $\mathcal{G}$ that produces an output y_u for any given input u. For this system, we can define a function g:ℝ → ℝ, which is the output produced by the plant input response system $\mathcal{G}$ in response to a unit impulse signal δ(t) where t ∈ ℝ is the time elapsed since the start of the signal. This function g is also called the unit impulse response of $\mathcal{G}$ and we have

$\begin{array}{l}\begin{array}{l}{y}_u(t)=g(t)*u(t):={\displaystyle {\int }_{-\infty }^{+\infty }g(t^{\prime })u(t-t^{\prime })d{t}^{\prime }.}\hfill \end{array}\end{array}$

with time t and * denotes the convolution over time. It leads to the total plant output

$\begin{array}{ll}y(t)=y_0(t)+y_u(t)=y_0(t)+g(t)*u(t).& (1)\end{array}$

With this choice of model, the contribution of the neurostimulation response to the total output is purely additive, allowing us to focus the analysis on $\mathcal{G}$, which represents the neurostimulation response part of the plant system. We also see that y₀, the resting state activity, contains the stochastic part of the output, while y_u can be predicted for any known input signal u if we have a model for the system $\mathcal{G}$. A method to estimate the plant input response model $\mathcal{G}$ is presented in Section 2.4.

2.3. Closed-loop control

In this section, we suppose that the function g is known. The estimation of g will be the aim of Section 2.4.

To close the loop, we generate the plant input signal u as the output of a linear controller $\mathcal{K}$ in response to the plant output y

$\begin{array}{l}u(t)=k(t)*y(t),\end{array}$

where k:ℝ → ℝ is the unit impulse response of the controller $\mathcal{K}$. We can now rewrite Eq. (1) as

$\begin{array}{ll}y(t)=y_0(t)+g(t)*k(t)*y(t).& (2)\end{array}$

Here, we assume that no delay between observation and stimulation application is present. We will relax this condition in Section 2.5. To solve Eq. (2), we apply the Laplace transform defined for each time-dependent function x : ℝ → ℝ by

$\begin{array}{ll}X(s)=\mathcal{L}\{x(t)\}(s):={\displaystyle {\int }_{0^{-}}^{+\infty }x}(t)e^{-st}dt,& (3)\end{array}$

Thus, we define Y : ℂ → ℂ, Y₀ : ℂ → ℂ, G : ℂ → ℂ and K : ℂ → ℂ as the Laplace transforms of respectively y, y₀, g and k, allowing us to write Eq. (2) as

$\begin{array}{l}Y(s)=Y_0(s)+G(s)K(s)Y(s).\end{array}$

Hence

$\begin{array}{ll}Y(s)=\frac{1}{1-G(s)K(s)}{Y}_0(s).& (4)\end{array}$

We now have an equation for the closed-loop output in function of the resting state activity. A block diagram of the closed-loop circuit is shown in Figure 1B.

Hence to design the frequency distribution of y we tune the frequency distribution of the transfer function K of the controller $\mathcal{K}$.

2.3.1. Controller synthesis

Our closed-loop setup aims to tune the observation power spectrum, or equivalently, the choice of Y(s) subjected to the resting state Y₀(s). To this end, we define a linear filter $\mathcal{H}$ with transfer function H : ℂ → ℂ and

$\begin{array}{ll}Y(s)=Y_0(s)+H(s)Y_0(s).& (5)\end{array}$

Specifically, we intend to restore the physiological state of the brain, e.g., of a schizophrenic patient as our motivation, with an observed EEG presenting low alpha activity and high gamma activity. The chosen filter $\mathcal{H}$ is a weighted double bandpass filter with positive weight in the α-frequency band to increase α-power and negative weights in the γ-band to decrease the systems γ-activity. The filter's transfer function is defined as

$\begin{array}{l}H(s)=c_1\frac{2\pi B_{1}s}{s^2+2\pi B_{1}s+{(2\pi f_1)}^2}+c_2\frac{2\pi B_{2}s}{s^2+2\pi B_{2}s+{(2\pi f_2)}^2}.\end{array}$

The exact parameters of $\mathcal{H}$ are shown in Table 1.

TABLE 1

Table 1. Parameter set of the filter $\mathcal{H}$.

We can synthesize the closed-loop controller $\mathcal{K}$, by combining equations (4) and (5) and solving for K as

$\begin{array}{l}\frac{1}{1-G(s)K(s)}{Y}_0(s)=Y_0(s)+H(s)Y_0(s)\\ K(s)=\frac{H(s)}{(1+H(s))G(s)}.& (6)\end{array}$

Therefore, if we know the plant input response transfer function G, we can find that desired controller transfer function K by Eq. (6). Once the transfer function is obtained, we can use it to find a corresponding state-space representation (Hespanha, 2018) for time domain simulations. The state-space system's ordinary differential equations can then be implemented by any device that can measure the brain activity and produce a custom neurostimulation signal in real-time.

2.4. Model estimation

The design of our closed-loop controller requires estimating the plant input response system $\mathcal{G}$, which in practice includes the brain dynamics, the neurostimulation device and the observation device. Our approach includes the estimation of $\mathcal{G}$ directly from observed brain activity, such as EEG of the patient. This ensure that the estimated plant model will be as close as possible to the real brain dynamics in the corresponding experimental conditions. To this end, we first need to find a way to measure the plant input response without also measuring the plant resting state activity. This is not trivial since the observed signal is the sum of the resting state activity and the stimulation response.

2.4.1. Signal extraction

Let us consider an open-loop setup with an arbitrary input u applied to the plant, which generates the output described by Eq. (1). In this equation, we only know u and y, and want to estimate the impulse response g. The problem is that we cannot observe y₀ only during the stimulation. Hence, based on previous data recordings, we need to find a way to predict the dynamics of y₀ during the stimulation.

First, we provide the following standard definitions that are important in the subsequent discussion. For any time domain signal x : ℝ → ℝ, we denote the Fourier transform by

$\begin{array}{ll}\widehat{x}(f)=\mathcal{F}\{x(t)\}(f):={\displaystyle {\int }_{-\infty }^{\infty }x}(t)e^{-2\pi ift}dt.& (7)\end{array}$

We define α₀ : ℝ → ℝ and α_u : ℝ → ℝ such as α₀(t) = y₀(t) − ȳ₀ and α_u(t) = y_u(t) − ȳ_u where ȳ, ȳ₀ and ȳ_u are respectively the ensemble means of y, y₀ and y_u.

We assume that y₀ is a wide-sense-stationary (WSS) random process, i.e. its mean and variance do not depend on time. According to the Wiener-Khinchin theorem (Khintchine, 1934; Gardiner, 2004), the autocorrelation function of a wide-sense-stationary random process has a spectral decomposition given by the power spectrum of that process

$\begin{array}{l}{S}_{yy}(f)=|\widehat{\alpha }(f){|}^2,\end{array}$

where $\widehat{\alpha }:ℝ\to ℂ$ is the Fourier transform of α(t) = y(t) − ȳ ∈ ℝ and $S_{yy}:ℝ\to {ℝ}^{+}$ is the spectral density of y.

Then, we can write Eq. (1) as

$\begin{array}{l}ȳ+\alpha (t)={ȳ}_0+{\alpha }_0(t)+{ȳ}_u+{\alpha }_u(t),\end{array}$

where ȳ = ȳ₀ + ȳ_u. The equation then simplifies to

$\begin{array}{l}\alpha (t)={\alpha }_0(t)+{\alpha }_u(t).\end{array}$

By application of the Fourier transform, we obtain

$\begin{array}{l}\widehat{\alpha }(f)={\widehat{\alpha }}_0(f)+{\widehat{\alpha }}_u(f)\end{array}$

and

$\begin{array}{l}|\widehat{\alpha }(f){|}^2=|{\widehat{\alpha }}_0(f){|}^2+|{\widehat{\alpha }}_u(f){|}^2+2\text{Re}\left[{\widehat{\alpha }}_0(f){\widehat{\alpha }}_u{(f)}^{*}\right].\end{array}$

In the following, we compute the ensemble average of each term of this equation. Since α and α_u are two independent processes sampled at different times and $\langle {\widehat{\alpha }}_0\rangle =\langle {\widehat{\alpha }}_u\rangle =0$.

Hence

$\begin{array}{l}\text{〈}2\text{Re}({\widehat{\alpha }}_0(f){\widehat{\alpha }}_u{(f)}^{*})\text{〉}=2\text{Re}\left[\text{〈}{\widehat{\alpha }}_0(f){\widehat{\alpha }}_u{(f)}^{*}\text{〉}\right]=0.\end{array}$

Here and in the following, 〈·〉 denotes the ensemble average. We point out that although Eq. (8) does hold when considering the ensemble average of the signals, fluctuations around 0 still remain in Eq. (8) for finite ensemble number of finite time signals.

We define the amplitude ratio (AR) as

$\begin{array}{ll}\text{AR}=\langle \frac{|â(f){|}^2}{|{â}_0(f){|}^2}\rangle ,& (8)\end{array}$

which quantifies the gain of amplitude between the resting-state output and the stimulated output. The fluctuations mentioned above can also be reduced by increasing the input current which will lead to a higher AR and therefore, a higher contribution to the total signal of the â_u term (which is known) compared to the other terms.

Nevertheless, this yields

$\begin{array}{lll}\langle |{\widehat{\alpha }}_u(f){|}^2\rangle & =\langle |\widehat{\alpha }(f){|}^2\rangle -\langle |{\widehat{\alpha }}_0(f){|}^2\rangle .& (9)\end{array}$

Using Eq. (1), we can express ${\widehat{\alpha }}_u$ in terms of the input impulse response g and the input u

$\begin{array}{l}{\widehat{\alpha }}_u(f)=\mathcal{F}\{y_u(t)-{\overline{y}}_u\}(f)\\ =\mathcal{F}\{g(t)*[u(t)-\overline{u}]\}(f)\\ =\widehat{g}(f)\mathcal{F}\{u(t)-\overline{u}\}(f) .& (10)\end{array}$

This equation permits to estimate the transfer function ĝ, see Section 3.

To express the transfer function ĝ in Laplace space, we use the fact that a unit impulse response function is non-zero only for positive time values t. Hence, based on equations (3) and (7), for s = 2πif, we can write the Laplace transform G as

$\begin{array}{l}G(2\pi if)={\displaystyle {\int }_{0^{-}}^{+\infty }g}(t)e^{-2\pi ift}dt={\displaystyle {\int }_{-\infty }^{+\infty }g}(t)e^{-2\pi ift}dt=\widehat{g}(f),\end{array}$

where the unit impulse response function g directly relates the output y to the resting-state output y₀ and the stimulation signal u, cf. Eq. (1).

We now need a method to generate a LTI system with a transfer function that matches the magnitude data computed with the formula. This is achieved by the magnitude vector fitting algorithm.

2.4.2. Magnitude vector fitting

Our goal is now to find a transfer function G corresponding the magnitude data |ĝ(f)|². For this purpose, we use a variant of the vector fitting algorithm design to work even with only the magnitude data. This method is called magnitude vector fitting (De Tommasi et al., 2010).

It allows to fit a passive LTI system to data by fitting the model transfer function. The system is synthesized such that the mean square error between the magnitude data sample and the transfer function evaluated at the same frequency points is minimized. De Tommasi et al. (2010) show that the transfer function of the fitted model reproduces both the magnitude and the phase shift of the original transfer function, although the fitting has been performed using sampled magnitude data only.

By minimizing the mean square error, the algorithm ensures that the transfer function of the fitted model accurately matches the original model as represented by the reconstructed gain data. Furthermore, to assess the accuracy of the reconstruction, we also compare the fitted model to the transfer function of the linearized brain model used for the simulation. This allows to double-check the validity of the reconstructed magnitude and also to verify if the reconstructed phase fits the phase of the original model as closely as possible (cf. Figures 3C, D).

We define the root-mean-square error (RMSE) as

$\begin{array}{ll}\text{RMSE}=\sqrt{\langle |\frac{\tilde{G}(2\pi if)-G(2\pi if)}{G(2\pi if)}{|}^2\rangle },& (11)\end{array}$

where $\tilde{G}$ is the fitted model's transfer function, G is the original model's transfer function and f ∈ ℝ⁺ are the frequency points used for the fitting. This allows to quantify the accuracy of the fitting step.

2.5. Delay compensation

Realistic feedback loops exhibit conduction delays between the moment of observation and feedback stimulation. Reasons for such delays are finite conduction speeds in cables, electronic switches, interfaces and delays caused by the controller device to compute numerically adapted stimuli. In systems with large time scales, such as controlled mechanical devices on the centimeter or larger scale, such delays may be negligible. Conversely biological systems such as the brain evolve on a millisecond scale and conduction delays may play an important role. Preliminary estimation of input and output devices of desktop computers have revealed an approximate delay of ~10ms. By virtue of such delays, it is important to take them into account in the closed-loop between the moment of observation and stimulation.

The different sources of delay can be represented as plant input and output delays. Since the controller $\mathcal{K}$ is LTI, the input and output delays can be concatenated into one single plant input delay. Hence, in our setup, we model the delay as an input delay τ in the system $\mathcal{G}$, modifying y(t) = g(t)*u(t) in Eq. (1) to y(t) = g(t)*u(t − τ). The Smith predictor (Smith, 1959; Morari and Zafiriou, 1989) is a known method to compensate such delay times. However, in the present problem, this approach allows controlling a limited frequency band only (Figures 7A, B). Consequently, it was necessary to invent another method. Since the plant input u is generated by the controller $\mathcal{K}$, we modify the controller to compensate the delay. To this end, the new controller $\mathcal{K}$ is chosen to estimate the future value of u instead of the present value. The new proposed method to apply this controller modification is presented in the Results Section 3.2.

2.6. Comparison to the state of the art

Our method is tested against two main control schemes commonly used in adaptive closed-loop neurostimulation. These control schemes are closed-loop control with a PI controller (Westover et al., 2015; Bolus et al., 2018; Su et al., 2019; Zhu et al., 2021) and Linear Quadratic Gaussian (LQG) control (Yang et al., 2018, 2019; Bolus et al., 2021) which refers to the combination of a Linear Quadratic Estimator (LQE) (or Kalman filter) with a Linear Quadratic Regulator (LQR) (Åström, 2012). For both these methods, the tracked reference signal is generated from pre-recorded pathological resting state activities to which we apply a filter restoring the target α- and γ-activities. In order to prevent closed-loop destabilization and regulate high frequency disturbances, a Smith predictor (Smith, 1959) is used for the PI and LQG controller to compensate the delay, while we use our own delay compensation method for our controller.

2.7. Brain models

Our closed-loop control method works for any LTI brain model. Furthermore, we want to show that it also produces good results on non-linear brain models, for which the neurostimulation input response behaves closely to an LTI system, when the input is sufficiently small. To this end, we present two models used to test our method. The first one is a linear neural population model of cortical activity, and the second one is a non-linear cortico-thalamic neural population model with cortico-thalamic delay.

2.7.1. Linear brain model

We describe neural population activity with a noise-driven linear model (Hutt, 2013). The model is composed of two pairs of interacting excitatory and inhibitory populations. Here we have $V_{e,i}^{(1,2)}:ℝ\to ℝ$, representing the mean activity of the associated population, where $V_e^{(1,2)}$ and $V_i^{(1,2)}$ correspond respectively to excitatory and inhibitory populations. Each population is driven by noise ξ_{1, 2} : ℝ → ℝ and the external input u : ℝ → ℝ, according to the following differential equations:

$\begin{array}{ll}\{\begin{array}{l}{\tau }_{e,1}\frac{d{V}_e{}^{(1)}(t)}{dt}\text{ = }(-1+N_{11})V_e^{(1)}(t)-N_{11}{V}_i^{(1)}(t)+b_{1}u(t)+{\xi }_1(t),\\ {\tau }_{i,1}\frac{d{V}_i{}^{(1)}(t)}{dt}=N_{21}{V}_e^{(1)}(t)+(-1-N_{21})V_i^{(1)}(t)+b_{2}u(t),\\ {\tau }_{e,2}\frac{d{V}_e{}^{(2)}(t)}{dt}=(-1+N_{12})V_e^{(2)}(t)-N_{12}{V}_i^{(2)}(t)+b_{3}u(t)+{\xi }_2(t),\\ {\tau }_{i,2}\frac{d{V}_i{}^{(2)}(t)}{dt}=N_{22}{V}_e^{(2)}(t)+(-1-N_{22})V_i^{(2)}(t)+b_{4}u(t),\end{array}& (12)\end{array}$

where the noise ξ_{1, 2} is uncorrelated Gaussian distributed with zero mean and variance ${\kappa }_{1,2}^2=10^{-7}$, and the stimulation u is weighted by the coupling constants b_i > 0 of the corresponding population. In addition, τ_{(e, i), (1, 2)} are the synaptic time constants of the populations, and constants N_ij > 0 are interaction gains of the respective population. Table 2 provides the parameters employed in subsequent simulations.

TABLE 2

Table 2. Parameter set of model (12).

The observed output

$\begin{array}{l}y(t)=V_e^{(1)}(t)-V_i^{(1)}(t)+V_e^{(2)}(t)-V_i^{(2)}(t)\end{array}$

is a sum of the effective field potential $V_e^{(j)}-V_i^{(j)}$ of both populations j = 1, 2, (cf. Figure 2A).

FIGURE 2

Figure 2. Healthy and pathological resting state activity of the linear and cortico-thalamic brain models. The pathological state is characterized by a decreased α-activity and an increased γ-activity. (A) Comparison of the healthy and pathological state of the linear brain model. The first row shows the last 500 ms of the simulated time series. The second row shows the power spectral densities estimated from the time series and the last row shows the estimated α- and γ-activities computed from the spectral densities and averaged over 50 simulations. (B) Same as (A) but using the cortico-thalamic brain model. “***” corresponds to a p-value less than 0.0005 using Welch's t-test.

The simulation of the linear brain model in time domain is done using the library control of python. The numerical integration is computed thanks to matrix exponential (Van Loan, 1978), with a simulation sampling time of 1ms. The resting state activity of the linear brain model in shown in Figure 2A.

2.7.2. Cortico-thalamic brain model

A different model considers the cortico-thalamic feedback circuit (Riedinger and Hutt, 2022). It describes the cortex layers I-III and the cortico-thalamic loop between cortical layers IV-VI, the thalamic relay cell population and the reticular structure. The cortical layer I-III exhibits mean activity of excitatory cells v and inhibitory cells w. Similarly, layer IV-VIs exhibits the mean activity V_e and V_i and thalamic relay cell populations the mean activity V_th,e and V_th,i. Moreover, the reticular structure has the mean activity V_ret. The fibers between the cortex and thalamus and the cortex and reticular structure exhibit a finite conduction delay τ (Hashemi et al., 2015; Riedinger and Hutt, 2022). The 7-dimensional dynamical system of the brain state $\mathbf{\text{x}}=(v,w,V_e,V_i,V_{th,e},V_{th,i},V_{ret})\in {ℝ}^7$ obeys

$\begin{array}{ll}\{\begin{array}{l}\stackrel{·}{\text{x}}(t)=\text{F}(\text{x}(t),\text{x}(t-\tau ))+\xi (t)+\text{B}u(t),\\ y(t)=\text{Cx}(t),\end{array}& (13)\end{array}$

where the superscript t denotes transposition, F ∈ ℝ⁷ is a nonlinear vector function, B ∈ ℝ^{7 × 1} is the input coupling matrix and C ∈ ℝ^{1 × 7} is the observation matrix. We mention that $\mathrm{\text{B}}={(b_1,b_2,b_3,b_4,0,0,0)}^t,b_i>0$, i.e. only the cortical layers are stimulated with weights b_i. The observation y captures the activity of the cortical excitatory populations (Nunez and Srinivasan, 2006; Riedinger and Hutt, 2022) with C = (c₁, 0, c₃, 0, 0, 0, 0), c_i > 0. For more details, please see the Supplementary material.

The time domain simulations of the cortico-thalamic model is done by numerical integration using the fourth-order Runge-Kutta method implemented by the scipy library in python with a maximum simulation time step of 1 ms. The resting state activity of the cortico-thalamic brain model is shown in Figure 2B.

2.8. Measuring brain activity

The general activity in EEG measurement is measured by first estimating the power spectral density of the signal using frequency bins of 1 Hz and then summing all the frequency bins up to the Nyquist frequency. In practice, the sum will be mostly determined by the activity in low frequencies and more precisely, near the α- and γ-activity peaks. To measure the results of our method, we define the α- and γ-activities as the sum of 1Hz frequency bins only in their respective frequency bands, i.e., respectively 8–12 Hz and 25–55 Hz. Meanwhile, the total activity of the neurostimulation signal, that we call the mean amplitude of u is computed from 0 Hz to the Nyquist frequency.

3. Results

The present work addresses two major problems in closed-loop control: the correct model choice of the systems dynamics and the present conduction delay. The subsequent sections propose solutions for both problems and illustrate them in some detail by applying them to the linear brain activity model from Section 2.7.1. The final section demonstrates the closed feedback loop for the cortico-thalamic brain model from Section 2.7.2.

3.1. Model estimation

Equations (9) and 10 permit to express the magnitude of ĝ(f) in terms of the spectral densities of observable signals

$\begin{array}{l}|\widehat{g}(f){|}^2|\mathcal{F}\{u(t)-\overline{u}\}(f){|}^2=|\widehat{\alpha }(f){|}^2-|{\widehat{\alpha }}_0(f){|}^2\\ |\widehat{g}(f){|}^2{S}_{uu}(f)=S_{yy}(f)-S_{y_0{y}_0}(f)\\ |\widehat{g}(f){|}^2=\frac{S_{yy}(f)-S_{y_0{y}_0}(f)}{S_{uu}(f)}.& (14)\end{array}$

The spectral density functions S_y0y0 and S_yy may be estimated numerically from output data before and during a stimulation with a known chosen stimulation function u. The estimation may be performed by applying conventional methods, such as the Welch method (Welch, 1967). These estimations provide the magnitude of the transfer function |ĝ| by utilizing Eq. (14). In detail, at first, we considered the linear model (12) and injected a white noise current into the plant gaining the system's response signal together with the resting state activity, (cf. Figure 3A). The subsequent estimation of S_yy(f), S_y0y0(f) and S_uu(f) (Figure 3B) from the data permitted to compute the brain input response model ĝ(f) by Eq. (14). We observe a very good accordance of the original model response function and its estimation in magnitude (Figure 3C) and phase (Figure 3D).

FIGURE 3

Figure 3. The magnitude vector fitting algorithm successfully reconstructs the transfer function G from magnitude-only data. (A) Time series of the resting state activity (blue), the input signal (green) and the stimulation response (red). (B) Spectral densities of the simulated input signal (green), the resting state activity (blue) and the stimulation response (red). The input signal is a white noise with chosen standard deviation 0.005. (C) Reconstructed gain |ĝ| of the plant input response. The fitted model (dashed cyan) accurately matches the original model (black) with a root mean square error (RMSE) of 5.4%±2.2% (Confidence Interval (CI) 95%). The RMSE represents the error percentage between the fitted model's transfer function and the original model's transfer function averaged over 50 trials. The pink curve is the raw data used for fitting, computed from the spectral density data in (A) using Eq. (14). (D) Reconstructed phase of the plant input response ĝ.

3.1.1. Robustness

The remaining error in the estimated model compared to the original model depends on the amplitude ratio between the stimulated output y and the resting state output y₀, (cf. Figure 4). Low stimulation currents or high driving noise can also cause the magnitude vector fitting algorithm not to converge, leading to a non-minimal mean-square error between the fitted and the original models when evaluated at the frequency sample points used for the algorithm.

FIGURE 4

Figure 4. The magnitude vector fitting algorithm's performances depend on the amplitude ratio of the stimulation current and the driving noise. (A) Reconstructed magnitude and phase data for an amplitude ratio (AR) of 4.488 ± 0.065. The AR is computed as the ratio between the mean amplitude of y and the mean amplitude of y₀ averaged over 50 simulations (CI 95%). The fitted model has a corresponding RMSE of 2.4% ± 1.1% (CI 95%). (B) Same as (A), but with an AR of 2.406 ± 0.042, the presented data are the same as in Figures 3C, D. (C) Same as (A), but with an AR of 1.479 ± 0.026 resulting in a RMSE of 15.6% ± 10.2% for the fitted model. We see that the noise level in magnitude data is higher for smaller AR, which leads to higher RMSE (CI 95%) between the fitted model's transfer function and the original model's transfer function. The fitted model is coded in dashed cyan and deviates more from the original model for higher noise levels. All data have been computed from 30s long simulated time series.

This problem can be solved by increasing the amplitude of the input current u that we inject in the plant, which decreases the contribution of the resting state driving noise ξ to the output signal relative to the input current. Although the remaining dominant input current is also noisy, its value at any time or frequency is known, meaning that it is canceled out in the ratio $\frac{S_{yy}}{S_{uu}}$ in Eq. (14). This effectively leads to lower noise in the transfer function magnitude data extracted with Eq. (14). The limitation is then set by the maximum amplitude of the current we are allowed to inject into the brain in a given neurostimulation setup. Indeed, the amplitude of the current is limited both for safety reasons that are beyond the scope of this paper and because of the assumption of linearity on which our method is based and which requires small currents. On the other hand, we can also decrease the noise in the spectral density data by increasing the stimulation time, and hence increasing the amount of data which decreases the contribution of the noise in the power spectral density estimation. Therefore, the accuracy of the model fitting step can be optimized by finding a trade-off between the maximum amplitude of the stimulation current, and the maximum duration of the stimulation.

The accuracy of the fitting is generally easily assessed by computing the root mean square error between the data and the fitted model's transfer function. In practise, this could be used as an indicator to evaluate if sufficient stimulation amplitude and time were chosen and then possibly reiterate this step with different parameters. The root mean square error could also be used to directly quantify the error between the transfer function of the estimated model and the original model, which can be an important parameter regarding the stability and the robustness of the closed-loop while performing delay compensation as discussed in the next section.

3.2. Delay compensation

Delay compensation is achieved by adding another LTI system at the output of the controller $\mathcal{K}$ (cf. Figure 5), whose purpose is to reproduce the transfer function of a negative delay. We call this system the predictor ϕ.

FIGURE 5

Figure 5. Closed-loop neurostimulation circuit with predictor.

However, perfectly reproducing the transfer function of a negative delay would be impossible since the associated time-domain system would then be a perfect predictor, which is a non-causal, i.e. un-physical, system. Nonetheless, we can build a causal and stable system that behaves almost like a perfect predictor, however only in the frequency ranges of interest.

The numerical implementation of the controller necessitates discretization in time. Consequently, it is reasonable to choose the predictor design as a discrete-time system, meaning that for any input signal at x_t : ℝ → ℝ at an instant t ∈ ℝ, it approximately predicts the future signal x_t+Δt where Δt ∈ ℝ is the sampling time chosen when building the predictor. Since x is a discrete sequence, its transfer function is obtained using the Z-transform, defined as

$\begin{array}{l}X(z)=Z\{x_{n\text{Δ}t}\}(z):={\displaystyle \sum _{n=0}^{\infty }}{x}_{n\text{Δ}t}{z}^{-n},\end{array}$

with z ∈ ℂ and X : ℂ → ℂ. Then the transfer function Φ : ℂ → ℂ of a negative delay of one step Δt applied to x would simply be Φ(z) = z, the Z-transform of a one-step delay. However, this choice would be non-causal, which is not implementable numerically in time. Nevertheless, to obtain a stable and implementable system with a transfer function as close as possible to z, we chose the ansatz

$\begin{array}{ll}\text{Φ}(z_0)=\frac{b_0{z}_0+b_1}{z_0-a}=z_0,& (15)\end{array}$

for a fixed value z = z₀ and where a ∈ ℝ is the pole of the system and b₀ ∈ ℝ and b₁ ∈ ℝ are the polynomial coefficients of the numerator of Φ. This equation corresponds to the transfer function of a discrete LTI system with exactly one pole and one zero, which is the closest form of a proper rational function to the identity function of z in the sense that it has only one more pole. We add the additional constraints that |a| < 1, since this is the necessary and sufficient condition for the discrete predictor $\mathcal{\varphi }$ to be stable.

We choose to reformulate this problem by setting a as a free parameter. This way, we can select any a between −1 and 1, and the remaining parameters are found by solving the linear equation b₀z₀ + b₁ = z₀(z₀ − a), where z ∈ ℂ is a chosen complex frequency point at which we want this equation to hold. Since there are two unknowns, we can write a second equation in which we want the derivative of each side of the equation also to be equal, yielding b₀ = 2z₀ − a. By replacing b₀ in the first equation, we obtain

$\begin{array}{l}{z}_0(2z_0-a)+b_1=z_0(z_0-a)\\ b_1\text{ =}-z_0^2.\end{array}$

In the z-domain, the zero frequency corresponds to z₀ = 1. We choose to solve this equation for this point, hence we can replace a, b₀ and b₁ in Eq. (15) which yields

$\begin{array}{ll}\text{Φ}(z)=\frac{(2-a)z-1}{z-a}.& (16)\end{array}$

This transfer function can then be converted to an associated state-space representation and used for time domain simulations with a sampling time Δt. The output of this system will then be y_t ≈ u_t+Δt for any input signal u_t. Simulating delays greater than the system sampling time is simply achieved by concatenating multiple times this predictor system. Here the delay has to be a multiple of the sampling time. This predictor can then be appended to the output of the digital controller $\mathcal{K}$.

To avoid closed-loop instability, we must limit the amplitude of the feedback signal computed from the controller input signal. This amplitude is determined by the three systems $\mathcal{G}$, $\mathcal{H}$ and $\mathcal{K}$. Since $\mathcal{G}$ is defined by the system under study and $\mathcal{H}$ is the chosen filter defining the desired modifications in the frequency distribution of the observed signal, $\mathcal{\varphi }$ (or equivalently parameter a) is the only degree of freedom. Figure 6 shows the region of closed-loop stability as a function of the predictor pole a and the delay.

FIGURE 6

Figure 6. The predictor pole location affects the closed-loop stability. The magnitude of the pole with the highest magnitude in the closed-loop transfer function parameterizes the stability of the closed-loop. Indeed, if this value is less than 0 dB, then all the poles of the closed-loop transfer function have a magnitude less than 0 dB, meaning that the system is stable. The system is unstable otherwise. Here the full curve, the dashed curve and the dotted curve correspond to predictors for delays of 3 ms, 5 ms and 10 ms, respectively. The higher the delay is, the lower is the size of the region of closed-loop stability for a.

Because the predictor has a gain that is still slightly greater than one in the frequency ranges of interest, we reduce the weights of the filter $\mathcal{H}$ to compensate for the excess gain at the α and γ-peaks. To do this, we simply divide the weight of each band by the magnitude of the predictor system evaluated at the band's natural frequency. This reduces the errors in the closed-loop transfer function in the α and γ-ranges.

Figure 7C shows results combining the model estimation by vector fitting and the delay compensation. The proposed closed-loop control yields an increase in α-power and a decrease in γ-power according to the employed target filter $\mathcal{H}$. The application of PI and LQG control with a Smith predictor for delay compensation (Figures 7A, B) has poor performances for higher γ-frequency activity.

FIGURE 7

Figure 7. Our method successfully decreases α- and γ-activities in the presence of a 5ms delay, while maintaining a low stimulation current. (A) Simulation time series (left panel) and power spectral densities (right panel) of the PI control loop with Smith predictor. (B) Same as (A) for the LQG control loop with Smith predictor. (C) Same as (A) for our method, including delay compensation. (D) α-, γ-activities and mean amplitude of u for each method, averaged over 50 simulations. “***” corresponds to a p-value less than 0.05 with Welch's t-test, while “n.s.” correspond to a p-value higher than 0.05. Our method provides both the closest match to the target α- and γ-activities and the lowest stimulation current (u) amplitude. The parameters for all controllers have been chosen to match the target activities as closely as possible without destabilizing the closed-loop. The activities are computed by averaging the spectral densities in their corresponding ranges while the u-amplitude correspond to the average spectral density of u from 0Hz to the Nyquist frequency.

3.2.1. Accuracy

The error between the achieved closed-loop output activity levels and the target activity levels is highly affected by the delay (cf. Figure 8). This is caused by the phase shift between the input and the output signal, which changes the effect of the control signal on the output in a frequency dependent manner. Our frequency range of interest is limited to frequencies below 55 Hz, which is the higher limit we use for the γ-range. In this case, the effect of delays of the order of milliseconds will be more visible for higher frequencies and higher delays (cf. Figure 8). However, using our predictor design, we significantly mitigate the effects of the delay in the frequency ranges of interest (cf. Figure 8) red curve. Nonetheless, these effect are still present, creating a limit to the maximum delay our predictor is able to compensate, which in our case is situated around 10ms.

FIGURE 8

Figure 8. Delay affects the stability and robustness of the closed-loop transfer function. (A) Stability margin and α- and γ-activities errors in function of the delay with no model errors. The stability margin is computed as the shortest distance between the Nyquist plot and the point −1. The higher the stability margin, the more robust the closed-loop is to model uncertainties. The α- and γ-errors represent the error percentage between the measured activity and the target activity. (B–D) Same as row (A) but with model errors (quantified as the RMSE of the model fitting step) corresponding respectively to the model fits shown in Figures 4A–C. The vertical bars represent the standard deviation of their corresponding point.

3.2.2. Stability and robustness

As discussed earlier, delay compensation can destabilize the closed-loop system depending on the parameters of its components. However, if the correct predictor pole is chosen based on Figure 6, the closed-loop will remain stable. These values are computed under the assumption that there are no model estimation errors. If we take into account the inaccuracies in the fitted brain model compared to the original brain model, extra gain can add up in the feedback signal, introducing again the risk of destabilizing the closed-loop (cf. Figure 8). The solution is either to simply reduce the amplitude of the spectral density modification that we want to apply by reducing the amplitude of the transfer function of filter $\mathcal{H}$, or to reduce the amplitude of the predictor $\mathcal{\text{Φ}}$ reducing its accuracy and possibly increasing delay errors. There is then a trade-off between how much we want to change the gain of the closed-loop transfer function while also compensating delay errors and how much we want to avoid closed-loop destabilization caused by model uncertainties. In any case, the inaccuracies in the estimated brain model create errors in the closed-loop transfer function regardless of the delay, which makes them the main determinant of the performance limits of our method in any given setup. Nevertheless, our method produces smaller α- γ-activities errors than LQG control and produces the smallest error for γ-activities for delays of 1 to 12ms (cf. Figure 8). Comparisons lead to p-values less than 0.0005 using Welch's t-test, except for LQG control in row D where the variance of the data was to high to find any significant difference with this test.

3.3. Application to cortico-thalamic circuit model

To extend the analysis to a biologically more realistic model, we employed a nonlinear cortico-thalamic brain model (cf. Section 2.7.2). Fitting a linear transfer function to the brain model activity as described above, we found a good accordance of fitted and original model as can be seen in Figures 9A, B. Small deviations in the gain and the phase resulted from the internal delay in the brain model and its non-linearity. Indeed, the magnitude vector fitting algorithm does not reproduce this delay but instead synthesizes a linear system that has no delay but still approximates well the transfer function of the original model. Nonetheless, the non-linearity of this model can also decrease the accuracy of the fitting, as we are trying to represent a non-linear input response model by a linear one. However, this effect is only seen when the current is large enough for the non-linear part of the response to be significant.

FIGURE 9

Figure 9. Fitted model-based control using the cortico-thalamic brain model successfully reproduces the target transfer function in the frequency domains of interest. (A) Magnitude of the fitted brain model transfer function (dashed cyan) obtained from the power spectral density data (pink) compared to the magnitude of the original cortico-thalamic brain model transfer function (black). (B) Phase shift of the fitted transfer function (dashed cyan) compared to the phase shift of the original transfer function (black). (C) Spectral densities of the resting state activity signal (blue), the stimulated brain output (red) and the stimulation signal (green). (D) α- and γ-activities of the closed-loop output averaged over 50 trials for which the fitting step was repeated each time. “***” corresponds to a p-value less than 0.0005 with Welch's t-test while “*” corresponds to a p-value less than 0.05. The α- and γ-activities are respectively increased and decreased after application of the closed-loop.

In fact, the model-based control enhances α-activity and diminishes γ-activity in good accordance to the imposed filter $\mathcal{H}$ (Figure 9C). The closed-loop transfer function deviates from the target transfer function for large frequencies beyond the γ-frequency range. This results from the employed conduction delay.

To elucidate better the functions of the different elements of the proposed method, we applied a second closed-loop setup, where the neurostimulation input was applied to the first three layers of the cortex modeled by u and v and to the reticulum modeled by V_ret (Figure 10). In this setting, the response in the high-frequency ranges are mainly produced by the cortex, while the response in low-frequency ranges originates mainly from the reticulum and the thalamic relay structure, with a gap approximately between 10Hz and 20Hz. The weak response between 10 and 20 Hz observable (cf. Figure 10A) is compensated by the controller, which produces a high magnitude stimulation in the closed-loop for these frequencies (cf. Figure 10C). The second consequence is the inaccuracy of the closed-loop output in the low-frequency ranges, this is caused by the rather long cortico-thalamic internal delay. This delay yields a larger phase shift at low-frequencies and originates from the fact that we observe signals in the cortex, but stimulate in the reticulum.

FIGURE 10

Figure 10. Reticulum stimulation yields incorrect closed-loop gain in low-frequency ranges. (A) Magnitude of the fitted brain model transfer function (dashed cyan) compared to the magnitude of the original cortico-thalamic brain model transfer function (black). (B) Phase shift of the fitted transfer function (dashed cyan) compared to the magnitude of the original transfer function (black). (C) Spectral densities of the resting state activity signal (blue), the stimulated brain output (red) and the stimulation signal (green). (D) Closed-loop transfer function (dashed red), compared to the target transfer function 1 + H(s) (black).

4. Discussion

The goal of the proposed method was to design a delayed closed-loop control method to apply defined modifications to the spectral distribution of an observed signal, such as EEG or LFP. Under the assumption of linear brain stimulation response, the presented work explicitly describes all the steps needed to build a delayed closed-loop neurostimulation setup to restore the physiological brain state of a patient (Hebb et al., 2014). Since the controller is modeled as a linear time-invariant system, its implementation is lightweight, straightforward, and easily applicable in most embedded systems. Applications to a simple neural populations model (Figure 7) and to a biologically plausible cortico-thalamic feedback system (Figures 9, 10) demonstrate its elements and their impact on the control performance.

4.1. Main contributions

Our method allows to precisely specify the desired frequency-domain modifications we want to apply to the brain activity. The resulting closed-loop controller can then synthesize in real-time the required closed-loop neurostimulation signal necessary to reach the desired output, without the need to track a pre-defined reference signal. This make the method more flexible since it requires to specify relative rather that absolute signal modifications which is preferable considering the intra- and inter-patients variability of the EEG spectrum. Furthermore, specifying a reference signal which is a stochastic signal uncorrelated with the noise of the current plant output introduces additional noise in the feedback signal, since the controller needs to compensate for both the mean difference between the frequency distributions of the two signal and the difference between the driving noises of the two signals. Therefore, our method is able to track a target frequency distribution for the brain output with a lower current amplitude than classical methods (cf. Figure 7).

4.1.1. Model estimation

We assume resting state activity signal driven by noise, when no neurostimulation is applied. Injecting a stimulation creates an additional response that adds to the resting state. Consequently, both the resting state signal and response signal can be observed separately in experimental practice and they serve to estimate a linear state-space model as outlined in Section 3.1. This approach is successful for both simplified linear models (cf. Figures 3, 4) and neurophsysiological realistic nonlinear models (cf. Figure 9). This approximation is suitable for nonlinear systems whose dynamics evolve close to a stationary state. Several studies have already exposed evidence confirming that the measured brain dynamics behave mostly linearly at macroscopic scales (Popivanov et al., 1996; Liu et al., 2010). Moreover, in the case of the brain response to small neurostimulation input, our assumption of the linear brain response is supported by results of Kim and Ching (2016). The authors of this study measured the controllability Gramian of their brain model with nonlinear sigmoid transfer function, similar to the cortico-thalamic brain model (Riedinger and Hutt, 2022) used in this paper. If the system exhibits nonlinear dynamics far from any linear approximation, such as bistable dynamics and chaotic evolution, the proposed vector fitting technique may yield a too large model error and thus instability of the closed-loop feedback. The hypothesis of macroscopically linear dynamics has also recently been tested against various nonlinear models (Nozari et al., 2020). While that work included fitting methods for both linear and nonlinear brain models, our work chose the paradigm of purely frequency domain model fitting with the magnitude vector fitting algorithm (De Tommasi et al., 2010) and applied it to the brain input response system, which we could isolate thanks to a simple open-loop neurostimulation setup. While models have already been studied in application to neurostimulation (Modolo et al., 2011; Wagner et al., 2014), we propose a straightforward black box modeling approach that is directly usable for adaptive closed-loop neurostimulation, and is technically applicable easily for each individual patients before any closed-loop neurostimulation sessions.

4.1.2. Delay compensation

Conduction delays of a few milliseconds in the transmission between observation and stimulation may be negligible in systems evolving on time scales of seconds or longer, but may play an important role in neural systems. Our study demonstrates that such feedback delays may introduce control errors and we show how these errors can be avoided by a novel delay compensation method (Section 3.2). Application to the linear model (7) demonstrated its superior performance compared to conventional delay compensation methods. Delay compensating systems have already been described in other work (Guo et al., 2004; Hosseini et al., 2019). However, we used a design primarily focused on the correction of a gain error in the closed-loop transfer function, whereas the majority of the current research is based on time domain criterion and stability enforcement (Sönmez and Ayasun, 2015; Ledva et al., 2017). The methods performance, i.e. how well the total gain function fits to the pre-defined transfer function, is good for low-frequencies but weakens for frequencies exceeding a limit frequency. Note that frequency domain compensation has also already been achieved, notably via delay equalizers (Podilchak et al., 2009). However, this would restrict the frequency range in which the delay is compensated, and create additional errors in the surrounding frequencies. Other designs include filters with negative group delays, however their applications are limited to band limited input signals (Bukhman and Bukhman, 2004; Voss, 2017). The predictor design we presented also relies on negative group delay, enabling delay compensation in a large frequency band, while still being applicable to the brain EEG, which is inherently not band limited, because of the noise. Nonetheless, while our predictor design allows to significantly decreases the delay errors in the closed-loop transfer function, the delay still imposes a limit on the controllable frequency range. The larger the delay, the smaller is this limit frequency. Low performance may induce instability in the feedback loop (Mirkin and Palmor, 2005) and thus should be avoided. A corresponding stability criteria has been proposed, (cf. Figure 6). Better predictor designs could allow better performance of the closed-loop system for larger delays. The improvement of the accuracy of our closed-loop neurostimulation setup by building more efficient predictor designs is in progress and we refer the reader to future work.

4.2. Limits of our methodology

4.2.1. Experimental stimulation parameters and safety

Experimental stimulation protocols have to ensure the subjects safety (Ko, 2021) and thus avoide stimulus-induced health risks and complications. For instance, tDCS may be administered for a duration of 60 minutes and a maximum current of 4 mA without yielding health risks. However, parameters beyond these limits may yield adverse effects in subjects, such as skin lesions similar to burns and mania or hypomania in patients with depression (Matsumoto and Ugawa, 2017). The proposed method does not limit the stimulation duration per se, but of course the duration can be chosen accordingly without constrating the method. The method adapts the systems brain rhythms to the target rhythms very rapidly on a time scale of less than a second and hence permits rather short stimulation duration longer than a second.

Moreover, the proposed method does not specify absolute stimulation current magnitude applied. The impact of stimulation at certain magnitudes depends heavily on the stimulation type. In tDCS, anodal stimulation with positive currents have a different impact as cathodal stimulation with negative currents. In addition, currents are thought to have to pass a certain threshold to yield a measurable effect. In tACS (Moliadze et al., 2012), stimulating in the α-frequency range large and small magnitudes yield excitation and inhibition, respectively, while intermediate magnitudes yield weak effects. Stimulating with a range of frequencies, as in tRNS (Potok et al., 2022), a 1mA peak-to-peak amplitude for 10 minutes stimulation duration does not yield adverse effects. We conclude that it is not straight-forward to decide which stimulation magnitude applied in the presented method would be safe for human subjects, since the stimulation signal is neither constant, single frequency oscillation nor random noise. In sum, we argue that a maximum peak-to-peak amplitude of 1mA for few tens of minutes may not yield adverse effects, but still may evoke a measurable impact on observations and the brain state. Of course, future experimental studies will gain deeper insights.

4.2.2. Model internal delay

The internal delay in the brain is not reproducible by the magnitude vector fitting algorithm, which relies on the time invariance of the signals. Hence, this will cause errors in the transfer function of the fitted model (cf. Figure 9) that are larger for higher contribution of the delay in the output (cf. Figure 10). To limit this effect, we must minimize the delay between the application of the neurostimulation input and the measurement of the response to this input as much as possible by taking into account the delay between the different brain regions.

4.2.3. Estimating the closed-loop delay

For delay compensation, in this paper, we assumed that we know the conduction delay in the closed-loop. However, although it is a single constant parameter, we would need a method to measure it for a real closed-loop neurostimulation setup. A straightforward way to do this would be to inject any current into the plant and measure the time lag between the moment at which we inject the input current and the moment at which we measure the output signal. This estimated delay would then correspond to the total closed-loop delay except for the computation delay of the digital controller $\mathcal{K}$. This computation delay can be easily measured with the same software used for computation, as it corresponds to the delay needed to perform constant-size matrix multiplications. Moreover, several methods have already been developed to estimate the conduction delays in linear systems (Schier, 1997; Dudarenko et al., 2014).

4.2.4. Direct input current measurements

One of the main challenges to solve for closed-loop neurostimulation is the elimination of direct transmission artifacts from the measured EEG signal (Iturrate et al., 2018). Indeed, when measuring the plant output signal, a portion of the measured signal might be a direct measurement of the input current without any influence from the brain dynamics. In the ideal case, one intends to minimize the contribution of the stimulation input to the observed signal since it would mean that the measured EEG signal does not fully correspond to the brain activity. Hence, reading the EEG of the patient would be more difficult for the user of our closed-loop setup, and the contribution of the brain dynamics to the closed-loop would be smaller. A simple solution to this problem is discussed further below.

4.3. Perspectives

The control proposed allows to perform accurate frequency shaping of the systems' activity spectral distribution. However, this approach is limited to linear models of the brain stimulation response. This may be disadvantageous if the systems dynamics exhibit nonlinear behavior (see e.g., Hutt and beim Graben, 2017) as we want to represent the brain dynamics realistically. Furthermore, in real-case scenarios, we would also have to take into account the noise in the acquisition of the signal by the sensor and in the application of the input signal by the actuator.

4.3.1. Filtering out direct input current measurements

Filtering out the direct input current measurements is achievable with our setup removing the strictly proper system requirement while using the magnitude vector fitting algorithm to measure the brain input response. In other words, while fitting the brain input response system, we want the fitted model to be able to contain a direct transmission term corresponding to the direct current measurement. Hence, if the real plant input response contains a significant direct transmission term, it will be identified by the magnitude vector fitting algorithm when synthesizing the estimated plant input response. The second step is them simply to substract the feedtrough term multiplied by the input current to the plant output signal. Thus, the remaining part of the signal would only correspond to the brain dynamics.

4.3.2. Application to multiple inputs multiple outputs plants

For now, we only focused on plant with a signal input signal and a single output signal. However, in a real setup, the EEG measurement is typically composed of multiple channels corresponding to different electrodes. This can also be true for the neurostimulation device. For example, with electric current stimulation, we can inject multiple signals using multiple electrodes. This can be simply solved by feeding a single input to each input channel and summing each output to a single output channel. However, when we separate the different channels, we can have more control over each individuals output channels. When we have multiple inputs and output, the plant is then a Multiple-Inputs Multiple-Outputs (MIMO) system. Everything developed in this paper is generalizable to MIMO systems, with one caveat: when solving E.q., (6), a unique solution only exists if the system has as more outputs than it has inputs. The user can always ensure this, by using as many neurostimulation input channels than there are EEG output channels. In this generalized setup, we can also define the filter $\mathcal{H}$ to apply different modifications to each output channel.

4.3.3. Neurostimulation effects on larger time scales

Our method relies only on the short term dynamics of the brain, using signal feedback and delay compensation to produce an adaptive stimulation current and obtain the desired EEG frequency distribution. However, more traditional neurostimulation techniques rely on the long term dynamics of neural plasticity, which is not modeled in the brain models we use in this paper. Long term brain adaptation to neurostimulation could cause the EEG frequency distribution to diverge from the desired frequency distribution after several minutes of stimulation. This effect could be compensated either by reiterating the model identification step and performing neurostimulation again, or by adjusting the weight of the filter $\mathcal{H}$ according to the observed changes in real-time. Incorporating the effect of neural plasticity in the brain models would allow our method to produce predictable and durable modification to the EEG frequency distribution, even after we stop the stimulation.

Data availability statement

The source code used for the simulation results can be freely accessed here: https://github.com/Thomas-Wahl/neuroclodec.

Author contributions

TW, AH, and MD contributed to the development of the methods presented in this study. TW produced the source code used for the simulations. JR and AH wrote the introduction section. TW and AH wrote the other sections of the manuscript. All authors read and approved the submitted version.

Funding

This research was funded by Inria in the Action Exploratoire project A/D Drugs.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2023.1183670/full#supplementary-material

References

Åström, K. J. (2012). Introduction to Stochastic Control Theory. Chelmsford, MA: Courier Corporation.

Google Scholar

Basar, E. (2013). Brain oscillations in neuropsychiatric disease. Dialogues Clin. Neurosci. 15, 291–300. doi: 10.31887/DCNS.2013.15.3/ebasar

CrossRef Full Text | Google Scholar

Bennabi, D., and Haffen, E. (2018). Transcranial direct current stimulation (tdcs): a promising treatment for major depressive disorder? Brain Sci. 8, 81. doi: 10.3390/brainsci8050081

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolus, M. F., Willats, A. A., Rozell, C. J., and Stanley, G. B. (2021). State-space optimal feedback control of optogenetically driven neural activity. J. Neural Eng. 18, 036006. doi: 10.1088/1741-2552/abb89c

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolus, M. F., Willats, A. A., Whitmire, C. J., Rozell, C. J., and Stanley, G. B. (2018). Design strategies for dynamic closed-loop optogenetic neurocontrol in vivo. J. Neural Eng. 15, 026011. doi: 10.1088/1741-2552/aaa506

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunelin, J., Mondino, M., Gassab, L., Haesebaert, F., Gaha, L., Suaud-Chagny, M., et al. (2012). Examining transcranial direct-current stimulation (tDCS) as a treatment for hallucinations in schizophrenia. Am. J. Psychiatry. 169, 719–724. doi: 10.1176/appi.ajp.2012.11071091

PubMed Abstract | CrossRef Full Text | Google Scholar

Bukhman, N., and Bukhman, S. (2004). On the negative delay time of a narrow-band signal as it passes through the resonant filter of absorption. Radiophys. Quantum Electron. 47, 68–76. doi: 10.1023/B:RAQE.0000031672.70934.3a

CrossRef Full Text | Google Scholar

Chen, Z. S., Kulkarni, P. P., Galatzer-Levy, I. R., Bigio, B., Nasca, C., and Zhang, Y. (2022). Modern views of machine learning for precision psychiatry. Patterns 3,100602. doi: 10.1016/j.patter.2022.100602

PubMed Abstract | CrossRef Full Text | Google Scholar

De Tommasi, L., Gustavsen, B., and Dhaene, T. (2010). Robust transfer function identification via an enhanced magnitude vector fitting algorithm. IET Control Theory & Applicat. 4, 1169–1178. doi: 10.1049/iet-cta.2009.0025

CrossRef Full Text | Google Scholar

Dudarenko, N., Polinova, N., and Ushakov, A. (2014). Fundamental matrix of linear continuous system in the problem of estimating its transport delay. Nauchno-Tekhnicheskii Vestnik Informatsionnykh Tekhnologii, Mekhaniki i Optiki. 14, 5.

Google Scholar

Edel, Y., and Caroli, F. (1987). Histoire de l'électrochoc : des traitements électriques à la convulsivothérapie en psychiatrie. Bulletin d'histoire de l'électricité. 9, 87–114. doi: 10.3406/helec.1987.1012

CrossRef Full Text | Google Scholar

Fang, H., and Yang, Y. (2022). Designing and validating a robust adaptive neuromodulation algorithm for closed-loop control of brain states. J. Neural Eng. 19, 036018. doi: 10.1088/1741-2552/ac7005

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, H., and Yang, Y. (2023). Predictive neuromodulation of cingulo-frontal neural dynamics in major depressive disorder using a brain-computer interface system: A simulation study. Front. Comput. Neurosci. 17, 1119685. doi: 10.3389/fncom.2023.1119685

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleming, J., Dunn, E., and Lowery, M. (2020). Simulation of closed-loop deep brain stimulation control schemes for suppression of pathological beta oscillations in parkinson's disease. Front. Neurosci. 14, 166. doi: 10.3389/fnins.2020.00166

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardiner, C. (2004). Handbook of Stochastic Methods. Berlin: Springer. doi: 10.1007/978-3-662-05389-8

CrossRef Full Text | Google Scholar

Guerrero Moreno, J., Biazoli, J.r. C, Fontes Baptista, A., and Remoaldo Trambaiolli, L. (2021). Closed-loop neurostimulation for affective symptoms and disorders: an overview. Biol. Psychol., 161:108081. doi: 10.1016/j.biopsycho.2021.108081

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, L., Cardullo, F., Houck, J., Kelley, L., and Wolters, T. (2004). “New predictive filters for compensating the transport delay on a flight simulator,” in AIAA Modeling and Simulation Technologies Conference and Exhibit (American Institute of Aeronautics and Astronautics). p. 5441. doi: 10.2514/6.2004-5441

CrossRef Full Text | Google Scholar

Haeusermann, T., Lechner, C., Fong, K., Sideman, A., Jaworska, A., Chiong, W., et al. (2023). Closed-loop neuromodulation and self-perception in clinical treatment of refractory epilepsy. AJOB Neurosci. 14, 32–44. doi: 10.1080/21507740.2021.1958100

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartshorn, A., and Jobst, B. (2018). Responsive brain stimulation in epilepsy. Ther. Adv. Chronic Dis. 9, 135–142. doi: 10.1177/2040622318774173

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashemi, M., Hutt, A., and Sleigh, J. (2015). How the cortico-thalamic feedback affects the EEG power spectrum over frontal and occipital regions during propofol-induced anaesthetic sedation. J. Comput. Neurosci. 39, 155. doi: 10.1007/s10827-015-0569-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hebb, A. O., Zhang, J. J., Mahoor, M. H., Tsiokos, C., Matlack, C., Chizeck, H. J., et al. (2014). Creating the feedback loop: closed-loop neurostimulation. Neurosurg. Clini. 25, 187–204. doi: 10.1016/j.nec.2013.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Hespanha, J. P. (2018). Linear Systems Theory. Princeton: Princeton University Press. doi: 10.23943/9781400890088

CrossRef Full Text | Google Scholar

Hirano, Y., Orbie, N., Kanba, S., Onitsuka, T., Nestor, P., and Spencer, K. (2015). Spontaneous gamma activity in schizophrenia. JAMA Psychiat. 72, 813–821. doi: 10.1001/jamapsychiatry.2014.2642

PubMed Abstract | CrossRef Full Text | Google Scholar

Holtzheimer, P., and Mayberg, H. (2011). Deep brain stimulation for psychiatric disorders. Annu. Rev. Neurosci. 34, 2890307. doi: 10.1146/annurev-neuro-061010-113638

PubMed Abstract | CrossRef Full Text | Google Scholar

Hosain, M., Kouzani, A., and Tye, S. (2014). Closed loop deep brain stimulation: an evolving technology. Australas. Phys. Eng. Sci. Med. 37, 619–634. doi: 10.1007/s13246-014-0297-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hosseini, S. A., Toulabi, M., Dobakhshari, A. S., Ashouri-Zadeh, A., and Ranjbar, A. M. (2019). Delay compensation of demand response and adaptive disturbance rejection applied to power system frequency control. IEEE Transac. Power Syst. 35, 2037–2046. doi: 10.1109/TPWRS.2019.2957125

CrossRef Full Text | Google Scholar

Howells, F., Temmingh, H., Hsieh, J., van Dijen, A., Baldwin, D., and Stein, D. (2018). Electroencephalographic delta/alpha frequency activity differentiates psychotic disorders: a study of schizophrenia, bipolar disorder and methamphetamine-induced psychotic disorder. Transl. Psychiatry. 8:75. doi: 10.1038/s41398-018-0105-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutt, A. (2013). The anesthetic propofol shifts the frequency of maximum spectral power in eeg during general anesthesia: analytical insights from a linear model. Front. Comput. Neurosci. 7, 2. doi: 10.3389/fncom.2013.00002

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutt, A., and beim Graben, P. (2017). Sequences by metastable attractors: interweaving dynamical systems and experimental data. Front. Appl. DYn. Syst. Stat. 3, 11. doi: 10.3389/fams.2017.00011

CrossRef Full Text | Google Scholar

Iturrate, I., Pereira, M., and Millán, J. R. (2018). Closed-loop electrical neurostimulation: challenges and opportunities. Curr. Opini. Biomed. Eng. 8, 28–37. doi: 10.1016/j.cobme.2018.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Khintchine, A. (1934). Korrelationstheorie der stationären stochastischen prozesse. Mathematische Annalen. 109, 604–615. doi: 10.1007/BF01449156

CrossRef Full Text | Google Scholar

Kim, S. A., and Ching, S. (2016). “Quasilinearization-based controllability analysis of neuronal rate networks,” in 2016 American Control Conference (ACC). Boston, MA: IEEE. p. 7371–7376.

Google Scholar

Ko, M. (2021). Safety of transcranial direct current stimulation in neurorehabilitation. Brain Neurorehabil. 14, e9. doi: 10.12786/bn.2021.14.e9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kühn, A. A., Kempf, F., Brücke, C., Gaynor Doyle, L., Martinez-Torres, I., Pogosyan, A., et al. (2008). High-frequency stimulation of the subthalamic nucleus suppresses oscillatory β activity in patients with parkinson's disease in parallel with improvement in motor performance. J. Neurosci. 28, 6165–6173. doi: 10.1523/JNEUROSCI.0282-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ledva, G. S., Vrettos, E., Mastellone, S., Andersson, G., and Mathieu, J. L. (2017). Managing communication delays and model error in demand response for frequency regulation. IEEE Transact. Power Syst. 33, 1299–1308. doi: 10.1109/TPWRS.2017.2725834

CrossRef Full Text | Google Scholar

Leicht, G., Vauth, S., Polomac, N., Andreou, C., Rauh, J., Mußmann, M., et al. (2015). EEG-Informed fMRI reveals a disturbed gamma-band–specific network in subjects at high risk for psychosis. Schizophr. Bull. 42, 239–249. doi: 10.1093/schbul/sbv092

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Rios, C., Zhang, N., Yang, L., Chen, W., and He, B. (2010). Linear and nonlinear relationships between visual stimuli, eeg and bold fmri signals. Neuroimage. 50, 1054–1066. doi: 10.1016/j.neuroimage.2010.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, S., Iturrate, I., Chavarriaga, R., Leeb, R., Sobolewski, A., Li, A., et al. (2018). Differential contributions of subthalamic beta rhythms and 1/f broadband activity to motor symptoms in parkinson's disease. NPJ Parkinson's Dis. 4, 32. doi: 10.1038/s41531-018-0068-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsumoto, H., and Ugawa, Y. (2017). Adverse events of tdcs and tacs: a review. Clini. Neurophysiol. Pract. 2, 19–25. doi: 10.1016/j.cnp.2016.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirkin, L., and Palmor, Z. J. (2005). Control Issues in Systems with Loop Delays, chapter 59. Boston, MA: Birkhäuser Boston. p. 627–648.

Google Scholar

Modolo, J., Legros, A., Thomas, A. W., and Beuter, A. (2011). Model-driven therapeutic treatment of neurological disorders: reshaping brain rhythms with neuromodulation. Interface Focus. 1, 61–74. doi: 10.1098/rsfs.2010.0509

PubMed Abstract | CrossRef Full Text | Google Scholar

Moliadze, V., Atalay, D., Antal, A., and Paulus, W. (2012). Close to threshold transcranial electrical stimulation preferentially activates inhibitory networks before switching to excitation with higher intensities. Brain Stimul. 5, 505–511. doi: 10.1016/j.brs.2011.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Morari, M., and Zafiriou, E. (1989). Robust Process Control. Madison, WI: Morari.

Google Scholar

Nasr, K., Haslacher, D., Dayan, E., Censor, N., Cohen, L. G., and Soekadar, S. R. (2022). Breaking the boundaries of interacting with the human brain using adaptive closed-loop stimulation. Prog. Neurobiol. 216, 102311. doi: 10.1016/j.pneurobio.2022.102311

PubMed Abstract | CrossRef Full Text | Google Scholar

Nejati, V., Salehinejad, M., Nitsche, M., Najian, A., and Javadi, A. (2020). Transcranial direct current stimulation improves executive dysfunctions in adhd: Implications for inhibitory control, interference control, working memory, and cognitive flexibility. J. Atten. Disord. 24, 1928–1943. doi: 10.1177/1087054717730611

PubMed Abstract | CrossRef Full Text | Google Scholar

Nozari, E., Bertolero, M. A., Stiso, J., Caciagli, L., Cornblath, E. J., He, X., et al. (2020). Is the brain macroscopically linear? a system identification of resting state dynamics. arXiv. [preprint]. doi: 10.1101/2020.12.21.423856

CrossRef Full Text | Google Scholar

Nunez, P., and Srinivasan, R. (2006). Electric Fields of the Brain: The Neurophysics of EEG. New York – Oxford: Oxford University Press.

Google Scholar

Paulus, W. (2011). Transcranial electrical stimulation (tes – tdcs; trns, tacs) methods. Neuropsychol.Rehabilitat. 21, 602–617. doi: 10.1080/09602011.2011.557292

PubMed Abstract | CrossRef Full Text | Google Scholar

Podilchak, S. K., Frank, B. M., Freundorfer, A. P., and Antar, Y. M. (2009). “High speed metamaterial-inspired negative group delay circuits in cmos for delay equalization,” in 2009 2nd Microsystems and Nanoelectronics Research Conference. Ottawa, ON: IEEE. p. 9–12.

Google Scholar

Popivanov, D., Dushanova, J., Mineva, A., and Krekule, I. (1996). “Detection of successive changes in dynamics of eeg time series: linear and nonlinear approach,” in Proceedings of 18th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. Amsterdam: IEEE. p. 1590–1591.

Google Scholar

Potok, W., van der Groen, O., Bächinger, M., Edwards, D., and Wenderoth, N. (2022). Transcranial random noise stimulation modulates neural processing of sensory and motor circuits, from potential cellular mechanisms to behavior: a scoping review. eNeuro. 9, ENEURO.0248-21.2021. doi: 10.1523/ENEURO.0248-21.2021

PubMed Abstract | CrossRef Full Text | Google Scholar

Prosky, J., Cagle, J., Sellers, K., Gilron, R., de Hemptinne, C., Schmitgen, A., et al. (2021). Practical closed-loop strategies for deep brain stimulation: lessons from chronic pain. Front. Neurosci. 15, 762097. doi: 10.3389/fnins.2021.762097

PubMed Abstract | CrossRef Full Text | Google Scholar

Riedinger, J., and Hutt, A. (2022). Mathematical model insights into eeg origin under transcranial direct current stimulation (tdcs) in the context of psychosis. J. Clin. Med. 11, 1845. doi: 10.3390/jcm11071845

PubMed Abstract | CrossRef Full Text | Google Scholar

Scangos, K., Khambhati, A., Daly, P., Makhoul, G., Sugrue, L., Zamanian, H., et al. (2021). Closed-loop neuromodulation in an individual with treatment-resistant depression. Nat. Med. 27, 1696–1700. doi: 10.1038/s41591-021-01480-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Schier, J. (1997). Estimation of transport delay using parallel recursive modified gram-schmidt algorithm. Int. J. Adapt. Control Signal Proc. 11, 431–442. doi: 10.1002/(SICI)1099-1115(199708)11:5<431::AID-ACS417>3.0.CO;2-Q

CrossRef Full Text | Google Scholar

Schulman, J. J., Cancro, R., Lowe, S., Lu, F., Walton, K. D., and Llinás, R. R. (2011). Imaging of thalamocortical dysrhythmia in neuropsychiatry. Front. Hum. Neurosci. 5, 69. doi: 10.3389/fnhum.2011.00069

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, O. J. (1959). A controller to overcome dead time. ISA J. 6, 28–33.

Google Scholar

Sönmez, S., and Ayasun, S. (2015). Stability region in the parameter space of pi controller for a single-area load frequency control system with time delay. IEEE Transact. Power Syst. 31, 829–830. doi: 10.1109/TPWRS.2015.2412678

CrossRef Full Text | Google Scholar

Stagg, C., Antal, A., and Nitsche, M. (2018). Physiology of transcranial direct current stimulation. J. ECT. 34, 144–152. doi: 10.1097/YCT.0000000000000510

PubMed Abstract | CrossRef Full Text | Google Scholar

Stanslaski, S., Farooqi, H., Sanabria, D., and Netoff, T. (2022). Fully closed loop test environment for adaptive implantable neural stimulators using computational models. J. Med. Device. 16, 034501. doi: 10.1115/1.4054083

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, F., Kumaravelu, K., Wang, J., and Grill, W. M. (2019). Model-based evaluation of closed-loop deep brain stimulation controller to adapt to dynamic changes in reference signal. Front. Neurosci. 13, 956. doi: 10.3389/fnins.2019.00956

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, F., and Morrell, M. (2014). Closed-loop neurostimulation: the clinical experience. Neurotherapeutics 11, 553–563. doi: 10.1007/s13311-014-0280-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Tervo, A. E., Nieminen, J. O., Lioumis, P., Metsomaa, J., Souza, V. H., Sinisalo, H., et al. (2022). Closed-loop optimization of transcranial magnetic stimulation with electroencephalography feedback. Brain Stimul. 15, 523–531. doi: 10.1016/j.brs.2022.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Loan, C. (1978). Computing integrals involving the matrix exponential. IEEE Trans. Automat. Contr. 23, 395–404. doi: 10.1109/TAC.1978.1101743

CrossRef Full Text | Google Scholar

Voss, H. U. (2017). A universal negative group delay filter for the prediction of band-limited signals. arXiv. [preprint].

PubMed Abstract | Google Scholar

Wagner, T., Eden, U., Rushmore, J., Russo, C. J., Dipietro, L., Fregni, F., et al. (2014). Impact of brain tissue filtering on neurostimulation fields: a modeling study. Neuroimage. 85, 1048–1057. doi: 10.1016/j.neuroimage.2013.06.079

PubMed Abstract | CrossRef Full Text | Google Scholar

Welch, P. (1967). The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodogram. Trans. Audio Electroacoustics. AU-15, 70–73. doi: 10.1109/TAU.1967.1161901

CrossRef Full Text | Google Scholar

Westover, M. B., Kim, S.-E., Ching, S., Purdon, P. L., and Brown, E. N. (2015). Robust control of burst suppression for medical coma. J. Neural Eng. 12, 046004. doi: 10.1088/1741-2560/12/4/046004

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Connolly, A. T., and Shanechi, M. M. (2018). A control-theoretic system identification framework and a real-time closed-loop clinical simulation testbed for electrical brain stimulation. J. Neural Eng. 15, 066007. doi: 10.1088/1741-2552/aad1a8

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Lee, J. T., Guidera, J. A., Vlasov, K. Y., Pei, J., Brown, E. N., et al. (2019). Developing a personalized closed-loop controller of medically-induced coma in a rodent model. J. Neural Eng. 16, 036022. doi: 10.1088/1741-2552/ab0ea4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Y., Wang, J., Li, H., Liu, C., and Grill, W. M. (2021). Adaptive parameter modulation of deep brain stimulation based on improved supervisory algorithm. Front. Neurosci. 15, 750806. doi: 10.3389/fnins.2021.750806

PubMed Abstract | CrossRef Full Text | Google Scholar