Temperature-dependent bursting pattern analysis by modified Plant model

Many electrophysiological properties of neuron including firing rates and rhythmical oscillation change in response to a temperature variation, but the mechanism underlying these correlations remains unverified. In this study, we analyzed various action potential (AP) parameters of bursting pacemaker neurons in the abdominal ganglion of Aplysia juliana to examine whether or not bursting patterns are altered in response to temperature change. Here we found that the inter-burst interval, burst duration, and number of spike during burst decreased as temperature increased. On the other hand, the numbers of bursts per minute and numbers of spikes per minute increased and then decreased, but interspike interval during burst firstly decreased and then increased. We also tested the reproducibility of temperature-dependent changes in bursting patterns and AP parameters. Finally we performed computational simulations of these phenomena by using a modified Plant model composed of equations with temperature-dependent scaling factors to mathematically clarify the temperature-dependent changes of bursting patterns in burst-firing neurons. Taken together, we found that the modified Plant model could trace the ionic mechanism underlying the temperature-dependent change in bursting pattern from experiments with bursting pacemaker neurons in the abdominal ganglia of Aplysia juliana.


Introduction
To date, researchers have investigated the effect of temperature on the electrical activity and firing patterns in neurons from many animals, including Aplysia juliana, crabs, crayfish, frogs, lobsters, locusts, snails, and squids [1][2][3][4][5][6][7][8][9]. Especially the functions and properties of bursting pacemaker neuron R15 in the abdominal ganglion of Aplysia have been extensively studied [10][11][12] and mathematical simulation of bursting activity has been successfully performed [13][14][15][16][17][18][19]. However, few studies on the temperature dependence of action potential (AP) parameters in the R15 bursting pacemaker neuron have been reported [20,21]. They analyzed typical changes of AP parameters in burst-firing neurons by investigating the effect of heat on R15 bursting pacemaker neuron activity as the temperature increased and reported temperaturedependent changes in inter-burst interval, burst duration, number of spike per burst, intra-burst spike broadening and spike height [21]. However, the reproducible properties and mechanism of temperature-dependent changes of AP parameters remain unknown yet.
Changes in temperature can produce numerous effects on the neural tissue of most organisms. Indeed it has been reported that hyperthermic temperature may induce depolarization and spontaneous firing of pyramidal neurons leading to enhanced excitability of hippocampus [22,23]. On the other hand, a body of evidence indicating a therapeutic effect of hypothermia has accumulated in several conditions. Orlowski and colleagues published that in refractory status epilepticus unresponsive to conventional treatment, systemic hypothermia (30-31°C) was an effective therapeutic method leading to burst suppression on electroencephalography (EEG) [24]. Despite a lot of works aiming at demonstrating the neural effect of temperature, the mechanism of temperaturedependent electrophysiological change containing alteration of bursting and firing pattern in neuron remains unclear. Here, we examined the effects of temperature changes on the neuronal activity and bursting patterns during several consecutive heating − cooling cycles by using bursting pacemaker neurons which are a proper specimen with capability of long-lasting recording for mathematical modeling. Next, we sought to identify the mechanism underlying temperature-dependent bursting patterns of these neurons by analyzing and comparing the experimental data to computational simulation data calculated by modified Plant equations with temperaturedependent scaling factors, ρ(T) and ϕ(T).

Animals and dissection
Animals (A. juliana) were collected locally in Seogwipo, South Korea and all experimental procedures were approved by the Jeju National University Animal Care and Use Committee. These animals were dissected to observe the temperature-dependency of bursting patterns generated by pacemaker neurons R15 in the abdominal ganglia of A. juliana as described in previously published paper [1]. Briefly, the animals were anesthetized with an injection of 0.38 M MgCl 2 amounting to half of each specimen's weight before the abdominal ganglia were removed. Each abdominal ganglion was incubated at 34°C for 40 minute in a solution containing equal volumes of artificial sea water (ASW; in mM: 460 NaCl, 10 KCl, 11 CaCl 2 , 55 MgCl 2 , and 10 HEPES; pH 7.6) and isotonic Leibovitz's L-15 media (Cat. No. L-4386; Sigma) containing 1% protease (type IX; Sigma) (hereafter, ASW:L-15 (1:1)). They were then washed several times with ASW and were placed in a low-temperature incubator (VS-1203PIN; Hanback Co., Daejeon, Korea) at 18°C. At last they were removed onto Petri dishes (50x9 mm), pinned down on Sylgard plates (Dow Corning, USA) filled with ASW:L-15(1:1), and desheathed.

