Extracting thermodynamic properties from van ’t Hoff plots with emphasis on temperature-sensing ion channels

ABSTRACT

Transient receptor potential (TRP) ion channels are among the most well-studied classes of temperature-sensing molecules. Yet, the molecular mechanism and thermodynamic basis for the temperature sensitivity of TRP channels remains to this day poorly understood. One hypothesis is that the temperature-sensing mechanism can simply be described by a difference in heat capacity between the closed and open channel states. While such a two-state model may be simplistic it nonetheless has descriptive value, in the sense that it can be used to compare overall temperature sensitivity between different channels and mutants. Here, we introduce a mathematical framework based on the two-state model to reliably extract temperature-dependent thermodynamic potentials and heat capacities from measurements of equilibrium constants at different temperatures. Our framework is implemented in an open-source data analysis package that provides a straightforward way to fit both linear and nonlinear van ’t Hoff plots, thus avoiding some of the previous, potentially erroneous, assumptions when extracting thermodynamic variables from TRP channel electrophysiology data.

KEYWORDS:

• TRP channels
• temperature sensors
• van ’t Hoff plots
• thermodynamics
• spline fitting

Introduction

An organism’s ability to sense its environment is crucial to its survival. One of the most well-studied families of biological temperature sensors in humans and other eukaryotes is the transient receptor potential (TRP) family of ion channels [1]. Members of this family have temperature sensitivity across the biologically relevant range of temperatures, but the most well-known are the heat and capsaicin-sensitive TRPV1 [2] in the TRPV subfamily, and the cold and menthol sensitive TRPM8 [3,4] in the TRPM subfamily. Hypotheses about the principles guiding the temperature-sensitivity of TRP channels were already being postulated within a few years of their discovery, with proposed mechanisms relating to phenomena from voltage-sensing to elongations of open channel burst times [5,6]. However, while these molecules have been identified as intrinsically sensitive to temperature [7] and playing a critical role as temperature sensors in our nervous system [8–10], we still do not understand the molecular and thermodynamic mechanism(s) that dictates their temperature-dependent activation.

One characteristic of TRP ion channels that seems clear are the large positive enthalpy differences between states for heat-sensitive TRPs like TRPV1 [5,11] and large negative enthalpy differences for cold-sensitive TRPs like TRPM8 [11,12]. Entropy and enthalpy differences between the open and closed states of a channel can be extracted from linear fits to the logarithm of the equilibrium constant $K_{\mathrm{eq}}$ as a function of the reciprocal temperature $1/\mathit{T}$ if said thermodynamic potentials are independent of temperature. However, as is well known in the literature of physical biochemistry, large conformational changes in proteins are usually accompanied by changes in their heat capacities, which leads to temperature-dependent enthalpies and entropies [13]. This is the premise of a model-free framework proposed by Clapham and Miller [11], which can explain both cold and heat-sensitive changes in the equilibrium constant. Recently [14], it was demonstrated that a generalized version of the framework by Clapham and Miller is able to capture a broad range of experimentally-consistent channel behaviors, which were previously only attributed to elaborate multi-state models. While it remains debated whether temperature-dependent gating in channels is also accompanied by observable changes in heat capacity [15], this is the main mechanism to induce temperature dependence in the relevant thermodynamic potentials.

Here, we embrace the approach of Clapham and Miller [11], and introduce a procedure to reliably extract temperature-dependent thermodynamic potentials and heat capacities from equilibrium constant measurements performed at different temperatures. We thereby assume that a TRP channel can, to a first approximation, be described as a two-state system, which may not provide the same mechanistic insight as more involved models [16], but has the benefit of being universally applicable and allows for a direct comparison of thermodynamic variables obtained for different ion channels or the same channel at differing experimental conditions. Our theory is implemented in an open-source data analysis package [17] written in Julia [18], and should provide practitioners a straightforward way to fit linear and nonlinear van ’t Hoff plots, thus avoiding previous potentially false assumptions about the nature of temperature sensors. To facilitate the use of our data analysis script, we offer a link on our GitHub page [17] to an interactive cloud environment, where the code can be run without any prior installation. Users can upload their own data and analyze it within this cloud environment, where all uploaded data gets purged as soon as the session is closed.

The paper is structured as follows. At the beginning of the Theory section, we list the thermodynamic relations relevant to our discussion, and briefly review their common use in the literature of temperature sensors. A subsection titled “Spline fitting of discrete data points” introduces cubic splines as continuously differentiable functions used to fit discrete measurements of $ln(K_{\mathrm{eq}})$. Under the assumption of a two-state model, the differentiability of splines allows us to calculate robust estimates for the thermodynamic potentials $\mathrm{\Delta }\mathit{H}(\mathit{T})$ and $\mathrm{\Delta S}(\mathit{T})$, and the associated heat capacity difference $\mathrm{\Delta }{C}_p(\mathit{T})$. To avoid overfitting, we rely on a Bayesian information criterion (BIC) [19] to penalize splines with many degrees of freedom, as described in the subsection “Model selection”. For illustrative purposes, we apply our data analysis package to two distinct data sets in the Results and discussion, and the Conclusions provide a summary of our results.

Theory

One of the simplest ways to model a TRP channel is to treat it as a two-state system. Every channel in an ensemble of channels can then either be in the open or closed state, such that the composition of the ensemble is encoded in the equilibrium constant $K_{\mathrm{eq}}$:

$\mathrm{open}\mathrm{state}\stackrel{{\mathit{K}}_{\mathrm{eq}}}{\rightleftharpoons }\mathrm{closed}\mathrm{state}.$

In electrophysiological experiments the charge current through a single channel or a collection of channels is measured at different temperatures, which can be used to calculate the so-called “open probability” $\mathit{P}$, i.e. the probability of finding a channel in the open state. The equilibrium constant $K_{\mathrm{eq}}$ and the open probability $\mathit{P}$ are related via

