In the previous study we analyzed the temperature-dependent change of 14 AP parameters(such as the AP amplitude, membrane potential at the positive peak, ISI, first half andlast half of the temperature rising phase, first half and last half of the temperaturefalling phase, absolute value of the membrane potential at negative peak, absolute valueof the maximum slope of the AP during the temperature rising and falling phases, andspiking frequency) in *Aplysia* neurons [[1]]. With these findings, in this study we tried to examine thetemperature-dependent change of bursting patterns in pacemaker neuron at abdominalganglion of *Aplysia juliana.* Furthermore, we attempted to identify themechanism underlying the temperature dependence of busting by developing an equationmodel predicting neuronal electrophysiological activity with a function oftemperature.

### Temperature-dependent change of bursting patterns

Firstly we conducted real experiments with bursting pacemaker neurons in order toexplore the functional relation between temperature and bursting patterns as well asAP parameters. We examined the temperature-dependent activity of bursting pacemakerneurons during several consecutive heating − cooling cycles. Wetested the temperature dependencies of AP parameters in bursting pacemaker neurons inthe temperature range of 16–30°C. The temperature-dependent properties ofAP 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 burstdecreased as the temperature increased, and it was possible to confirm functionalproperties of temperature dependencies in six burst parameters. When compared to datareported by other research group, Fletcher and Ram analyzed changes of AP parametersof burst-firing neurons during heating by investigating *Aplysia* R15pacemaker neuron activity in the temperature range of 23 – 37°C; theyfound that interburst interval, burst duration, and number of spike per burst alldecreased within the temperature range of 23 to 32°C [[21]]. We found that our data are consistent with results reported by Fletcherand Ram [[21]] within the temperature range of 23 to 30°C. Furthermore, the presentstudy revealed that the values of ISI and Δt_{r1} decreased and thenincreased, 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 wasraised. These results were similar to those shown in previously published paper [[1]]. Fletcher and Ram also showed that intra-burst spike broadening and spikeheight decreased as temperature increased [[21]]. These results were also in good accordance with our experimental resultswithin the temperature range of 23 to 30°C. Because Δt _{1/2} andA_{AP,ave} decreased as temperature increased in the temperature range of23 to 30°C, as shown in Figure 4D and H: twoparameters, spike-broadening and spike height, could be compared with two APparameters, Δt_{AP,1/2} and A_{AP,ave}, respectively. Moreoverthe present study generated quantitative data regarding the temperature-dependentproperties of the other three AP parameters of burst-firing neurons (ISI duringburst, number of burst per minute, and number of spike per minute) and nine other APparameters (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 themodel

Numerous revised model of bursting in R15 [[17]-[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 staticvariables; 10 currents (the fast Na^{+}-current I_{Na}, the fastCa^{2+}-current I_{Ca,} the delayed rectifierK^{+}-current I_{K}, the slow inward Ca^{2+}-currentI_{SI}, the nonspecific cation current I_{NS}, the anomalousrectifier current I_{R}, the leak current I_{L},K^{+}-current I_{K}, the Na^{+}-Ca^{2+} exchangercurrent I_{NaCa}, the Na^{+}-K^{+} pump currentI_{NaK}, and the Ca^{2+} pump current I_{CaP}), and 11static variables (membrane potential V, intracellular concentration ofCa^{2+} [Ca]_{i}, the occupancy of the intracellularCa^{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 competingprocesses, the slow voltage-dependent activation of the slow inwardCa^{2+}-current (I_{SI}), and the still more slower calcium-dependentinactivation of this same current. Butera et al. [[19]] also developed the models of bursting in R15 with 10 currents and 12static variables; 10 currents (the fast Na^{+}-current I_{Na}, thefast Ca^{2+}-current I_{Ca,} the delayed rectifierK^{+}-current I_{K}, the slow inward C^{2+}-currentI_{SI}, the nonspecific cation current I_{NS}, the anomalousrectifier current I_{R}, the leak current I_{L},K^{+}-current I_{K}, the Na^{+}-Ca^{2+} exchangercurrent I_{NaCa}, the Na^{+}-K^{+} pump currentI_{NaK}, and the Ca^{2+} pump current I_{CaP}), and 12static variables (membrane potential V, intracellular concentration ofCa^{2+} [Ca^{2+}]_{i}, intracellular concentration of cAMP[cAMP], the occupancy of the intracellular Ca^{2+}-buffer O_{c}, andeight voltage-dependent state-variables (m, h, d, f, n, l, s, b)). The slow inwardCa^{2+} current (I_{SI}) in this model is the key currentresponsible for leading to bursting phenomena.

In addition to these models, Bertram [[18],[29]] developed a mathematical model of bursting neuron R15 including 8currents and 11 static variables; 8 currents (excitatory sodium currentI_{Na}, excitatory calcium current I_{Ca}, inhibitory potassiumcurrent I_{K1} and I_{K2} , small leakage current I_{L},calcium current which initiate the burst I_{NSR}, and cation-nonspecificcurrent I_{D}, and a potassium current which is activated only near thepotassium equilibrium potential I_{R'}) and 11 static variables (membranepotential V, intracellular concentration of Ca^{2+} c, and ninevoltage-dependent state-variables (m, h, n, j, q, y, μ, r, x)). Two subthresholdburst currents, the calcium current which initiate the burst I_{NSR} and thecation-nonspecific current I_{D}, induce the busting oscillation.

To mathematically clarify the temperature-dependent changes of bursting patterns inburst-firing neurons, it was suggested to perform computational simulations of thesephenomena by using a proper model composed of equations with temperature-dependentscaling factors.

We have tried to conduct computational simulations by using modified models: inmodifying the models that were developed by Canavier et al. [[17]] and Butera et al. [[19]], the maximal conductances of currents {\overline{\mathrm{g}}}_{\mathrm{i}} were multiplied by temperature-like scaling factor for currents \mathrm{\rho}\left(\mathrm{T}\right)\phantom{\rule{0.12em}{0ex}}({\overline{\mathrm{g}}}_{\mathrm{i}}\times \mathrm{\rho}\left(\mathrm{T}\right)={\overline{\mathrm{g}}}_{\mathrm{i}}\times 1{.3}^{\frac{{\mathrm{T}\u2010\mathrm{T}}_{0}}{10\xb0\mathrm{C}}}, i = Na, Ca, K, SI, NS), and the relaxation timeconstant T_{j} were multiplied by temperature-like scaling factor for ionickinetics {\mathrm{T}}_{\mathrm{j}}\left(\left(1/{\mathrm{T}}_{\mathrm{j}}\right)\times \mathrm{\varphi}\left(\mathrm{T}\right)\right.=\left(1/{\mathrm{\tau}}_{\mathrm{j}}\right)\times {3}^{\frac{{\mathrm{T}\u2010\mathrm{T}}_{0}}{10\xb0\mathrm{C}}}, j = m, h, d, f, n, l, s, b). In order to obtaincomputer simulation of temperature-dependent bursting patterns, it was necessary toselect 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 datafiles were expressed by symbols L1, M1 and H1, respectively, shown on the panels atFigure 8A and B. After temperature values shown onthese panels should be substituted for T in the temperature-dependent scalingfactors, ρ(T) and ϕ(T) involved in these modified models(T_{0} = 23°C), we could get these figures drawn in redsolid lines and these were shown on these panels. In the results of simulations withthese modified models, it could be found that burst duration and number of spikeduring burst increased as temperature increased. But these results were notconsistent with our experimental results. So, although the model presented byCanavier et al. [[17]] and Butera et al. [[19]] were based on experimental data to a great extent, it might be thoughtthat these models are not appropriate for investigation about the ionic mechanismunderlying the temperature-dependent change in bursting pattern from experiments withbursting 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 spikeduring burst that were obtained from simulation by using this modified modeldecreased as temperature increased was consistent with our experimental results in asmall temperature range, it was very difficult to simulate the temperature-dependentchanges of bursting patterns in a large temperature range; from 18.1°C to27.9°C. So we tried to simulate temperature-dependent bursting patterns in asmall temperature range from 22.2°C to 25.8°C and the results were shown onFigure 8D. But these results were also not consistentwith our experimental results. Namely, the values of interspike interval during burstobtained from simulations were larger than those got from experiments, and thesevalues could not be reduced: the values of the other bursting parameters might begiven rise to inconsistencies with our experimental results. Thus, it might bethought 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 bursttermination in *Aplysia* bursting neurons Kramer and Zucker [[30]]. Even though calcium-activated potassium channels are located in the somaof R15 neuron, it was not necessarily proved that this channel was important forbursting pattern. Next, we have to discuss on the following question: are there anyother experimental evidences suggesting that I_{K(Ca)} is involved inbursting pattern of *Aplysia* neurons, particularly R15 neuron? Although wedidn’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 thispaper could be given us any chance to run back over the experimental evidence onit.

Now, it might be necessary to compare the temperature-dependent changes in burstfiring patterns of subthreshold currents (or burst currents) involved in Plant modelwith those included in the other models. Time series of temperature-dependentsubthreshold current activities (ρ(T)∙I_{SI}) simulated by usingthe equation of slow inward calcium current based on the model developed by Canavieret al. [[17]] at temperature values 18°C, 23°C, and 28°C are showntogether in blue lines (L), black lines (M), and red lines (H), respectively on panelA in Figure 9. Those (ρ(T)∙I_{SI}) and(ρ(T)∙(I_{Ca} + I_{K(Ca)})) obtained fromsimulation 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 inFigure 9, respectively. Those(I_{D} + ρ(T)∙I_{NSR}) taken fromsimulation data calculated by using the model made by Bertram [[18]] at temperature values 22.3°C, 23°C, and 23.7°C are showntogether on panel C in Figure 9. Next, it might be wantedto calculate the maximum values of temperature-dependent burst currents underlyingthe hyperpolarization of the inter-burst intervals, and then compare the values thatwere 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 values18°C, 23°C, and 28°C were 0.28 μA, 0.12 μA, and0.081 μA, respectively. The maximum values calculated by using Canavier etal. (Butera et al.) developed model at the same temperature values were−1.62 nA, −0.9 nA, and −0.75 nA (−1.5 nA,−1.1 nA, and −0.77 nA), respectively. The maximum valuescalculated by using the model that Bertram developed at the temperature values22.3°C, 23°C, and 23.7°C were −0.18 nA,−0.38 nA, and −0.58 nA, respectively. These absolute maximumvalues calculated by using the models that Plant, Canavier, and Butera developeddecreased as temperature increased, but those calculated by Bertram developed modelincreased as temperature increased. But more detailed analysis of these data remainedas a challenge for future study.

Anyway, these facts might show that it was hard to clarify mathematically thetemperature-dependent changes of bursting patterns at burst-firing neurons of*Aplysia juliana* with these three models involving many currents and alarge number of static variables. So, it was necessary to investigate much simplermodel with small number of currents and static variables such as Plant model.

At this stage, it might be useful to remind that temperature-dependent impulsepatterns of mammalian cold receptors could be well simulated by using nonlineardifferential equations involving I_{K(Ca)}. Different types of impulsepatterns of mammalian cold receptors can be observed as a function of skintemperature: irregular and less frequent burst discharges, regular and frequentbursting patterns, and irregular single spike patterns are observed from low to hightemperatures. These patterns could be simulated by Braun et al., and the Huber-Brauncold receptor model has been described in detail with 5 currents and 5 staticvariables (membrane potential and four voltage-dependent state-variables) [[26],[27],[31]-[36]]. This model consisted of two minimal sets of ionic conductances operatingat different voltage levels with different delays and the leakage current. The firsttwo voltage-dependent currents that generate the action potentials mean depolarizingcurrent I_{Na} and repolarizing current I_{K}. The next twovoltage-dependent slow currents that generate subthreshold potential oscillationswere slow depolarizing noninactivating Na^{+}-current I_{Nap} andslow repolarizing Ca^{2+}-dependent K^{+}-current I_{K(Ca)}with voltage dependent activation of Ca^{2+}-current. The temperaturedependences were given by temperature-like scaling factors for the maximumconductances and the time constants, with reference temperatureT_{0} = 25°C. Q_{10} of 3.0 for the activationvariable and Q_{10} of 1.3 for temperature dependences of maximumconductances were chosen.

Now, it would be necessary to look into another simple model to investigate themechanisms of temperature-dependent changes of bursting patterns which share a fewsimilarities with mammalian cold receptors. Plant model consisted of 5 currents and 6static variables: 5 currents (fast sodium current I_{Na}, fast potassiumcurrent I_{K}, inward calcium current I_{Ca}, calcium-activatedpotassium current I_{K(Ca)}, and leakage current, I_{L}) and 6 staticvariables (membrane potential, intracellular concentration of Ca^{2+}, andfour voltage-dependent state-variables). There are two fast and slow processes inthis model. The fast process had three components: the activation and inactivationvariables for Na^{+} channels “m” and “h”,respectively, and the activation of K^{+} channels n. The slow process hadtwo components: a slow conductance for Ca^{2+}-current “X”, andintracellular free calcium concentration “Ca”. Not only Plant model had asimple structure, but also theoretical analyses of Plant’s model were performedalready [[16],[37],[38]].

However, it is not easy to justify the reason of our choice of Plant model, becausesome papers excluded I_{K(Ca)} as a key bursting current in R15 and focusedon the inward current with slow activating and inactivating components. Plant modelwas based on the works of Gorman and Thomas [[39]], who demonstrated that I_{K(Ca)} was linearly dependent onincreasing concentration of intracellular calcium ions (Ca^{2+}) injectedinto 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+}-activatedK^{+} channel: intracellular Ca^{2+} concentration increasesgradually 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 terminationin bursting neurons of *Aplysia*, Adams and Levitan [[41]] did not exclude a role for I_{K(Ca)} in the repolarization ofaction 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 ofcalcium ions that entered into the cell during the activity of action potential wasinsufficient to generate it. Even if Canavier et al. [[17]] developed the model without I_{K(Ca),} it was suggested as alimitations of this model that more experimental data of I_{K(Ca)} wererequired for further investigation. Bertram [[18]] revised the model with I_{K(Ca)} conductance depended only on theconstant calcium concentration of the soma. Besides this, it has been known that thecalcium activated potassium channels are located in R_{15} soma [[42]] and electrode was inserted into the soma of the cell in the abdominalganglion of *A. juliana* to measure the membrane potentials in ourexperiments. Thus, we cannot underestimate the impact of I_{K(Ca),} but wedid not here want to claim the quantitative application of Plant model, rather, weused it fr understanding the mechanism underlying the temperature dependence ofbusting patterns by computer simulations: it might be challenging to work later witha Rinzel and Lee’s model based on the hypothesis that it is necessary for burstgeneration 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 resultscalculated with modified Plant equations, it was suggested that the mechanismunderlying temperature-dependent bursting patterns of bursting pacemaker neurons atabdominal ganglion of *Aplysia juliana* might be derived from temperature-likescaling factors ρ(T) and ϕ(T) for the maximum conductances and the timeconstants, respectively; together with reference temperatureT_{0} = 23°C, Q_{10} of 3.0 for the activationvariable, and Q_{10} of 1.3 for temperature dependences of maximumconductances.

Taken together, it was suggested that a modified Plant model could be used tosimulate the temperature-dependent bursting activity of bursting pacemaker neurons inthe abdominal ganglia of *Aplysia juliana* and to unravel the mechanism oftemperature 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}

{\overline{\mathrm{g}}}_{\mathrm{Ca}}, maximal conductance for calcium ion current

{\overline{\mathrm{g}}}_{\mathrm{K}}, maximal conductance for potassium ion current

{\overline{\mathrm{g}}}_{\mathrm{L}}, maximal conductance for leak current

{\overline{\mathrm{g}}}_{\mathrm{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 form_{∞}

μ_{h}, constant value involved in the defining equation forh_{∞}

μ_{n}, *v*_{n}, constant values involved in the defining equation forn_{∞}

α, β, γ, δ, constant values involved in the defining equation forsteady state values of activation (inactivation) variables and relaxationconstants

1/λ, maximal relaxation time constant of h and n

{\overline{\mathrm{\tau}}}_{\mathrm{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; \rho \left(\mathrm{T}\right)\equiv {1.3}^{\frac{{\mathrm{T}\u2010\mathrm{T}}_{0}}{{10}^{\xb0\mathrm{C}}}}

ϕ(T), temperature-dependent scaling factor; \varphi \left(\mathrm{T}\right)\equiv {3}^{\frac{{\mathrm{T}\u2010\mathrm{T}}_{0}}{{10}^{\xb0\mathrm{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 potassiumequilibrium potential

I_{Nap,} slow depolarizing noninactivating Na^{+}-current