Fitting of TC model according to key parameters affecting Parkinson's state based on improved particle swarm optimization algorithm

Biophysical models contain a large number of parameters, while the spiking characteristics of neurons are related to a few key parameters. For thalamic neurons, relay reliability is an important characteristic that affects Parkinson's state. This paper proposes a method to fit key parameters of the model based on the spiking characteristics of neurons, and improves the traditional particle swarm optimization algorithm. That is, a nonlinear concave function and a Logistic chaotic mapping are combined to adjust the inertia weight of particles to avoid the particle falling into a local optimum in the search process or appearing premature convergence. In this paper, three parameters that play an important role in Parkinson's state of the thalamic cell model are selected and fitted by the improved particle swarm optimization algorithm. Using the fitted parameters to reconstruct the neuron model can predict the spiking trajectories well, which verifies the effectiveness of the fitting method. By comparing the fitting results with other particle swarm optimization algorithms, it is shown that the proposed particle swarm optimization algorithm can better avoid local optima and converge to the optimal values quickly.

In recent years, neurological diseases caused by abnormal brain function have become a class of diseases that seriously affect human physical and mental health 1 . However, due to the limitations of human experiments, modelbased research is more common 2 . The researches on neurological diseases have given birth to the demand for neuron calculation and modeling 3 . Generally speaking, computational models can verify the proposed hypotheses, reflect individual differences through the internal structure and parameters of the model, and link the basic neural circuits with brain functions, thereby providing new tools for understanding neurological diseases. Since the specific characteristics of the model depend on the neuron, it is necessary to fit the computational model to the electrophysiological record 4 .
The conductance-based biophysical model can predict various spiking patterns with high accuracy, but the large network simulation of this type of model is computationally expensive 5 . Therefore, scientists have proposed various simple phenomenological spiking models as an alternative, such as linear integrate-and-fire model and various extensions to it-square integrate-and-fire model 6 , exponential integrate-and-fire model 7 , and Poisson generalized integrate-and-fire model 8 , etc. In recent years, most of the fitting and prediction based on neuron models have adopted simple spiking models rather than detailed conductance models describing cell biophysics.
Hertäg used standard f-I curves obtained by in vitro electrophysiologists to fit the adaptive integrate-and-fire neuron model and the adaptive exponential integrate-and-fire model 9 . Ghosh et al. implemented an adaptive exponential integrate-and-fire model on the NEURON simulation platform, estimated its parameters to describe the reference neuron, and verified the initiation dynamics of the model 4 . Pozzorini proposed a protocol combined with automatic patch clamp recording technology, which can efficiently transform a large amount of in vitro electrophysiological data into a spiking neuron model 10 . As large-scale network simulation has gradually become a hot research topic, more realistic models describing neuron characteristics are needed 11 . However, the complex structure and numerous parameters of the biophysical model reduce its computational efficiency and increase the computational difficulty. And serious overfitting is likely to occur. While in most biophysical models, only a few key ion channel parameters affect the intrinsic spiking characteristics of neurons 12 .
Parkinson's disease is a common degenerative disease of the nervous system, which is more common in the elderly. The clinical manifestations are static tremor, muscle rigidity, bradykinesia and postural instability 13 . The main pathological changes of Parkinson's disease are in the basal ganglia and thalamus neurons of the brain. When more than 70% of the substantia nigra dopaminergic neurons in the basal ganglia are destroyed, a significant decrease in striatal dopamine is caused and the normal function of the basal ganglia is hindered, inhibiting the relay function from thalamus neurons to cortex, thus leading to the symptoms of parkinson's disease 14,15 . The hindrance of the relay function of thalamic neurons is mainly due to the impact of the rebound cluster spiking behavior of the basal ganglia. Researchers have found that this rebound cluster spiking behavior is mainly related to the low threshold calcium ion current 15,16 . Therefore, the important parameters are fitted based on the significant spiking characteristics affecting Parkinson's state in this paper.
Since the signal transmission of the nervous system contains complex nonlinear processes 17 , it is difficult to obtain an effective system model using traditional linear fitting methods. Mitra et al. used the AugMAT model to predict the experimental discharge data, and adopted the simplex method and the gradient-based optimization method to fit the model parameters 18 . These algorithms are local optimization algorithms, that is, the search effect depends on the initial values of the parameters. In recent years, some evolutionary algorithms have been introduced into the parameter fitting process, such as genetic algorithm 19 , differential evolution algorithm 20 and particle swarm optimization algorithm 21 . These algorithms optimize the optimal solution through multi-dimensional parallel work. And the search process has higher efficiency and better results than traditional algorithms. The particle swarm optimization algorithm is a multi-dimensional and large-scale intelligent optimization algorithm proposed by Kenney and Eberhart 22 . Bangyal et al. gave a comprehensive survey of the various particle swarm optimization and DE initialization approaches based on the family of quasirandom sequences such as Sobol sequence, Halton sequence, and uniform random distribution 23 . They show that the low-discrepancy sequences-based initialization performed exceptionally better than a uniform random number. Pervaiz et al. presented a thorough overview of different particle swarm optimization initialization approaches which were dependent on semi-arbitrary successions systems and the best in class in the populace instatement was uncovered 24 . It was found that the proposed version of the algorithm is more effective in finding the global minimum of unimodal and multimodal test functions. Although there are many advantages of the particle swarm optimization algorithm, particles are prone to fall into local optima during the search process, or the case of premature convergence occurs 12 . Using a single method to improve particle swarm optimization algorithm may still fall into a local optimum. Therefore, a combination of nonlinear concave function and Logistic chaotic mapping is adopted to improve the particle swarm optimization algorithm in this paper.
In summary, we propose a method to fit the key parameters in the thalamic cell (TC) model and predict spiking trajectories of the Parkinson's state and the normal state. The particle swarm optimization algorithm is improved, and the fitness function is selected based on the spiking characteristics of the thalamic neurons in the Parkinson's state. After the key parameters that affect the Parkinson's state are fitted, the TC neuron model is reconstructed using the fitted parameters. The spiking trajectories of the reconstructed neuron model are compared with that of the original model to verify the effectiveness of the improved algorithm. Finally, four improved particle swarm optimization algorithms are compared with to prove the advantages of the proposed algorithm.

