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.


. 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 obsessivecompulsive 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 nonpatient 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 closedloop, 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(Yang et al., , 2019Bolus et al., 2021) or model predictive control (MPC) 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  (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.
. Materials and methods

. . 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 : R → R and y : R → R 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 0 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 P as the system that takes u as its input and generates an output y which is equal to y 0 when no input is applied. By modeling the dynamics of 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.

. . 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 G that produces an output y u for any given input u. For this system, we can define a function g : R → R, which is the output produced by the plant input response system G in response to a unit impulse signal δ(t) where t ∈ R is the time elapsed since the start of the signal. This function g is also called the unit impulse response of G and we have with time t and * denotes the convolution over time. It leads to the total plant output y(t) = y 0 (t) + y u (t) = y 0 (t) + g(t) * u(t). (1) 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 G, which represents the neurostimulation response part of the plant system. We also see that y 0 , 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 G. A method to estimate the plant input response model G is presented in Section 2.4.

. . 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 K in response to the plant output y where k : R → R is the unit impulse response of the controller K. We can now rewrite Eq. (1) as y(t) = y 0 (t) + g(t) * k(t) * y(t). (2) . /fnins. .

Closed-loop neurostimulation setup. (A)
The di erent 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 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.
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 : R → R by Thus, we define Y : C → C, Y 0 : C → C, G : C → C and K : C → C as the Laplace transforms of respectively y, y 0 , g and k, allowing us to write Eq. (2) as Hence 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 K.

. . . 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 0 (s). To this end, we define a linear filter H with transfer function H : C → C and 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 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 The exact parameters of H are shown in Table 1. We can synthesize the closed-loop controller K, by combining equations (4) and (5) and solving for K as 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.
. /fnins. . Parameter Description Value f 1 α-band natural frequency 10ms The frequency parameters are chosen based on the alpha frequency range (8-12Hz) and the gamma frequency range (25-55 Hz) in an EEG. The weighting parameters c1 and c2, respectively positive and negative, corresponding to the choice to increase the alpha activity and decrease the gamma activity.

. . Model estimation
The design of our closed-loop controller requires estimating the plant input response system G, which in practice includes the brain dynamics, the neurostimulation device and the observation device. Our approach includes the estimation of 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.

. . . 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 0 only during the stimulation. Hence, based on previous data recordings, we need to find a way to predict the dynamics of y 0 during the stimulation.
First, we provide the following standard definitions that are important in the subsequent discussion. For any time domain signal x : R → R, we denote the Fourier transform bŷ We define α 0 : R → R and α u : R → R such as α 0 (t) = y 0 (t) −ȳ 0 and α u (t) = y u (t) −ȳ u whereȳ,ȳ 0 andȳ u are respectively the ensemble means of y, y 0 and y u .
We assume that y 0 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-sensestationary random process has a spectral decomposition given by the power spectrum of that process S yy (f ) = |α(f )| 2 , whereα : R → C is the Fourier transform of α(t) = y(t) −ȳ ∈ R and S yy : R → R + is the spectral density of y.
Then, we can write Eq. (1) as whereȳ =ȳ 0 +ȳ u . The equation then simplifies to By application of the Fourier transform, we obtain and 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 α 0 = α u = 0. Hence 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 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 Using Eq.
(1), we can expressα u in terms of the input impulse response g and the input û 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 where the unit impulse response function g directly relates the output y to the resting-state output y 0 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.

. . . Magnitude vector fitting
Our goal is now to find a transfer function G corresponding the magnitude data |ĝ(f )| 2 . 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 whereG is the fitted model's transfer function, G is the original model's transfer function and f ∈ R + are the frequency points used for the fitting. This allows to quantify the accuracy of the fitting step.

. . 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 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 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 K, we modify the controller to compensate the delay. To this end, the new controller 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.

. . 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(Yang et al., , 2019Bolus 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.

. . 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.

. . . Linear brain model
We describe neural population activity with a noise-driven linear model (Hutt, 2013 The choice of parameter is partially based on the paper in which it was developed (see Hutt, 2013). The difference between the healthy and pathological state is modeled by changes in the amplitude of the finite size fluctuations of each population.
noise ξ 1,2 : R → R and the external input u : R → R, according to the following differential equations: where the noise ξ 1,2 is uncorrelated Gaussian distributed with zero mean and variance κ 2 1,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.
The observed output is a sum of the effective field potential V i of both populations j = 1, 2, (cf. Figure 2A).
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.
. . . 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 where the superscript t denotes transposition, F ∈ R 7 is a nonlinear vector function, B ∈ R 7×1 is the input coupling matrix and C ∈ R 1×7 is the observation matrix. We mention that B 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 1 , 0, c 3 , 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.

. . 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.

. Results
The present work addresses two major problems in closedloop 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.

FIGURE
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 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 simulations. (B) Same as (A) but using the cortico-thalamic brain model. "***" corresponds to a p-value less than .

. . Model estimation
Equations (9) and 10 permit to express the magnitude ofĝ(f ) in terms of the spectral densities of observable signals The spectral density functions S y 0 y 0 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 y 0 y 0 (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).

. . . 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 0 , (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.
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 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 Frontiers in Neuroscience frontiersin.org . /fnins. .  ± .
. The AR is computed as the ratio between the mean amplitude of y and the mean amplitude of y averaged over simulations (CI %). The fitted model has a corresponding RMSE of . % ± . % (CI %). (B) Same as (A), but with an AR of .
± . , the presented data are the same as in Figures C, D. (C) Same as (A), but with an AR of . ± . resulting in a RMSE of . % ± . % for the fitted model. We see that the noise level in magnitude data is higher for smaller AR, which leads to higher RMSE (CI %) 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 s long simulated time series.
Frontiers in Neuroscience frontiersin.org . /fnins. . 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 closedloop while performing delay compensation as discussed in the next section.

. . Delay compensation
Delay compensation is achieved by adding another LTI system at the output of the controller K (cf. Figure 5), whose purpose is to reproduce the transfer function of a negative delay. We call this system the predictor φ.
However, perfectly reproducing the transfer function of a negative delay would be impossible since the associated timedomain system would then be a perfect predictor, which is a noncausal, 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 : R → R at an instant t ∈ R, it approximately predicts the future signal x t+ t where t ∈ R 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 with z ∈ C and X : C → C. Then the transfer function : C → C 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 for a fixed value z = z 0 and where a ∈ R is the pole of the system and b 0 ∈ R and b 1 ∈ R 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 φ 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 0 z 0 + b 1 = z 0 (z 0 − a), where z ∈ C 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 0 = 2z 0 − a. By replacing b 0 in the first equation, we obtain In the z-domain, the zero frequency corresponds to z 0 = 1. We choose to solve this equation for this point, hence we can replace a, b 0 and b 1 in Eq. (15) which yields This transfer function can then be converted to an associated statespace 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 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 G, H and K. Since G is defined by the system under study and H is the chosen filter defining the desired modifications in the frequency distribution of the observed signal, φ (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.
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 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 closedloop control yields an increase in α-power and a decrease in γ-power according to the employed target filter 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
The predictor pole location a ects 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 dB, then all the poles of the closed-loop transfer function have a magnitude less than 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 ms, ms and ms, respectively. The higher the delay is, the lower is the size of the region of closed-loop stability for a. simulations. "***" corresponds to a p-value less than . with Welch's t-test, while "n.s." correspond to a p-value higher than . . 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 Hz to the Nyquist frequency.

. . . 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 .
/fnins. . 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.

. . . 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 H, or to reduce the amplitude of the predictor reducing its accuracy and possibly increasing delay errors. There is then a tradeoff between how much we want to change the gain of the closedloop 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.
. . 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.
In fact, the model-based control enhances α-activity and diminishes γ-activity in good accordance to the imposed filter 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.

. 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.

. . Main contributions
Our method allows to precisely specify the desired frequencydomain 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).

. . . 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 closedloop 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 .
/fnins. . (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 trials for which the fitting step was repeated each time. "***" corresponds to a p-value less than .
with Welch's t-test while "*" corresponds to a p-value less than . . The αand γ-activities are respectively increased and decreased after application of the closed-loop. 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.

. . . 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.
. . Limits of our methodology . . . 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.

. . . 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.

. . . 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 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).

. . . Direct input current measurements
One of the main challenges to solve for closed-loop neurostimulation is the elimination of direct transmission Frontiers in Neuroscience frontiersin.org . /fnins. . 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.

. . 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.

. . . 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.

. . . 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 H to apply different modifications to each output channel.

. . . Neurostimulation e ects 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 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.
. /fnins. . 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