Abstract
A closed-loop neuromodulation automatically adjusts stimuli according to brain response in real time. It is viewed as a promising way to control medically intractable epilepsy. A suitable closed-loop modulation strategy, which is robust enough to unknown nonlinearities, dynamics, and disturbances, is in great need in the clinic. For the specialization of epilepsy, the Jansen’s neural mass model is utilized to simulate the undesired high amplitudes epileptic activities, and active disturbance rejection control is designed to suppress the high amplitudes of epileptiform discharges. With the help of active disturbance rejection control, closed-loop roots of the system are far from the imaginary axis. Time domain response shows that active disturbance rejection control is able to control seizure no matter whether disturbances exist or not. At the same time, frequency domain response presents that enough stability margins and a broader range of tunable controller parameters can be obtained. Stable regions have also been presented to provide guidance to choose the parameters of active disturbance rejection control. Numerical results show that, compared with proportional-integral control, more accurate modulation with less energy can be achieved by active disturbance rejection control. It confirms that the active disturbance rejection control-based neuromodulation solution is able to achieve a desired performance. It is a promising closed-loop neuromodulation strategy in seizure control.
Keywords
Introduction
Epilepsy is a typical chronic neurological disorder resulting from the development of synchronous firing in a massive group of neurons. It afflicts nearly 1% of the people worldwide. The newest research reports illuminate that almost 1/3 patients never achieved seizure freedom by antiepileptic drug therapy or even experience unexpected side effects. 1 For those drug-resistant patients, only a small fraction qualifies for surgical resection of the epileptogenic zone. Patients, who are not suitable for surgery and who are still suffering seizures after surgeries, attract much attention not only from the physicians and surgeons but also from the social societies and governments. Since epilepsy imposes significant psychological, physical, and financial burdens on those patients and their families. In order to address it, numerous new and effective treatments have been proposed. In the last decades, neuromodulation technologies, such as vagus nerve stimulation, deep brain stimulation, and responsive neurostimulation, have been utilized to alleviate the intractable epilepsy. 2 It has been confirmed that neuromodulation is one of the most impressive achievements in epilepsy treatments. Additionally, invasive, especially, noninvasive neurostimulation techniques will become more effective and popular in the next couple of years. 3
For the stimulation of nerve tissues, there are two modes. One is the open-loop way, that is, stimulus signals are generated from the device and delivered to the brain. Such mode provides an invariant train of stimulatory pulses in a fixed frequency. It is independent of ongoing and varying neuronal activities. In contrast to the scheduled stimulation, the other is the closed-loop mode. It delivers stimulatory pulses upon the detection of seizure activities, that is, it generates the stimulatory pulses according to the real-time neural activities. Such a system is a ‘clever’ system which is flexible and individually adaptable to patients’ specific ictal pattern. However, by comparison with the open-loop mode, it is technically more challenging.
For the closed-loop mode, the challenge is the ‘on-demand’ stimulation strategy. It is known that system control theory is just the method generating the control signals based on the real-time system outputs. From this point of view, system control theory is the proper approach to address such a challenge. In addition, based on the closed-loop system, high efficacy of neurostimulation, a lower number of stimuli, and fewer side effects can be obtained. 4 Therefore, the closed-loop mode may be a promising method to suppress epileptic seizures. However, for improvements in seizure suppression, a ‘clever’ stimulation strategy, that is, an advanced control algorithm, which is able to automatically and adaptively generate stimuli with appropriate amplitudes and frequencies, is in great need. In order to acquire improvements, much effort has been made. 1,5 –11 Proportional–integral (PI) control 5 and proportional–derivative (PD) control 6 have been designed to suppress epileptic activities simulated by a neural mass model (NMM). Relationship between parameters of PI/PD control and parameters of an NMM has been analyzed to provide insights into the mechanism of suppressing epileptic activities. Considering nonlinearities and uncertainties of a neural system, a fuzzy inference system is taken to acquire the parameters of proportional–integral–derivative (PID) control, that is, fuzzy PID control is designed to restrain the epileptiform spikes. 7 For minimizing control energy and estimating system states from noisy measurements, unscented Kalman filter (UKF) 8,9 has been employed, and the UKF control strategy has been proposed for the NMM and the Wilson–Cowan equations. Based on the Pinsky–Rinzel model, input–output linearization (IOL) approach is designed to terminate epilepsy. 10 Reported approaches are mostly designed to suppress the epileptic activities simulated by those computational epileptic models. For an epileptic animal model, electroencephalogram (EEG) features are used to detect seizures, and stimulation for seizure elimination is applied in a predetermined mode with a fixed frequency and amplitude. 11 –14 It does greatly reduce the rate of ictal events, 4 however, such mode is not able to adjust stimulations in response to the specific electrical activities. Therefore, it cannot provide a personalized solution.
If a stimulation protocol is adaptive, and the closed-loop control algorithm is able to adjust its amplitudes and frequencies based on the real-time electrical activities, better control of seizure activities is able to achieve. Proportional feedback simulation has been tried in rats, 15 mean reduction in amplitude variance is significant. It shows that proportional feedback stimulation is a promising way to suppress seizure activities, and more complicated control algorithms for generating feedback stimulation may provide further improvements in seizure suppression. 4,15 However, for improving the seizure suppression, larger feedback gains have been tested and it has been found that higher gains are unable to maintain the stability of a system. As a result, the animals used were killed. Obviously, stimulation parameters are determined by physicians or surgeons empirically, and it cannot guarantee the safety of animals or even the patients. It may be effective for some patients, but it still needs some improvements. Firstly, if stable range of a closed-loop neuromodulation system can be provided in advance, one can guarantee animals’ or even the patients’ safety. Additionally, in a stable range, one may get a group of proper stimulation parameters so as to get better modulation. However, such requirements cannot be satisfied in animal experiments or even clinic trials on epileptic patients. Solid guidance for both closed-loop neuromodulation and determining a group of better stimulation parameters can be obtained by numerical simulations with a computational model of epilepsy.
In this article, a computational model is also utilized, and closed-loop neuromodulation by an advanced control algorithm is designed. Jansen’s NMM is characterized by interactions of the interlinked excitatory and inhibitory feedback loops. It is able to successfully simulate epileptic activities that are similar to those observed in physiological experiments. Therefore, it has been commonly utilized in most of the studies. 5,6,16 –23 PI/PD control 5,6 is able to get its control signals when tracking errors exist, however, it is a passive control approach. Such an approach depends on tracking errors to generate the control signals, it is susceptive to nonlinearities and uncertainties, and it needs more energy to modulate epileptic activities. Fuzzy PID 7 is robust, but fuzzy rules depend largely on experiences. Additionally, fuzzy rules lower system performance more or less. For the complexities and specialties of epilepsy, it is difficult to set a faithful model. Additionally, uncertainties and nonlinearities are ubiquitous in neuromodulation systems. Therefore, rather than PID, UKF, or IOL approaches, an advanced control algorithm, which is robust enough to those undesired factors, is able to deliver stimulation signals adaptively to real-time seizure activities, reduce the average number of stimulations and the adverse effects, and prolong the battery life, is in great need. On the other hand, closed-loop neurostimulation holds the greatest therapeutic potential, 24 but there still exist numerous obstacles in realizing the closed-loop stimulations. A very conspicuous one is the contradiction between the computational complexity of a closed-loop neuromodulation algorithm and its realization. Current systems have limited processing abilities. It is not practical to run computationally complex algorithms smoothly, although those algorithms are promising. In other words, they are too complex to be realized in the real-time modulation. 24
From the analysis above, a clear picture can be drawn. Closed-loop neuromodulation attracts much attention and it is increasingly utilized in the United States, 24 however, neuromodulation with a fixed frequency and amplitude or realized via simple closed-loop control approaches need to be improved. Additionally, neuromodulation realized by complex algorithms is not practical for clinical applications. Therefore, a control algorithm, which can achieve the desired performance and relatively simple to be realized on a chipset, is of great importance.
In recent years, disturbance estimation and rejection approaches attract much more attention. For example, disturbance observer-based sliding model control has been designed for a class of underactuated robotic system
25
and mobile wheeled inverted pendulum systems.
26
A flexible spacecraft
27
and polymer electrolyte membrane fuel cell systems
28
have also employed disturbance observer-based control techniques. Active disturbance rejection control (ADRC), a special kind of disturbance estimation and rejection approaches proposed by Han in the 1990s, is a promising way to deal with nonlinearities and uncertainties. By virtue of an extended state observer (ESO), one can get system output, the derivatives of the system output, and the generalized disturbance. Then, a control signal can be constructed. Undesired factors, which affect system performance, are estimated and cancelled before they corrupt the system output. Therefore, ADRC is an active control approach and it is capable of providing proper neuromodulation stimuli according to real-time EEG signals. Moreover, it can be seen later that ADRC is simple enough to be realized. Therefore, in recent years, ADRC has been accepted and achieve satisfactory performance in numerous fields.
29
–32
From the engineering point of view, ADRC is the very suitable closed-loop neuromodulation solution for its practical design and satisfied performance. Main contributions are ADRC has been firstly introduced in suppressing epileptic activities, closed-loop performance has been confirmed via both time domain and frequency domain when excitatory and inhibitory gains of an NMM vary in a certain range, stable region has been presented to provide a solid guidance to determine the parameters of ADRC in seizure control, and advantages of ADRC over PI in seizure control have been confirmed from the time domain indexes, including root mean square error (RMSE), mean absolute error (MAE), integral of time-multiplied absolute-value of error (ITAE), and the energy expenditure in the neuromodulation.
In this article, the underlying reasons for epilepsy are analyzed by frequency analysis of the NMM. Then, ADRC is taken to suppress epileptiform activities generated by an NMM. Parameters and performance are discussed. The following parts are organized like this. The NMM and its performance are presented in the second section. ADRC is introduced in the third section. Seizure control by ADRC, root loci, stable regions, frequency domain performance, and numerical simulations are presented in the fourth section. Finally, concluding remarks are drawn in the fifth section.
The neural mass model and its performance
The neural mass model
Epileptiform activities, resulting from the lack of balance between the excitation and inhibition, can be described by an NMM. Its structure can be illustrated in Figure 1. 6,16

