1. Introduction

Sensitivities for environmental temperature and body temperature are vital functions for creature as well as human beings. Thermosentive receptors or neurons, which seem to be broadly distributed on the skin and in the brain, are classified into two categories, those that sense warmth and those that sense cold[1],[2]. In particular, cold receptors are usually excited in the presence of a cold stimulus. Its activities paradoxically respond when the temperature is relatively high (for example, above 45℃ in the cold receptor on slowly conducting fibers[3]). Such activities were modelled with a Boltzmann description of voltage-gated membrane channels[4].

Since the Boltzmann model does not differ functionally from standard Hodgkin-Huxley (HH) formulations[5], the Boltzman as well as HH models can be interpreted as models of detection of body temperature in the brain area. In fact, the HH model, which begins repetitive firings with zero frequency via a saddle-node bifurcation occurring on an invariant cycle (called the SNIC), shows the aforementioned spike activities dependent on temperature variations[6]. In[7], a frequency increase with a temperature Copyright © 2012 Scientific & Academic Publishing. All Rights Reservedreducti-on was calculated even with the Morris-Lecar (ML) model[8],[9], which has a common structure with two dimensional excitability models such as the (V, n)-reduced HH model[10]. This indicates that such an increase may be caused by nonlinearity of the main current dynamics. We thus suppose it meaningful to understand mechanisms on such temperature-dependent neuronal properties.

It is also important to examine whether or not neurons sensing body temperature variations exist in the nervous system, for the sake of homeostasis or stabilities in brain functions. The brain functions have been well known to be affected by small changes in brain temperature (as little as 2-3℃)[11-13]. Effects of physiological body temperature on neural activities in the hippocampal slice have been more intensively studied since 2007, those being, regulation of resting membrane potentials[14], increased synchronized network activity[15] and increased frequency of γ-oscillation[16]. These experiments give us one of indications that cold-receptor-like neuronal activities may even be observed in the brain.

The concrete target of this work is to analyze a frequency-current (f-I) curve that is V-shaped with a temperature (T) increase, recorded from pyramid cells in the hippocampal CA3 slices, indicating that the pyramidal cells may functionally act as cold neurons. For this target, the HH model exhibiting the SNIC bifurcation with an I increase is employed with the parameter sets as referred to[17]. We simulate that the HH model behaves like the cold receptor activation as mentioned above, calculating f at different T and different I values.

Next, making phase descriptions for small perturbations of ΔI and ΔT within the framework of the phase reduction method[18-20], we derive the theoretical formulation of frequency gradients for I and T. The frequency gradient formulation gives us an explanation of how the theoretical curve for the HH model should best-fit to the empirical plots. In the most appropriate fittings, a temperature coefficient Q₁₀ in the HH model is rescaled to the value that is in a good agreement with the empirical Q₁₀ calculated by the experimental data. Finally, discussion and conclusion will be given.

2. Methods and Materials

2.1. Slice culture preparation

Entorhinal-hippocampal organotypic slices were prepared from 7-d-old Wistar/ST rats. Rat pups were anesthetized by hypothermia and decapitated[21]. The brains were removed and placed in aerated ice-cold Gey’s balanced salt solution supplemented with 25 mM glucose. Horizontal entorhinal- hippocampal slices were cut at a thickness of 300 μm using a vibratome. The slices were placed on Omnipore membrane filters and incubated in 5% CO₂ at 37℃. The culture medium, which was composed of 50% minimal essential medium, 25% Hanks’ balanced salt solution, 25% horse serum, and antibiotics, was changed every 3.5 d. The slices were cultured for 10-15 days in vitro and used for electrophysiological experiments.

2.2. Electrophysiology

An entorhinal-hippocampal slice was placed in a recording chamber and perfused at 25, 30, and 35℃ at a rate of 3-4 ml/min with artificial cerebrospinal fluid (aCSF) consisting of 127 mM NaCl, 26 mM NaHCO₃, 3.5 mM KCl, 1.24 mM KH₂PO₄, 1.0 mM MgSO₄, 1.8-2.0 mM CaCl₂, 10 mM glucose, 10 μM CNQX, and 50 μM AP5[22]. The temperature was controlled using an in-line heater (SC-20l, Warner instruments). For whole-cell recordings borosilicate glass pipettes (4-6 MO) were filled with solution consisting of 130 mM K-gluconate, 10 mM KCl, 10 mM HEPES, 10 mM phosphocreatine, 4 mM MgATP, 0.3 mM NaGTP. Data were obtained using a MultiClamp 700B amplifier and a Digidata 1440A digitizer controlled by pCLAMP 10 software (Molecular Devices). Signals were low-pass filtered at 2 kHz and digitized at 20 kHz.