(1) $K_{\mathrm{eq}}=\frac{P}{1-\mathit{P}}$(1)

for a two-state system. Note that macroscopic ionic currents are subject to a multitude of additional sources of variability that can predominate at the temperature extremes, where channel activity is either very low or near-maximal. Because the quality of the associated $\mathit{P}$ and $K_{\mathrm{eq}}$-estimates is directly affected, we recommend users to carefully select the temperature range of the data to be fitted to avoid contributions from sources of signal variability unrelated to channel gating.

Thermodynamic description

Thermodynamics tells us, on the one hand, that the differentials of enthalpy $\mathit{H}$ and entropy $\mathit{S}$ are related via

$\mathrm{d}\mathit{H}=\mathit{T}\mathrm{d}\mathit{S}$

for systems at constant pressure, where $\mathit{T}$ denotes the absolute temperature. On the other hand, it can be shown that the heat capacity $C_p$ at constant pressure satisfies

(2) $C_p={\left(\frac{\mathrm{\partial H}}{\mathrm{\partial T}}\right)}_{\mathit{p}},$(2)

where $\mathrm{\partial }/\mathrm{\partial }\mathit{T}$ denotes a partial derivative with respect to $\mathit{T}$ and the index of the bracket reveals which quantity is being held constant (in this case it is the pressure $\mathit{p}$). We therefore conclude that the enthalpy and entropy differences between two metastable states must be integral functions of the heat capacity difference $\mathrm{\Delta }{C}_p=C_p^{\mathrm{s}\mathrm{t}\mathrm{ate}2}-C_p^{\mathrm{s}\mathrm{t}\mathrm{ate}1}$ between the states, i.e.

(3) $\begin{array}{cc}\mathrm{\Delta }\mathit{H}(\mathit{T})& =\mathrm{\Delta H}(T_0)+{\int }_{T_0}^{\mathit{T}}\mathrm{d}\mathit{T}{ }^{\mathrm{\prime }}\mathrm{\Delta }{C}_p(\mathit{T}{ }^{\mathrm{\prime }}),\\ \mathrm{\Delta }\mathit{S}(\mathit{T})& =\mathrm{\Delta S}(T_0)+{\int }_{T_0}^{\mathit{T}}\mathrm{d}\mathit{T}{ }^{\mathrm{\prime }}\frac{\mathrm{\Delta }{C}_p(\mathit{T}{ }^{\mathrm{\prime }})}{\mathit{T}{ }^{\mathrm{\prime }}},\end{array}$(3)

for some arbitrary reference temperature $T_0$. Note that $\mathrm{\Delta }\mathit{H}(\mathit{T})$ and $\mathrm{\Delta S}(\mathit{T})$ are independent of $\mathit{T}$ whenever $\mathrm{\Delta }{C}_p(\mathit{T})=0$, e.g. for bistable systems whose states have the same heat capacity. However, in the case of protein folding, we know that large heat capacity differences exist between their folded and unfolded state [13].

Another thermodynamic potential of interest is the Gibbs free energy, which is given by

(4) $\mathrm{\Delta G}(\mathit{T})=\mathrm{\Delta H}(\mathit{T})-\mathit{T}\Delta \mathit{S}(\mathit{T}).$(4)

It can be related to the equilibrium constant $K_{\mathrm{eq}}$ of the two-state system via the fundamental relation of chemical thermodynamics, namely

(5) $\mathrm{\Delta }\mathit{G}(\mathit{T})=-\mathit{RT}ln(K_{\mathrm{eq}}),$(5)

where $\mathit{R}=8.31446261815324$ J mol⁻¹ K⁻¹ denotes the molar gas constant. The logarithm of the equilibrium constant and its derivative with respect to $\mathit{T}$ therefore have the form

(6) $ln(K_{\mathrm{eq}})=-\frac{\mathrm{\Delta }\mathit{H}(\mathit{T})}{\mathit{RT}}+\frac{\mathrm{\Delta S}(\mathit{T})}{R},$(6)

(7) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta H}(\mathit{T})}{R{T}^2},$(7)

where the latter is the well-known van ’t Hoff equation, which is sometimes also written as follows:

(8) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}(1/\mathit{T})}=-\frac{\mathrm{\Delta H}(\mathit{T})}{R}.$(8)

Equations (6)(6) $ln(K_{\mathrm{eq}})=-\frac{\mathrm{\Delta }\mathit{H}(\mathit{T})}{\mathit{RT}}+\frac{\mathrm{\Delta S}(\mathit{T})}{R},$(6) and (8(8) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}(1/\mathit{T})}=-\frac{\mathrm{\Delta H}(\mathit{T})}{R}.$(8) ) reveal that a so-called van ’t Hoff plot, where $ln(K_{\mathrm{eq}})$ is plotted against the reciprocal of the absolute temperature $\mathit{T}$, will be linear whenever $\mathrm{\Delta H}$ and $\mathrm{\Delta S}$ are constant with respect to temperature. The thermodynamic potentials can then be read off the slope and intercept of $ln(K_{\mathrm{eq}})$, respectively. This convenient fact often seems to guide the decision of practitioners to fit their data to straight lines, even when the van ’t Hoff plot is highly nonlinear (see, e.g. Refs. [12] and [20]), which can be an indication for temperature-dependent behavior. Overall, it is important to note that Eqs. (7)(7) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta H}(\mathit{T})}{R{T}^2},$(7) and (8(8) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}(1/\mathit{T})}=-\frac{\mathrm{\Delta H}(\mathit{T})}{R}.$(8) ) are valid for all functions $\mathrm{\Delta }\mathit{H}(\mathit{T})$ and $\mathrm{\Delta S}(\mathit{T})$ of the form given in Eq. (3)(3) $\begin{array}{cc}\mathrm{\Delta }\mathit{H}(\mathit{T})& =\mathrm{\Delta H}(T_0)+{\int }_{T_0}^{\mathit{T}}\mathrm{d}\mathit{T}{ }^{\mathrm{\prime }}\mathrm{\Delta }{C}_p(\mathit{T}{ }^{\mathrm{\prime }}),\\ \mathrm{\Delta }\mathit{S}(\mathit{T})& =\mathrm{\Delta S}(T_0)+{\int }_{T_0}^{\mathit{T}}\mathrm{d}\mathit{T}{ }^{\mathrm{\prime }}\frac{\mathrm{\Delta }{C}_p(\mathit{T}{ }^{\mathrm{\prime }})}{\mathit{T}{ }^{\mathrm{\prime }}},\end{array}$(3) , and not just constant thermodynamic potentials.

