Nonlinear computational models of dynamical coding patterns in depression and normal rats: from electrophysiology to energy consumption

Major depressive disorder (MDD) is one of the most serious neuropsychiatric disorders. Exploring the pathogenesis and dynamical coding patterns of MDD can provide new targets for clinical drug treatment and new ideas for the research of other neuropsychiatric and neurodegenerative diseases. We selected the medium spiny neuron (MSN) of nucleus accumbens (NAc) as the research objective. NAc is located in the dopaminergic pathway, regulating rewards, emotions and other behaviors. Abnormalities in these behaviors are considered as the main clinical symptoms of MDD. We simulated the different spike patterns of MSNs in MDD group and control group by dynamical Hodgkin–Huxley model. The simulated results can match the electrophysiological experiments, which occurred due to following reasons: (1) The external stimulus current of MDD group was amplified by the local neural microcircuit; (2) the selective permeability to sodium was abnormally decreased; and (3) the dopamine D2 receptor signaling pathway was abnormal in the MDD group. Furthermore, we proposed a dynamical energy model, and the energy results demonstrated that the energy cost in MDD group was lower, which led to persistent depression in patients with MDD. Simultaneously, the negative-to-total energy ratio of MSN in MDD group was higher than that in control group, and the delay time of the power peak and the potential peak in MDD group was shorter than that in the control group. The results showed that the abnormal firing patterns were the direct cause of abnormal behaviors of MDD and indicated that subthreshold activities of MDD group were more intense.