Figure 1. Frequency (f Hz) is temperature-dependent. (a) f-I plots at 25, 30 and 35℃. (b) The average frequency over all I steps, , is plotted at each temperature. P†< 0.01, compared with T = 25℃ sampling data, and P* < 0.001 (n = 6), compared with T = 30℃, ANOVA. (c) log(f)-T line results in the temperature coefficient Q10 = 1.22 ± 0.38

Figure 2. Bifurcation diagram for the SNIC HH model when increasing temperature from 10, 25 to 41.6℃ for I = 0.161. The broken lines represent unstable fixed points (FP). One is the saddle while the other is focus. The solid lines represent stable focus and stable limit cycle (LC)

Figure 3. Frequency f versus temperature T on simulation with the HH model. The frequency passed through a minimum at 38℃and then increased at the higher temperature (Red full circles), on which a certain piecewise line (green line) is put to fit experimental data plots to theoretical curves for the frequency gradient, which will be derived below

We used one slice for increased temperature experiments, raising the temperature from 25℃, 30℃ to 35℃ and two slices for a decreased temperature from 35℃, 30℃ to 25℃. Two pyramidal cells picked up in each hippocampal CA3 slice, are for whole-cell clamped records at each temperature. In the whole-cell clamped records, the current I increases from 0 pA to 150 pA with a current step of 5 pA. Record duration is 2 sec. and the latency before and after the recording is 1 sec in total. We have thus found the f-I curve at 25, 30 and 35℃ respectively. The curve at each temperature is drawn with the onset of the repetitive firing of the 0-frequency as the I increases. They can straightforwardly be expected to be the class I f-I curves. Interestingly enough, we have observed as follows: the entire frequency level has decreased at 30℃, despite the highest frequency-curve at 35℃ (Figure 1(a)). We did a significant test of the frequency level at each temperature, in order to confirm that the frequency level at 30℃ is consistently lower than the one at 25℃. As shown in Figure 1(b), we examined P^† < 0.01 (n = 6), compared with that of T = 25℃, and P^* < 0.001 (n = 6), and compared with that of T = 30℃, ANOVA[23]. The temperature coefficient is also calculated as Q₁₀ = 1.22 ± 0.38, obtaining the monotonic increased line based on sampling data for the frequency of the pyramidal cell (see Figure 1(c)).

3. Neuron Model and Frequency Gradient Curves

In Figure 1, we indicate that still-unknown nonlinearities exist in the neuronal firing property. To understand why a firing frequency can be recovered again once it has decreased, and what the mechanism is, the Hodgkin-Huxley (HH) model is employed. The HH model is rewritten with small perturbations of temperature and current, ΔT and ΔI, in the general forms:

(1)

where x = (V, y)^t ∈ R⁴ and ^t denotes the transpose. V is the membrane potential variable and y = (m, n, h)^t is the non- dimensional gating variable. While G(x) represents a perturbation vector for ΔT and ΔI (whose detailed forms will be described below), F(x) is a baseline vector field:

(2)

where C_m (= 1 μF/cm²) is the membrane capacity. E_Na, E_K, and E_L are the reversal potentials of Na⁺, K⁺ and leak currents respectively, while g_Na, g_K, and g_L are the conductances. μ(T) = (Q₁₀)^(T-Te)/10 (Q₁₀ = 3 in this work) rescales time courses of the ion-channel activations with temperature T (^◦C) and T_e (= 25^◦C in this work). α_m(V) = -0.1(V + V₁) / (exp(-0.1(V + V₁)) -1), β_m(V) = 4exp(-(V + V₂) / 18), α_h(V) = 0.07exp(-(V + V₃) / 20), β_h(V) = 1 / (exp(-0.1(V + V₄)) + 1), α_n(V) = 0.01(V + V₅)/(exp(-0.1(V + V₅)) - 1) and β_n(V) = 0.125exp(- (V + V₆) / 18). We suppose that g_Na = 35 mS/cm², E_Na = 55 mV, g_K = 9 mS/cm², E_K = -90 mV, g_L = 0.1 mS/cm², E_L = -65 mV, V₁ = 35, V₂ = 60, V₃ = 58, V₄ = 28, V₅ = 34 and V₆ = 44 (see[17]) so that the HH model represents an oscillatory system exhibiting repetitive firings via the SNIC bifurcation mechanism (sometimes called the SNIC HH model). It is noticed that the firing frequency gradually decreases and increases again, correspondingly the amplitude of the membrane potential becomes smaller at higher temperature T = 41.6℃. Figure 2 shows the emergence of a small limit cycle (LC) attractor via the SNIC bifurcation around I = 0.161. The LC size varies with an increase in temperature. The bi-stability of the stable equilibrium and the stable LC attractor exists at T = 10℃ (Analysis on the SNIC bifurcation, referred to in Appendix A). In the frequency-current (f-I) curve for T = 10℃, the HH model shows class II excitable, at the same time, class I spiking behaviour. Here the class II means the emergence of oscillation with a nonzero frequency while the class I is a system that begins oscillations with 0-frequency[24],[25]. The SNIC also shows a system that is simultaneously class I excitable and class I spiking at T = 25℃. The SNIC at T = 41.6℃ shows again a system that simultaneously class II excitable and class I spiking. As the results of the f-I curve, we can depict the frequency-temperature (f-T) curve as shown in Figure 3, at I = 0.161. In an increase of T, the frequency decreases to the minimum at 38 ^◦C through its peak at 13℃. Immediately afterwards, the frequency slope becomes extraordinary sharp. The f-T curve plays a crucial role in fitting to the V-shaped f-T plot recorded from the hippocampal CA3 pyramid cells. The detailed best-fitting method will be mentioned in the next section.