A popular empirical measure of temperature sensitivity is $Q_{10}$, which is used to characterize temperature sensitivity in electrophysiological experiments on TRP channels [21,22]. It is defined as the ratio of $K_{\mathrm{eq}}$ measured at two temperatures that are $10\mathrm{K}$ apart, i.e.

$Q_{10}=\frac{{\mathit{K}}_{\mathrm{eq}}\left(\mathit{T}+10\mathrm{K}\right)}{K_{\mathrm{eq}}\left(\mathit{T}\right)}.$

The reason for its wide-spread use is the fact that on a logarithmic scale it approximately reproduces the van ’t Hoff equation, i.e.

(9) $\begin{array}{cc}\frac{ln\left({\mathit{Q}}_{10}\right)}{10\mathrm{K}}& =\frac{ln\left({\mathit{K}}_{\mathrm{eq}}\left(\mathit{T}+10\mathrm{K}\right)\right)-ln\left(K_{\mathrm{eq}}\left(\mathit{T}\right)\right)}{\left(\mathit{T}+10\mathrm{K}\right)-\mathrm{T}}\\ \mathit{\hspace{0.17em}}& \approx \frac{\mathrm{d}ln\left({\mathit{K}}_{\mathrm{eq}}\right)}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta }H\left(T\right)}{R{T}^2},\end{array}$(9)

and can therefore be used to estimate $\mathrm{\Delta }\mathit{H}(\mathit{T})$. Again, it is common practice to assume that $\mathrm{\Delta H}$ is temperature-independent and $ln(Q_{10})$ is therefore only evaluated at a single temperature $\mathit{T}$, which can lead to arbitrary and skewed results whenever $\mathrm{\Delta H}(\mathit{T})=\mathrm{const}.$ is not satisfied.

Here, we call for a different approach to extract thermodynamic information from measured $K_{\mathrm{eq}}$-values, without any ad-hoc assumptions. Instead of performing a linear fit of $ln(K_{\mathrm{eq}})$ plotted against $T^{-1}$, we propose to fit the data points using piecewise continuous polynomials and extracting the corresponding potentials $\mathrm{\Delta }\mathit{H}(\mathit{T})$, $\mathrm{\Delta S}(\mathit{T})$ and $\mathrm{\Delta G}(\mathit{T})$ via Eqs. (4)(4) $\mathrm{\Delta G}(\mathit{T})=\mathrm{\Delta H}(\mathit{T})-\mathit{T}\Delta \mathit{S}(\mathit{T}).$(4) , (5(5) $\mathrm{\Delta }\mathit{G}(\mathit{T})=-\mathit{RT}ln(K_{\mathrm{eq}}),$(5) ) and (7(7) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta H}(\mathit{T})}{R{T}^2},$(7) ) [or, equivalently, Eq. (8)(8) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}(1/\mathit{T})}=-\frac{\mathrm{\Delta H}(\mathit{T})}{R}.$(8) ], as discussed in the next subsection. In cases, where $\mathrm{\Delta H}(\mathit{T})$ and $\mathrm{\Delta S}(\mathit{T})$ vary strongly with the temperature, this novel approach also allows us to estimate the change in heat capacity $\mathrm{\Delta }{C}_p$.

Spline fitting of discrete data points

At a fixed temperature $\mathit{T}$, a measured value of $ln(K_{\mathrm{eq}})$ can be used to calculate the corresponding free-energy difference $\mathrm{\Delta }\mathit{G}(\mathit{T})$ via Eq. (5)(5) $\mathrm{\Delta }\mathit{G}(\mathit{T})=-\mathit{RT}ln(K_{\mathrm{eq}}),$(5) . The associated potentials $\mathrm{\Delta H}(\mathit{T})$ and $\mathrm{\Delta S}(\mathit{T})$ then follow from Eqs. (7)(7) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta H}(\mathit{T})}{R{T}^2},$(7) and (4(4) $\mathrm{\Delta G}(\mathit{T})=\mathrm{\Delta H}(\mathit{T})-\mathit{T}\Delta \mathit{S}(\mathit{T}).$(4) ), respectively, and $\mathrm{\Delta }{C}_p$ is given by Eq. (2)(2) $C_p={\left(\frac{\mathrm{\partial H}}{\mathrm{\partial T}}\right)}_{\mathit{p}},$(2) . However, Eqs. (2)(2) $C_p={\left(\frac{\mathrm{\partial H}}{\mathrm{\partial T}}\right)}_{\mathit{p}},$(2) and (7(7) $\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta H}(\mathit{T})}{R{T}^2},$(7) ) require us to calculate the first and second derivative of $ln(K_{\mathrm{eq}})$ with respect to $\mathit{T}$, and this can be somewhat tricky for discrete points, because direct numerical differentiation amplifies the noise in the data. To circumvent this problem, we propose to fit a so-called cubic spline [23] to the data, which is a twice continuously differentiable function made up of third-order polynomials. Our choice of a fit function is motivated by the fact that splines are generally preferred over ordinary polynomials when fitting nonlinear data, because the latter have many undesirable properties, e.g. the fit in one region can influence the behavior in other regions, leading to extreme oscillations at high polynomial orders. Ordinary polynomials also cannot fit certain functional forms, such as logarithmic functions and functions with sharp slope changes. The order of the splines should be cubic to guarantee a continuous second derivative, while keeping the number of fit parameters to a minimum, where the former is needed for the proper estimation of $\mathrm{\Delta }{C}_p$.