Data acquisition and analysis
A PT100 platinum resistance temperature sensor connected to a digital thermometer (TRM-006, Toho, Japan) was placed near the abdominal ganglia of Aplysia soaked in ASW:L-15 (1:1, v/v) to measure the medium temperature, which could be increased or decreased by activating the temperature controller system (HMN 3940, Acetec Co., Korea). We used a low flow rate peristaltic pump (BJ100-2 J; Baoding Longer Precision Pump Co., China) to maintain a viable solution with ASW:L-15 (1:1) media. The flow rate and the speed were 0.14 mL/min and 1.0 rpm, respectively. Intracellular recordings of bursting pacemaker neurons were performed to measure the membrane potential. The glass intracellular electrodes were filled with 3 M KCl and the membrane potentials recorded were simultaneously saved using a DAQ card (NI PCI-6221, National Instruments) and the Labview program (National Instruments). The electrical signals of the APs were identified using a digital oscilloscope (54622A; Agilent, Colorado Springs, Colo., USA). These data were recorded at a rate of 3 kHz to reproduce the APs on the computer because the bursting cells generated APs at a rate of 0-4 Hz; each component of this data is composed of 180,000 pairs of temperature and membrane potential values with unique file names, and the dataset was composed of 1,713 files as shown in Table 1. With this selected dataset, the average values of each of the AP parameters of the burst-firing neuron could be calculated using Origin 6.0 (Microcal Software, Inc.) and a computer program we designed. For each experiment, the bursting neuron was held at room temperature for one hour as shown in Table 2 and then changed with continuous ramp from about 16°C to 30°C. Then, the data recorded by increments (or decrements) of 2°C from about 16°C to 30°C (16-18-20-22-24-26-

Definition of AP parameters in burst-firing neuron
The graph in Figure 1A represents the typical bursting patterns of electrical signals shown in Aplysia bursting pacemaker neuron. The intraburst interspike interval, interburst interval, and burst duration are represented by symbols (1), (2), and (3), respectively. The membrane potentials at the positive peak, V pp (mV), and negative peak, V np (mV), are defined as the values of membrane potentials at points T and B shown in Figure 1B, respectively. The definition of the first half of the rising phase of AP, Δt r1 (ms), and the following AP parameters are very similar to those defined in our previous study [1].
The last half of the rising phase of AP, Δt r2 (ms) and the first half of the falling phase of AP, Δt f1 (ms) are time intervals shown in Figure 1B. The last half of the falling phase Δt f2 (ms) is defined as the values of the time interval between two points H2 and B shown in Figure 1B. The AP half-width duration, Δt AP, 1/2 (ms), is defined as Δt r2 + Δt f1 (ms). The interspike interval, ISI (ms), is defined as Δt r1 + Δt AP, 1/2 + Δt f2 . The spontaneous firing frequency, simply referred to as Frequency (s −1 ), is defined as ISI −1 (s −1 ). For convenience, the average angles θ 1 and θ 2 are defined to compare the slope of the second half of the rising phase of the AP with that of the first half of the falling phase. The angle θ 1 is defined as the inverse tangent of the ratio of half A AP to Δt r2 , and the angle θ 2 is also defined as the inverse tangent of the ratio of half A AP to Δt f1 : i . Although angles are usually defined  as dimensionless quantities, these angles are defined differently for convenience. The AP amplitude, A AP (mV), is defined as the absolute value of difference between the two values of membrane potentials at points T and B, respectively. More detailed analytical techniques regarding the definition and computational analysis of AP parameters and others were described previously [1,25]. Data were processed using the scientific data analysis tools, C ++ , Origin 6.0, and Mathematica 5.1, and the resulting data were presented as mean ± SE (standard error).

Structure of the selected dataset for AP parameter analysis
Eight experiments (A − H) were conducted using A. juliana specimens, and Table 1 shows the animal weights, selected data for analysis, total number of spike, and total recording time.
Animal weights were between 130 and 310 g, and the average weight was 226 g. The average total number of spike and average total recording time of each experiment were 15,039 spikes and 607 min, respectively. It was necessary to select the middle portion of each recorded dataset for analysis. The average values of time intervals for analysis, numbers of cycles of temperature change, and number of spike selected for analysis of these selected data were 214 min, 2 cycles, and 5,433 spikes, respectively. The dataset gathered for analysis was composed of these selected data from each experiment. The total time interval for analysis and the total number of spike of this dataset were 1,713 minutes and 43,468 spikes, respectively. The numbers written in parentheses in the middle column of Table 1 represent the total number of cycles of temperature change for each experiment.
In each experiment, the temperature was maintained at room temperature for the first one hour. As shown in Table 2, average values of temperature and four AP parameters of burst-firing neurons were calculated using the selected dataset saved by intracellular recordings during maintaining at room temperature for 30 minutes; all data were selected after 5 to 10 minutes from temperature change onset. Because the standard errors of the values of these parameters were low, the bursting patterns of these experiments might be considered as regular while the temperatures were held constant at room temperature.

Temperature dependence of AP parameters in burst-firing neurons
We selected burst trains in a continuous time series from the intracellular recording data of experiment A during the fourth falling phase (from 440 to 487 min) ( Figure 2) out of the eight falling phases in temperature recorded during eight consecutive heating − cooling cycles: it was composed of 48 panels. Temperature and membrane potentials are represented as the upper and lower traces in Figure 2, respectively. These figures demonstrated that interburst interval, burst duration, and number of spike during burst decreased as temperature increased. Next, we extracted six burst signals (6 panels) among many burst signals to examine the temporal change of instantaneous ISI during burst at specific temperature; 30.0°C (at 440 min), 26.9°C (at 454 min), 23°C (at 464 min), 21.3°C (at 468 min), 19.1°C (at 474 min), and 16.4°C (487 min). As shown in Figure 3, ISI during burst versus time displayed parabolic pattern, suggesting that temperature-dependent bursting patterns are parabolic bursts characterized by lower spike frequency at the beginning and end of the burst. The analyzing dataset was composed of 1,713 files corresponding to the sum of each selected time interval for analysis as shown in the 4th column of the Table 1.
With the selected data shown in Table 1, the temperature-dependent properties of AP parameters were analyzed by using techniques similar to those described in Hyun et al. [1]. As shown in Figure 4A and C, all values of ISI and Δt r1 decreased and then increased as the temperature was raised, but frequency shown in Figure 4B rose and then fell as the temperature increased. Figure 4D, E and G show that the parameters Δt f2 , Δt 1/2 , Δt r2 , Δt f1 , and V np,ave decreased as the temperature was raised with small standard errors for each value. On the other hand, the values of angles shown in Figure 4F increased as the temperature was raised, and the values of angle θ 1 were larger than those of θ 2. In addition, both A AP, max and A AP, ave shown in Figure 4H increased and then decreased as the temperature was raised. These temperature dependencies of AP parameters (except A AP, max , A AP, ave , and angles) were similar to the case of those analyzed using a dataset obtained from beating cells: in that case, A AP,ave decreased as the temperature was raised between 16°C and 28°C, but A AP, ave , and angles were not shown [1].

Temperature dependence of bursting patterns in burst-firing neurons
The values of six temperature-dependent burst parameters were averaged by using these selected dataset shown in Table 1. As shown in Figure 5A, B, and C respectively,  interburst interval, burst duration, and number of spike during burst decreased as temperature increased. ISI during burst shown in Figure 5D firstly decreased and then increased as temperature increased. In contrast, the numbers of burst per minute and of spike per minute shown in Figure 5E and F respectively increased and then decreased as the temperature was raised.
However, since these six burst parameters were not independent parameters, some equations representing an interrelationship between these parameters of burst-firing neurons could be driven; number of spike during burst = (burst duration)/(inter spike interval during burst), number of spike per minute = (number of spike during burst) × (number of burst per minute), and number of burst per minute = 60/(interburst interval + burst duration).
The values of these six temperature-dependent burst parameters during lowering the temperature (•) and raising the temperature (○) are shown in Figure 6. In the cases of interburst interval, burst duration, number of spike during burst, and ISI during burst shown in Figure 6A, B, C, and D respectively, there were no significant differences in these parameter values between the cooling-off and heat periods, which could be considered that there were reproducible properties of temperature dependencies in these burst parameters. Moreover, the number of burst per minute and spike per minute during heating shown in Figure 6E and F were slightly larger than those during cooling between 19°C and 25°C. However, it might be thought that the reproducible property of temperaturedependent changes in bursting patterns could not be undermined by these. It is because the values of the experimental percentage error of interspike interval during burst were between 0.1% and 14.1%; experimental percentage error ≡ 100% × | value of bursting parameter during heating (cooling) -average value of bursting parameter |/ average value of bursting parameter.
Thus, each mean value of functional forms of reproducible temperature-dependent bursting parameters in Figures 5 and 6 could be calculated by using these dataset and the computer program we designed. It was shown that these six bursting parameters changed in response to temperature variation, but the mechanism underlying these correlations remains unverified. So, it was suggested that performing computational simulations of these phenomena by using a modified Plant model which were composed of equations with temperature-dependent scaling factors was necessary to mathematically clarify the temperature-dependent changes of bursting patterns in burst-firing neurons.

Simulation of temperature-dependent bursting patterns
In order to simulate various temperature-dependent spiking patterns of bursting pacemaker neuron, we set up nonlinear differential equations by modifying the Plant model [15] with temperature-dependent scaling factors [25][26][27]. The Plant model studied by Rinzel and Lee was designed to assess parabolic bursters by analyzing the fast and slow processes to show how spike activities were generated by their mutual interaction [16].
The cell membrane of modified Plant model contains sodium channels carrying a fast sodium current, I Na , and potassium channels underlying a fast potassium current, I K . The slow processes include the slow inward calcium current, I Ca , and the slow changes in intracellular free calcium concentration, Ca. The accumulation of calcium turns on outward calcium-activated potassium current, I K(Ca) , and undermines the slow inward calcium current, and finally brings to an end of bursting activity. This system also contains a leak current, I L .
Where ; Where g Na , g Ca , g K , and g L are maximal conductances for the Na + , Ca 2+ , K + , and Cl − currents, respectively, and V Na , V Ca , V K , and V L are the reversal potentials for the respective currents: g K Ca ð Þ is maximal conductance for calcium-activated potassium current. Here, voltage, V indicates membrane potential (mV).
The voltage-dependent activation and inactivation variables for sodium channels are m and h respectively, and the activation variable for potassium channels is n. A slowly activating conductance for calcium current is denoted by χ, and a slow change in intracellular free calcium concentration Ca, is treated as parameters. The

Comparison of experimental and simulation results
In order to obtain good computer simulation of temperature-dependent bursting patterns generated respectively by these eight bursting pacemaker neurons, it was necessary to select three good data files of each experiment. The figures drawn by using data files selected during the same rising (or falling) phase of temperature change from dataset of experiment A at temperatures below 20°C, between 20°C and 25°C, above 25°C, were expressed by symbols L1, M1, H1, respectively, shown on Figure 7. Similarly, the figures shown by using data files selected from datasets of experiments B -H at the same temperature ranges were symbolized by L2-L8, M2 -M8, and H2 -H8, respectively, shown on Figure 7.
Then, it was necessary to describe how to draw three figures shown on the panels L1, M1, and H1 in Figure 7. Firstly, temperature values shown on these figures should be substituted for T in the temperature-dependent scaling factors, ρ(T) and ϕ(T) involved in modified Plant model. Secondly, computer simulation should be carried out until percentage errors of all parameters had to be calculated below 50%: the results of these calculations were shown at first three rows in Table 3. However, percentage error could be defined with the following formula shown below: Finally, numerical values of each of the 24 parameters involved in equations (1) to (4) would be fixed: C m = 1 μF/cm 2 , g KCa ¼ 0:018 mmho=cm 2 , g Na ¼ 4:0 mmho=cm 2 , g Ca ¼ 0:007 mmho=cm 2 , g K ¼ 0:60 mmho=cm 2  To simulate the other figures shown in Figure 7, it was necessary to fix the same numerical values of each of the 22 parameters as those used in experiment A, except the values of ρ and τ x . And then, computer simulation should be carried out by changing parameter values ρ and τ x until percentage errors of all parameters had to be calculated below 50%: interburst interval, burst duration, and number of spike per burst were increased (number of burst per minute decreased) as ρ, and τ x increased. The results of these calculations were shown from 4th row to the last row in Table 3: the values of all calculated percentage errors were between 0 and 47.7%.
Here, the results about the number of spike per minute among burst parameters was excluded from Table 3, because all percentage errors in the number of burst per minute were 0.0 and then the values of number of spike during burst would be the same as the values of number of spike per minute.
The results of computer simulation with the values of temperature shown on each panel in Figure 7 are represented as the time series of bursting activity. Although the results of simulation (red solid line) of AP amplitude were not completely optimal, simulation data for number of burst per minute, number of spike during burst, inter-burst interval, burst duration, and interspike interval during burst were well reflected by modified Plant model when these simulation results were compared to the experimental results (blue dotted line).
Therefore, analysis and simulation of the experimental data by our equation model might be helpful for understanding the mechanisms underlying the change of temperature-dependent bursting activity in neurons.

Discussion
In the previous study we analyzed the temperaturedependent change of 14 AP parameters (such as the AP amplitude, membrane potential at the positive peak, ISI, first half and last half of the temperature rising phase, first half and last half of the temperature falling phase, absolute value of the membrane potential at negative peak, absolute value of the maximum slope of the AP during the temperature rising and falling phases, and spiking frequency) in Aplysia neurons [1]. With these findings, in this study we tried to examine the temperature-dependent change of bursting patterns in pacemaker neuron at abdominal ganglion of Aplysia juliana. Furthermore, we attempted to identify the mechanism underlying the temperature dependence of busting by developing an equation model predicting neuronal electrophysiological activity with a function of temperature.
Percentage error ≡ the average value of the simulated data -the average value of the experimental data j j the average value of the experimental data Â 100 % Figure 6 Reproducible temperature-dependent bursting patterns of burst-firing neurons while decreasing (•) and increasing the temperature (○). There is no remarkable difference in each AP parameters when compared the values according to temperature between the cooling-off and heating period which means the reproducible properties of temperature-dependent burst parameters (A to F).

Temperature-dependent change of bursting patterns
Firstly we conducted real experiments with bursting pacemaker neurons in order to explore the functional relation between temperature and bursting patterns as well as AP parameters. We examined the temperature-dependent activity of bursting pacemaker neurons during several consecutive heating − cooling cycles. We tested the temperature dependencies of AP parameters in bursting pacemaker neurons in the temperature range of 16-30°C. The temperature-dependent properties of AP parameters were analyzed with the selected data shown in Table 1 by using techniques similar to those described in Hyun et al. [1]. The values of inter-burst interval, burst duration, and number of spike during burst decreased as the temperature increased, and it was possible to confirm functional properties of temperature dependencies in six burst parameters. When compared to data reported by other research group, Fletcher and Ram analyzed changes of AP parameters of burst-firing neurons during heating by investigating Aplysia R15 pacemaker neuron activity in the temperature range of 23 -37°C; they found that interburst interval, burst duration, and number of spike per burst all decreased within the temperature range of 23 to 32°C [21]. We found that our data are consistent with results reported by Fletcher and Ram [21] within the temperature range of 23 to 30°C. Furthermore, the present study revealed that the values of ISI and Δt r1 decreased and then increased, frequency increased and then decreased, and the parameters Δt f2 , Δt 1/2 , Δt r2 , Δtf 1 , and V np all decreased as the temperature was raised. These results were similar to those shown in previously published paper [1]. Fletcher and Ram also showed that intra-burst spike broadening and spike height decreased as temperature increased [21]. These results were also in good accordance with our experimental results within the temperature range of 23 to 30°C. Because Δt 1/2 and A AP,ave decreased as temperature increased in the temperature range of 23 to 30°C, as shown in Figure 4D and H: two parameters, spike-broadening and spike height, could be compared with two AP parameters, Δt AP,1/2 and A AP,ave , respectively. Moreover the present study generated quantitative data regarding the temperature-dependent properties of the other three AP parameters of burst-firing neurons  (ISI during burst, number of burst per minute, and number of spike per minute) and nine other AP parameters (ISI, Δt r1 , Δt r2 , Δt f1 , Δt f2 , θ 1 , θ2, V np , and A AP, max ).
Mechanism underlying the temperature dependence of busting: choice of the model Numerous revised model of bursting in R15 [17][18][19]28] have been developed since the mid-1980s. Canavier et al. [17] developed the model of bursting in R15 with 10 currents and 11 static variables; 10 currents (the fast Na + -current I Na , the fast Ca 2+ -current I Ca, the delayed rectifier K + -current I K , the slow inward Ca 2+ -current I SI , the nonspecific cation current I NS , the anomalous rectifier current I R , the leak current I L , K + -current I K , the Na + -Ca 2+ exchanger current I NaCa , the Na + -K + pump current I NaK , and the Ca 2+ pump current I CaP ), and 11 static variables (membrane potential V, intracellular concentration of Ca 2+ [Ca] i , the occupancy of the intracellular Ca 2+ -buffer O c , and eight voltage-dependent state-variables (m, h, d, f, n, l, s, b)). The parabolic bursting pattern results from two competing processes, the slow voltage-dependent activation of the slow inward Ca 2+ -current (I SI ), and the still more slower calciumdependent inactivation of this same current. Butera et al. [19] also developed the models of bursting in R15 with 10 currents and 12 static variables; 10 currents (the fast Na + -current I Na , the fast Ca 2+ -current I Ca, the delayed rectifier K + -current I K , the slow inward C 2+ -current I SI , the nonspecific cation current I NS , the anomalous rectifier current I R , the leak current I L , K + -current I K , the Na + -Ca 2+ exchanger current I NaCa , the Na + -K + pump current I NaK , and the Ca 2+ pump current I CaP ), and 12 static variables (membrane potential V, intracellular concentration of Ca 2+ [Ca 2+ ] i , intracellular concentration of cAMP [cAMP], the occupancy of the intracellular Ca 2+ -buffer O c , and eight voltage-dependent statevariables (m, h, d, f, n, l, s, b)). The slow inward Ca 2+ current (I SI ) in this model is the key current responsible for leading to bursting phenomena.
In addition to these models, Bertram [18,29] developed a mathematical model of bursting neuron R15 including 8 currents and 11 static variables; 8 currents (excitatory sodium current I Na , excitatory calcium current I Ca , inhibitory potassium current I K1 and I K2 , small leakage current I L , calcium current which initiate the burst I NSR , and cation-nonspecific current I D , and a potassium current which is activated only near the potassium equilibrium potential I R' ) and 11 static variables (membrane potential V, intracellular concentration of Ca 2+ c, and nine voltage-dependent state-variables (m, h, n, j, q, y, μ, r, x)). Two subthreshold burst currents, the calcium current which initiate the burst I NSR and the cation-nonspecific current I D , induce the busting oscillation.
To mathematically clarify the temperature-dependent changes of bursting patterns in burst-firing neurons, it was suggested to perform computational simulations of these phenomena by using a proper model composed of equations with temperature-dependent scaling factors.
We have tried to conduct computational simulations by using modified models: in modifying the models that were developed by Canavier et al. [17] and Butera et al. [19], the maximal conductances of currents g i were multiplied by temperature-like scaling factor for currents T-T 0 10 C , i = Na, Ca, K, SI, NS), and the relaxation time constant Tj were multiplied by temperature-like scaling factor for ionic kinetics h, d, f, n, l, s, b). In order to obtain computer simulation of temperaturedependent bursting patterns, it was necessary to select three good data files obtained from experiment A at different temperatures; 18.1°C, 22.1°C, and 27.9°C. The figures drawn by using these data files were expressed by symbols L1, M1 and H1, respectively, shown on the panels at Figure 8A and B. After temperature values shown on these panels should be substituted for T in the temperature-dependent scaling factors, ρ(T) and ϕ (T) involved in these modified models (T 0 = 23°C), we could get these figures drawn in red solid lines and these were shown on these panels. In the results of simulations with these modified models, it could be found that burst duration and number of spike during burst increased as temperature increased. But these results were not consistent with our experimental results. So, although the model presented by Canavier et al. [17] and Butera et al. [19] were based on experimental data to a great extent, it might be thought that these models are not appropriate for investigation about the ionic mechanism underlying the temperature-dependent change in bursting pattern from experiments with bursting pacemaker neurons in the abdominal ganglia of Aplysia juliana.
In a similar way, we could get three figures shown on Figure 8C by simulation using the modified models that were developed by Bertram [18]. Although the fact that the values of burst duration and number of spike during burst that were obtained from simulation by using this modified model decreased as temperature increased was consistent with our experimental results in a small temperature range, it was very difficult to simulate the temperature-dependent changes of bursting patterns in a large temperature range; from 18.1°C to 27.9°C. So we tried to simulate temperature-dependent bursting patterns in a small temperature range from 22.2°C to 25.8°C and the results were shown on Figure 8D. But these results were also not consistent with our experimental results. Namely, the values of interspike interval during burst obtained from simulations were larger than those got from experiments, and these values could not be reduced: the values of the other bursting parameters might be given rise to inconsistencies with our experimental results. Thus, it might be thought that this model was not also appropriate for investigation about that. Any way, it was shown that that I K(Ca) is not involved in burst termination in Aplysia bursting neurons Kramer and Zucker [30]. Even though calcium-activated potassium channels are located in the soma of R15 neuron, it was not necessarily proved that this channel was important for bursting pattern. Next, we have to discuss on the following question: are there any other experimental evidences suggesting that I K(Ca) is involved in bursting pattern of Aplysia neurons, particularly R15 neuron? Although we didn't yet find any one of research papers that were published after year 1985, it might be necessary to discuss whether the figures that we have suggested in this paper could be given us any chance to run back over the experimental evidence on it. Now, it might be necessary to compare the temperaturedependent changes in burst firing patterns of subthreshold currents (or burst currents) involved in Plant model with those included in the other models. Time series of temperature-dependent subthreshold current activities (ρ(T)•I SI ) simulated by using the equation of slow inward calcium current based on the model developed by Canavier et al. [17] at temperature values 18°C, 23°C, and 28°C are shown together in blue lines (L), black lines (M), and red lines (H), respectively on panel A in Figure 9. Those (ρ(T)•I SI ) and (ρ(T)•(I Ca + I K(Ca) )) obtained from simulation data by using the model made by Butera et al. [19] and Plant [15] at the same temperature values are shown together on panel B and D in Figure 9, respectively. Those (I D + ρ(T)•I NSR ) taken from simulation data calculated by using the model made by Bertram [18] at temperature values 22.3°C, 23°C , and 23.7°C are shown together on panel C in Figure 9. Next, it might be wanted to calculate the maximum values of temperature-dependent burst currents underlying the hyperpolarization of the inter-burst intervals, and then compare the values that were calculated by using Plant developed model with those by using the other models. The maximum values of these calculated by using Plant model at temperature values 18°C, 23°C, and 28°C were 0.28 μA, 0.12 μA, and 0.081 μA, respectively. The maximum values calculated by using Canavier et al. (Butera et al.) developed model at the same Anyway, these facts might show that it was hard to clarify mathematically the temperature-dependent changes of bursting patterns at burst-firing neurons of Aplysia juliana with these three models involving many currents and a large number of static variables. So, it was necessary to investigate much simpler model with small number of currents and static variables such as Plant model.
At this stage, it might be useful to remind that temperature-dependent impulse patterns of mammalian cold receptors could be well simulated by using nonlinear differential equations involving I K(Ca) . Different types of impulse patterns of mammalian cold receptors can be observed as a function of skin temperature: irregular and less frequent burst discharges, regular and frequent bursting patterns, and irregular single spike patterns are observed from low to high temperatures. These patterns could be simulated by Braun et al., and the Huber-Braun cold receptor model has been described in detail with 5 currents and 5 static variables (membrane potential and four voltage-dependent state-variables) [26,27,[31][32][33][34][35][36]. This model consisted of two minimal sets of ionic conductances operating at different voltage levels with different delays and the leakage current. The first two voltage-dependent currents that generate the action potentials mean depolarizing current I Na and repolarizing current I K . The next two voltage-dependent slow currents that generate subthreshold potential oscillations were slow depolarizing noninactivating Na + -current I Nap and slow repolarizing Ca 2+ -dependent K + -current I K(Ca) with voltage dependent activation of Ca 2+ -current. The temperature dependences were given by temperature-like scaling factors for the maximum conductances and the time constants, with reference temperature T 0 = 25°C. Q 10 of 3.0 for the activation variable and Q 10 of 1.3 for temperature dependences of maximum conductances were chosen. Now, it would be necessary to look into another simple model to investigate the mechanisms of temperaturedependent changes of bursting patterns which share a few similarities with mammalian cold receptors. Plant model consisted of 5 currents and 6 static variables: 5 currents (fast sodium current I Na , fast potassium current I K , inward calcium current I Ca , calcium-activated potassium current I K(Ca) , and leakage current, I L ) and 6 static variables (membrane potential, intracellular concentration of Ca 2+ , and four voltage-dependent state-variables). There are two fast and slow processes in this model. The fast process had three components: the activation and inactivation variables for Na + channels "m" and "h", respectively, and the  [18] are shown together. (D) Those (ρ(T)•(I Ca + I K(Ca) )) obtained from simulation data at temperature values 18°C, 23°C and 28°C using the model made by Plant [15] are shown together.
activation of K + channels n. The slow process had two components: a slow conductance for Ca 2+ -current "X", and intracellular free calcium concentration "Ca". Not only Plant model had a simple structure, but also theoretical analyses of Plant's model were performed already [16,37,38].
However, it is not easy to justify the reason of our choice of Plant model, because some papers excluded I K (Ca) as a key bursting current in R15 and focused on the inward current with slow activating and inactivating components. Plant model was based on the works of Gorman and Thomas [39], who demonstrated that I K (Ca) was linearly dependent on increasing concentration of intracellular calcium ions (Ca 2+ ) injected into the cytoplasm and expected to be activated during a burst. Chay [40] constructed mathematical model applicable to the Aplysia bursting neurons with involving the properties of Ca 2+ -activated K + channel: intracellular Ca 2+ concentration increases gradually during activity to levels where it activates the I K(Ca) [10]. Although Kramer and Zucker [30] concluded that I K(Ca) did not play a role in burst termination in bursting neurons of Aplysia, Adams and Levitan [41] did not exclude a role for I K(Ca) in the repolarization of action potential. But they asked the question: why a long-lasting I K(Ca) did not be activated by action potential? The plausible answer was the amount of calcium ions that entered into the cell during the activity of action potential was insufficient to generate it. Even if Canavier et al. [17] developed the model without I K(Ca), it was suggested as a limitations of this model that more experimental data of I K(Ca) were required for further investigation. Bertram [18] revised the model with I K(Ca) conductance depended only on the constant calcium concentration of the soma. Besides this, it has been known that the calcium activated potassium channels are located in R 15 soma [42] and electrode was inserted into the soma of the cell in the abdominal ganglion of A. juliana to measure the membrane potentials in our experiments. Thus, we cannot underestimate the impact of I K(Ca), but we did not here want to claim the quantitative application of Plant model, rather, we used it fr understanding the mechanism underlying the temperature dependence of busting patterns by computer simulations: it might be challenging to work later with a Rinzel and Lee's model based on the hypothesis that it is necessary for burst generation to consider an inward current which slowly inactivates as Ca 2+ accumulates during the burst to the exclusion of I K(Ca) [16]. By making a comparison between experimental data and simulation results calculated with modified Plant equations, it was suggested that the mechanism underlying temperature-dependent bursting patterns of bursting pacemaker neurons at abdominal ganglion of Aplysia juliana might be derived from temperature-like scaling factors ρ(T) and ϕ(T) for the maximum conductances and the time constants, respectively; together with reference temperature T 0 = 23°C, Q 10 of 3.0 for the activation variable, and Q 10 of 1.3 for temperature dependences of maximum conductances.
Taken together, it was suggested that a modified Plant model could be used to simulate the temperaturedependent bursting activity of bursting pacemaker neurons in the abdominal ganglia of Aplysia juliana and to unravel the mechanism of temperature dependence in bursting patterns.

Nomenclature
A AP , action potential amplitude A AP, ave , averaged action potential amplitude A AP, max , maximum action potential amplitude Ca, slow change in intracellular free calcium concentration Frequency, spontaneous firing frequency (=1/ISI) ISI, interspike interval (=Δt r1 + Δt AP, 1/2 + Δt f2 ) I Ca , slow inward calcium current I K , fast potassium current I K(Ca) , calcium-activated potassium current I L , leak current I Na , fast sodium current T 0 , reference temperature (=23°C) V, membrane potential V Ca , reversal potential for calcium current V K , reversal potential for potassium current V L , reversal potential for leak current V Na , reversal potential for sodium current V np , membrane potential at the negative peak V pp , membrane potential at the positive peak Δt r1 , first half of the rising phase of AP Δt r2 , last half of the rising phase of AP Δt f1 , first half of the falling phase of AP Δt f2 , last half of the falling phase of AP Δt AP, 1/2 , AP half-width duration (=Δt r2 + Δt f1 ) θ 1, inverse tangent of the ratio of half A AP to Δt r2 θ 2 , inverse tangent of the ratio of half A AP to Δt f1 g Ca , maximal conductance for calcium ion current g K , maximal conductance for potassium ion current g L , maximal conductance for leak current g Na , maximal conductance for sodium ion current h, voltage-dependent inactivation variable for sodium channel; h ∞ = steady state value of inactivation h m, voltage-dependent activation variable for sodium channel; m ∞ = steady state value of activation m n, activation variable for potassium channel; n ∞ = steady state value of activation n μ m , constant value involved in the defining equation for m ∞ μ h , constant value involved in the defining equation for h ∞ μ n , v n , constant values involved in the defining equation for n ∞ α, β, γ, δ, constant values involved in the defining equation for steady state values of activation (inactivation) variables and relaxation constants 1/λ, maximal relaxation time constant of h and n τ n , maximal relaxation time constant of n τ h, relaxation constant of h τ χ, relaxation constant of χ τ n, relaxation constant of n χ, slowly activating conductance for calcium current; χ ∞ = steady state value of activation χ ρ −1 , estimation for the time-constant of the Ca-equation ρ(T), temperature-dependent scaling factor; ρ T ð Þ≡1:3 T-T 0 10 C ϕ(T), temperature-dependent scaling factor; ϕ T ð Þ≡3 T-T 0 10 C I SI, slow inward Ca 2+ -current I NS, nonspecific cation current I R, anomalous rectifier current I NaCa, Na + -Ca 2+ exchanger current I NaK, Na + -K + pump current I CaP, Ca 2+ pump current I K1, inhibitory potassium current I NSR, calcium current which initiate the burst I D, cation-nonspecific current I R, potassium current which is activated only near the potassium equilibrium potential I Nap, slow depolarizing noninactivating Na + -current