Introduction
Major depressive disorder (MDD) is one of the most serious neuropsychiatric disorders accompanied by long-time depression, prolonged response time, lack of pleasure and other symptoms [1]. Patients with serious MDD have a strong suicidal tendency. According to the statistics, about 20% of the world's population are suffering (or have suffered) from depression [2]. At present, there is no drug or other treatment that can cure MDD; therefore, it is necessary to explore the pathogenesis of MDD.
In the past, the research on the pathogenesis of MDD mainly focused on the following four aspects: (1) exploring the pathogenesis of the abnormal behaviors in patients with MDD based on the psychological experiments [3,4]; (2) studying the abnormal functional connectivity of brain regions in patients with MDD by functional magnetic resonance imaging (fMRI) [5,6]; (3) exploring the abnormal firing patterns and synaptic transmission process in brain nervous system of patients with MDD based on the typical pharmacological, neurobiological and electrophysiological methods [7][8][9]; and (4) exploring more effective therapies by analyzing the clinical data (containing EEG data and behavior statistical data, etc.) of depressive patients during the medical treatments [10][11][12]. Meanwhile, since the incidents of MDD involve the pathological changes in multiple neural circuits and regions, researchers usually focus on some regions rather than on the whole brain. For example, the prefrontal cortex (PFC) is usually considered as the part of the brain responsible for executive function, the ventral tegmental area (VTA), nucleus accumbens (NAc) and amygdala are linked with emotional and reward-punishment decisions, and the hippocampus is linked with memory research.
However, the main locus and pathogenesis of MDD have not yet been fully elucidated by these studies. One of the most fundamental reasons is that there are few studies on the dynamical coding and decoding processes of MDD at this stage, and hence, the corresponding relation between the experimental data and behavioral performance has not yet been found. Besides, since the clinical evaluation of MDD mainly relies on behavioral and psychological measurement, and the measurement mainly relies on the psychological scale, the distinction between MDD and other psychological disorders (such as anxiety disorder and bipo-lar disorder) is not clear. Thus, with the evidence-based medicine (EBM) approach, doctors' clinical experience could be of crucial importance when treating patients.
The development of computational neuroscience has greatly advanced the research of MDD [13,14]. In 2016, Nature Neuroscience launched a special issue 'Focus on neural computation and theory' in order to summarize the progress of computational neuroscience, and it called for the establishment and improvement of theoretical neuroscience [15]. With these dynamical theoretical models, different types of the complex nervous brain activity can be explained [16][17][18]. In addition, neuropsychiatric disorders have received a special attention in Nature [19]. Although the experimental research has provided many useful data and there are already many methods of machine learning and data mining [20], the way to analyze the data and build a theoretical model remains a major difficulty. Therefore, Huys et al. put forward a new concept of computational psychiatry, which is to study the pathogenesis and mechanism of neuropsychiatric diseases and evaluate the degree of patients' disease and the treatment effect through quantitative computational models [19]. According to this approach, new drugs as well as treatments can be discovered. Moreover, these researches can provide a new perspective for the study of the pathogenesis of some psychiatric disorders and neurodegenerative diseases such as depression, Parkinson's disease and Alzheimer's disease [21][22][23][24].
According to the reductionism view of cognitive neuroscience, since the brain is a complex nonlinear dynamical network coupled by simple neurons, the underlying causes of MDD can only fall in one of the following two categories: First, one is a single neuron disease, such as abnormal pathological changes of the ion channels. The other is the neural network (or neural circuit) disease, such as abnormal changes of neuron types and number of neurons or synaptic transmission process (including the concentration and activity of neurotransmitter, receptor number, etc.). Therefore, the study of MDD with the methods of nonlinear dynamics can include everything from the abnormal changes of ion channels to neural networks in brain regions corresponding to some behaviors.
Since the abnormal behaviors such as rewards, punishments and emotions are considered as the main clinical manifestations of MDD, and these behaviors are mainly regulated by dopamine neurotransmitters, the NAc from the dopamine circuit was chosen as the main objective in this paper [25]. The medium spiny neurons (MSNs), which are the most important neurons in NAc and striatum, can affect the expression of dopamine receptors (D1-type and D2-type) and mediate the above listed behaviors. MSN is a special type of GABAergic inhibitory neuron [26]. It has the complex dendritic structure, which explains significant morphological differences among individuals [27][28][29]. A previous study has shown that the abnormal firing patterns of MSNs in NAc lead to abnormal emotional and reward behaviors [30]. In the above reported study, Liu et al. have successfully established a MDD rat model (Wistar Kyoto rats). Applying the electrophysiological approach, MSNs of brain slices of Wistar Kyoto rats and control rats were experimented before and after drug treatment (quinpirole and nomifensine). According to the results, MSN was abnormal in MDD group before and after the drugs administration, which was a strong prove that the pathological changes in neurons and their local microcircuits were the main cause of the dysfunction. In this paper, using nonlinear Hodgkin-Huxley model (H-H model) based on the above experimental data [30,31], we built a model that can describe the firing patterns of MSN in MDD group and control group of rat to discuss the abnormalities of its ion channels.
Furthermore, we analyzed the energy cost of depressive behavior for the following reasons: • Is depressive disorder related to different levels of brain energy cost? In recent years, the relation between neurological diseases and energy metabolism is becoming popular. In particular, abnormalities in the energy supply of glial cells to neurons are thought to be one of the causes of neurodegenerative diseases [32]. Although such significant results had not been obtained in studies of depression, it will be certain that depression must also be related to brain energy metabolism since Hu et.al found that Astroglial Kir4.1 in the lateral habenula drove neuronal bursts in depression [33]. Moreover, 'subthreshold depression' has been proposed as a subcategory of depressive disorders, which suggests that patients may have depressive symptoms but not achieve clinical diagnostic standards. Some researchers believe that it may be a prodromal period of depression, and that it may exist in the state of the intermission or incomplete recovery of MDD [34][35][36]. We would like to analyze this phenomenon with the methods of dynamical computational models and explore whether the different state of depression is related to energy. • The onset of depressive disorder involves different brain regions, circuit and even small-scale neural networks. Thus, it is difficult to find the main origin to focus on. The energy analysis of dynamical coding patterns can provide new ideas and new treatment methods [37,38]. For example, we can analyze the dynamical energy cost in different levels of nervous system, and the abnormal energy cost may imply pathological changes in this location [39]. In addition, the neural system also follows energy saving mechanisms [40], and we are curious if there are abnormalities in MDD. Such analyses can provide a theoretical basis for EEG and fMRI's results [41]. • Regarding the simplified calculation, since H-H model is a group of nonlinear differential equation models, if large-scale neural network model is established, the numerical calculation becomes very complicated, and the calculation speed is very slow. A series of recent studies have shown that the method of neural energy coding is more effective than the traditional method of neural coding in the study of local and global nerve activity in the brain. What's more, the energy method has many advantages such as no information loss [33,42,43,45]. Since energy is expressed as a scalar quantity, if the energy model is used to calculate the energy, the calculation can be greatly simplified.
This was a new perspective for the study of computational models of MDD. Previously, computational neuroscience focused on the study of fMRI and EEG, and there had been no computational model for MDD that matched those experimental data at the molecular and cellular level. Although the studies of fMRI and EEG were meaningful to explore the condition of brain regions of depression, the fundamental mechanism cannot be known. Therefore, the dynamical computational models given in this paper were of great significance: On the one hand, the models established were based on the previous electrophysiological experiments, according to which can provide guidance for future experiments; on the other hand, further analysis of the models can explore some new phenomena that had not been found in the experiments; for example, this paper can calculate the dynamical energy con-sumption of MSN, which was a conclusion difficult to verify in electrophysiological experiments, but can be easily calculated in the computational model.

