Amplitude Adaptive Modulation of Neural Oscillations Over Long-Term Dynamic Conditions: A Computational Study

Closed-loop deep brain stimulation (DBS) shows great potential for precise neuromodulation of various neurological disorders, particularly Parkinson’s disease (PD). However, substantial challenges remain in clinical translation due to the complex programming procedure of closed-loop DBS parameters. In this study, we proposed an online optimized amplitude adaptive strategy based on the particle swarm optimization (PSO) and proportional–integral–differential (PID) controller for modulation of the beta oscillation in a PD mean field model over long-term dynamic conditions. The strategy aimed to calculate the stimulation amplitude adapting to the fluctuations caused by circadian rhythm, medication rhythm, and stochasticity in the basal ganglia–thalamus–cortical circuit. The PID gains were optimized online using PSO, based on modulation accuracy, mean stimulation amplitude, and stimulation variation. The results showed that the proposed strategy optimized the stimulation amplitude and achieved beta power modulation under the influence of circadian rhythm, medication rhythm, and stochasticity of beta oscillations. This work offers a novel approach for precise neuromodulation with the potential for clinical translation.


I. INTRODUCTION
C LOSED-LOOP deep brain stimulation (DBS) has been used in clinical research to modulate neural activity and improve the symptoms of neurological disorders, including Parkinson's disease (PD), tremor, and dystonia [1], [2], [3].Several closed-loop DBS strategies have been used in clinical practice, including the on-off strategy [4], the dual-threshold strategy [5], and linear modulation strategy, these strategies could modulate beta power of PD patients and improve symptoms [6], [7].Although these strategies were effective at modulating the beta power of PD patients and improving symptoms, they still require optimization to achieve precise neuromodulation, improve symptoms comprehensively, and reduce side effects.Excessive beta power is a biomarker of PD, and beta power amplitude has been associated with symptom severity [8].Beta power may be suppressed by levodopa, with the degree of suppression correlating positively with symptom improvement [9], [10].Beta power also exhibits diurnal fluctuations, with increases during the day and reductions at night [11], [12], [13].Fluctuations due to medication metabolism and circadian rhythm introduce uncertainty and nonlinearity into the system, which may lead to suboptimal modulation performance.Therefore, neuromodulation to manage neural activity dynamics remains a substantial challenge.
Both model-free and model-based closed-loop DBS strategies can modulate neural activity to a preset target [14], [15], [16], [17], [18].The model-free proportional-integraldifferential (PID)-based closed-loop DBS strategy balances modulation performance and the ability to manage neural dynamics [16], [17].It avoids the error of modeling the closed-loop DBS system and provides the advantage of robustness and computational efficiency [19], thus offering great potential for clinical translation.However, the gains tuning is time-consuming and nonadaptive to neural activity dynamics, resulting in sub-optimal modulation performance [20].The backpropagation neural network was trained offline with eight different PD states to find PID gains which were also effective for an additional three PD states [21].Although the intelligent algorithm has a strong ability to solve system uncertainty and nonlinearity, its data dependency made real-time and online performances to be optimized.The self-tuning proportional controller was also designed to modulate beta-band oscillations in both a static firing-rate model and a static physiological cortical basal ganglia model [22].However, neural activity evolves over time due to various factors, such as circadian rhythm and medication rhythm.The control targets of beta power suppression and treatment requirements also vary between individuals, such as the electrode implantation location and constraints on stimulation increment, which may require the tuning of all three PID gains for better modulation performances, rather than only proportional control [16], [23], [24].Therefore, the online optimization of the amplitude adaptive neuromodulation strategy, based on transient and long-term dynamics of neural activities, remains unsolved.
To address this problem, we propose a particle swarm optimization (PSO) based amplitude adaptive neuromodulation strategy [25].PSO offers several advantages, including self-optimization, online tuning, and computational efficiency, which is suitable for closed-loop neuromodulation and has the potential for clinical translation.PSO has been used across many domains that apply optimization techniques in PID gain tuning, such as mechanical engineering, power systems, and industrial systems, proving to be suitable for realization and capable of optimization [26], [27].The remaining sections of this paper are organized as follows: Section II establishes a dynamic PD basal ganglia-thalamus-cortical circuit meanfield model, which can simulate the PD beta fluctuations in the presence of circadian rhythm and medication metabolism.Section III presents a PSO-PID based amplitude adaptive modulation strategy that is designed to modulate beta power in long-term dynamic conditions.The proposed strategy was tested with the stochasticity of particles and different dynamic degrees of beta power, validating the strategy's robustness and potential for clinical translation.Finally, Section IV provides a discussion of the results and conclusions.