A cubic spline is a piecewise continuous function

(10) $\mathcal{S}(\mathit{x})=\left(\begin{array}{cc}{\mathcal{S}}_1(\mathit{x})& s_0\le \mathit{x}\le s_1\\ {\mathcal{S}}_2(\mathit{x})& s_1\le \mathit{x}\le s_2\\ ⋮& \mathit{\hspace{0.17em}}\\ {\mathcal{S}}_i(\mathit{x})& s_{\mathit{i}-1}\le \mathit{x}\le s_i\\ ⋮& \mathit{\hspace{0.17em}}\\ {\mathcal{S}}_N(\mathit{x})& s_{\mathit{N}-1}\le \mathit{x}\le s_N\end{array}\right.$(10)

made up of third-order polynomials of the form

${\mathcal{S}}_i(s_{\mathit{i}-1}\le \mathit{x}\le s_i)=\sum _{\mathit{n}=0}^3{a}_n^{(i)}(\mathit{x}-s_{\mathit{i}-1}{)}^n,$

which are joined together in the spline knots $s_i$. It satisfies the continuity conditions

(11) $\begin{array}{cc}{\mathcal{S}}_i(s_i)& ={\mathcal{S}}_{\mathit{i}+1}(s_i),\\ {\mathcal{S}{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{i}}(s_i)& ={\mathcal{S}{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{i}+1}(s_i),\\ {\mathcal{S}{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{i}}(s_i)& ={\mathcal{S}{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{i}+1}(s_i),\end{array}$(11)

where the notation $\mathcal{S}{ }^{\mathrm{\prime }}(\mathit{x})=\mathrm{d}\mathcal{S}(\mathit{x})/\mathrm{d}\mathit{x}$ and $\mathcal{S}{ }^{\mathrm{\prime \prime }}(\mathit{x})={\mathrm{d}}^2\mathcal{S}(\mathit{x})/\mathrm{d}{x}^2$ was introduced to abbreviate the expressions. We also require some appropriate boundary conditions, e.g. the natural boundary conditions

(12) $\mathcal{S}{{\mathit{ }}^{\mathrm{\prime \prime }}}_1(s_0)=\mathcal{S}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{N}}(s_N)=0,$(12)

which are a popular choice, because they reduce the number of free parameters by two and lead to a less divergent behavior of the spline beyond the edge knots $s_0$ and $s_N$. Equations (11)(11) $\begin{array}{cc}{\mathcal{S}}_i(s_i)& ={\mathcal{S}}_{\mathit{i}+1}(s_i),\\ {\mathcal{S}{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{i}}(s_i)& ={\mathcal{S}{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{i}+1}(s_i),\\ {\mathcal{S}{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{i}}(s_i)& ={\mathcal{S}{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{i}+1}(s_i),\end{array}$(11) and (12(12) $\mathcal{S}{{\mathit{ }}^{\mathrm{\prime \prime }}}_1(s_0)=\mathcal{S}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{N}}(s_N)=0,$(12) ) constrain the values of the spline coefficients $a_n^{(i)}$, such that only $\mathit{N}+1$ of them can be varied independently.

In our fitting procedure, the edge knots are held fixed to define a finite interval $[s_0,s_N]$, on which the “inner” knots $\{s_1,\dots ,s_{\mathit{N}-1}\}$ are allowed to vary. We also vary the values $\{\mathcal{S}(s_i){\}}_{i=0}^N$ of the splines at the knots (see Figure 1). The best fit of $\mathcal{S}(\mathit{x})$ to the data $\{x_m,y_m{\}}_{m=0}^M$ minimizes the sum of squared residuals between the data points and the spline, i.e.

Figure 1. Visualizing the principle of spline fitting. When fitting a cubic spline $\mathcal{S}(\mathit{x})$ to some data points $\{x_m,y_m{\}}_{m=0}^M$, we exploit the fact that for every set of knot coordinates $(s_0,\mathcal{S}(s_0)),(s_1,\mathcal{S}(s_1)),\dots ,(s_N,\mathcal{S}(s_N))$ (blue circles) there exists a unique cubic spline (blue solid lines) satisfying the boundary conditions in Eq. (12)(12) $\mathcal{S}{{\mathit{ }}^{\mathrm{\prime \prime }}}_1(s_0)=\mathcal{S}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{N}}(s_N)=0,$(12) . Thus, by varying the knot coordinates (black arrows), we can change the shape of the spline to minimize ${\chi }^2$ in Eq. (13)(13) ${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{y_m-\mathcal{S}(x_m)}{{\sigma }_m}\right)}^2,$(13) . Note that the edge knots at $\mathit{x}=s_0$ and $\mathit{x}=s_N$ can only be varied in $\mathit{y}$-direction, whereas the “inner” knots are allowed to take arbitrary $\mathit{x}$-values within the interval $[s_0,s_N]$.

(13) ${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{y_m-\mathcal{S}(x_m)}{{\sigma }_m}\right)}^2,$(13)

where ${\sigma }_m$ denotes the standard error of $y_m$. Here, we consider $y_m=ln(K_{\mathrm{eq}}(x_m))$ for either a linear ($x_m=T_m$) or reciprocal temperature scale ($x_m=1/T_m$), and set the values of $s_0$ and $s_N$ equal to the lowest and highest values of $x_m$ found in the data set, respectively. The reason why we consider both scales is because one cannot distinguish between a linear and reciprocal temperature dependence for the temperature ranges realized in electrophysiological experiments (see also Figure 2).

Figure 2. Distinguishing between linear and reciprocal functions on physiologically relevant temperature scales is impossible. (a) A linear function $\mathit{f}(\mathit{T})=\mathit{a}+\mathit{bT}$, evaluated at discrete $\mathit{T}$-values and perturbed by small noise (blue circles), is plotted next to a reciprocal noisy function $\mathit{g}(\mathit{T})=\mathit{c}-\mathit{d}/\mathit{T}$ (red squares) on a temperature scale ranging from 0 to 100. Both functions appear linear, because the absolute temperature is not varied by orders of magnitude to reveal the nonlinearity of $\mathit{g}(\mathit{T})$. (b) Same data as in (a) plotted on a reciprocal temperature scale. Again, both functions seem linear in $1/\mathit{T}$, although only $\mathit{g}(\mathit{T})$ truly is. Insets: same data as in (a) and (b) plotted on a wider temperature scale to visualize the linear and reciprocal trends of $\mathit{f}(\mathit{T})$ (blue solid line) and $\mathit{g}(\mathit{T})$ (red dashed line), respectively.

Evaluating $\mathcal{S}(\mathit{x}(\mathit{T}))$ for the parameters that minimize Eq. (13)(13) ${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{y_m-\mathcal{S}(x_m)}{{\sigma }_m}\right)}^2,$(13) therefore gives the best estimate of $ln(K_{\mathrm{eq}}(\mathit{T}))$, which can be used to extract the heat capacity difference $\mathrm{\Delta }{C}_p(\mathit{T})$ and the thermodynamic potentials $\mathrm{\Delta G}(\mathit{T})$, $\mathrm{\Delta H}(\mathit{T})$, and $\mathrm{\Delta S}(\mathit{T})$ as follows:

(14) $\mathrm{\Delta G}(\mathit{T})=-\mathit{RT}\mathcal{S}(\mathit{x}(\mathit{T})),$(14)

(15) $\mathrm{\Delta H}(\mathit{T})=\mathit{R}{T}^2\frac{\mathrm{d}x(T)}{\mathrm{d}T}\mathcal{S}{ }^{\mathrm{\prime }}(\mathit{x}(\mathit{T})),$(15)

(16) $\mathrm{\Delta S}(\mathit{T})=\frac{\mathrm{\Delta H}(\mathit{T})-\mathrm{\Delta G}(\mathit{T})}{T},$(16)

$\mathrm{\Delta }{C}_p(\mathit{T})=\mathit{R}\frac{\mathrm{d}}{\mathrm{d}T}\left[T^2\frac{\mathrm{d}x(T)}{\mathrm{d}T}\right]\mathcal{S}{ }^{\mathrm{\prime }}(\mathit{x}(\mathit{T}))$

(17) $+\mathit{R}{T}^2{\left[\frac{\mathrm{d}x(T)}{\mathrm{d}T}\right]}^2\mathcal{S}{ }^{\mathrm{\prime \prime }}(\mathit{x}(\mathit{T})).$(17)

Note that for $\mathit{x}=\mathit{T}$ we have $\mathrm{d}\mathit{x}/\mathrm{d}\mathit{T}=1$, whereas the reciprocal relation $\mathit{x}=1/\mathit{T}$ gives $\mathrm{d}\mathit{x}/\mathrm{d}\mathit{T}=-T^{-2}$ for which the first term of Eq. (17)(17) $+\mathit{R}{T}^2{\left[\frac{\mathrm{d}x(T)}{\mathrm{d}T}\right]}^2\mathcal{S}{ }^{\mathrm{\prime \prime }}(\mathit{x}(\mathit{T})).$(17) vanishes.

Model selection

The choice between a linear and a reciprocal fit, as well as the number of spline knots $\mathit{N}+1$, gives rise to a multitude of models that fit the data set to varying degree. For a fixed $\mathit{N}$, one can distinguish between the qualities of a linear and a reciprocal fit by comparing their corresponding ${\chi }^2$ values, but if $\mathit{N}$ is allowed to vary then models with $\mathit{N}\gg 1$ will always be preferred. We therefore propose the use of an information criterion [19] to penalize models with too many fit parameters. By interpreting Eq. (13)(13) ${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{y_m-\mathcal{S}(x_m)}{{\sigma }_m}\right)}^2,$(13) as a negative log-likelihood for Gaussian distributed residuals, we obtain the following BIC:

(18) $\mathrm{BIC}\left(\mathit{N}\right)={\chi }_{min}^2+2\mathit{N}ln\left(\mathit{M}+1\right),$(18)

where $\mathit{M}$ is the number of data points as defined in the previous section. Equation (18)(18) $\mathrm{BIC}\left(\mathit{N}\right)={\chi }_{min}^2+2\mathit{N}ln\left(\mathit{M}+1\right),$(18) is evaluated using the optimal values for the $2\mathit{N}$ spline parameters $\{s_i{\}}_{i=1}^{N-1}$ and $\{\mathcal{S}(s_i){\}}_{i=0}^N$ that minimize ${\chi }^2$ [Eq. (13)(13) ${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{y_m-\mathcal{S}(x_m)}{{\sigma }_m}\right)}^2,$(13) ], resulting in the minimum value ${\chi }_{min}^2$. The model that best fits the data, while avoiding overfitting, minimizes $\mathrm{BIC}(\mathit{N})$ with respect to $\mathit{N}$.

Our data analysis package automatically varies $\mathit{N}$, finds the associated optimal parameter values that minimize ${\chi }^2$, and subsequently calculates the corresponding BIC-value. It finally returns the model and associated parameter values that best fit the data at hand.

Extension to multi-state models

In principle, our data-fitting approach can be extended to (and used to generalize) models with multiple states, such as the ones presented in Ref. [16], by replacing Eq. (1)(1) $K_{\mathrm{eq}}=\frac{P}{1-\mathit{P}}$(1) with an expression $\mathit{P}(\mathit{T})$ for the open probability involving multiple spline functions ${\mathcal{S}}^{(\mathit{k})}(\mathit{x}(\mathit{T}))$ with $\mathit{k}=1,2,\dots$. The fitting must then be performed on the level of $\mathit{P}$, instead of $ln(K_{\mathrm{eq}})$, which implies that Eq. (13)(13) ${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{y_m-\mathcal{S}(x_m)}{{\sigma }_m}\right)}^2,$(13) gets replaced with