In vivo recording
According to the experimental results reported in previous study [30], this paper studied the abnormal firing patterns of MDD from three perspectives: firing frequency, delay time before the action potential and stimulated current threshold. Using the H-H model above, the membrane potential of MDD group and control group MSN at different times was computed by adjusting the parameters of ion channel and stimulated current mode, thus matching the firing patterns that were measured in real neurons.

H-H model
The H-H model is used to describe the firing pattern of neurons, which can transform the complex nerve signals (mainly including the electrical signals generated by dendrites receiving AMPA/NMDA/GABA-A ligand-gated neurotransmitters and the intracellular expression of the second messenger mediated by Gprotein-coupled dopamine D1/D2 receptors) into the properties of neuron ion channels and local microcircuits. Because the individual differences in morphological structure of MSN are great, the calculation for MSNs would be difficult, and thus, it is necessary to properly simplify the model [46][47][48]. In the current study, MSN is considered as a soma-only neuron, since those above-mentioned nerve signal processing is almost entirely carried out in the nucleus of neurons. According to Rall model [49], electrical signals on dendrites are transmitted to soma through the mode of cable transmission. The simplified H-H model is as follows: where C m is the capacitance of neural cell membranes, which was taken as 1 µF/cm 2 in this paper. V m is the membrane potential in mV (V m = V in − V out ). I ion is the ion current, and I ext is the stimulated current. The unit of these currents is µA/(cm) 2 . The circuit diagram is shown in Fig. 1. Fig. 1 The circuit diagram of H-H model. Note that although I M is a special kind of potassium ion channel, its corresponding Nernst potential is still E k . We can get the dynamical membrane potentials by solving this nonlinear circuit model using numerical calculation methods Compared with the original H-H model, a noninactivating K + current (M-type K + ) was added in this paper. This current exists in the soma and dendrites of neurons and is responsible for spike frequency adaptation: When the external stimulated current is initially applied, the interval of spikes is relatively short; afterward, it gradually prolongs and reaches the stability (Supplementary Fig. S1) [31]. Therefore, there are four types of ion current in the model investigated in this paper: Na + current I Na , K + current I K , M-type K + current I M and leak current I L : Ion currents can be described by ion conductance, which is determined by the ion channels of neuron membranes. The general expressions of different kinds of ion currents are as follows: where i denotes the different types of ion currents. g i is the maximum conductivity in µS/(cm) 2 . m, n and h denote different types of voltage-gated variables, which are used to describe the gated characteristics of ion channels. M i , N i and H i are exponential coefficients. E i is the Nernst potential of the ion channels.
In this model, because the Cl − channel property was generally stable and the M-type K + current was only used to adjust the interval of spikes, g L and g M were taken as constants. g K , g Na and I ext are adjustable parameters. The detailed calculation method and parameter setting of the model are shown in Supplementary Table S1.