In this work, we select the perturbation G(x) as

(3)

Referring to the phase reduction method[Appendix B], we make a phase description with respect to small perturbation (ΔI, ΔT)^t:

(4)

where f = f(I₀, T₀) =[1/T_p].

Next, the frequency at (I₀, T₀) can be expanded in (ΔI, ΔT). The quadratic Taylor series is obtained

(5)

(4) and (5) gives us the relational expressions

(6)

(7)

where ρ=log(Q₁₀)/10, but in this work, let ρ be 1. We then obtain

(8)

The normalization condition described in Appendix B is used. (8) is thus represented as

(9)

Since[∂f / ∂T] in the H(I, T) does not depend on any width of temperature step, (9) is considered to rewritten as:

(10)

where

Let experimental (f, I, T) data be applied into the H-function. ΔT is the temperature difference between two of 25, 30 and 35℃, namely, ΔT = ± (30 - 25), ± (35 - 30) and ± (35 - 25). ΔI is ±50 pA, which represents the difference of one sampling current to its nearest neighbours.

4. Numerical Calculations

Next, we calculate the H-function for experimental data on (f, I, T) (called the empirical plot). The H-function is also given by the simulation data of (f, I, T) with the HH neuron (called the theoretical curve). We here investigate how such theoretical curves can fit to the empirical plots. Under the most appropriate condition for the best fitting, the temperature coefficient in the physiological experiment is predicted by rescaling the model temperature in the HH model to the empirical temperature.

For the empirical plot, the H(I) is calculated for ΔT = ± (30 - 25), ± (35 - 30) and ± (35 - 25). Figures 4(a) and 4(b) represent respectively the H plots for ΔT > 0 and ΔT < 0 as the functions of I. In analogy, we obtain theoretical curves on simulations with the HH neuron model (solid lines in Figure 4). For this, T_i (i = 0, 1, 2) is sampled on the f-T curve at each current. T₀ and T₂ are relatively high temperatures while T₁ is close to the minimum (see Figure 3). We then select the most appropriate value of T_i (T₀ = 15, T₁ = 40 and T₂ = 41.7) that the theoretical curve passes through the respective plot distribution.

In order to confirm that the empirical plot best-fits to the corresponding theoretical curve, we compute confidence that the theoretical curve can be shared with the empirical plot distribution, under an assumption to obey the normal distribution. This confidence is computed as referred to Appendix C. The standard error of the mean is calculated as H₁ ± 1.96σ /√n where σ is the standard deviation of the sampling distribution and n = 6. In Figure 5, almost all samplings in plot distribution are of confidence as the corresponding theoretical curve.

Finally, we will give additional supporting information to the best fitting by rescaling the temperature coefficient Q₁₀ in the HH neuron. For this, we rescale the temperature coefficient Q₁₀ = 3 in the HH neuron to be in the empirical range of 1.22 ± 0.38 that is calculated in Figure 1(c). In Figures 4 and 5, the temperature ranges for obtaining empirical plots are respectively 5℃ and 10℃. The temperature ranges for obtaining theoretical curves are correspondingly 25℃ and 26.7℃. The rescaled temperature coefficient is thus calculated in the range of 1.24 to 1.51. The range value presumably is satisfied with Q₁₀ = 1.22 ± 0.38 in the experiment.