${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{{\bar{P}}_{\mathit{m}}-\mathit{P}(T_m)}{{\tilde{\sigma }}_{\mathit{m}}}\right)}^2,$

where ${\bar{P}}_{\mathit{m}}$ and ${\tilde{\sigma }}_{\mathit{m}}$ are the mean and associated standard error of the measured open probability at temperature $T_m$. The corresponding BIC takes the form

$\mathrm{BIC}(\{N^{(\mathit{k})}\})={\chi }_{min}^2+2\sum _k{N}^{(\mathit{k})}ln(\mathit{M}+1),$

where $N^{(\mathit{k})}+1$ is the number of spline knots in ${\mathcal{S}}^{(\mathit{k})}(\mathit{x})$.

For a concrete example, consider the four-state model in Ref. [16], where the open probability is given by

$\mathit{P}(\mathit{T})=\frac{\mathit{L}(1+\mathit{CJ})}{1+\mathit{J}+\mathit{L}(1+\mathit{CJ})}$

and the coefficients $\mathit{L}$, $\mathit{C}$, and $\mathit{J}$ are related to the equilibrium constants between the two open (“O”) and two closed (“C”) states as follows:

$K_{\mathrm{eq}}^{{\mathrm{O}}_0{\mathrm{C}}_0}=\mathit{L},K_{\mathrm{eq}}^{{\mathrm{C}}_0{\mathrm{C}}_1}=\mathit{J},$