Energy model
Because the ion current and the membrane potential at every period in H-H model above are available, the power of neurons can be calculated, and the energy cost can be obtained by integrating the power with time. Therefore, this section further explores the energy cost characteristics of rats with depression by studying the energy coding of a single MSN. On the basis of H-H model reported above, the power was calculated according to the following formula [50]: where the power of Na + was negative. This is because the directions of Na + current I Na and Nernst voltage E Na are opposite to that of K + and Cl − . Thus, the calculation of Na + power needed to be negative. Generally, in the calculation of ion power, the selection of symbols is based on the direction of K + current: If the direction of ion current is the same as that of K + current, the ion power is positive; conversely, if the direction of ion current is opposite to that of K + current, the ion power is negative. Relevant theoretical deduction can be seen in some research of Zhu et al. [50]. By integrating the above reported power according to time, the energy cost of neurons in any time interval [a, b] was obtained through the following formula: 3 Results

Firing frequency
Firstly, we explored the reasons for the abnormal increase in firing frequency in MDD rats. Liu et al. [30] have shown that the firing frequency of MSN in MDD group rats was not less compared to control group rats when stimulated with the same external clamp current. The result of this electrophysiological experiment is consistent with the clinical symptoms of MDD. According to neuropsychiatric and neuropsychological literature, there are abnormal enhancements of 'bottom-up processing' and weakening of 'top-down control' in patients with MDD to negative affective stimuli. 'Bottom-up processing' refers to the processes in which external stimuli are received by sensory organs (such as the retina, skin and olfactory mucosa) and transmitted through a series of nerves to the advanced central nervous system (such as the cerebral cortex and thalamus). In this process, stimulus information is transformed into electrical and chemical information, which are transmitted and processed through synapses and neurotransmitters. 'Top-down control' refers to the process in which the cortex receives and analyzes the above-mentioned information and generates this information (such as vision, olfaction and emotion). Patients with MDD usually magnify negative information, especially negative affective stimuli, and it is more difficult to regulate their negative mood. In the vicious circle of 'processing' and 'control,' the condition gradually deteriorates [51]. In fact, these phenomena also exist at the cellular level. If we consider the neural microcircuit composed of a single MSN cell and other neurons connecting the cell where synaptic connections exist, the neurons are regulated by neurotransmitters, which directly or indirectly alter the properties of ion channels on the cell membranes [52]: (1) Ligand-gated ion channel receptors (such as AMPA and NMDA receptors, GABA A receptors and N-ACh receptors) directly regulate the ion channel properties of postsynaptic membranes and alter firing patterns; (2) G-protein coupled receptors (such as GABA B receptors and dopamine D1 and D2 receptors) indirectly alter the firing patterns by regulating the intracellular protein synthesis through signal transduction [53]. That is to say that the local microcircuit amplifies the external stimulus current due to the pathological changes in the network structure of MDD, which leads to the increase in the firing frequency of a single MSN. This occurs due to neuroplasticity. By studying the abnormal regulation of these microcircuits, we can explain the abnormal enhancement of 'bottom-up processing' in MDD. In our model, we established the H-H stimulated current models, respectively: • Regardless of the magnification and weakening of MSN stimulated current in the control group, the stimulated current model was as follows: where I ext is the external stimulated current in H-H model of control group, t 0 is the moment of stimuli, and I Clamp is the stimulated current given by current clamp.
• As for the rats with MDD, the microcircuit of MSN amplifies the stimulated current. In this paper, we called this current the processed current in H-H model, and we used it as the external stimuli. The processed current model was as follows: when t ∈ [t 0 , t 1 ) k · I Clamp , when t > t 1 (8) where I ext is the processed current in H-H model of MDD group, and I Clamp is the stimulated current given by current clamp. t 0 is the moment of stimuli, and t 1 was the moment when the microcircuit amplified the stimulated current. Amplification of the current was simulated by first-order function model, where k is the magnification (k > 1). According to these H-H stimulated current models, the time-current image is shown in Fig. 2.
The neuroelectrophysiological experiment of MSN demonstrated that: (1) The firing frequency of a single MSN corresponding to the stronger current was also higher; (2) when the intensities of external stimuli were the same, the firing frequency of MDD group was higher than that of the control group [30]. This phenomenon was simulated by H-H model. For example, Fig. 3 shows that the firing count of MSN from control group and MDD group was stimulated by different current. Detailed parameter settings for Fig. 2 are shown in Supplementary Tables S3 and S4: The reasons for the g Na and g K are detailed in the following passage, and the stimulated current was adjusted according to the previously reported experiment [30]. The simulation results were in great agreement with the neuroelectrophysiological experiment.  2) The results of ERP have demonstrated that the latency of N2 and P3 components in P300 of depressive patients could be longer than that of control group when stress stimulation, reward events and emotional information are produced, or even longer when the condition is more serious [54]. The behavioral explanation for this phenomenon is that depressive patients are duller and more insensitive in the face of external stimuli, and the delay time is positively correlated with the severity of depressive disorder. Previous neuroelectrophysiological experiments [30] have shown that the delay time of MDD group was longer showed that the firing count of MDD group was higher than the control group, which was similar to the experimental results, and detailed parameter settings are shown in Supplementary Table S4 than that of control group. The following modeling was based on that experiment.
According to the Nernst equation and H-H model, the membranes of neurons in rest state are selectively permeable to K + ; however, the selective permeability to Na + is poor, which leads to positive charges outside the cell membrane and negative charges inside. Therefore, the potential outside membrane is higher than the inside one. When the neuron is stimulated by external stimuli to reach the threshold, the selective permeability of Na + rapidly increases. Under the high concentration gradient and electric field, a large quantity of Na + flows into neuron. Consequently, the intramembrane potential of neuron increases rapidly and the membrane potential decreases, resulting in action potential. In this process, the selective permeability of Na + is the most important factor depending on the speed of Na + flowing into neurons. If the Na + permeability is high, the delay time will be short.
In H-H model, g Na was the parameter corresponding to Na + . By adjusting g Na , the speed of Na + flowing into neurons can be simulated. It was demonstrated that when g Na increases in the theoretical range, the delay time shortens. Hence, we simulated the characteristic of 'prolonged response time' of MDD by reducing g Na , and thus increasing the delay time. The simulation results are shown in Fig. 4.