Results
System description. The TC model parameter fitting process based on the improved particle swarm optimization algorithm is shown in Fig. 1. Firstly, by inputting signals to the TC model with given parameters, the output firing data is obtained as standard data. Then give the same input signals to the reconstructed model with unknown key parameters, and compare its output data with the standard output data through the selected fitness function. Finally use the improved particle swarm optimization algorithm proposed in this paper to fit the unknown key parameters.
Fitness function of the TC model. The TC model in normal state can relay the input signal without losing the spiking information. While the TC model in Parkinson's state cannot relay the input signal well, that is, the output signal cannot match the input signal well, resulting in the spiking information loss. Also the reliability of the relay is greatly reduced. Therefore, relay reliability is an important spiking characteristic of the TC model. In order to quantify relay reliability of the TC model, we define: where error index is the error rate of the relay, including the spiking information of missed relay and the spiking information with relay error. Obviously, the value range of RI is between [0, 1]. 0 means that the output signal does not match the input signal at all, and 1 means that the output completely matches the input signal. Generally speaking, there is no obvious boundary between the TC model in normal state and in Parkinson's state. But it is generally believed that the TC model is in normal state if RI > 0.9 , and it is in Parkinson's state if RI < 0.5.
In addition to the reliability of the relay, we also select the total number of spikes, the average value of the spike peak and the average value of the subthreshold voltage together as the spiking characteristics of the TC www.nature.com/scientificreports/ model for parameter fitting to ensure the uniqueness of the spiking trajectory. Therefore, the fitness function for parameter fitting of TC model is selected as follows: where RI i , N i , PK i and ST i ( i = rec, ori) respectively represent the relay reliability, the number of spikes, the average spike peak value and the average subthreshold voltage value of the reconstructed TC model and the TC model with given parameters. w 1T , w 2T , w 3T and w 4T are the weight of each part respectively.
Fitting process of the TC model. Given the input signal I SM of the TC model as a square wave pulse signal: where the amplitude is A SM = 5 pA/μm 2 , the period is ρ SM = 25 ms, the duration is δ SM = 5 ms, and H is the Heaviside step function: In the TC model with given parameters, the key parameters of Parkinson's state are given as 3,120], which refers to the data in Ref. 26 and comes from biophysical experiments. In the TC model to be fitted, except for the three parameters I Gi→Th , g T and E T , the remaining parameters are set as constants. The improved particle swarm optimization algorithm are used to fit three important parameters. The fitting process of the three parameters is shown in Fig. 2. The red, pink and blue curves represent the parameters I Gi→Th , g T and E T to be fitted respectively. It can be seen from the figure that as the number of iterations increases, the three parameters all tend to the given value from the initial value. After a period of oscillation, it finally stops after the 31st iteration. The change of the fitness function with the number of iterations is shown in Fig. 3. It is shown in the figure that as the number of iterations increases, the fitness function gradually converges. When the iteration termination condition is met, the iteration stops. The values of the three parameters are the final iteration results, that is Effects of the fitting method. Define the parameter error of the TC model as: where e I Gi→Th , e g T and e E T are the errors between the optimal solution and the given value of the three parameters I Gi→Th , g T and E T . It can be seen from Fig. 4 that as the number of iterations increases, the parameter error to be fitted gradually decreases. And the final iteration error satisfies ln e T < −2.
(3) www.nature.com/scientificreports/ In order to explore the relationship between the fitness function of the optimal particle and the parameter error, we perform multiple simulation experiments on the TC model, and statistically analyze the results. The logarithmic relationship between the fitness function q T and the parameter error function e T of multiple simulation results is shown in Fig. 5. The results show that the relationship between the two is almost linear. Therefore, it can be ensured that when the fitness function is small enough, the error between the optimal solution and the true value is small enough.
To test the effectiveness of the fitted parameters, we use the fitted parameters to reconstruct the TC model and check whether the spiking trajectories of the reconstructed model and the original model are consistent. Figure 6a shows the spiking trajectories of the TC model with given parameters and the reconstructed model with fitted parameters in Parkinson's state. The black curve is the original model, the red curve is the reconstructed model, and the blue curve is the Poisson pulse input. It can be seen that the two curves overlap perfectly. Figure 6b shows the spiking trajectories of the given parameter model and the reconstructed model in normal state of the TC model. The black curve is the original model, the red curve is the reconstructed model, and the blue curve is  www.nature.com/scientificreports/ the impulse input. The curves can also overlap well. It shows that whether it is in Parkinson's state or in normal state, the method proposed in this paper can fit the key parameters of the TC model effectively.