Figure 4. H-I plots for experimental data and its curves for the HH model. (a) ΔT > 0. (b) ΔT < 0

Figure 5. Theoretical curves best-fit the empirical plots. (a) and (b) are ΔT > 0. (c) and (d) are ΔT < 0

5. Discussion and Conclusions

In this work, we have shown that V-shaped frequency -current responses, recorded from the hippocampal CA3 pyramidal cells, are analytically consistent with the ones simulated with the HH model. However, we cannot declare that the pyramidal cell of the hippocampal CA3 behaves like activations of cold-receptor-like neurons. Rather, in order to ascertain this, we have to examine more in vitro frequency versus current relationships with fine temperature variations of 1-2^◦C and then observe the frequency -temperature curve as shown in Figure 3. Also, we have to give physiological reasons of why the frequency increases at 30 - 35 ^◦C. Such a frequency increase may result from activations of transient receptor potential (TRP) family of channels such as TRPV3[26-28] and TRPV4[14], which seem to be widely distributed in a brain, are sensitive to body temperature. Both the TRP channels have the temperature threshold of 30 - 35 ^◦C for activation, and particularly the TRPV4 channels exist even in the neuronal circuit in the hippocampus.

In the theoretical aspects, we have to figure out detailed transition mechanisms of low to high frequencies at T = 41 ^◦C. To figure out such transition mechanisms, we analyze the difference of the f-I curves at high and low temperatures. The two dimensional reduced Hindmarsh-Rose model (2DHR)[29] is suitable for clarifying the detailed transition mechanisms or the difference of the two f-I curve, in terms of (9) and (10). The 2DHR model has less complexity than the HH model and retains the analytic tractability of the FitzHugh- Nagumo model[30]. All bifurcations of the Hopf, the SNIC and the SN are computed with appropriate parameter sets in the 2DHR model. Drawing the whole diagram of these bifurcations, the corresponding[∂f / ∂T]-I curve will also be calculated. We thus expect to find differences of the two f-I curves with the same class if the[∂f / ∂T]-I curves have obvious differences.

Figure 1(b) showed a decreased frequency with the temperature increase in individual neurons. However in Figure 1(c), the statistic property of the frequency was the monotonic increase, although it was very slight. This indicates the increased activities of the assemblies, in which the underlying individual cell paradoxically exhibits the frequency descent with the temperature increase. Such a paradoxical mechanism should be clarified. One of considerable solutions may be to confirm whether or not pyramidal cells in the hippocampal CA3 are sensitive to current fluctuations on the f-I relationships, because it is known that f-I responses are significantly influenced by the Gaussian distributed current[31]. In addition, under a statistical assumption of empirical plots obeying the Gaussian distribution, the linearity of the log(f)-T is obtained as shown in Figure 1(c). This implies that the Gaussian distributed current contributes to finding a positive slope of the log(f)-T.

In order to confirm the statistical assumption, we should do simulations on the HH models with the current fluctuated by Gaussian distribution, analogous to the work of[31] for investigating effects of the variation on f-I responses. Also, the H-I curve for the HH model obeying the Gaussian distribution will be calculated for finding the most suitable variance for best-fit to the empirical plots. Doing such studies becomes more valuable, in comparison with the other cases such as the Gamma, Poisson and exponential interspike interval distributions. This is because we have to show which distribution is of significant essence for understanding the mechanisms of neural information processing[32].

We have to discuss if it is necessary to visualize even spike patterns of pyramidal cells in the hippocampal CA3 slice, because it is shown that the spike patterns of neocortical layer 2/3 and 5/6 pyramidal cells change with temperature. The change in spike pattern is due to temperature-dependent changes in the intrinsic properties of the pyramidal cells[33],[34]. To understand such physiological mechanisms on the spike pattern generations of regular spikes, slow frequency adaptation and bursting, we will have to study and model the complex dynamics of various types of the ion channels, referring to[25],[35]. This gives arise to an additional issue if we understand rigorously the detailed dynamics of the ion channels for the sake of the modeling.

To the additional issue, we may answer as follows: spike pattern observed in the recording duration as well as detailed channel structures on the current dynamics is not so important in fitting analysis on the temperature controlled frequency-current curve. In our theory, we assume a stable periodic solution for unperturbed firing dynamics, to derive the phase equation and the frequency gradient equation. The stable periodic solution is regardless of any spike pattern such as regular spiking and bursting. Moreover, we can assume the recording duration of 2 sec in the experiment as the one periodic solution. Such assumptions enable us to obtain the frequency gradient equation corresponding to the experimental data.