Stimulated current threshold
In addition to the behavioral symptoms reported above, patients with MDD often have a long-term low mental state, such as losing concentration or interest in anything. Hence, they may lose confidence to have the ability to overcome difficulties. Clinical analysis has suggested that one of the main causes of MDD is the abnormality of the dopaminergic pathway. Because dopamine can directly affect mood and regulate emotion, the weak release of dopamine associated with the abnormal pathway could reduce the patient's sense of pleasure and interest.
In our model, these behavioral symptoms were described by stimulated current thresholds, which were defined as the smallest external stimulated current causing the neurons (except spontaneous firing neurons) to generate an action potential.
Since MSN studied in this paper was located on the dopamine mesolimbic system, we postulated that it was the abnormality of the dopaminergic pathway in patients with MDD that changed the property of ion channel of MSN, thus affecting the stimulated current threshold. This hypothesizes was confirmed in previous study [30]: (1) When the dopamine D2 receptor agonist drug quinpirole was given externally, the stimulated current threshold of MSN in both MDD group and control group increased, but the threshold of MDD group Fig. 4 The relation between Na + channel and delay time. A The change of firing delay time after changing g Na . B Delay time corresponding to different g Na , the results showed that when the g Na rose, the delay time will decrease, and detailed parameter settings are shown in Supplementary  Table S5 was significantly lower than that of control group; (2) when the dopamine reuptake inhibitor drug nomifensine (the pharmacological effect of nomifensine is to increase the concentration of dopamine by inhibiting the dopamine reuptake activity) was given to the MDD group, the firing patterns of MSN in MDD group returned to normal. This above experiment demonstrated that MSNs in dopamine D2 receptor signaling pathway could be abnormal in the excitatory regulation.
Some researchers have shown that the activation of dopamine D2 receptors indirectly regulates the activity of K + channel on the membrane of MSN by intracellular metabolism [55,56]. The relation between the K + channel parameter g K and the stimulated current threshold in our model was simulated through the control variable, and the result in Fig. 5 showed that when the g K was increased in the theoretical range, the stimulus current threshold simultaneously increased.
This result has a physiological basis. When quinpirole is given to MSN at rest state, the activity of K + channel increases, resulting in the faster speed of Na + flowing into neurons. Therefore, if action potential is generated, the required stimulus current is higher, i.e., the stimulated current threshold is increased.
Hence, the results of neuroelectrophysiological experiment could be simulated by the following methods: (1) For the control group, if the dopamine D2 receptor signaling pathway was normal, after giving quinpirole, g K corresponding to the activity of K + channel rapidly increased, and the stimulated current threshold increased; (2) for the MDD group, if the dopamine D2 receptor signaling pathway was abnormal, after giving quinpirole, although the activity of K + channel was increased, it was lower than that of the control group.
In our model, we let the maximum conductance of K + in control group before quinpirole administration is g K (control, before), and that after quinpirole administration is g K (control, after). And the maximum conductance of K + in MDD group after quinpirole administration is g K (MDD, after). Since the MSN of control group before quinpirole administration is normal, g K (control, before) is standard value, and thus, we let g K (control, before) = 100; then, the following conclusion is obvious. 100 = g K (control, before) < g K (MDD, after) < g K (control, after)