$K_{\mathrm{eq}}^{{\mathrm{C}}_1{\mathrm{O}}_1}=\mathit{CL},K_{\mathrm{eq}}^{{\mathrm{O}}_0{\mathrm{O}}_1}=\mathit{CJ}.$

This model can be generalized by replacing $\mathit{L}$, $\mathit{C}$, and $\mathit{J}$ with $exp({\mathcal{S}}^{(\mathit{k})}(\mathit{x})){|}_{\mathit{k}=1,2,3}$, respectively, if all equilibrium constants are assumed to be temperature dependent. After fitting the data in analogy to the two-state case, Eqs. (14)(14) $\mathrm{\Delta G}(\mathit{T})=-\mathit{RT}\mathcal{S}(\mathit{x}(\mathit{T})),$(14) to (17(17) $+\mathit{R}{T}^2{\left[\frac{\mathrm{d}x(T)}{\mathrm{d}T}\right]}^2\mathcal{S}{ }^{\mathrm{\prime \prime }}(\mathit{x}(\mathit{T})).$(17) ) can then be evaluated by replacing $\mathcal{S}(\mathit{x})$ with ${\mathcal{S}}^{(1)}(\mathit{x})$, ${\mathcal{S}}^{(3)}(\mathit{x})$, ${\mathcal{S}}^{(1)}(\mathit{x})+{\mathcal{S}}^{(2)}(\mathit{x})$, and ${\mathcal{S}}^{(2)}(\mathit{x})+{\mathcal{S}}^{(3)}(\mathit{x})$ to extract the thermodynamic potentials and heat capacity differences related to all the different equilibrium constants.

Even though the generalization to multiple states is fairly straight-forward, our data analysis package currently only supports a two-state description.

Results and discussion

For illustrative purposes, we applied the data analysis package to two previously published data sets, one for the warm-sensitive TRPV3 channel [20], and another for the heat and capsaicin-sensitive TRPV1 channel [5]. Each data set was analyzed by performing a van ’t Hoff fit of measured values of $ln(K_{\mathrm{eq}})$ for different temperatures $\mathit{T}$ to extract heat capacity differences and thermodynamic potentials, as described in the subsections “Spline fitting of discrete data points” and ”Model selection” above. Here, we deliberately avoid a direct comparison with the results of the associated publications, as it is not our intention to question their conclusions, but to demonstrate how our data analysis package works in practice.

In the case of the TRPV3 channel, we considered the measured open probabilities $\mathit{P}$ that are tabulated in the source data associated with the extended data Figure 1 in Ref. [20]. For each temperature $T_m$, we calculated the sample mean and variance of $P_m$, i.e.

(19) ${\bar{P}}_{\mathit{m}}=\frac{1}{N_m}\sum _{\mathit{n}=1}^{N_m}{P}_m^{(n)},$(19)

(20) ${\tilde{\sigma }}_{\mathit{m}}^2=\frac{1}{N_m-1}\sum _{\mathit{n}=1}^{N_m}(P_m^{(n)}-{\bar{P}}_{\mathit{m}}{)}^2,$(20)

where $P_m^{(n)}$ denotes the $\mathit{n}$th measurement (of $N_m$ in total) of the open probability $P_m$ at temperature $T_m$. The equilibrium constant $K_{\mathrm{eq}}$ can be calculated via Eq. (1)(1) $K_{\mathrm{eq}}=\frac{P}{1-\mathit{P}}$(1) and according to the variance formula of error propagation one has

$\mathrm{var}(ln(K_{\mathrm{eq}}))\approx {\left(\frac{\mathrm{d}ln({\mathit{K}}_{\mathrm{eq}})}{\mathrm{d}\mathit{P}}\right)}^2\mathrm{var}(\mathit{P}).$