In conclusion, with respect to the temperature sensitive activities on neurons, the fact that the frequency increases with temperature reduction and paradoxical response when the temperature is increased can be regarded to be of importance for information processing in the brain. We then predict that such temperature sensitive activities are potentially observed even in the hippocampal CA3 pyramidal cells. To support this prediction, we propose an experimental data fitting method of a frequency gradient curve for current and temperature simulated with the SNIC HH oscillator. In this method, the theoretical curves were established by expanding the HH oscillator with small perturbations of I and T while the empirical plot was also obtained with f-I responses recorded from the hippocampal CA3 pyramidal cells at different temperatures. The empirical plot was best fitted to the theoretical curve under a certain temperature range on simulations with the HH oscillator.

The next task in the future is to confirm that the experimental data analysis proposed in this work is a useful and powerful method. For this, we will have to study more the experimental data fitting method by using the ML model, the 2DHR model and the bursting model, which undergo the SNIC bifurcation mechanism with or without slow spike frequency adaptation. In addition, we suppose it very interesting to simulate the HH model with Gaussian noise for getting the higher reproducibility of the empirical plot. This is because we can easily expect that theoretical curves computed with the Gaussian distributions cover the empirical data plot when the noise intensity is larger.

ACKNOWLEDGEMENTS

The authors thank Natsume at Kyushu Institute of Technology and Aihara at The University of Tokyo for fruitful and active discussions. Y.I. was supported by the Funding Program for Next Generation World-Leading Researchers (no. LS023).

Appendix A. Analysis on Stabilities of Stationary Solutions

The stability analysis of the differential equations was briefly reviewed. The detailed analysis is referred to[36],[37]. We analyzed the character of the equilibrium points x₀ of[dx/dτ] = 0. This was done by linearizing the HH model without the perturbation term, close to x₀ and by solving the characteristic equation

(A.1)

where I denotes the identity matrix, and

. (A.2)

All elements, which are functions of the equilibrium solutions, contain the T- and I-parameters. λ is the eigenvalue of the Jacobian matrix and is of principal relevance for determining the stability of the equilibrium points. It is given as a fourth order equation:

(A.3)

a₀ is the determinant. a₁ is a sum of products of the three elements from different rows and columns. a₂ is a sum of products of the two elements. a₃ is a sum of the diagonal elements. We investigated the cases of node and saddle points[Table A.1].

Table A.1. Stability of equilibrium points and the relation to the eigenvalues of the Jacobian matrix DxF character λ0 λ1 λ2 λ3 Stable node Re < 0 Re < 0 Re < 0 Re < 0 Unstable node Re > 0 Re > 0 Re > 0 Re > 0 Saddle point Re > 0 Re < 0 Re < 0 Re < 0

Appendix B. Phase Reduction Method

The detailed derivation is referring to[18],[19]. Supposing a small perturbation in the HH oscillator, the phase reduction method should be applied to:

(B.1)

Let x_p(τ) denote a stable T_p-periodic solution for the unperturbed dynamics,

(B.2)

where x_p(τ + T_p) = x_p(τ). Then a stable solution for (B.1) is approximated as

(B.3)

where η(τ) means a small perturbation in the phase direction on the periodic orbit. u(τ + η(τ)) denotes the orbital deviation to the periodic orbit x_p(τ). Substituting (B.3) into (B.1) and expanding the both-hand sides into a Taylor series

(B.4)

D_x =[∂/∂x] and D_xF is the Jacobian matrix. In the perturbed oscillator, the orbital deviation vector u(τ) = x(τ) - x_p(τ) evolves as

(B.5)

The vector Z(τ) tangent to the periodic orbit x_p(τ), which is the unique solution to

(B.6)

where the normalization condition

(B.7)

is satisfied for every τ. (B.1) is reduced to the evolution equation for η:

(B.8)

Introducing phase variables defined by θ = (τ+η) / T_p (∈[0, 1)), and then averaging to obtain the closed form of the evolution equation for θ(τ), (B.8) is rewritten as

(B.9)

Appendix C. Averaged Frequency gradient equations

Let the sampling index be i (= 0, ..., n - 1), that is, the number of pyramid cells employed for recording in the experiment, referring to[23]. The ith neuron’s equation for frequency gradient is

(C.1)

The sampling average of the above equation is

(C.2)

where. Then, both the sides are divided over averaged frequency:

(C.3)

(C.2) gives us the variance

(C.4)

where σ is the standard deviation. The standard error of (C.3) is rewritten as σ / (f √n).