Scheme of the NMM. NMM: neural mass model.
As shown in Figure 1, an NMM model is composed of three coupling subpopulations, including the main subpopulation (in the middle of the structure), the excitatory feedback subpopulation (on the top of the structure), and the inhibitory feedback subpopulation (at the bottom of the structure). The main subpopulation comprises principal cells that receive feedback signals from both the excitatory and the inhibitory subpopulations.
The excitatory and inhibitory subpopulations transform the average presynaptic firing rates into the average postsynaptic membrane potentials. Each of them can be represented by following impulse responses, 6,16 that is
where He
and Hi
are excitatory and inhibitory average synaptic gains, and
Based on impulse responses given in equation (1), dynamics of the excitatory and inhibitory subpopulations can be obtained and they can be described by two second-order ordinary differential equations
where
Nonlinear static function S presented in Figure 1 transforms the average membrane potential into the average firing rate, and it can be depicted by a sigmoid function
where e 0 determines the maximum firing rate, v 0 is the postsynaptic potential when 50% firing rate is achieved, r is the steepness of the sigmoid function, and v is the average presynaptic membrane potential.
Interactions among main cells, excitatory and inhibitory interneurons are characterized by four connectivity constants C 1 through C 4. It includes the total number of synapses established by interneurons. An excitatory input is represented as p(t), which is utilized to simulate influences from neighboring areas and modeled by the white noise. 5 Output of an NMM, that is, y(t), is the average synaptic activities of the pyramidal cells, which can be interpreted as the EEG signals. 5,6
Nominal values of those parameters can be found in Table 1. 6,16
Interpretation and values of the parameters in an NMM.
NMM: neural mass model.
Characteristic of the NMM
In Figure 1, only sigmoid function S is nonlinear. In order to analyze the NMM by linear system theory, nonlinear function S is linearized. According to small perturbation around the equilibrium point, dynamics of S depends on its slope at the equilibrium point. It can be defined as the sigmoid ‘gain’ K s
Remark 1
Linearization is utilized so as to present the root loci and frequency responses in the following sections. It does not affect the analysis of ADRC, since ADRC is less dependent on the faithful system model in its design and analysis. Additionally, the local linearization model also coincides with the properties of the nonlinear NMM model in terms of the influence of He s and Hi s.
Obviously, the local linearization approach is much simpler, and it is enough in the analysis. 5,6 After the Laplace transformation, the transfer function of equation (2) can be obtained
where