Simulation
According to the methods above, the parameters for the model were set in relation to the following four aspects: (1) When the external stimulated current increased, the number of spike increased, and the firing frequency of MSN increased; (2) depression amplified the stimulated current through microcircuit, resulting in abnormal increase in the firing frequency of MSN; (3) the activity of Na + channel was abnormal in MDD group, and the selective permeability of Na + in cell membrane was weakened, thus leading to a longer delay time; and (4) the dopamine D2 receptor signaling pathway was abnormal in MDD group, when dopamine inhibitor was given to MSN, and the threshold of stimulated current was different from that of control group. The firing patterns of MSN from MDD group and control group were, respectively, simulated, and the symptoms were matched with the experimental results. Thus, the model proposed in this section described the physiological Fig. 5 The relation between K + channel and the stimulated current threshold (g Na = 200). A Firing when g K = 100 and the stimulated current is 2.60. B Not firing when g K = 200 and the stimulated current is 2.60; C firing when g K = 200 and the stimulated current is 2.62. The results demonstrated that when g K increased, the stimulated current threshold will also increase. Detailed parameter settings are shown in Supplementary Table  S6 mechanism of MDD in a single neuron, indicating that the model has a physiological and behavioral basis. The simulation results are shown in Fig. 6. Furthermore, the power and energy of a single MSN in MDD group and control group without quinpirole were calculated, and the simulation results are shown in Fig. 7.

