April 2021: Fitting sigmoids with exponential functions.
Published:
Some thoughts on how we fit exponentials to sigmoidal waveforms.
Table of Contents
Introduction
One of the main goals of ion channel electrophysiologists is to describe the kinetics of ion channel currents - that is, how currents develop over time. For example, to characterize a voltage-gated ion channel, an experimenter would record a series of protocols that close, open, or inactivate the channel at multiple voltages. The resulting data would describe the voltage-dependence of channel deactivation, activation, and inactivation, respectively. However, analyzing and comparing current traces directly is often cumbersome, impractical, and unintuitive. Therefore, electrophysiologists have come up with ways to summarize their data, one of the most widely-used being fits of exponential functions.
Why exponential functions?
Before explaining why exponential functions are used, let’s first see why we can’t understand underlying kinetics by directly comparing raw current traces. Some basic reasons include experimental or technical variability, such as differences in
- transfection efficiency, expression, or membrane trafficking
- the size (conductance) of endogenous leak currents, pipette and seal properties (membrane resistance, membrane capacitance, series resistance, etc.)
These factors modify how the currents are conducted and/or recorded, so it doesn’t make sense to directly compare current recordings (in most cases).
Besides comparing currents directly, we could potentially normalize current amplitudes, or come up normalized statistics such as the time between the onset of activation and half-maximal activation at a given voltage. Actually, I think these methods are quite reasonable, depending on the context. However, they are not used because they are unintuitive. Simply put, fits of exponential functions provide a more mechanistic view of current kinetics. And, as with many things in electrophysiology, a lot of it has to do with Hodgkin and Huxley.
For those unfamiliar with Hodgkin and Huxley’s work, the relevant part that I’ll reference here is very simple. Hodgkin and Huxley were recording currents from preparations of the squid giant axon (so-called giant axons from a normal squid, not an axon from a giant squid - a common misconception), and used a mathematical model to describe how these currents are generated.
Insights from Hodgkin & Huxley
Part of their model assumed that the currents they recorded were generated from a process with multiple independent gating subunits. Let’s just consider the case where a channel exists in only two states - closed or open.
\[C \overset{\alpha(V)}{\underset{\beta(V)}\rightleftarrows} O\]Where, $C$ and $O$ represent closed and open states, and $\alpha(V)$ and $\beta(V)$ are voltage-dependent rates of opening and closing, respectively. The rate of change in open occupancy can be written as:
\[\frac{dO}{dt} = \alpha C - \beta O = \alpha (1 - O) - \beta O = \alpha - (\alpha + \beta) O\]Where, we have dropped the voltage-dependency of $\alpha$ and $\beta$ for clarity. This Ordinary Differential Equation (ODE) is separable, meaning we can put the variables in the differential - all the $O$ terms and $t$ terms - on separate sides of the equations:
\[\frac{1}{\alpha + \beta}\frac{dO}{dt} = \frac{\alpha}{\alpha + \beta} - O\] \[(\alpha + \beta) \ dt = \left(\frac{\alpha}{\alpha + \beta} - O \right)^{-1} dO\]Integrating both sides,
\[(\alpha + \beta) \ t + C_t = - \ln \left(\frac{\alpha}{\alpha + \beta} - O \right) + C_O\] \[-(\alpha + \beta) \ t + C = \ln \left(\frac{\alpha}{\alpha + \beta} - O \right)\]Where, we’ve combined the two constants of integration into one (since it’s an arbitrary number anyhow). Next, we exponentiate both sides to remove the logarithm:
\[\exp (-(\alpha + \beta)t) \exp(C) = \frac{\alpha}{\alpha + \beta} - O(t)\] \[O(t) = \frac{\alpha}{\alpha + \beta} - C\exp (-(\alpha + \beta)t)\]Where, we’ve let $C = \exp(C)$, which is possible because it’s an arbitrary number. Let’s assume we know the initial and steady-state values $O_0 = O(0)$ and $O_{\infty} = O(\infty)$, respectively. Then, we have what’s called an Initial Value Problem, and we can solve for $C$ as follows
\[O_\infty = O(\infty) = \frac{\alpha}{\alpha + \beta}\] \[O_0 = O(0) = \frac{\alpha}{\alpha + \beta} - C \implies C = \frac{\alpha}{\alpha + \beta} - O_0 = O(\infty) - O_0\]Therefore, we can rewrite $O(t)$ as
\[O(t) = O_\infty - (O(\infty) - O_0) \exp (-t/\tau)\]Where we have defined $\tau = 1/(\alpha + \beta)$, which is known as the time constant for exponential functions. Exponential functions of this type appear throughout nature, including various growth and decay processes, in part because they are solutions to processes whose rates of change depend linearly on their present values. For instance, we have shown that this is the case for this simple model of ion channel gating.
Thus, because currents depend on the occupancy of open states, it is natural to use exponential functions to describe how currents develop. In particular, we see that, for a two-state model, current kinetics are entirely determined by the time constant $\tau$, whereas the other parameters control the amplitude. As such, time constants are often used to describe current kinetics. In summary,
- Exponentials naturally arise from mathematical modelling of how channels open and close.
- The kinetics of exponentials, and thus of current and open occupancy, are determined by their time constants $\tau$.
- Current kinetics can are often summarized by showing time constants derived from fitting exponential functions.
Considerations for Sigmoid Currents
In general, the following approach works quite well: record current activation at multiple voltages, fit an exponential, then plot the time constants with respect to, for example, voltage to visualize the voltage-dependence of activation kinetics. There are several cases where this approach is inadequate. One is when a single exponential,
\[I(t) = A \exp (-t/\tau) + C\]doesn’t fit the current well. The equation above is the general form of a single-exponential function, where $A$ and $C$ are the amplitude of the exponential and a steady-state constant, respectively.
When a single exponential component is insufficient, the general approach is to use more:
\[I(t) = \sum^{N}_{i=1} \left( A_i \exp(-t/\tau_i) \right) + C\]Thus, an equation with $N$ exponential components, also sometimes referred to as the sum of $N$ exponentials, involves $2N + 1$ parameters. Since $N$ generally never exceeds three or four, this poses little issue computationally with modern resources. However, the important question is, when and how do we fit higher numbers of exponential components?
Sigmoid currents are particularly ill-fit by a single exponential, due to having distinct regions of positive and negative concavity/convexity. In this case, even two exponentials may be insufficient. A common workaround is to fit a single or double exponential function after an initial delay. often, the delay has pronounced concavity of opposite sign to the rest of the current trace, and is sometimes discussed as evidence of intermediate opening steps, implying mechanisms with more than just one closed and open state.
But, I won’t delve too far into what the so-called “delay” could represent mechanistically. Instead, I wnat to talk about: (1) how to quantify/measure the delay, and (2) how to fit exponentials when a delay is present.
Quantifying delays in sigmoid currents.
While the extent of delay usually appears inversely correlated with faster current kinetics, one challenge with quantifying the delay is that the degree of sigmoidicity can vary greatly between individual cells expressing the same ion channel. Due to this variability, currents with limited sigmoidicity can sometimes be well-fit by double exponentials, or even single exponentials, e.g. under conditions that yield maximal kinetics.
What I do is fit exponential functions while including the delay as an additional parameter. A cost function that considers the number of data points used is crucial, otherwise the delay can be extended to minimize the residual, simply by reducing the amount of data used for fitting. To inform the fit, it is useful to set reasonable initial values and bounds, if using algorithms that support or require these inputs. An initial value can usually be reasonably estimated as the inflection point of the current, e.g. as the zero-crossing of the second time derivative. Thereafter, bounds can be set around the initial value to maintain a consistent window for parameter searching while automatically accounting for how the extent of delay would change between traces, e.g. voltages.
In my experience, I’ve also found it useful to test different algorithms, as is often done for more complex optimization problems. Personally, some methods, such as BFGS tend to perform relatively poorly compared to others, such as Nelder-Mead or AMPGO.
Traditionally, however, it’s uncommon to ‘fit’ the delay as I describe above. Instead, a manual ‘fitting’ is performed by determining when in a current trace a single (or double) exponential results in an adequate fit. This method can lack precision, and should be somewhat dependent on the algorithm(s) used, but in general, I haven’t seen much evidence suggesting that it’s any less accurate than the previous method.
A Final Note.
When fitting with multiple exponential components, it can be unclear how to interpret the resulting parameters. With two exponential components, most people choose to show both time constants, as well as a proportion of amplitudes, e.g. $\frac{A_1}{A_1 + A_2}$. Distinct exponential components are sometimes described as representing distinct populations, or states, so that the proportions reflect relative occupancies of these states.
Conclusion
Thanks for reading this far! Here’s my Twitter handle, if you’d like to share your thoughts: @ydhaganeno