Scheme of the linearized NMM. NMM: neural mass model.
According to Figure 2, the transfer function of the NMM is
where Y(s) and P(s) are the Laplace transformation of y(t) and p(t), respectively.
Substitute
where
As we know, the unbalance between the excitation and inhibition results in epileptic activities. In the NMM, homeostasis can be achieved by regulating the excitation and inhibition. Abnormal values of average excitatory synaptic gain He or average inhibitory synaptic gain Hi generate the epilepsy. Here, by increasing He or decreasing Hi , we can see how He and Hi influence the dynamics of the NMM vividly from the root loci.
Note that the abnormal value of the excitatory gain is larger than the nominal value, and the abnormal value of the inhibitory gain is smaller than the nominal one. Therefore, with an attempt to see the influence of He
and Hi
clearly, firstly, let the excitatory gain He
range from 3.25 to 10, and the inhibitory gain Hi
still chooses the nominal value (i.e.

Root loci of the NMM. NMM: neural mass model.
Since the NMM is a six-order system, there are six branches. When He and Hi are abnormal, Figure 3(a) and (b) presents unstable branches. Obviously, abnormal excitatory gains or over excitations result in roots that are far from the imaginary axis in the right half s plane (see Figure 3(a)). It means larger excitation does lead to epileptic activities with higher probability. At the same time, decreasing inhibition also makes the roots of the NMM be still far from the imaginary axis (see Figure 3(b)). It also signifies that smaller inhibitory gains also result in epilepsy with higher probability.
For epilepsy, it results from the loss of homeostasis between the excitation and inhibition. Therefore, an approach, which is able to make the system be less sensitive to the excitatory gains and inhibitory gains, is in great need. Next, ADRC is introduced, and the neuromodulation will be checked.
Active disturbance rejection control
ADRC is a special approach, which is independent of concrete model of a system. It is able to deal with both tracking and disturbance rejection problem by its controller and ESO. For its satisfying performance, an increasing number of applications have been emerging. The basic structure of a closed-loop system by ADRC is presented in Figure 4.

Structure of a closed-loop system by ADRC. ADRC: active disturbance rejection control.
In order to illustrate ADRC, following system is considered
Here, y is the system output, and
z1, z2, and z3 are outputs of the ESO, which estimate a system output, derivative of a system output, and the generalized disturbance, respectively. An ESO can be designed as
where yr
is the set-point value;
After substituting equation (10) into (9), we have
where
Let estimation errors be
Obviously, when ESO performs well, that is,
Accordingly, the desired closed-loop transfer function is
By choosing proper kp and kd , desired closed-loop dynamics can be obtained. However, there are six parameters needed to be determined. To address it, let
then the bandwidth-parameterization based tuning approach can be obtained, 33 that is
Therefore, tunable parameters become
Taking Laplace transformation for both sides of equation (11), we have
Substituting equation (17) into (16), one has
Seizure control by active disturbance rejection approach
Control scheme
Seizure control by ADRC is designed. Based on equation (18), the equivalent structure of Figure 4 can be shown in Figure 5.

Equivalent structure of a closed-loop system by ADRC. ADRC: active disturbance rejection control.
Here, seizure control is discussed. Epileptic activities, usually, are characterized by high-amplitude oscillations. 34 –36 To suppress the high-amplitude epileptic activities, similarly to literatures., 1,5,6,8,23 the desired output is set to be zero, that is, yr = 0. Then, the structure of seizure control by ADRC can be shown in Figure 6.

Seizure control by ADRC. ADRC: active disturbance rejection control.
Remark 2
Here, ADRC is designed to overcome the influence of abnormal He s and Hi s. By an active control, the high amplitudes of epileptiform wave are expected to be suppressed. However, it is worth pointing out that the controlled EEG signals will not be zero, since EEG signals will not be zero even if He s and Hi s are normal. Setting value for seizure control is chosen to be zero so as to reduce the high amplitudes, and zero setting value can also be found in many other literatures (see literatures 1,5,6,8,23 ).
Performance analysis
Root loci analysis
From Figure 6, the open-loop transfer function can be obtained.
Obviously, if parameters of ADRC are fixed, root loci can be obtained. Here, let

Root loci of the closed-loop system by ADRC. ADRC: active disturbance rejection control.
Comparing with Figure 3, we can see that ADRC is able to correct the influence of the excitatory synaptic gain He and the inhibitory synaptic gain Hi . By ADRC, system dynamics will be more insensitive to He and Hi . It is the very aim of neuromodulation.
Stable region analysis
In order to determine the stable region of the closed-loop ADRC–NMM system in (b 0, ω c) parameter space, the Routh stability criterion is utilized. Based on the open-loop transfer function (20), we have the closed-loop characteristic equation
equation (21) is a ninth-order equation. It is somewhat complex. Let

Stable region of ADRC in suppressing epileptic activities generated by an NMM. ADRC: active disturbance rejection control; NMM: neural mass model.
For different abnormal He
s and Hi
s, stable regions have been discussed to check their influence. Let

Stable regions of ADRC for different abnormal excitatory and inhibitory gains. ADRC: active disturbance rejection control.
Frequency analysis
From the sections ‘Root loci analysis’ and ‘Stable region analysis’, we can see that the closed-loop system is less sensitive to He and Hi by ADRC. In this section, an explanation for the root loci will be presented in the section ‘Root loci analysis’ from the point of frequency domain will be presented. Then, the influence of the tunable ADRC parameters will also be discussed.
Figure 10(a) and (c) shows that the closed-loop neuromodulation system is insensitive to the varying parameters He
s. Gain margins and phase margins presented in Table 2 also confirm that enough stability margins can be guaranteed even if the excitatory synaptic gain He
varies. On the other hand, when inhibitory synaptic gain

Frequency response for a group of fixed parameters of ADRC. ADRC: active disturbance rejection control.
Frequency indexes for a group of fixed ADRC parameters.
ADRC: active disturbance rejection control.
However, is just a group of parameters effective? Here, we also check the closed-loop neuromodulation performance when choosing different groups of parameters. Figure 11 gives out the frequency responses corresponding to groups of controller bandwidths and observer bandwidths. Data presented in Table 3 also confirm that not only a group of fixed parameters given in the section ‘Root loci analysis’ is effective in suppressing epilepsy, but also groups of tunable ADRC parameters take effect. Especially, gain margins and phase margins are very close to each other. It signifies that ADRC is capable of taking effect in a wider range. It does help ADRC to apply in clinical practice.

Frequency response for groups of parameters of ADRC. ADRC: active disturbance rejection control.
Frequency indexes for groups of ADRC parameters.
ADRC: active disturbance rejection control.
Steady-state performance analysis
In this section, steady-state error of the closed-loop system is discussed. According to Figures 5 and 6, the tracking error transfer function from the input signal to tracking error can be derived, that is
where
When the set-point value
According to final value theorem, the steady-state error is
where
Since
Simulations
In this section, time domain responses are presented to illustrate that ADRC is able to suppress epileptic activities effectively. A group of hyper-excitatory parameters (i.e. He = 7, Hi = 22) is employed to simulate the epileptic behaviors. ADRC is designed to modulate the undesired activities. PI, a commonly used controller in control engineering, is also taken to make a comparison. The control action is activated at 8 s, and simulations last for 16 s. Two scenarios are considered. The first is the nominal case, and the second is the case when sinusoidal disturbance exists. Control structure is shown in Figure 12.

Seizure control structure.
In both scenarios, the tunable parameters of PI and ADRC are selected from Table 4.
Tunable parameters of PI and ADRC.
ADRC: active disturbance rejection control; PI: proportional-integral.
Scenario I
Here, a nominal case, that is, d = 0 in Figure 12, is considered. Both PI and ADRC are designed to control seizure. Figure 13 shows the comparisons of time domain responses and control signals between PI and ADRC.

Time responses and control inputs of PI and ADRC. PI: proportional-integral; ADRC: active disturbance rejection control.
Comparisons of the indexes are listed in Table 5.
Comparisons of indexes between PI and ADRC.a
ADRC: active disturbance rejection control; PI: proportional-integral; RMSE: root mean square error; MAE: mean absolute error; ITAE: integral of time-multiplied absolute-value of error.
a
Figure 13(a) shows that ADRC is able to control seizure with much shorter transient time. In other words, EEG signals modulated by ADRC are better than the ones modulated by PI. Moreover, Figure 13 signifies that, with less oscillations of the control signal, better modulation can be obtained by ADRC. It means that, in seizure control, less energy is needed by ADRC. RMSE, MAE, ITAE, and control energy listed in Table 5 also show that, with less energy, ADRC is able to get much smaller RMSE, MAE, and ITAE. It signifies that ADRC can control seizure more effectively.
Additionally, it is worth pointing out that the energy expenditure of the stimulation signal directly affects the life of the stimulation battery. Therefore, energy expenditure is a very important factor, and it is also considered in the comparison. Both Figure 13 and Table 5 tell us that, in a nominal scenario, ADRC is better than PI in suppressing epileptic activities.
Scenario II
With an attempt to check the robustness of the neuromodulation, sinusoidal signal, that is, d = 120sin(πt), is introduced at the 12th second to simulate the external disturbance. Both PI and ADRC are designed, and the tunable parameters are also selected from Table 4. Figure 14 shows the time domain responses and control signals of PI and ADRC, respectively.

Time responses and control inputs of PI and ADRC (sinusoidal disturbance exists). PI: proportional-integral; ADRC: active disturbance rejection control.
For the sinusoidal disturbance given in Figure 15, ADRC is capable of estimating and cancelling it effectively. Therefore, from Figure 14(a), we can see clearly that the time response of ADRC is much better than the one of PI. RMSE, MAE, and ITAE values of PI and ADRC given in Table 6 also confirm the fact mentioned above.

Sinusoidal disturbance.
Comparisons of the indexes between ADRC and PI (sinusoidal disturbance exists).
ADRC: active disturbance rejection control; PI: proportional-integral; RMSE: root mean square error; MAE: mean absolute error; ITAE: integral of time-multiplied absolute-value of error.
Additionally, from Figure 14(b) and the last column data shown in Table 6, it is obvious that the desired performance can also be achieved by ADRC with less control energy. It is the very target that we want in the neuromodulation of epilepsy.
So far, the results of both scenarios tell us that ADRC provides much better performance in seizure control. Therefore, it may be a more promising approach to suppress epileptic activities.
Discussion
An NMM provides a straightforward way to describe the activities of the populations of neurons. It cannot describe exactly how neurons interact, however, it is suitable to depict the role of exhibition and inhibition, and the balance/unbalance between them from the macroscopic view. 19 Based on such a computational model, which is also employed in most studies, 5,6,16 –23 an active disturbance rejection modulation approach is discussed. Here, in order to verify the robustness of the closed-loop system, excitatory gains are varied within ±5% of its nominal value. Both PI and ADRC are taken the same parameters listed in Table 4. Monte Carlo tests have been performed, and the results are presented in Figure 16.

Monte Carlo test results in the presence of He varies.
From Figure 16, one can see that, even if He varies in a certain range, RMSE values of ADRC are much less than the ones of PI.
Considering that an epileptic model is somewhat conservative, stable regions and frequency responses are also discussed when excitatory and inhibitory gains vary. System frequency domain performance is verified when different groups of parameters of ADRC are taken. Additionally, the time domain responses, in the presence and absence of external disturbance, are also discussed. Both the frequency responses and the time domain responses confirm that ADRC is robust enough to the nonlinearities, uncertainties, and external disturbances. In other words, there are enough stability margins for ADRC to deal with the conservative of a computational model. As a guidance for clinical applications, it may be a suitable closed-loop neuromodulation approach in the near future, if it is confirmed by animal experiments and clinical trials.
Conclusion
In this article, a closed-loop neuromodulation approach is discussed. Jansen’s NMM is taken to be the test platform, and the active disturbance rejection approach is employed to suppress the epileptic activities. Both time domain responses and frequency domain responses confirm that ADRC is capable of dealing with undesired factors, such as nonlinearities, uncertainties, and external disturbances, so as to guarantee the satisfied modulation. It is worth pointing out that ADRC is independent of a concrete model. Thus, it can be naturally developed to other neural models, animal models, or even clinical trials. From this point of view, the results presented in this article are able to provide a helpful reference, guidance, or even solution to the closed-loop neuromodulation in suppressing epilepsy. The future work should be online seizure control to standard animal epileptic models by ADRC or even applying ADRC to alleviate epilepsy patients as a long-term goal.
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the National Natural Science Foundation of China (grant numbers 61873005, 81771395), Support Project of High-level Teachers in Beijing Municipal Universities in the Period of 13th Five-year Plan (grant number CIT&TCD201704044), and Key Program of Beijing Municipal Education Commission (grant number KZ201810011012).