Result analysis
Based on previously reported neuroelectrophysiological experimental data [30], H-H model and energy model were, respectively, established in the current study, and these two dynamical models were used to describe the different firing patterns of a single MSN in MDD group and control group. The following discussion demonstrates that the two coding models can match the experimental results. Besides, the results of energy model can explain why depressive patients suffer from long-time depression.
Firstly, we considered the correspondence of membrane potential, power and energy in a single MSN, such as the sixth group (I Clamp = 5.22) of MDD group without quinpirole and with enlarged membrane potential, power and energy figures. The data are shown in Fig. 8.
After the enlargement of the above-mentioned membrane potential simulation results, the action potential was observed, which matched the physiological experiments. In other words, the simulation results could describe the five states of the spike process: (1) After being stimulated by external current, the selective permeability to Na + on the membrane begun to increase, Na + flew slowly into the cell, and the membrane potential gradually approached the firing threshold from the subliminal level; (2) the selective permeability to Na + continuously increased, and the membrane potential rapidly increased over the threshold, resulting in an action potential; (3) the selective permeability to Na + begun to decrease, conversely, the selective permeability to K + increased; K + begun to flow slowly into the cell, and the membrane potential started to decline after reaching its peak; (4) the selective permeability to K + continuously increased, and Na + flew out of the cell rapidly. As a result, the membrane potential declined lower than that of the rest state, which is also known as 'hyperpolarization'; (5) the selective permeability to K + decreased, and the membrane potential returned to a resting state; the next spike begun after certain period of time (refractory period).
From the time-power simulation results, we drew the following conclusions: (1) In the process of firing, MSN needs energy, and the longer it is stimulated externally, the more energy it consumes; (2) there is a negative power before the firing of MSN. In addition, we also found that when MSN starts to fire, it has a negative power, but when the membrane potential reaches its peak, the power changes to positive, and the ampli-  Table S7 tude of energy cost rapidly increases, indicating that the peak of neuron energy cost lags behind the peak of its membrane potential. Wang et al. have named this phenomenon 'negative energy' [33,42,43,45,50]. This phenomenon was explained with the fact that at the beginning of the spike, neurons may absorb a small amount of energy and oxygen from the microarteries and microvessels connected with them, thus showing a downtrend in the time-energy figure; also, when neurons go from repolarization to hyperpolarization, a large deal of biological energy needs to be consumed, so the time-energy figure has upward tendency. This phenomenon also has physiological significance, which can explain why the oxygen consumption is different from the blood flow volume: When the neurons are firing, the blood flow volume increases by about 31%, while the oxygen consumption increases by only 6% [39,57].
H-H model calculation and energy model calculation are used to analyze the mechanism of MDD from different perspectives. The energy model is calculated by H-H model, so the energy results are essentially the same to the membrane potential results. Since energy is scalar, the energy cost of a single MSN can describe that of neural networks partially, which is also the view of reductionism in cognitive neuroscience. As shown in Fig. 9, when MSN is stimulated, the action potential matches the energy cost, and the number of spikes is the same to that of energy cost activity. Additionally, by analyzing the energy cost of MSN in MDD group and Fig. 7 The power and energy simulation of a single MSN in MDD group and control group without quinpirole. A The power and energy of MDD group with different stimuli current (A1 power results, A2 energy results). B The power and energy of control group with different stimuli current (B1 power results, B2 energy results). In the energy simulation image, the vertical coordinate means the total biological energy cost before the corresponding time control group, the following three abnormalities were found: (1) The energy cost of MSN in MDD group was lower than that in control group; (2) the negative-tototal energy ratio of MSN in MDD group was higher than that in control group; and (3) the delay time of the power peak and the potential peak in MDD group were shorter than that in control group. Considering that the negative energy is supplied by local microarteries and microvessels, we can conclude that in rest state, the MSN of MDD group has higher energy levels than control group, thus suggesting that patients with MDD have more intense subthreshold activity in their brain [58]. Also, when the MSN fires, the energy cost of MDD group is lower than that of control group, which indicates that the dopamine production ability decreases in MDD group. Furthermore, the abnormal increase of the firing frequency in MDD group also indicates that the depression state can persist for long periods of time. Hence, the energy analysis is based on behavioral theory, and it can reflect the energy cost in depressive patients.