The data points and standard errors entering Eq. (13)(13) ${\chi }^2=\sum _{\mathit{m}=0}^M{\left(\frac{y_m-\mathcal{S}(x_m)}{{\sigma }_m}\right)}^2,$(13) are therefore given by

(21) $y_m=ln\left(\frac{{\bar{P}}_{\mathit{m}}}{1-{\bar{P}}_{\mathit{m}}}\right),$(21)

(22) ${\sigma }_m=\sqrt{\frac{1}{N_m}{\left(\frac{1}{{\bar{P}}_{\mathit{m}}-{\bar{P}}_{\mathit{m}}^2}\right)}^2{\tilde{\sigma }}_{\mathit{m}}^2}.$(22)

Our results for the TRPV3 data are shown in Figure 3. The model that best fits the data is reciprocal in the temperature ($\mathit{x}=1/\mathit{T}$) and contains $\mathit{N}+1=3$ spline knots. The model predicts a temperature-dependent heat capacity difference $\mathrm{\Delta }{C}_p(\mathit{T})$ that decreases monotonically beyond $\mathit{T}\approx 300$ [Figure 3(c)]. The resulting enthalpy and entropy-temperature product differences, $\mathrm{\Delta H}(\mathit{T})$ and $\mathit{T}\Delta \mathit{S}(\mathit{T})$, are therefore nonconstant and vary between $-100$ and $+300$ [Figure 3(b)]. However, they mostly cancel each other out and give rise to a moderate free-energy difference, as seen in Figure 3(a).

Figure 3. Thermodynamic potentials predicted from a van ’t Hoff fit of TRPV3 channel data. The data are best fitted by a model that is reciprocal in the temperature with 3 spline knots. (a) Gibbs free energy $\mathrm{\Delta G}$ as a function of $\mathit{T}$, calculated from data (points) and compared to model prediction (solid line). Shaded areas (gray) mark temperature intervals, where the trend of the cubic spline is no longer constrained by data points and can therefore not be trusted. (b) Enthalpy difference $\mathrm{\Delta H}$ (blue solid line) and entropy-temperature product $\mathit{T}\Delta \mathit{S}$ (red dashed line) as functions of $\mathit{T}$. (c) Heat capacity difference $\mathrm{\Delta }{C}_p$ as function of $\mathit{T}$. The temperature-dependence of $\mathrm{\Delta H}$ and $\mathrm{\Delta S}$ emerges from a nonzero $\mathrm{\Delta }{C}_p$ predicted by the best-fitting model.

Figure 4 displays our results for the TRPV1 channel, where the data points $\{T_m,y_m\}$ and standard errors ${\sigma }_m$ were read off Figure 2c of the original publication. The model that best fits the data is reciprocal in the temperature ($\mathit{x}=1/\mathit{T}$) and has no inner knots, i.e. $\mathit{N}+1=2$. At first this may seem somewhat surprising, considering the fact that the data are not perfectly linear in $T^{-1}$, but is essentially a good example of how our data analysis package avoids overfitting. Apparently, one does not gain sufficiently large improvements in the ${\chi }_{min}^2$ term of Eq. (18)(18) $\mathrm{BIC}\left(\mathit{N}\right)={\chi }_{min}^2+2\mathit{N}ln\left(\mathit{M}+1\right),$(18) to warrant a more complex model than one with $\mathrm{\Delta }{C}_p=0$ and therefore constant thermodynamic potentials $\mathrm{\Delta }\mathit{H}$ and $\mathrm{\Delta S}$. While we are here only illustrating the use of our data analysis package and want to refrain from making scientific assessment of the results at this time, we would like to make note of the narrow temperature range in this particular data set, so as the readers do not conclude that we are definitively claiming that TRPV1 has a vanishing $\mathrm{\Delta }{C}_p$ across the physiological temperature range.

Figure 4. Thermodynamic potentials predicted from a van ’t Hoff fit of TRPV1 channel data. The data are best fitted by a reciprocal model with 2 spline knots. (a-c) Same as in Figure 3, except that here a model is preferred with vanishing heat capacity difference $\mathrm{\Delta }{C}_p$, which results in $\mathrm{\Delta }\mathit{H},\mathrm{\Delta S}=\mathrm{const}.$.

We can now compare our results to the output of alternative data analysis methods, such as the thermal coefficient $Q_{10}$. In Figure 5 we plot the enthalpies of Figures 3(b) and 4(b) next to predictions that arise when Eq. (9)(9) $\begin{array}{cc}\frac{ln\left({\mathit{Q}}_{10}\right)}{10\mathrm{K}}& =\frac{ln\left({\mathit{K}}_{\mathrm{eq}}\left(\mathit{T}+10\mathrm{K}\right)\right)-ln\left(K_{\mathrm{eq}}\left(\mathit{T}\right)\right)}{\left(\mathit{T}+10\mathrm{K}\right)-\mathrm{T}}\\ \mathit{\hspace{0.17em}}& \approx \frac{\mathrm{d}ln\left({\mathit{K}}_{\mathrm{eq}}\right)}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta }H\left(T\right)}{R{T}^2},\end{array}$(9) is solved for $\mathrm{\Delta }\mathit{H}(\mathit{T})$. The latter was evaluated using $ln(Q_{10})=\mathcal{S}(\mathit{T}+10)-\mathcal{S}(\mathit{T})$, where $\mathcal{S}(\mathit{T})\equiv \mathcal{S}(\mathit{x}=1/\mathit{T})$ because both data sets were fitted via reciprocal models. Figure 5(b) demonstrates that $Q_{10}$ gives a decent estimate for the enthalpy whenever $\mathrm{\Delta H}$ is independent of temperature. If this is not the case [Figure 5(a)], then the differences can become arbitrarily large, as can be seen in Figure 5(a). Note that the discrepancy between the spline-fitting estimate and the $Q_{10}$-estimate for $\mathrm{\Delta H}(\mathit{T})$ vanishes when the temperature difference entering the definition of $Q_{10}$ goes to zero, i.e. when the finite-difference approximation in Eq. (9)(9) $\begin{array}{cc}\frac{ln\left({\mathit{Q}}_{10}\right)}{10\mathrm{K}}& =\frac{ln\left({\mathit{K}}_{\mathrm{eq}}\left(\mathit{T}+10\mathrm{K}\right)\right)-ln\left(K_{\mathrm{eq}}\left(\mathit{T}\right)\right)}{\left(\mathit{T}+10\mathrm{K}\right)-\mathrm{T}}\\ \mathit{\hspace{0.17em}}& \approx \frac{\mathrm{d}ln\left({\mathit{K}}_{\mathrm{eq}}\right)}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta }H\left(T\right)}{R{T}^2},\end{array}$(9) becomes exact.

