Temperature-dependent bursting pattern analysis by modified Plant model

Abstract

Many electrophysiological properties of neuron including firing rates and rhythmicaloscillation change in response to a temperature variation, but the mechanismunderlying these correlations remains unverified. In this study, we analyzed variousaction potential (AP) parameters of bursting pacemaker neurons in the abdominalganglion of Aplysia juliana to examine whether or not bursting patterns arealtered in response to temperature change. Here we found that the inter-burstinterval, burst duration, and number of spike during burst decreased as temperatureincreased. On the other hand, the numbers of bursts per minute and numbers of spikesper minute increased and then decreased, but interspike interval during burst firstlydecreased and then increased. We also tested the reproducibility oftemperature-dependent changes in bursting patterns and AP parameters. Finally weperformed computational simulations of these phenomena by using a modified Plantmodel composed of equations with temperature-dependent scaling factors tomathematically clarify the temperature-dependent changes of bursting patterns inburst-firing neurons. Taken together, we found that the modified Plant model couldtrace the ionic mechanism underlying the temperature-dependent change in burstingpattern from experiments with bursting pacemaker neurons in the abdominal ganglia ofAplysia juliana.

Introduction

To date, researchers have investigated the effect of temperature on the electricalactivity and firing patterns in neurons from many animals, including Aplysiajuliana, crabs, crayfish, frogs, lobsters, locusts, snails, and squids [[1]-[9]]. Especially the functions and properties of bursting pacemaker neuron R15 inthe abdominal ganglion of Aplysia have been extensively studied [[10]-[12]] and mathematical simulation of bursting activity has been successfullyperformed [[13]-[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 byinvestigating the effect of heat on R15 bursting pacemaker neuron activity as thetemperature increased and reported temperature-dependent changes in inter-burstinterval, burst duration, number of spike per burst, intra-burst spike broadening andspike height [[21]]. However, the reproducible properties and mechanism of temperature-dependentchanges of AP parameters remain unknown yet.

Changes in temperature can produce numerous effects on the neural tissue of mostorganisms. Indeed it has been reported that hyperthermic temperature may inducedepolarization and spontaneous firing of pyramidal neurons leading to enhancedexcitability of hippocampus [[22],[23]]. On the other hand, a body of evidence indicating a therapeutic effect ofhypothermia has accumulated in several conditions. Orlowski and colleagues publishedthat in refractory status epilepticus unresponsive to conventional treatment, systemichypothermia (30-31°C) was an effective therapeutic method leading to burstsuppression on electroencephalography (EEG) [[24]]. Despite a lot of works aiming at demonstrating the neural effect oftemperature, the mechanism of temperature-dependent electrophysiological changecontaining alteration of bursting and firing pattern in neuron remains unclear. Here, weexamined the effects of temperature changes on the neuronal activity and burstingpatterns during several consecutive heating − cooling cycles by usingbursting pacemaker neurons which are a proper specimen with capability of long-lastingrecording for mathematical modeling. Next, we sought to identify the mechanismunderlying temperature-dependent bursting patterns of these neurons by analyzing andcomparing the experimental data to computational simulation data calculated by modifiedPlant equations with temperature-dependent scaling factors, ρ(T) and ϕ(T).

Materials and methods

Animals and dissection

Animals (A. juliana) were collected locally in Seogwipo, South Korea and allexperimental procedures were approved by the Jeju National University Animal Care andUse Committee. These animals were dissected to observe the temperature-dependency ofbursting patterns generated by pacemaker neurons R15 in the abdominal ganglia ofA. juliana as described in previously published paper [[1]]. Briefly, the animals were anesthetized with an injection of 0.38 MMgCl₂ amounting to half of each specimen’s weight before theabdominal ganglia were removed. Each abdominal ganglion was incubated at 34°Cfor 40 minute in a solution containing equal volumes of artificial sea water(ASW; in mM: 460 NaCl, 10 KCl, 11 CaCl₂ , 55 MgCl₂ , and 10HEPES; 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 werethen 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 removedonto Petri dishes (50x9 mm), pinned down on Sylgard plates (Dow Corning, USA) filledwith 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 Aplysiasoaked in ASW:L-15 (1:1, v/v) to measure the medium temperature, which could beincreased or decreased by activating the temperature controller system (HMN 3940,Acetec Co., Korea). We used a low flow rate peristaltic pump (BJ100-2 J; BaodingLonger 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 performedto measure the membrane potential. The glass intracellular electrodes were filledwith 3 M KCl and the membrane potentials recorded were simultaneously savedusing a DAQ card (NI PCI-6221, National Instruments) and the Labview program(National Instruments). The electrical signals of the APs were identified using adigital oscilloscope (54622A; Agilent, Colorado Springs, Colo., USA). These data wererecorded at a rate of 3 kHz to reproduce the APs on the computer because thebursting cells generated APs at a rate of 0–4 Hz; each component of thisdata is composed of 180,000 pairs of temperature and membrane potential values withunique file names, and the dataset was composed of 1,713 files as shown inTable 1. With this selected dataset, the average valuesof each of the AP parameters of the burst-firing neuron could be calculated usingOrigin 6.0 (Microcal Software, Inc.) and a computer program we designed. For eachexperiment, the bursting neuron was held at room temperature for one hour as shown inTable 2 and then changed with continuous ramp fromabout 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-28-30°C) wereaveraged over. To reduce experimental variability, in each experiment the burstingneuron was incubated at least 5 min to adapt to a new steady-state temperaturebefore recording rhythmic change.

Table 1 Animal weights and electrophysiological properties of signals generated byeight R15 bursting pacemaker neurons in A. juliana specimens

Table 2 Temperature and selected AP parameters in bursting pacemaker neuronsmaintained for 30 min

Results

Definition of AP parameters in burst-firing neuron

The graph in Figure 1A represents the typical burstingpatterns of electrical signals shown in Aplysia bursting pacemaker neuron.The intraburst interspike interval, interburst interval, and burst duration arerepresented by symbols (1), (2), and (3), respectively. The membrane potentials atthe positive peak, V_pp (mV), and negative peak, V_np (mV), aredefined as the values of membrane potentials at points T and B shown inFigure 1B, respectively. The definition of the firsthalf of the rising phase of AP, Δt_r1 (ms), and the following APparameters are very similar to those defined in our previous study [[1]]. The last half of the rising phase of AP, Δt_r2 (ms) andthe first half of the falling phase of AP, Δt_f1 (ms) are timeintervals shown in Figure 1B. The last half of the fallingphase Δt_f2 (ms) is defined as the values of the time interval betweentwo points H2 and B shown in Figure 1B. The AP half-widthduration, Δt_{AP, 1/2} (ms), is defined asΔt_r2 + Δt_f1 (ms). The interspikeinterval, ISI (ms), is defined as Δt_r1 + Δt_AP,1/2 + Δt_f2. The spontaneous firing frequency,simply referred to as Frequency (s⁻¹), is defined asISI⁻¹ (s⁻¹). For convenience, the average anglesθ₁ and θ₂ are defined to compare the slope of thesecond half of the rising phase of the AP with that of the first half of the fallingphase. The angle θ₁ is defined as the inverse tangent of the ratio ofhalf A_AP to Δt_r2, and the angle θ₂ is alsodefined as the inverse tangent of the ratio of half A_AP toΔt_f1: ${\mathrm{\theta }}_1\equiv \mathrm{Arctan}\left[\frac{\mathrm{V}\left({\mathrm{t}}_{\mathrm{T}}\right)‐\mathrm{V}\left({\mathrm{t}}_{\mathrm{H}1}\right)}{{\mathrm{t}}_{\mathrm{T}}‐{\mathrm{t}}_{\mathrm{H}1}}\right]$ and ${\mathrm{\theta }}_2\equiv \mathrm{Arctan}\left[\frac{\mathrm{V}\left({\mathrm{t}}_{\mathrm{T}}\right)‐\mathrm{V}\left({\mathrm{t}}_{\mathrm{H}2}\right)}{{\mathrm{t}}_{\mathrm{H}2}‐{\mathrm{t}}_{\mathrm{T}}}\right]$. 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 thetwo values of membrane potentials at points T and B, respectively. More detailedanalytical techniques regarding the definition and computational analysis of APparameters 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 werepresented as mean ± SE (standard error).

Figure 1

Typical bursting patterns and AP parameters of bursting pacemakerneuron. (A) Bursting patterns of bursting pacemaker neuron; membrane potentialvs. time. The electrophysiological parameters related to burst firing containthe intraburst interspike interval, the duration between the two points ofsequential V_pp within one burst firing event [A (1)];interburst interval, the duration between the two points of the lastV_pp of one burst firing event and the first V_pp ofnext burst firing event [A (2)]; and burst duration, the durationbetween the two points of the first V_pp and the last V_ppwithin one burst firing event [A (3)]. (B) Definition of APparameters. The maximum point T on the upper part of AP means a positive peakmembrane potential, V_pp and the minimum point B on the lower part ofAP means a negative peak membrane potential, V_np. H1 and H2 meansthe half point of potential difference during the rising phase and the fallingphase of AP, respectively. The last half of the rising phase,Δt_r2 designate the time intervals between the two point H1and T. The first half of the falling phase of AP, Δt_f1 isdefined as the time interval between point T and point H2. The angleθ₁ shows the inverse tangent of the ratio of halfA_AP to Δt_r2 and angle θ2 means the inversetangent of the ratio of half A_AP to Δt_f1.

Structure of the selected dataset for AP parameter analysis

Eight experiments (A − H) were conducted using A. julianaspecimens, and Table 1 shows the animal weights, selecteddata for analysis, total number of spike, and total recording time.

Animal weights were between 130 and 310 g, and the average weight was226 g. The average total number of spike and average total recording time ofeach experiment were 15,039 spikes and 607 min, respectively. It was necessaryto select the middle portion of each recorded dataset for analysis. The averagevalues of time intervals for analysis, numbers of cycles of temperature change, andnumber of spike selected for analysis of these selected data were 214 min,2 cycles, and 5,433 spikes, respectively. The dataset gathered for analysis wascomposed of these selected data from each experiment. The total time interval foranalysis and the total number of spike of this dataset were 1,713 minutes and43,468 spikes, respectively. The numbers written in parentheses in the middle columnof Table 1 represent the total number of cycles oftemperature change for each experiment.

In each experiment, the temperature was maintained at room temperature for the firstone hour. As shown in Table 2, average values oftemperature and four AP parameters of burst-firing neurons were calculated using theselected dataset saved by intracellular recordings during maintaining at roomtemperature for 30 minutes; all data were selected after 5 to 10 minutesfrom temperature change onset. Because the standard errors of the values of theseparameters were low, the bursting patterns of these experiments might be consideredas 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 recordingdata of experiment A during the fourth falling phase (from 440 to 487 min)(Figure 2) out of the eight falling phases intemperature recorded during eight consecutive heating − coolingcycles: it was composed of 48 panels. Temperature and membrane potentials arerepresented 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, weextracted six burst signals (6 panels) among many burst signals to examine thetemporal change of instantaneous ISI during burst at specific temperature;30.0°C (at 440 min), 26.9°C (at 454 min), 23°C (at464 min), 21.3°C (at 468 min), 19.1°C (at 474 min), and16.4°C (487 min). As shown in Figure 3, ISIduring burst versus time displayed parabolic pattern, suggesting thattemperature-dependent bursting patterns are parabolic bursts characterized by lowerspike frequency at the beginning and end of the burst. The analyzing dataset wascomposed of 1,713 files corresponding to the sum of each selected time interval foranalysis as shown in the 4th column of the Table 1.

Figure 2

Continuous bursting patterns obtained from the dataset of intracellularrecordings of experiment A during the 4th falling phase (from 440 to487 minutes) out of eight recorded falling phases. The temperaturevalues of each panel is as follows; 30.1 ~ 27.4°C (top),27.4 ~ 23°C (2nd from the top), 23 ~ 18.6°C(3rd from the top), and 18.6 ~ 16.4°C (bottom). As temperaturedecreased, the interburst interval, burst duration, and number of spike duringburst increased.

Figure 3

Plot of interspike intervals during a burst versus time. (A) Figuresdrawn by using data calculated from the burst shown on panels inFigure 2: parabolic shape. The values oftemperature corresponding to the figures from the left to the right are30.0°C, 26.9°C, 23.0°C, 21.3°C, 19.1°C, and16.4°C, sequentially. (B) The values of each figure are the same asthose shown in (A): all the figures were shown with the same start time.

With the selected data shown in Table 1, thetemperature-dependent properties of AP parameters were analyzed by using techniquessimilar to those described in Hyun et al. [[1]]. As shown in Figure 4A and C, all values ofISI and Δt_r1 decreased and then increased as the temperature wasraised, but frequency shown in Figure 4B rose and thenfell as the temperature increased. Figure 4D, E and G showthat the parameters Δt_f2, Δt_1/2,Δt_r2, Δt_f1, and V_np,ave decreased as thetemperature was raised with small standard errors for each value. On the other hand,the values of angles shown in Figure 4F increased as thetemperature was raised, and the values of angle θ₁ were larger thanthose of θ_2. In addition, both A_{AP, max} and A_AP,ave shown in Figure 4H increased and thendecreased as the temperature was raised. These temperature dependencies of APparameters (except A_{AP, max}, A_{AP, ave}, and angles) were similarto the case of those analyzed using a dataset obtained from beating cells: in thatcase, A_AP,ave decreased as the temperature was raised between 16°Cand 28°C, but A_{AP, ave}, and angles were not shown [[1]].

Figure 4

Temperature-dependent AP parameters of burst-firing neurons. The valuesof ISI (A) and Δt_r1(C) decreased and then increased as temperature increase, but frequency(B) and the values of A_{AP, max} and A_{AP, ave}(H) showed the reversed pattern according to the change of temperature.Δt_f2 and Δt_1/2(D)_, Δt_r2 and Δt_f1(E), V_np,ave(G) decreased as temperature increase and the values of angleθ₁ and θ₂(F) showed the reversed manner.

Temperature dependence of bursting patterns in burst-firing neurons

The values of six temperature-dependent burst parameters were averaged by using theseselected dataset shown in Table 1. As shown inFigure 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 andthen increased as temperature increased. In contrast, the numbers of burst per minuteand of spike per minute shown in Figure 5E and Frespectively increased and then decreased as the temperature was raised.

Figure 5

Temperature-dependent bursting patterns of burst-firing neurons.Interburst interval (A), burst duration (B), number of spikeduring burst (C) decreased as temperature increased. ISI during burst(D) showed a parabolic pattern which means that it decreased and thenincreased as temperature increased. Number of burst per minute (E) andnumber of spike per minute (F) increased and then decreased showing thereversed pattern of ISI during burst as temperature increased. Theinstantaneous ISI during a burst versus time firstly decreased and thenincreased showing the parabolic bursting pattern at each temperature rangingfrom 30°C to 16°C (G).

However, since these six burst parameters were not independent parameters, someequations representing an interrelationship between these parameters of burst-firingneurons could be driven; number of spike during burst = (burstduration)/(inter spike interval during burst), number of spike perminute = (number of spike during burst) × (number of burst perminute), and number of burst per minute = 60/(interburstinterval + burst duration).

The values of these six temperature-dependent burst parameters during lowering thetemperature (●) and raising the temperature (○) are shown inFigure 6. In the cases of interburst interval, burstduration, number of spike during burst, and ISI during burst shown inFigure 6A, B, C, and D respectively, there were nosignificant differences in these parameter values between the cooling-off and heatperiods, which could be considered that there were reproducible properties oftemperature dependencies in these burst parameters. Moreover, the number of burst perminute and spike per minute during heating shown in Figure 6E and F were slightly larger than those during cooling between 19°Cand 25°C. However, it might be thought that the reproducible property oftemperature-dependent changes in bursting patterns could not be undermined by these.It is because the values of the experimental percentage error of interspike intervalduring burst were between 0.1% and 14.1%; experimental percentageerror ≡ 100% × | value of bursting parameter during heating(cooling) - average value of bursting parameter |/ average value of burstingparameter.

Figure 6

Reproducible temperature-dependent bursting patterns of burst-firing neuronswhile decreasing (●) and increasing the temperature (○).There is no remarkable difference in each AP parameters when compared thevalues according to temperature between the cooling-off and heating periodwhich means the reproducible properties of temperature-dependent burstparameters (A to F).

Thus, each mean value of functional forms of reproducible temperature-dependentbursting parameters in Figures 5 and 6 could be calculated by using these dataset and the computer program wedesigned. It was shown that these six bursting parameters changed in response totemperature variation, but the mechanism underlying these correlations remainsunverified. So, it was suggested that performing computational simulations of thesephenomena by using a modified Plant model which were composed of equations withtemperature-dependent scaling factors was necessary to mathematically clarify thetemperature-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 burstingpacemaker neuron, we set up nonlinear differential equations by modifying the Plantmodel [[15]] with temperature-dependent scaling factors [[25]-[27]]. The Plant model studied by Rinzel and Lee was designed to assessparabolic bursters by analyzing the fast and slow processes to show how spikeactivities were generated by their mutual interaction [[16]].

The cell membrane of modified Plant model contains sodium channels carrying a fastsodium current, I_Na, and potassium channels underlying a fast potassiumcurrent, 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 finallybrings to an end of bursting activity. This system also contains a leak current,I_L.

$$\begin{array}{c}\hfill -{\mathrm{C}}_{\mathrm{m}}\frac{\mathrm{dV}}{\mathrm{dt}}=\mathrm{\rho }\left(\mathrm{T}\right)\left({\mathrm{I}}_{\mathrm{Na}}+{\mathrm{I}}_{\mathrm{Ca}}+{\mathrm{I}}_{\mathrm{K}}+{\mathrm{I}}_{\mathrm{K}\left(\mathrm{Ca}\right)}\right)+{\mathrm{I}}_{\mathrm{L}}\hfill \\ \hfill =\mathrm{\rho }\left(\mathrm{T}\right)\left\{{\overline{\mathrm{g}}}_{\mathrm{Na}}\cdot {\mathrm{m}}_{\infty }^3\left(\mathrm{V}\right)\cdot \mathrm{h}\cdot \left({\mathrm{V}‐\mathrm{V}}_{\mathrm{Na}}\right)+{\overline{\mathrm{g}}}_{\mathrm{Ca}}\cdot \mathrm{\chi }\cdot \left({\mathrm{V}‐\mathrm{V}}_{\mathrm{Ca}}\right)+\left({\overline{\mathrm{g}}}_{\mathrm{K}}\cdot {\mathrm{n}}^4+{\overline{\mathrm{g}}}_{\mathrm{K}‐\mathrm{Ca}}\frac{\mathrm{Ca}}{0.2+\mathrm{Ca}}\right)\cdot \left({\mathrm{V}‐\mathrm{V}}_{\mathrm{K}}\right)+{\overline{\mathrm{g}}}_{\mathrm{L}}\cdot \left({\mathrm{V}‐\mathrm{V}}_{\mathrm{L}}\right),\right.\hfill \end{array}$$

(1)

$$\begin{array}{c}\hfill \frac{\mathrm{dh}}{\mathrm{dt}}=\mathrm{\varphi }\left(\mathrm{T}\right)\cdot \mathrm{\lambda }\frac{\left[{\mathrm{h}}_{\infty }\left(\mathrm{V}\right)‐\mathrm{h}\right]}{{\mathrm{\tau }}_{\mathrm{h}}\left(\mathrm{V}\right)},\hfill \\ \hfill \frac{\mathrm{d\chi }}{\mathrm{dt}}=\mathrm{\varphi }\left(\mathrm{T}\right)\frac{{\mathrm{\chi }}_{\infty }\left(\mathrm{V}\right)-\mathrm{\chi }}{\mathrm{\tau \chi }},\hfill \\ \hfill \frac{\mathrm{dn}}{\mathrm{dt}}=\mathrm{\varphi }\left(\mathrm{T}\right)\cdot \mathrm{\lambda }\frac{\left[{\mathrm{n}}_{\infty }\left(\mathrm{V}\right)-\mathrm{n}\right]}{{\mathrm{\tau }}_{\mathrm{n}}\left(\mathrm{V}\right)},\hfill \\ \hfill \frac{\mathrm{dCa}}{\mathrm{dt}}=\mathrm{\rho }\cdot \left[{\mathrm{K}}_{\mathrm{c}}\cdot \mathrm{\chi }\cdot \left({\mathrm{V}}_{\mathrm{Ca}}-\mathrm{V}\right)-\mathrm{Ca}\right],\hfill \end{array}$$

(2)

Where

$$\begin{array}{c}\hfill {\mathrm{m}}_{\infty }\left(\mathrm{V}\right)=\frac{\mathrm{\mu m}\frac{50-\mathrm{\alpha V}-\mathrm{\beta }}{{\mathrm{e}}^{\left(\frac{50-\mathrm{\alpha V}-\mathrm{\beta }}{10}\right)}-1}}{\mathrm{\mu m}\frac{50-\mathrm{\alpha V}-\mathrm{\beta }}{{\mathrm{e}}^{\left(\frac{50-\mathrm{\alpha V}-\mathrm{\beta }}{10}\right)}-1}+4{\mathrm{e}}^{\left(\frac{25-\mathrm{\alpha V}-\mathrm{\beta }}{18}\right)}},\hfill \\ \hfill {\mathrm{h}}_{\infty }\left(\mathrm{V}\right)=\frac{{\mathrm{\mu }}_{\mathrm{h}}{\mathrm{e}}^{\left(\frac{25-\mathrm{\alpha V}-\mathrm{\beta }}{20}\right)}}{{\mathrm{\mu }}_{\mathrm{h}}{\mathrm{e}}^{\left(\frac{25-\mathrm{\alpha V}-\mathrm{\beta }}{20}\right)}+\frac{1}{{\mathrm{e}}^{\left(\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{10}\right)}}},\hfill \\ \hfill {\mathrm{\chi }}_{\infty }\left(\mathrm{V}\right)=\frac{1}{1+{\mathrm{e}}^{\mathrm{\gamma }\left(\mathrm{\delta }-\mathrm{V}\right)}},\hfill \end{array}$$

(3)

and

$$\begin{array}{c}\hfill {\mathrm{n}}_{\infty }\left(\mathrm{V}\right)=\frac{\mathrm{\mu n}\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{{\mathrm{e}}^{\left(\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{10}\right)}}}{\mathrm{\mu n}\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{{\mathrm{e}}^{\left(\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{10}\right)}-1}+{\mathrm{\nu }}_{\mathrm{n}}\mathrm{e}\left(\frac{45-\mathrm{\alpha V}-\mathrm{\beta }}{80}\right)},\hfill \\ \hfill {\mathrm{\tau }}_{\mathrm{h}}\left(\mathrm{V}\right)=\frac{1}{{\mathrm{\mu }}_{\mathrm{h}}{\mathrm{e}}^{\left(\frac{25-\mathrm{\alpha V}-\mathrm{\beta }}{20}\right)}+\frac{1}{{\mathrm{e}}^{\left(\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{10}\right)}+1}},\hfill \end{array}$$

(4)

and

$${\mathrm{\tau }}_{\mathrm{n}}\left(\mathrm{V}\right)=\frac{{\overline{\mathrm{\tau }}}_{\mathrm{n}}}{{\mathrm{\mu }}_{\mathrm{n}}\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{{\mathrm{e}}^{\left(\frac{55-\mathrm{\alpha V}-\mathrm{\beta }}{10}\right)}-1}+{\mathrm{\nu }}_{\mathrm{n}}{\mathrm{e}}^{\left(\frac{45-\mathrm{\alpha V}-\mathrm{\beta }}{80}\right)}}$$

(5)

Where ${\overline{\mathrm{g}}}_{\mathrm{Na}}$, ${\overline{\mathrm{g}}}_{\mathrm{Ca}}$, ${\overline{\mathrm{g}}}_{\mathrm{K}}$, and ${\overline{\mathrm{g}}}_{\mathrm{L}}$ are maximal conductances for the Na⁺, Ca²⁺,K⁺, and Cl⁻ currents, respectively, andV_Na, V_Ca, V_K, and V_L are the reversalpotentials for the respective currents: ${\overline{\mathrm{g}}}_{\mathrm{K}\left(\mathrm{Ca}\right)}$ 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 arem and h respectively, and the activation variable for potassiumchannels is n. A slowly activating conductance for calcium current isdenoted by χ, and a slow change in intracellular free calcium concentration Ca,is treated as parameters. The maximal relaxation time constants of h and n aredefined as $\frac{1}{\mathrm{\lambda }}$. It is taken that ρ⁻¹ is an estimate for thetime-constant of the Ca-equation. The temperature-dependent scaling factors,ρ(T) and ϕ(T), are defined as $\mathrm{\rho }\left(\mathrm{T}\right)\equiv {1.3}^{\frac{{\mathrm{T}‐\mathrm{T}}_0}{10°\mathrm{C}}}$ and $\mathrm{\varphi }\left(\mathrm{T}\right)\equiv 3^{\frac{\mathrm{T}-{\mathrm{T}}_0}{10°\mathrm{C}}}$, respectively [[25]]. The steady-state values of activation or inactivation variablesm_∞, h_∞, χ_∞, andn_∞ are functions of voltage. The relaxation time constants arerepresented by τ_h, τ_χ, andτ_n.

Comparison of experimental and simulation results

In order to obtain good computer simulation of temperature-dependent burstingpatterns generated respectively by these eight bursting pacemaker neurons, it wasnecessary to select three good data files of each experiment. The figures drawn byusing data files selected during the same rising (or falling) phase of temperaturechange from dataset of experiment A at temperatures below 20°C, between20°C and 25°C, above 25°C, were expressed by symbols L1, M1, H1,respectively, shown on Figure 7. Similarly, the figuresshown by using data files selected from datasets of experiments B – H at thesame temperature ranges were symbolized by L2- L8, M2 - M8, and H2 - H8,respectively, shown on Figure 7.

Figure 7

Time series of bursting activity drawn by using the data for the experimentsfrom A to H and computer simulations (A to H from top to bottom panel).The results of experiments (blue, dotted line) and computer simulation (red,solid line) are represented as the activity with the same correspondingtemperatures shown on each panel. The figures shown by using data filesselected from datasets of experiments A – H (from top tobottom panel) were symbolized by L1- L8, M1 - M8, and H1 - H8, at thetemperature ranges below 20°C, between 20°C and 25°C, above25°C, respectively. Computer simulations should be carried out untilpercentage errors of all parameters had to be calculated below 50%.

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 shownon these figures should be substituted for T in the temperature-dependent scalingfactors, ρ(T) and ϕ(T) involved in modified Plant model. Secondly, computersimulation should be carried out until percentage errors of all parameters had to becalculated below 50%: the results of these calculations were shown at first threerows in Table 3. However, percentage error could bedefined with the following formula shown below:

$$\begin{array}{l}\mathrm{Percentage}\mathrm{error}\equiv \\ \frac{\left|\mathrm{the}\mathrm{average}\mathrm{value}\mathrm{of}\mathrm{the}\mathrm{simulated}\mathrm{data}‐\mathrm{the}\mathrm{average}\mathrm{value}\mathrm{of}\mathrm{the}\mathrm{experimental}\mathrm{data}\right|}{\mathrm{the}\mathrm{average}\mathrm{value}\mathrm{of}\mathrm{the}\mathrm{experimental}\mathrm{data}}\times 100\%\end{array}$$

(6)

Table 3 Comparison of experimental and simulation results

Finally, numerical values of each of the 24 parameters involved in equations (1) to(4) would be fixed: C_m = 1 μF/cm², ${\overline{\mathrm{g}}}_{\mathrm{KCa}}=0.018\mathrm{mmho}/{\mathrm{cm}}^2$, ${\overline{\mathrm{g}}}_{\mathrm{Na}}=4.0\mathrm{mmho}/{\mathrm{cm}}^2$, ${\overline{\mathrm{g}}}_{\mathrm{Ca}}=0.007\mathrm{mmho}/{\mathrm{cm}}^2$, ${\overline{\mathrm{g}}}_{\mathrm{K}}=0.60\mathrm{mmho}/{\mathrm{cm}}^2$, ${\overline{\mathrm{g}}}_{\mathrm{L}}=0.017\mathrm{mmho}/{\mathrm{cm}}^2$, V_Na = 40 mV,V_Ca = 140 mV,V_K = −75 mV,V_L = −40 mV, λ = 0.18,ρ = 0.000074 ms⁻¹,τ_χ = 1, 500,K_c = 0.0275 mV⁻¹,α = 127/105, β = 8265/105,γ = 0.3, δ = −18,μ_m = 0.1, μ_h = 0.08,μ_n = 0.016, ν_n = 0.1, ${\overline{\mathrm{\tau }}}_{\mathrm{n}}=1$, T₀ = 23.0°C. Here, we took a value of23°C as the reference temperature T₀ because bursts are usuallyactivated from 22°C to 25°C, and this represented the middle value oftemperature range of the experiment A; from 16.0°C to 30.0°C.

To simulate the other figures shown in Figure 7, it wasnecessary to fix the same numerical values of each of the 22 parameters as those usedin 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 calculatedbelow 50%: interburst interval, burst duration, and number of spike per burst wereincreased (number of burst per minute decreased) as ρ, and τ_xincreased. The results of these calculations were shown from 4th row to the last rowin Table 3: the values of all calculated percentage errorswere between 0 and 47.7%. Here, the results about the number of spike per minuteamong burst parameters was excluded from Table 3, becauseall percentage errors in the number of burst per minute were 0.0 and then the valuesof number of spike during burst would be the same as the values of number of spikeper minute.

The results of computer simulation with the values of temperature shown on each panelin Figure 7 are represented as the time series of burstingactivity. Although the results of simulation (red solid line) of AP amplitude werenot completely optimal, simulation data for number of burst per minute, number ofspike during burst, inter-burst interval, burst duration, and interspike intervalduring burst were well reflected by modified Plant model when these simulationresults were compared to the experimental results (blue dotted line).

Therefore, analysis and simulation of the experimental data by our equation modelmight be helpful for understanding the mechanisms underlying the change oftemperature-dependent bursting activity in neurons.

Discussion

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₁, 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, θ₁, θ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²⁺-current I_Ca, the delayed rectifierK⁺-current I_K, the slow inward Ca²⁺-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²⁺ exchangercurrent I_NaCa, the Na⁺-K⁺ pump currentI_NaK, and the Ca²⁺ pump current I_CaP), and 11static variables (membrane potential V, intracellular concentration ofCa²⁺ [Ca]_i, the occupancy of the intracellularCa²⁺-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²⁺-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²⁺-current I_Ca, the delayed rectifierK⁺-current I_K, the slow inward C²⁺-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²⁺ exchangercurrent I_NaCa, the Na⁺-K⁺ pump currentI_NaK, and the Ca²⁺ pump current I_CaP), and 12static variables (membrane potential V, intracellular concentration ofCa²⁺ [Ca²⁺]_i, intracellular concentration of cAMP[cAMP], the occupancy of the intracellular Ca²⁺-buffer O_c, andeight voltage-dependent state-variables (m, h, d, f, n, l, s, b)). The slow inwardCa²⁺ 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²⁺ 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)({\overline{\mathrm{g}}}_{\mathrm{i}}\times \mathrm{\rho }\left(\mathrm{T}\right)={\overline{\mathrm{g}}}_{\mathrm{i}}\times 1{.3}^{\frac{{\mathrm{T}‐\mathrm{T}}_0}{10°\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}‐\mathrm{T}}_0}{10°\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₀ = 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.

Figure 8

Time series of bursting activity drawn by using the data for the experimentsA and computer simulations at different temperatures. Experimentallyrecorded spike trains (blue, dotted line) and the result from computersimulations of modeling (red, solid line) are compared. The results fromsimulation using the equations developed by Canavier et al. [[17]] and modified by us with temperature-dependent scaling factors atthe temperatures, 18.1°C, 22.1°C and 27.8 are shown on panels L1, M1,and H1 respectively, in the upper traces (A). Those using the modelsrevised by Butera et al. [[19]] and Bertram [[18]] and modified by us with the same temperature dependent scalingfactors are shown in the intermediate traces (B, C). The comparison of those using the equations revised by Bertram andmodified by us are shown in the lowest panels (D) at the temperatures22.2°C, 23.4°C and 25.8°C with small range of temperaturechanges.

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.

Figure 9

Comparison of temperature-dependent subthreshold current activities. (A)Time series of temperature-dependent subthreshold current activities(ρ(T)∙I_SI) simulated by using the results based on themodel developed by Canavier et al. [[17]] at temperature values 18°C (L, blue lines)), 23°C (M,black lines), and 28°C (H, red lines) are shown together. (B) Those(ρ(T)∙I_SI) obtained from simulation data at the sametemperature values using the model made by Butera et al. [[19]] are shown together. (C) Those(I_D + ρ(T)∙I_NSR) taken fromsimulation data at temperature values 22.3°C, 23°C, and 23.7°Cusing the model made by Bertram et al. [[18]] are shown together. (D) Those(ρ(T)∙(I_Ca + I_K(Ca))) obtainedfrom simulation data at temperature values 18°C, 23°C and 28°Cusing the model made by Plant [[15]] are shown together.

Anyway, these facts might show that it was hard to clarify mathematically thetemperature-dependent changes of bursting patterns at burst-firing neurons ofAplysia 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²⁺-dependent K⁺-current I_K(Ca)with voltage dependent activation of Ca²⁺-current. The temperaturedependences were given by temperature-like scaling factors for the maximumconductances and the time constants, with reference temperatureT₀ = 25°C. Q₁₀ of 3.0 for the activationvariable and Q₁₀ 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²⁺, 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²⁺-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²⁺) injectedinto the cytoplasm and expected to be activated during a burst. Chay [[40]] constructed mathematical model applicable to the Aplysiabursting neurons with involving the properties of Ca²⁺-activatedK⁺ channel: intracellular Ca²⁺ 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₁₅ 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²⁺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₀ = 23°C, Q₁₀ of 3.0 for the activationvariable, and Q₁₀ 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₀, 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

θ₂, 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 χ

ρ⁻¹, 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}‐\mathrm{T}}_0}{{10}^{°\mathrm{C}}}}$

ϕ(T), temperature-dependent scaling factor; $\varphi \left(\mathrm{T}\right)\equiv 3^{\frac{{\mathrm{T}‐\mathrm{T}}_0}{{10}^{°\mathrm{C}}}}$

I_SI, slow inward Ca²⁺-current

I_NS, nonspecific cation current

I_R, anomalous rectifier current

I_NaCa, Na⁺-Ca²⁺ exchanger current

I_NaK, Na⁺-K⁺ pump current

I_CaP, Ca²⁺ 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