A. Long-Term Dynamic Model of the Basal Ganglia-Thalamus-Cortical Circuit
To establish a PD pathological circuit with circadian and medication rhythm, as well as verify the effectiveness of the proposed neuromodulation strategy, we modified a mean-field model of the basal ganglia-thalamus-cortical circuit [28], [29], [30].An implementation of this model in Python can be downloaded from ModelDB (https://senselab.med.yale.edu/modeldb/).
In the mean field model of the basal ganglia-thalamuscortical circuit, the parameters included connection strengths and firing thresholds among the subcortical nucleus and cerebral cortex (Table I).The bold parameters in Table I are  [28], [29], [30].Based on the influence of medication and circadian rhythm on the neural circuits, we established the dynamic process of medication administration by changing the connection strength between the D1/D2 dopamine receptors and the cerebral cortex with a preset rate in (1), i.e. v d 1 e (k) and v d 2 e (k) [31], [32], [33].To simulate the 24-h diurnal cycle, other parameters related to PD changed as sine functions as (1), i.e. v p 2 p 2 (k), v ee (k), v ie (k), v ei (k), v ii (k), v p 2 d 2 (k), θ p 2 (k) and θ ς (k) [11], [12].The ranges of parameters were set between the healthy state and the PD med-off state.For example, the range of v d 1 e (k) was [0.5, 1], and the medication was fully effective when v d 1 e = 1.v p 2 p 2 changed from -0.1 to -0.07, representing the dynamic process during the 24-h diurnal cycle.Moreover, 10% noise was added to each parameter to represent stochasticity in the basal ganglia-thalamus-cortical circuit.
where par m m represents the mean value of parameters related to circadian rhythm, and par m r represents the range of related parameters.δ (k) represents the noise added to the related parameters, and η represents noise strength.
To simulate the effects of DBS within the mean field model, the stimulus was incorporated as a direct current injection into the target nucleus of STN.Combining all structures discussed above, the model could (1) simulate the healthy and dynamic PD states with circadian and medication rhythm; (2) provide LFP signals; (3) provide beta oscillation, the power of which changed with medication and circadian rhythm; and (4) respond to closed-loop DBS.

B. Diagram of the PSO-PID Based Amplitude Adaptive Neuromodulation
A PSO-PID based controller was designed to implement the amplitude adaptive neuromodulation strategy.First, a schematic diagram of beta power when taking medication over four cycles across 24 h was proposed, corresponding to the dynamic and stochastic basal ganglia-thalamus-cortical circuit mean field model.The dynamic model involved the med-off period during 00:00-06:00 (l 1 ), and the medication to take effect (l 2 ) and wear off (l 3 ) periods.l 2 and l 3 were repeated four times, representing four medication cycles during 06:00-23:00.After the last cycle, beta power followed the circadian rhythm during 23:00-24:00, with the med-off PD state (l 4 ) (Fig. 2(a)).The closed-loop DBS strategy included two modules: the PSO-based gains searching module and the PID-based closed-loop DBS module.After initializing PSO, the online gains searching process started and found a set of feasible PID gains K p , K i , and K d , denoted as the global best position.Based on the exported PID gains, adaptive stimulation amplitudes were calculated by the PID controller (Fig. 2(b)).Note that each trial in the PSO based gains searching module was a process of closed-loop DBS with a PID controller with different gains K p , K i , and K d (Fig. 2(c)).After each stimulation process, the fitness value J t i was calculated to evaluate the performance of the tested PID gains.
For each trial in the PSO-based gains searching module and the PID-based closed-loop DBS module, the connection strength and firing threshold of the mean field model varied dynamically.The frequency and pulse width of stimuli were all set at 130 Hz and 60 µs, respectively.Other than K p , K i , and K d , the parameters of the PID controllers were held constant.The sampling rate was 1,000 Hz, denoted as f s .The LFPs were sampled from the basal ganglia-thalamuscortical mean field model and set as the feedback signal for the PID controller.The LFPs were band-pass filtered to obtain beta oscillations (12-35 Hz) with a Chebyshev I filter.The cut-off frequency of the band-pass filter was 10 Hz and 37 Hz respectively.Beta power was calculated by applying a rectification filter and smooth filter.The window length of the smooth filter was 200 ms.Then, beta power was converted to a log scale by calculating 20lg[b(k)].This conversion served to amplify the value of the feedback signal, which made it more robust for the PID controller to track the dynamic control target without changing gains in a single trial.The amplitude of the stimuli could be calculated every few pulses, resulting in different control frequencies, denoted as f c .If the stimulation amplitude was adapted at each pulse, the control frequency was 130 Hz.
In this paper, a dynamic control target was calculated for the PID controller, which included short-term targets and a longterm target, modulating beta power with multiple time scales.Short-term targets changed at each sampling point, calculated by ( 3) and ( 4).The diagram of the dynamic control target was illustrated in supplementary materials (Fig. S1).
where b(k) is the beta power at sampling point k, and b(k) is the average of beta power with n points before k.N s is the total number of data points for each closed-loop stimulation period.b l was the long-term target, representing the desired average of beta power in the longest time scale.In the mean field model, the longest time scale was 48 min.In practice, the longest time scale can be minutes, hours, days, months, and years.p∈ [0, 1] represents the modulation strength determined by how much the beta power would be modulated at each step.After setting b l and p, the short-term target b s was calculated at each sampling point.Based on the dynamic target, beta power can be modulated with different time scales from milliseconds to minutes.Through setting different modulation strengths and long-term targets, it could also provide individual control targets for different patients.

C. PSO-PID Based Control Strategy and Evaluation
Since each trial in the PSO based parameter searching module was a process of closed-loop PID-DBS with different controller parameters, we first proposed the PID control strategy.The error at each sampling point was defined as where N c = f s / f c represents the number of sampling points between every two stimuli.In this study, f s was 1,000 Hz and f c was 130 Hz.When k was an integer multiple of N c , the stimulation amplitude was calculated by the PID controller.The structure of the incremental PID controller was designed as where u(k) represents the adaptive stimulation amplitude whose upper and lower bound are denoted as u max and u min .u (k) represents the change in stimulation amplitude, bounded by u max and − u max .For clinical translation, u max was set as 0.069 mA, which was the mean value of the changing speed for the stimulation amplitude.Details of calculating u max was presented in part B of supplementary materials.The constraints of u(k) and u (k) guarantees the safety and tolerance of the gains searching process in clinic.e(k) represents the mean error between the short-term target and the feedback beta power within N c points.K p , K i , and K d represent the controller gains.To simulate the actual condition when the DBS is turned on, the initial stimulation amplitude u(0) was set to zero.For a complete searching space, the upper bounds of K p , K i , and K d were all set as 0.2.
To be viable for clinical translation, the PID gains needed to be optimized online according to neural activity dynamics.We employed the PSO algorithm to tune PID gains.In the PSO algorithm, each particle possessed a velocity vector and a position vector.The position vector represented the PID gains, while the velocity vector indicated the direction and step length of the PID gains during each generation.The fitness value was calculated after each trial and used to evaluate the tested PID gains.The best local position of the ith particle with the dimension of D was denoted as pbest i =( pbest i1 , • • • , pbest i D ).Meanwhile, the best position of the whole swarm was denoted as gbest=(gbest 1 , • • • ,gbest D ) and termed as the global best position, which was the desired PID gains.During each generation, the position and velocity of each particle were updated by ( 8) and ( 9) towards the global best position.
where ω is the inertia weight of the last velocity, c 1 and c 2 are acceleration constants that determine the extent to which the particle was influenced by the history of its individual best position and by the global best position.In this study, ω was set to 0.9, and c 1 and c 2 were each set to 2, which are typical values for these parameters [34], [35].r 1 and r 2 represented random numbers in the range of [0, 1].G max was the maximum generation and t≤G max was the generation index.N p was the number of particles and i= 1, • • • ,N p was the index of particles.v t i and x t i were the velocity vector and the position vector of the ith particle at the tth generation, respectively.K p , K i , and K d corresponded to three dimensions of the position vector.
The performance of tested PID gains was evaluated using three indices, including the mean error, the mean stimulation amplitude, and the stimulation variation, calculated by ( 10)- (12).The mean error between the short-term target and the feedback beta power of the stable stimulation period was calculated as (10), measuring the neuromodulation accuracy.
where N s represents the total number of data points for each trial and N 1 represents the starting point for calculating the index.Furthermore, N s was 10,000 with the sampling point of 1,000 Hz.N 1 was 7,000, by choosing the last 30% of each trial for performance evaluation.The mean stimulation amplitude was calculated as (11).
The stimulation amplitude variation-i.e., the variance of u(k)-was calculated as (12).
To simultaneously consider modulation accuracy, stimulation amplitude, and ramping rate, the three indices were normalized to the same order of magnitude.Specifically, P/u max was calculated as the normalized mean stimulation amplitude with the order of 10 −1 , denoted as P ′ .Note that E is typically on the order of 10 −2 and V s is typically on the order of 10 −3 .E and V s were multiplied by 10 and 100, denoted as E ′ and V ′ s , respectively.Finally, a comprehensive index was calculated as (13), which was the weighted summation of these three indices.
In order to reflect the degree of demand for different indices, ω j1 , ω j2 , and ω j3 can be assigned different values according to individual clinical requirements.As an example, ω j1 , ω j2 , and ω j3 were set to 0.5, 0.25, and 0.25 in the present study, indicating that the modulation accuracy was guaranteed first.The objective of the particles was to minimize J.The global best fitness threshold was denoted as J min .The implementation procedure of the PSO algorithm for optimally tuning the PID parameters of the closed-loop DBS system is summarized as follows.For further explanation, the pseudocode is also presented (Supplementary Algorithm 1).
Step 1: Specify the PSO population size N p ; the maximum generation G max ; particle dimension D; the inertia weight factor ω; the acceleration constant c 1 and c 2 ; total sampling points of each trial N s ; the starting point for performance evaluation N 1 ; sampling points between each calculation of stimulation amplitudes N c ; long-term target of closed-loop DBS b l ; window length of the dynamic control target n; modulation strength of the dynamic control target p; the maximum stimulation amplitude u max ; the constraint of the stimulation amplitude increment u max ; the global best fitness threshold J min ; weights of the comprehensive evaluation index ω j1 , ω j2 , and ω j3 .To summarize, in order to solve the problem of online optimized amplitude adaptive neuromodulation based on neural activity dynamics, a closed-loop DBS strategy based on PSO-PID was designed.We applied PSO to realize the online gains tuning of PID-based closed-loop DBS, providing the potential for clinical translation.

A. Modeling Long-Term Dynamic and Stochastic Basal Ganglia-Thalamus-Cortical Circuits
The subthalamic LFPs showed dynamic characteristics from milliseconds to tens of minutes, following both the medication and circadian rhythm.Changing rates of v d 1 e and v d 2 e , that is, step d 1 e and step d 2 e in (1) were presented in Table II.In the off medication period, changing rates were set at zero.
The LFPs, the time-frequency analysis, and the beta power of the dynamic and stochastic mean-field model are presented in Fig3.Details of the variation in , and θ ς (k) are presented in supplementary materials (Fig. S2).Moreover, considering fluctuations in connection strengths and firing thresholds, 10% white noise was added to the ten parameters discussed above during the entire process.Yellow dashed lines in Fig. 3 mark the moment of med-off and med-on over four medication cycles.
The LFPs, time-frequency analysis, and beta power are also illustrated in Fig. 3.In the dynamic model, the beta power reached a high level before taking medication and decreased after taking medication.Beta power also followed the circadian rhythm, reaching the lowest level at 12 min and the peak level at 36 min.For the concern of computational time, the numerical simulation was implemented using a shrunken time scale of 1/30.The dynamic model represented stimulation Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.processes lasting for 48 min, involving 12 min for the medoff period (l 1 ), 1.5 min for the medication to take effect (l 2 ), and 7 min for the medication to wear off (l 3 ).l 2 and l 3 were repeated four times, representing four medication cycles.After the last cycle, beta power followed the circadian rhythm with the med-off PD state of 2 min (l 4 ).
The effect of different stimulation amplitudes on beta power of subthalamic LFPs was investigated by applying continuous DBS (cDBS) at the STN.The stimulation experiments were repeated five times for each condition.
The dots in Fig. 4 represent the mean values of the five experiments, with the shaded area representing variance.The beta power of STN in the healthy model was about −22.9 dBmV (green dashed line in Fig. 4).The stimulation amplitude increased from 0 mA to 6 mA, with a step of 1 mA.Additionally, the stimulation frequency and pulse width were 130 Hz and 60 µs, respectively, which are consistent with commonly used clinical stimulation parameters.The beta power was gradually suppressed when the stimulation amplitude increased in the dynamic model.Moreover, the beta power almost reached a healthy level when the stimulation amplitude was 5.8 mA.The maximum stimulation amplitude of the PID controller, i.e., u max in ( 6), was set as 5.8 mA.For complete testing, b l in (4) was selected in the range of −19 dBmV and −23 dBmV.

B. Amplitude Adaptive Neuromodulation Strategy
The PSO-based gains searching process and the long-term closed-loop DBS with the obtained PID gains are presented.An example with five particles and five generations is shown in Fig. 5. Magnification of colored areas are presented in the supplementary materials (Fig. S3).
During the gains searching process, the middle level of feasible beta power was chosen as the long-term target, meaning b l in (4) was set as −21.09dBmV, corresponding to the beta power when the stimulation amplitude was 3 mA in Fig. 4. The modulation strength p was set as 0.2, and the length for b(k) was 50 ms.PID-DBS was performed with each particle position for a duration of 10 s.The process of gains searching within 250 s is presented in Fig. 5(a), illustrating LFPs, stimulation amplitude, u(k), dynamic control target, and closed-loop beta power.The red shaded area marks the global best gains.The blue shaded area marks PID gains whose fitness value was at the middle level, denoted as the suboptimal gains.The yellow shaded area marks PID gains with the lowest fitness value, denoted as the inappropriate gains.In this example, the global best gains were K p = 0.01404247, K i = 0.00120627 and K d = 0.0148426; the suboptimal gains were K p = 0.09803478, K i = 0.00137931, and K d = 0.09767946; and the inappropriate gains were K p = 0.0375386, K i = 0.00014538 and K d = 0.2.
The performance of particle position-the fitness value calculated by ( 13)-converged with increasing generation index (Fig. 5(b)).With the global best gains, the average stimulation amplitude of PID-DBS was expected not to be higher than cDBS which suppressed beta power to the same long-term target.Moreover, the mean error and the stimulation amplitude variation were expected to be close to zero.Therefore, based on ( 13) and the normalization principle described in the Method section, 0.25× (3/5.8)= 0.13, the minimum bound of the fitness value J min was set as 0.15.The global best fitness value converged to 0.1446 at the fifth generation, which was lower than J min , and completed the gains searching process.The positions of each particle and the global best position of each generation are illustrated in Fig. 5(c).The colored dots represent the positions of particles at each generation.The larger dots represent an increasing number of generations.The red dots represent the global best position after each generation.The changes of PID gains K p , K i , and K d during five generations are presented in Fig. 5(d)-(f), respectively.In this example, the obtained PID gains was K p , K i , K d = [0.01404247,0.00120627, 0.0148426].
The performance of PID-DBS was dependent on the selection of PID gains.Fig. 6 compared the effectiveness of PID-DBS with different gains, showing the stimulation amplitude, the stimulation amplitude increments, and the closed-loop beta power.In this experiment, the long-term target was -21.09 dBmV, the modulation strength was 20%, and the length of b(k) was 50 ms.The red, blue and yellow solid lines represented the global best gains, the suboptimal gains, and the inappropriate gains respectively.
The stimulation amplitudes corresponding to the global best gains and the suboptimal gains both showed fluctuations consistent with circadian and medication rhythms, while that of the inappropriate gains did not (Fig. 6(a)).The variance of stimulation amplitude increments was the smallest with the global best gains, while the inappropriate gains corresponded to the largest variance (Fig. 6(b)), meaning the stimulation sequence was the smoothest with the global optimal gains.As the three horizontal solid lines shown in Fig. 6(c), the long-term beta power average of the global best gains, the suboptimal gains, and the iappropriate gains were −20.97 dBmV, −20.02 dBmV, and −19.10 dBmV, respectively.Unlike the global best gains, neither the suboptimal gains nor the inappropriate gains could modulate beta power to the long-term target.Additionally, the mean error between the closed-loop beta power and shor-term targets with the global best gains, the suboptimal gains, and the inappropriate gains was 0.01 dBmV, 0.19 dBmV, and 0.39 dBmV, indicating that the global best gains had the highest modulation accuracy.To summerize, neuromodulation with the global best gains had the highest modulation accuracy with the most stable stimulation, and better reduced the fluctuations of beta power caused by circadian and medication rhythms.
To verify the effectiveness and robustness of the obtained global best gains with different dynamic control targets in the dynamic and stochastic basal ganglia-thalamus-cortical model, different combinations of long-term targets and modulation strengths in (4) were investigated.The long-term target was set at −19 dBmV, −20 dBmV, −21 dBmV, −22 dBmV, and −23 dBmV, whereas the modulation strength was set as 20%, 40%, 60%, 80%, and 100%.For each parameter combination, five independent random experiments were repeated using the 48-min process of the dynamic model.Examples of closed-loop DBS with the global best PID gains presented in supplementary materials (Fig. S4).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.To summarize, PSO could search PID gains online with a self-optimization process and high computational efficiency within minutes.With the obtained PID gains, beta power could track dynamic control targets and be modulated with the time scale of milliseconds and minutes in the long-term dynamic basal ganglia-thalamus-cortical mean field model.

C. Closed-Loop Deep Brain Stimulation With Different Searching Moments
From the results presented above, the start time of gains searching was the beginning of the dynamic model, which included the circadian rhythm but not the medication rhythm.In order to verify the influence of different dynamic degrees of neural activitiy on the amplitude adaptive neuromodulation strategy, we investigated the effectiveness of the proposed strategy at different start times of the gains searching process.
The times 12 min, 24 min and the 36 min were taken as examples.Among them, the 12 min marked the lowest point of beta power in the circadian rhythm and the beginning of the first medication cycle; the 24 min was the middle moment when beta power increased with the circadian rhythm, occurring in the medication failure process of the second medication cycle; and the 36 min was the moment when beta power was at the highest point in the circadian rhythm, occurring in the medication failure process of the third medication cycle.In the following experiments, the bound of particle position vectors was 0.2, the number of particles was 5, and the maximum generation was 10.The dynamic control target was set with a long-term target of −21.09 dBmV, modulation strength of 0.2, and length for b(k) of 50 ms.The closed-loop stimulation was performed for 10 s at each particle position.
Random initialized position and velocity of particles resulted in different searching processes among trials.Due to the stochasticity of both PSO and the mean field model, we conducted five searches at each moment, and five stimulation experiments with every group of PID gains  (10).(c) Stimulation amplitude, P, calculated using (11).(d) Stimulation variation, variance of ∆u(k), calculated using (12).(e) Long-term mean error of the beta power average and the long-term target.The duration of calculating the dynamic quantity was 50ms.The modulation strength was 20%, and the long-term target was −21.09 dBmV.The colored bars represent the mean value of five repeat experiments, and the black dots represent the five repeat experiments for five gains at each searching moment.
obtained from each search.PID gains obtained through the five searches at three moments are shown in supplementary materials Table I.At different starting moments, PID gains that met the closed-loop neuromodulation requirements could be automatically searched online, based on the PSO algorithm.
Taking gains 1 at 12 min as an example, the LFPs, stimulation amplitude, u(k), dynamic control target, and beta power of three searching moments are illustrated in Fig. 8(a).Examples of 24 min and 36 min are presented in supplementary materials (Fig. S5).The stimulation amplitudes fluctuated with the circadian and medication rhythms.The beta oscillations tracked the short-term target and were suppressed to the long-term target.The modulation accuracy remained high with a mean error of less than 0.05 dBmV across all experiments (Fig. 8(b)).Furthermore, the average stimulation amplitude was not higher than cDBS which suppressed beta power to the same long-term target (Fig. 8(c)).The variance of u(k) was near zero, indicating the stability of closed-loop DBS (Fig. 8(d)).The long-term mean errors are also presented, which were less than 0.15 dBmV for all experiments (Fig. 8(e)).By observing all obtained gains, the balance of K p , K i , and K d determined the modulation performances, rather than any single element having a more significant influence.
To summarize, the proposed amplitude adaptive neuromodulation strategy was able to optimize the controller gains online and accomplish the long-term modulation of beta power, with different dynamic degrees of neural activity, demonstrating its potential for clinical translation.

IV. CONCLUSION AND DISCUSSION
An amplitude adaptive modulation strategy with PSO and PID controller was developed to modulate beta oscillations over long-term dynamic conditions.
First, the long-term dynamic model of the basal gangliathalamus-cortical circuit was established by continuously changing connection strengths and firing thresholds, as well as adding white noise.Considering the circadian and medication rhythms, the long-term dynamic model is closer to practical clinical situations than static models [9], [10], [11], [12].However, the effect of transitions between different sleep patterns on closed-loop DBS performances may be an issue that needs to be addressed in the future [36], [37].Abrupt state transitions may bring disturbances into LFPs.We verified the effectiveness of PID-DBS with the optimal gains searched by PSO under disturbances (Fig. S6).
Second, in order to precisely modulate dynamic beta oscillations, a PSO-PID based amplitude adaptive neuromodulation strategy was developed, including a gains searching module and a closed-loop DBS module.The proposed strategy tuned all three gains of PID controller, which may result in better performances than single gain [16].The strategy was designed by taking several factors related to clinical application into consideration.The PSO-PID based closed-loop DBS strategy provides smooth adjustment of stimulation parameters instead of rapid changes, which may cause side effects such as tingling sensations.The time required to tune the parameters of PSO-PID gains is usually around several minutes, allowing for optimization of stimulation across essential physiological and pathological rhythms, such as circadian rhythms and medication influence, while saving computational power.Safety measures include limiting the maximum stimulation to 5.8 mA and the maximum increment to 0.069 mA per pulse.
Third, both the gains searching module and the closedloop DBS module involved a dynamic control target.The beta power was modulated to track the dynamic control target with the obtained PID gains, which accomplished neuromodulation in the time scale of both milliseconds and minutes.The settings for the long-term target and the modulation strength verified the robustness of the PSO-PID based amplitude adaptive neuromodulation and also provided a personlized programming approach for closed-loop DBS [14], [38].Different patterns of the long-term target were also tested to show the effectiveness of the PSO-PID based strategy (Fig. S7).
Fourth, the influence of different dynamic degrees of beta power was explored by starting the gains searching process at different moments of the dynamic model.At the beginning, the lowest, the middle, and the highest level within single circadian cycles and four medication cycles, PID gains were obtained and the beta power was modulated to track dynamic targets with a modulation accuracy of less than 0.05 dBmV.
Due to the self-optimization and model-free characteristics of PSO and PID, the PSO-PID-based amplitude adaptive neuromodulation strategy has the potential to handle new tasks with computational efficacy, in various clinical settings, e.g., different control frequencies, stimulation amplitude increments, recording and stimulation sites, and types of neurological and psychiatric diseases.The strategy can additionally be transferred to other nuclei, such as the GPi, as well as other diseases, if feasible biomarkers are available to provide feedback signals.However, the individual variability can be a great challenge when applying PSO-PID based closedloop DBS strategy in clinic.Other data-driven methods may be helpful in optimizing the proposed strategy in the future.
The amplitude adaptive modulation strategy could meet the requirements of closed-loop precise neuromodulation of longterm conditions.However, the obtained PID gains may be infeasible due to various neural activity dynamics,such as the different time or action rate of medication intake, sleep-wake transitions, and the progression of disease over time scales of days, months, and years.In these situations, an evaluation module could be added to the closed-loop strategy to restart the PSO-based gains searching module and update the PID gains.For example, evaluation indice can be calculated every n seconds and setting a threshold for re-triggering the searching process.Considering the stochasticity and disturbances in practice, to prevent frequent triggers due to sporadic factors, a new search may be re-triggered only when the index exceeds the threshold for m consecutive evaluations.
Compared with other typical methods, such as the on-off strategy [4], the dual threshold strategy [5], and the linear modulation strategy [6], [7], PSO-PID based closed-loop DBS modulated beta power amplitude precisely to track dynamic control targets with automated tuned PID gains.Moreover, published research showed that results based on the computational model have the feasibility of translating to animal models and clinical [39].In the future, the closed-loop DBS strategy will be implemented in a research platform involving animal experimentation [40], which will further facilitate the translation of this approach to clinical practice.

Fig. 1 .
Fig. 1.Structure of the basal ganglia-thalamus-cortical mean field model.The letters in the brackets, i.e., r, s, d 1 , p 1 , etc., denote the corresponding nucleus.Black and red solid arrows represent constant excitatory and inhibitory connections.Black and red dotted arrows represent dynamic and stochastic excitatory and inhibitory connections respectively.Blue arrows represent the addition of Gaussian white noise to the output of each nucleus.Grey arrows represent the recording and stimulation at STN. Adapted from van Albada and Robinson, 2009, and Grado et al. 2018[28],[29],[30].

Fig. 2 .
Fig. 2. Diagrams of the dynamic model and the closed-loop neuromodulation strategy.(a) Schematic diagrams of the beta power when taking medication over four administration cycles in 24 h.Starting at 00:00, the beta power fluctuated following the circadian rhythm.The medication was taken at 6:00 a.m. for the first time.The beta power decreased to a med-on level after 0.75 h and increased to a med-off level after 3.5 h in each medication cycle.(b) Diagram of the amplitude adaptive neuromodulation strategy.The red square represents the PSO-based gains searching module, and the blue square represents the closed-loop DBS module with PID gains calculated by PSO.(c) Sequence diagram of the PSO-PID based amplitude adaptive neuromodulation strategy.

Fig. 3 .
Fig. 3.Modeling of the circadian rhythm and four cycles of levodopa administration in the basal ganglia-thalamus-cortical mean-field model.STN LFPs, time frequency analyses, and the beta power of the mean-field model over four cycles with the circadian rhythm.The beta power showed the dynamics corresponding with changes in parameters of the mean-field model.The yellow dashed lines represent the moments of med-on and med-off over four cycles.

Fig. 4 .
Fig. 4. Effect of stimulation amplitudes on the beta power of the STN.The black dotted line represented the beta power of the dynamic model.The green dashed line represented the beta power of healthy model.The beta power decreased as the stimulation amplitude increased from 0 mA to 6 mA.The shaded areas indicate the variance in five random experiments.

Fig. 5 .
Fig. 5. PID gains searching process with PSO and closed-loop DBS in the dynamic model.(a) Stimulation with each particle position during five generations.The red, blue, and yellow shaded area marks the global best gains, the suboptimal gains, and the inappropriate gains.(b) Fitness values of each particle and the global best fitness value over five generations.(c) Positions of five particles and the global best positions during five generations.The smallest dots represent the first generation, while the biggest dots represent the last generation.(d)-(f) PID gains, K p , K i , and K d of five particles and the global best positions during five generations.The blue lines represent five particles and the red line represents the global best position.

Fig. 6 .
Fig. 6.Comparison of closed-loop DBS with different PID gains in the dynamic model.(a) The stimulation amplitude of closed-loop DBS.(b) Increments of stimulation amplitudes.Magnification of 805 s to 815 s, showing detail (yellow shaded area).(c) Beta power of closedloop DBS.Red dotted lines represented closed-loop DBS with the global best gains, whose K p = 0.01404247, K i = 0.00120627 and K d = 0.0148426.Blue dotted lines represented closed-loop DBS with the suboptimal gains, whose K p = 0.09803478, K i = 0.00137931 and K d = 0.09767946.Yellow dotted lines represented closed-loop DBS with the inappropriate gains, whose K p = 0.0375386, K i = 0.00014538 and K d = 0.2.The red, blue and yellow horizontal solid lines represented the average of beta power respectively.The length of calculating the dynamic quantity was 50ms.The modulation strength of the dynamic target was 20%.The long-term target was −21.09 dBmV.The searching process started at 0 min.The yellow dashed lines represent the moments of med-on and med-off over four cycles.

Fig. 7 .
Fig. 7. Performance of closed-loop DBS in the dynamic model with the global best PID gains searched by PSO.(a) Modulation accuracy, E, calculated using (10).(b) Stimulation amplitude, P,calculated using (11).(c) Stimulation variation, variance of ∆u(k), calculated using (12).(d) Long-term mean error of the beta power average and the long-term target.The modulation strength was 20% and the long-term target was −21.09 dBmV.The length of calculating the dynamic quantity was 50ms.The searching process started at 0 min.The colored bars represent the mean value of five repeat experiments.The black dots represent the five repeat experiments for each combination of long-term targets and modulation strengths.

Fig. 8 .
Fig. 8. Influence of different searching moments.(a) Closed-loop DBS at 12 min.(b) Modulation accuracy, E, calculated using(10).(c) Stimulation amplitude, P, calculated using(11).(d) Stimulation variation, variance of ∆u(k), calculated using(12).(e) Long-term mean error of the beta power average and the long-term target.The duration of calculating the dynamic quantity was 50ms.The modulation strength was 20%, and the long-term target was −21.09 dBmV.The colored bars represent the mean value of five repeat experiments, and the black dots represent the five repeat experiments for five gains at each searching moment.

TABLE I PARAMETERS
OF THE BASAL GANGLIA-THALAMUS-CORTICAL CIRCUIT IN THE MEAN FIELD MODEL related to different PD states.For example, v d 1 e represents the connection strength between D1 receptors and the cerebral cortex.The bold parameters in the dynamic model shown in TableIvaried at each sampling point k, while they were held constant in the PD med-off model.To model the inherent stochasticity of the basal ganglia-thalamus-cortical circuit, white noise was added to the output of each nucleus.

TABLE II CHANGING
RATES OF v d 1 e AND v d 2 e OVER FOUR For the position K p ,K i , K d of each particle, calculate the fitness value.Step 4: Compare the obtained fitness values of single particles and save the position with the minimum value during its searching history as pbest.Meanwhile, compare fitness values of all particles and save the position with minimum value in the population as the global best position gbest.Step 5: Modify each particle position and velocity according to (8) and (9).Step 6: If (1) the number of generations reaches the maximum G max , or (1) reach the global best fitness threshold, go to Step 7; otherwise, go to Step 3. Step 7: Obtain the PID gains gbest.