Discussion
The improved particle swarm algorithm adopted in this paper is to adjust the inertia weight of particles through the combination of nonlinear concave function and Logistic chaotic mapping. In the particle swarm optimization algorithm, the inertia weight has a great influence on the convergence of the algorithm. In the classical particle swarm optimization algorithm, the inertia weight is usually chosen as a constant. While the usual improved algorithm generally selects the inertia weight as a certain function, which often causes the algorithm to fall into a local optimum and not to converge. In order to compare with the improved particle swarm optimization algorithm proposed in this paper, three key parameters of the TC model are fitted by selecting inertia weight as a constant, a single determined linear function and two sets of dynamic nonlinear concave functions, respectively. Figure 7 shows the fitting process of the three parameters I Gi→Th , g T and E T to be fitted in the TC model, in which: the inertia weight of Fig. 7a is a constant, and ω = 0.7 is selected. The inertia weight of Fig. 7b adopts a determined linear function:  www.nature.com/scientificreports/ A single dynamic nonlinear concave function is used to improve the inertia weight in Fig. 7c, which is shown in Eq. (13) on Methods part of this paper. Figure 7d uses another nonlinear concave function to improve the inertia weight: As the number of iterations increases, it can be seen from When the inertia weight is selected as a constant, a single determined linear function and two sets of dynamic nonlinear concave functions, the fitness functions change with the number of iterations as shown in Fig. 8. It can be seen from the figure that as the number of iterations increases, the fitness function decreases faster when the inertia weight is a constant, but the final value is large; the fitness function decreases slowly when the inertia weight is a linear function; and the fitness function fluctuates down when the inertia weight is a nonlinear concave function. The fitness function converges faster when the first concave function is selected; the fitness function converges slowly when the second concave function is selected. When the iteration is terminated, the final fitness functions of the four cases are much larger than that of the improved particle swarm optimization algorithm proposed in this paper. www.nature.com/scientificreports/ When the inertia weight is selected as a constant, a single determined linear function and two sets of dynamic nonlinear concave functions, the comparison chart of the total error of the three parameters I Gi→Th , g T and E T between the optimal solution and the given value is shown in Fig. 9. The natural logarithms of the total error in the four cases are [4.53, 3.79, 3.70, 4.11], which are much larger than that of the algorithm proposed in this paper.  www.nature.com/scientificreports/ It can be clearly seen from the comparisons that the method of combining dynamic nonlinear concave function and chaotic mapping adopted in this paper to construct inertial weight can avoid local optimum, and quickly converge to the optimal value, which can achieve good prediction characteristics.
The improved particle swarm optimization algorithm not only inherits the advantages of the particle swarm optimization algorithm, which has fast convergence speed and is easy to implement. It does not depend on the internal mechanism of the problem to be solved, which is very flexible in application. But it also avoids its disadvantages of tendency to fall into local optimum and appear premature convergence. The effectiveness of its search results only depends on the fitness function of the particles. So as long as the appropriate fitness function is selected, the data can be effectively fitted.
To conclude, a method of model parameter fitting based on neuron spiking characteristics is proposed in this paper. The traditional particle swarm optimization algorithm is improved, and the key parameters affecting neuron spiking characteristics are fitted using the improved particle swarm optimization algorithm. For thalamic neurons, three parameters that have a key impact on the Parkinson's state are fitted by using relay reliability as spiking characteristic index. Using the fitted parameters to reconstruct TC model, the spiking trajectories of the model can match well with the original model no matter in Parkinson's state or in normal state, which verifies the effectiveness of the fitting method. By comparing with the fitting effects of classical particle swarm optimization algorithm and several improved algorithms, it is shown that the proposed particle swarm optimization algorithm can better avoid local optimization and converge to the optimal value more quickly.