Figure 5. Comparison of enthalpy estimates obtained from van ’t Hoff fits (blue solid lines) and $Q_{10}$-based analysis (red dashed lines). (a) For the TRPV3 channel data analyzed in Figure 3, the finite-difference approximation of the derivative with respect to $\mathit{T}$ in Eq. (9)(9) $\begin{array}{cc}\frac{ln\left({\mathit{Q}}_{10}\right)}{10\mathrm{K}}& =\frac{ln\left({\mathit{K}}_{\mathrm{eq}}\left(\mathit{T}+10\mathrm{K}\right)\right)-ln\left(K_{\mathrm{eq}}\left(\mathit{T}\right)\right)}{\left(\mathit{T}+10\mathrm{K}\right)-\mathrm{T}}\\ \mathit{\hspace{0.17em}}& \approx \frac{\mathrm{d}ln\left({\mathit{K}}_{\mathrm{eq}}\right)}{\mathrm{d}\mathit{T}}=\frac{\mathrm{\Delta }H\left(T\right)}{R{T}^2},\end{array}$(9) results in a vastly different $\mathrm{\Delta H}(\mathit{T})$ estimate than obtained from our spline-fitting procedure. (b) The van ’t Hoff fit of the TRPV1 channel data predicted a temperature-independent enthalpy, for which $Q_{10}$ provides a decent estimate of $\mathrm{\Delta H}$ (in this case only 3% off).

Conclusions

We have developed an open-access data analysis package [17] to reliably extract thermodynamic potentials and heat capacities from empirical measurements of equilibrium constants at different temperatures. This is in line with most recent theoretical developments, which show that a nontrivial heat capacity difference $\mathrm{\Delta }{C}_p(\mathit{T})$ can be crucial for the correct modeling of the thermodynamics of temperature-dependent channel gating [14]. Our package accounts for the fact that on physiologically relevant temperature scales one cannot distinguish between a linear and reciprocal temperature dependence (see Figure 2), and therefore fits multiple models to the data, which vary in complexity (characterized by the number of parameters) and in the way they scale with temperature. A Bayesian information criterion [Eq. (18)(18) $\mathrm{BIC}\left(\mathit{N}\right)={\chi }_{min}^2+2\mathit{N}ln\left(\mathit{M}+1\right),$(18) ] is used to select the model that best fits the data, while minimizing the number of model parameters to avoid overfitting. Our software can therefore be used to fit nonlinear van ’t Hoff plots without any ad hoc assumptions and outperforms conventional methods, such as the thermal coefficient $Q_{10}$ (see Figure 5). Yet, we urge users to practice caution and not use our package to analyze data containing artifacts or unreasonably small error bars, because these can affect the resulting model selection and lead to faulty conclusions.

To demonstrate the use of the data analysis package, we applied it to measurements of equilibrium constants for the temperature-sensitive TRPV1 and TRPV3 channels, respectively. For both data sets, we found that models with the functional form $ln(K_{\mathrm{eq}})=\mathit{f}(1/\mathit{T})$, i.e. reciprocal in the temperature, were best suited to fit the data at hand, albeit with differing complexity. While the fit to the TRPV1 data predicted $\mathrm{\Delta }{C}_p=0$, and therefore a constant enthalpy and entropy (see Figure 4), an initially increasing and then monotonically decreasing $\mathrm{\Delta }{C}_p(\mathit{T})$ was required to capture the extremely nonlinear trends seen in the TRPV3 data (Figure 3). Note that our analysis relies on the common assumption that the TRPV1 and TRPV3 channels can, to a first approximation, be described as two-state systems. Even if this assumption is unlikely to hold for most (if not all) TRP channels, the thermodynamic variables that can be extracted from our approach provide a more reliable and accurate description of the temperature sensitivity of TRP channels than the results of a $Q_{10}$-analysis or linear fits to a van ’t Hoff plot. In principle, a more intricate analysis involving multiple states can also be conducted (see the subsection “Extension to multi-state models”), but has not been implemented in our data analysis package for the simple reason that every mechanistic model would have to be implemented separately.

It is our belief that our data analysis tool will not only benefit the community of electrophysiologists studying temperature-sensitive channels, but also help researchers in chemistry and biochemistry to rigorously analyze their van ’t Hoff plots. Temperature-dependent enthalpies and entropies open up exciting new possibilities in the theoretical modeling of the kinetics and dynamics of thermoresponsive systems, as the associated transition rates between the open and closed state intuitively must exhibit non-Arrhenius behavior. Whether such generalized models are applicable to TRP channels should be addressed in future research.

Acknowledgments

We thank Dr. Andrés Jara-Oseguera for fruitful discussions and critical comments on the manuscript. The Flatiron Institute is a division of the Simons Foundation.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Max-Planck-Gesellschaft (JTB) and the Simons Foundation (SMH).