Prediction
According to the view of reductionism in cognitive neuroscience, the pathogenesis of MDD must come from the abnormal pathological changes of single neuron (such as the abnormal properties of ion channel, the proliferation and apoptosis), or the abnormal neural network (such as abnormal neurotransmitter secretion and receptors). The H-H model mentioned above  The ratio in MDD group is higher. C The comparison of lag time between power peak and the potential peak: The lag time in MDD group is shorter simulated the firing patterns of MDD, and the energy model provided new indicators for the clinical diagnosis of depressive disorder, which can predict subthreshold depression and early stage of major depressive disorder. Therefore, the clinical symptoms can be predicted with the above models.

Integration
At present, the researches of depression mainly focused on different levels: molecular level, such as genetic and fluorescent protein experiments; cellular level, such as electrophysiological experiments; system level, such as EEG and fMRI experiments; individual level, such as clinical experiments; and population level, such as analyzing the big data of depressive patients. These experiment data were massive and messy, and it was necessary to integrate them. But the most difficult part of integrating these data was that these experiments were performed at different levels; for example, EEG and fMRI reflected the level of the whole brain, while genetic and electrophysiological experiments reflected the microlevel; due to the current technology, it was impossible to monitor the brain activities at microlevel and the whole brain level simultaneously. The models presented in this paper can solve this problem to a certain extent: Although the fundamental reason for the change of the properties of ion channels was the change of genes, the effects on the cell membrane were the firing properties, and thus, the H-H model unified the molecular level and the cellular level; the energy model was based on H-H model, and the energy was scalar, which can be summed to reflect the energy consumption of brain regions, and this was consistent with the core idea of fMRI, so the energy model integrated the cell level and the system level. Therefore, the models can reflect different research levels of patients with depression at the same time, which were of great significance to explore the pathogenesis of depression.

Conclusion
Based on the electrophysiological experiment, a nonlinear computational model of membrane potential for MDD was established for the first time. And the dynamical power and energy cost model in the firing process of MSN were calculated. The results demonstrated that there were three differences in energy cost of MSN between MDD group and control group: (1) The energy cost of MSN in MDD group was lower than that in the control group; (2) the negative-to-total energy ratio of MSN in MDD group was higher than that in the control group; and (3) the delay time of the power peak and the potential peak in MDD group were shorter than that in the control group. These were the causes identified in the depressive patients with persistent low mood and intense subthreshold activity. The dynamical action potential model and energy model proposed in this paper can predict somewhat the clinical symptoms in major depressive disorder, which are the basis for building largescale nonlinear complex networks.
Since we were more interested in analyzing the encoding patterns of depression than in analyzing complex neuronal structures, we simplified the MSN model. However, although the MSN model we established in this paper was a little bit simple, it was the first model which considered the modulation relations between the membrane potentials and ion channels in MDD. To some extent, because this nonlinear dynamical model we used can explain the dynamically kinetic changes of the different types of ion channels, it could provide some quantitative basis for the explanations of the physiological mechanisms in the electric activities of MDD. This will be a very new and useful research direction of depression, because decoding the electric activities and firing patterns could be the most important step in exploring the mechanism of MDD.
In future research, we will use some other mathematical and dynamical methods to consider the small-scale microcircuit composed of MSN, analyze the dynamic characteristics of receptors coupled with ion channels and metabotropic receptors under the influence of neurotransmitters, explore the coding mode of the network and gradually establish a large-scale model, so as to explain the pathological causes of MDD.

Supplementary Information
The supplementary materials showed detailed parameter settings in the simulation.