TC model and its key parameters. Signal relay related to Parkinson's disease is a very complex process.
In order to study the dynamic mechanism of Parkinson's disease, researchers have established a conductancebased neuron model that simulates the relay process of thalamic neurons 25 . For the convenience of research, the TC model can be simplified to a single-compartment model 26 . The dynamic equation is as follows: www.nature.com/scientificreports/ where V Th is the membrane potential of thalamic neurons. C m is the cell membrane capacitance. I L , I Na , I K and I T are leakage current, sodium ion current, potassium ion current and low-threshold calcium ion current, respectively. I Gi→Th is the inhibitory current flowing from the basal ganglia of the brain to the thalamic neurons, and I SM is the excitatory input current of the thalamic neurons. g L , g Na , g K , g T are the leakage conductance, the maximum conductance of sodium ion channel, the maximum conductance of potassium ion channel and the maximum conductance of low-threshold calcium ion channel, respectively. E L , E Na , E K , E T are the corresponding reversal potenitals. h Th and r Th are the gating variables of the ion channels. m ∞ (V Th ) and p ∞ (V Th ) are the steady state value of the sodium channel gating variable and low-threshold calcium channel gating variable. τ h (V Th ) and τ r (V Th ) are the time constants of the gating variables.
As the inhibitory current I Gi→Th from the basal ganglia to the thalamus neurons increases, the excitatory relay signal of the cortex is weakened, leading to Parkinson's state. And the low-threshold calcium current I T also has an important impact on the relay function of thalamic neurons. So we select I Gi→Th , g T and E T as the key parameters of the model in this paper. The values of other parameters are detailed in literature 26 . By fixing other parameter values, the TC model is abstracted into a model to be fitted with three variables. Improved particle swarm optimization algorithm. Particle swarm optimization algorithm is an intelligent optimization algorithm proposed by Kenney and Eberhart 22 . In order to avoid particles falling into local optimum during the search process, or appearing premature convergence, this paper adopts the method of combining nonlinear concave function and Logistic chaotic mapping to adjust the inertia weight of particles.
The fitness function of the particle is used as the criterion to judge the position of the particle. When the fitness function of a particle is smaller than the average fitness function of the particle population, that is, when the particle is evolving faster than other particles, we choose a certain nonlinear concave function to adjust the inertia weight of the particle to ensure that the convergence speed of the particle will not be too fast and lead to premature convergence. The concave function is selected because a large-scale global search is required at the beginning, so the function curve is relatively steep; later it needs to change from a large-scale global search to a small-scale fine search, so the function curve is relatively smooth. When the fitness function of the particle is larger than the average fitness function of the particle population, that is, when the evolution speed of the particle is slower than other particles, in order to avoid the particle from falling into the local optimum, we choose the non-deterministic Logistic chaotic mapping function 27 to adjust the inertia weight to force the particle to jump out of the local optimum and continue to search for the global optimum.
The steps of the improved particle swarm optimization algorithm proposed in this paper are as follows: I. Initialize the particle swarm; II. Calculate the fitness function q k id of each particle and the average fitness function q k avg of the particle swarm; III. If q k id < q k avg , meaning that the evolution speed of the particle is faster than other particles. We choose a nonlinear concave function to modify the inertia weight to make sure that the function curve is relatively www.nature.com/scientificreports/ steep at the beginning because large-scale global search is required, and relatively smooth at the end since it needs to change to a small-scale fine search: where k is the current iteration number and k max is the maximum iteration number. ω max and ω min are the maximum and minimum values of the inertia weight, respectively. If q k id > q k avg , meaning that the particle is evolving slower than other particles. We choose a nondeterministic Logistic chaotic mapping function to modify the inertia weight to make the particle to jump out of the local optimum: where α is the adjustment factor of inertia weight which value range is [0,1], α max and α min are the maximum and minimum values of α respectively. L map is the result of Logistic chaotic mapping. IV. Update the position and velocity of the particles according to the following formula: where i = 1, 2, ..., N , d = 1, 2, ..., D.x id , v id and p id respectively represent the position, velocity and historical optimal position of the d-dimensional search space of the i th particle, and p gd represents the historical optimal position of the entire particle population. c 1 and c 2 are constants, which are cognitive factor and social factor respectively, usually c 1 = c 2 = 2 . r 1 and r 2 are uniformly distributed random numbers in the range of 0 to 1. ω is the inertia weight. v max is the upper limit of the particle speed, when the particle speed |v id | > |v max | , then |v id | = |v max |. V. Determine whether the iteration termination condition is reached, that is, whether the maximum iteration number is reached or the fitness function of the particles is less than a certain threshold. If yes, the iteration ends; if not, repeat steps II to V until the iteration ends.