A thermodynamic description for physiological transmembrane transport

A general formulation for both passive and active transmembrane transport is derived from basic thermodynamical principles. The derivation takes into account the energy required for the motion of molecules across membranes and includes the possibility of modeling asymmetric flow. Transmembrane currents can then be described by the general model in the case of electrogenic flow. As it is desirable in new models, it is possible to derive other well-known expressions for transmembrane currents as particular cases of the general formulation. For instance, the conductance-based formulation for current turns out to be a linear approximation of the general formula for current. Also, under suitable assumptions, other formulas for current based on electrodiffusion, like the constant field approximation by Goldman, can be recovered from the general formulation. The applicability of the general formulations is illustrated first with fits to existing data, and after, with models of transmembrane potential dynamics for pacemaking cardiocytes and neurons. The general formulations presented here provide a common ground for the biophysical study of physiological phenomena that depend on transmembrane transport.


Introduction
One of the most important physiological mechanisms underlying communication within and between cells is the transport of molecules across membranes. Molecules can cross membranes either passively (Stein & Litman, 2014), or via active transport (Bennett, 1956). Molecules are passively transported across a membrane when they move along their (electro)chemical gradient. In contrast, active transport involves transmembrane motion of molecules against their electrochemical gradients. One important functional distinction between channels and pumps is that the rate of transport for channels is generally several orders of magnitude faster than the rate for pump-mediated transport (Gadsby, 2009;Ussing, 1949a). Such differences are reflected in the sizes of different transmembrane fluxes typically observed in excitable cells (Herrera-Valdez & Lega, 2011).
Passive transport may occur through transmembrane proteins (Hille, 1992; Stein & Litman, 2014) that may be selective for molecules of specific types (Almers & McCleskey, 1984;Doyle et al., 1998;Favre et al., 1996), typically mediating (electro)diffusion through them. Sometimes these proteins are also gated by conformational changes triggered by different signlling mechanisms. Passive transport has also been observed in pores spontaneously formed within synthetic lipid bilayers (Blicher & Heimburg, 2013), which could also occur in natural conditions (Gurtovenko & Anwar, 2007). One example of importance in the context of energy homeostasis is the transport of monosaccharides and other similar molecules through GLUT transporters. GLUT transporters first bind to their substrates, triggering a conformational change that allows the substrate to cross in the direction of its chemical gradient (e.g. GLUT5 is highly specific to fructose, Mueckler & Thorens (2013)). Channels are another important class of transmembrane proteins that typically mediate fast passive transport, often displaying selectivity for specific ion types, gated by changes in the membrane potential (Covarrubias et al., 1991;Peng & Wu, 2007) or the binding of a ligand molecule (e.g. AMPA-Kainate synaptic receptors, Bowie (2002)).
Active transport is mediated by transmembrane proteins commonly called pumps, or carriers, that mechanically translocate the molecules they transport (Bennett, 1956;Ussing, 1949b;Ussing, 1949c). The energy for primary active transport is usually obtained from biochemical reactions (e.g. ATPases, light-driven pumps). For instance, the energy to transport molecules against their (electro)chemical gradient via ATPases is obtained from hydrolysis of ATP (Chapman, 1973). In secondary and tertiary active transport, the electrochemical gradient of at least one the ion types provides the energy to transport other molecules against their (electro)chemical gradient (Skou, 1965). Two large classes of non-primary active transport pumps are the symporters and counterporters, carrying at least two types of molecules in the same, or in opposite directions, respectively, with at least one type against its electrochemical gradient. For instance Na-H exchangers carry Na + and H + in opposite directions, typically using the driving force from Na + . In contrast, K-Cl symporters carry K + and Cl − in the same direction, which means that one of the two ion types is carried against its concentration gradient. This is because the K + and Cl − concentration gradients are typically oriented in opposite directions. As a consequence, the movement of one of the two ions releases energy from its electrochemical gradient, enabling the transport of the other ion against its gradient.
Theoretical models of transmembrane transport play a critical role in developing our understanding of the function and mechanisms underlying electrical signalling and cellular excitability (Barr, 1965;Cole, 1965;DiFrancesco & Noble, 1985;Endresen et al., 2000;Gadsby, 2009;Goldman, 1943;Kell, 1979;Kimizuka & Koketsu, 1964;Läuger, 1973;Stevens & Tsien, 1979;Wiggins, 1985a;Wiggins, 1985b;Wiggins, 1985c), and some of its associated pathologies (Ashcroft, 2005;Marbán, 2002). The best known transmembrane transport models include the widely used conductance-based formulation from the seminal work of Hodgkin & Huxley (1952), the Goldman-Hodgkin-Katz equation (Goldman, 1943;Hodgkin & Katz, 1949;Pickard, 1976), and several other expressions for carrier and channel mediated transport with many different functional forms (DiFrancesco & Noble, 1985;Rasmusson et al., 1990a;Rasmusson et al., 1990b;Rosenberg & Wilbrandt, 1955). Other formulations for ionic transport across membranes derived from biophysical principles available in the literature include those in the work by Jacquez & Schultz (1974); Kimizuka & Koketsu (1964); Pickard (1969); Pickard (1976). See also Jacquez (1981) and similar work by Endresen et al. (2000), and those in the excellent book by Johnston et al. (1995). Such formulations describe the relationship between the activity and permeability of ions across membranes, and the transmembrane potential. However general models that describe physiological transport including passive and active transport of charged or non-charged molecules, with bidirectional but possibly asymmetric flows, are still missing. The work presented here builds upon the results previously mentioned by describing transport macroscopically in terms of the energy required to move molecules across a membrane. The result is a general formulation with a common functional form for both passive and

Amendments from Version 2
The manuscript has been edited to correct minor typos and rephrase sentences for clarity throughout the document and also, to include data and two sections that were previously in supplementary files. In more detail, formulas were corrected in one header and in the last column of Table 1; corrections have also been made in other sections where the reversal potential for the Na-Ca transport is mentioned. Also, amplitude descriptions in the third column of Table 3 and Table 4 were shortened, and some parameter values were adjusted in Table 3. The data in Figure 3 is now included in Table 2. The two sections added to the manuscript contain the derivation of the Goldman constant field approximation from the general formula derived in the article, and a model of neuronal membrane potential dynamics that uses the general formulation for transport.
Any further responses from the reviewers can be found at the end of the article REVISED active transport (Herrera-Valdez, 2014) that also includes a term that regulates the asymmetry in the flow (rectification). The details of the derivation and examples of fits to experimental data and features like asymmetric bidirectional flow can be found in the next section. An application of the general formulation is illustrated with models for the transmembrane potential dynamics in cardiac pacemaker cells and striatal fast spiking interneurons. A derivation and connection of the general formula with existing formulations like the Goldman expression for current can be in the paragraphs below.

General formulation for transmembrane flux and current
Work required for transmembrane molecular fluxes Consider a system consisting of a biological membrane surrounded by two aqueous compartments (e.g. extracellular and intracellular). Assume, to start with, that the compartments contain molecules of a single type s (e.g. Na + , K + , glucose), possibly in different concentrations. Let ∆G s be the energy required for the transport of the molecules across the membrane in a specific direction (e.g. inside to outside). To write an expression for ∆G s it is necessary to take the direction of motion of the s-molecules into account. To do so, label the extracellular and intracellular compartments as 0 and 1, respectively, and let c s and d s ∈ {0, 1} represent the source and the destination compartments for the transport of the s-molecules. The pair (c s , d s )=(0,1) represents inward transport and the pair (c s , d s )=(1,0) represents outward transport. The work required to transport n s molecules of type s from compartment c s to compartment d s can then be written as where v s is the Nernst potential for the s-molecules 1  If ∆G s < 0, then the molecules can be transported passively (e.g. electrodiffusion), decreasing the electrochemical gradient for s across the membrane. In contrast, if ∆G s > 0, the transmembrane transport of s from c s to d s is not thermodynamically favourable, which means the transport from c s to d s requires energy that is not available in the electrochemical gradient for s (active transport). As a consequence, active transport of s would increase the driving force for the motion of s across the membrane.
Joint transmembrane transport of different types of molecules To find an expression that describes a more general transport mechanism, assume that transport takes place as single events in which molecules of m different types move in parallel, or possibly sequentially (e.g. first Na + , then K + ), across the membrane. Let S be a set that represents the types of molecules that are jointly transported in a single event. For instance, for Na + -H + exchangers, S = {Na + , H + }, with m = 2. The energy required to transport these molecules is the sum of the energies required to transport each type of molecule in S. In other words, As before, transport is thermodynamically favourable when ∆G S ≤ 0. If not, extra energy is required. To distinguish between these two cases, define the total energy ∆G of the transport mechanism, possibly including an extra source of energy, as where δ Ext = 1 if ∆G S > 0, and 0 otherwise. If ∆G S > 0, then the energy from ATP hydrolysis or any other sources represented by δ Ext ∆G Ext should be negative and larger in size in comparison to ∆G S , so that ∆G ≤ 0, making the transport thermodynamically possible. In particular, for ATP-driven transport, the extra energy supplied by hydrolysis of ATP (Tanford, 1981;De Weer et al., 1988) is could be derived for active transport driven by light, or other sources of energy. The concentrations of ATP, ADP, and P i are assumed to be constant in most models presented here, but it should be noted that such concentrations are not necessarily constant, and in fact, may vary a lot in some cases, as it is the case for skeletal muscle (Wackerhage et al., 1998).
Flux due to transmembrane transport The formulation in Equation (5) can be combined with Equation (1) to derive a generalised expression for flux and model different known mechanisms of physiological transmembrane transport, possibly combining the transport of different molecules simultaneously (e.g. Na-H exchange). In this case, the forward direction of the transport would be described by the combined forward transport of each of the different molecules under consideration. For instance, the source and target compartments for Na + and Ca 2+ are different in Na-Ca exchangers. The stoichiometry for the transport mediated by Na-Ca exchangers in the forward direction involves three Na + molecules moving inward (along their electrochemical gradient) in exchange for one Ca 2+ molecule moving outward (against their electrochemical gradient) (Mullins, 1979;Venetucci et al., 2007).
Let α and β be the flux rates in the forward and backward directions, in units of molecules per ms per µm −2 . These rates depend, a priori, on the energy required for the transport of the molecules in S. The net flux rate associated to the net transmembrane transport, can then be written as How do α and β depend on ∆G? The steady state relationship between the energy ∆G and the the forward and backward flow rates, hereby represented by α and β, can be expressed by means of a Boltzmann distribution (Boltzmann, 1868) as Assuming that α and β are continuous functions, the rates α and β can be rewritten as where r and b represent the rate at which the transport takes place (molecules per ms per µm −2 ) and the bias of the transport in one of the two directions (b = 1/2 means the transport is symmetrical relative to the point at which ∆G=0). Note that the functional form of the α and β in Equations (9) are similar to those by Butler (1924); Erdey-Grúz & Volmer (1930). Also, notice that the steady state relationship between α and β in Equation (8) can be obtained from Equations (9), for any r and any b. However, it should be the case that r and b vary in specific ranges depending on the physicochemical characteristics of the pore through which molecules cross the membrane, and in general, on the transport mechanism. As already mentioned, the rate r should be larger for electrodiffusive transport mediated by ion channels, in comparison to the slower transport rates for facilitated diffusion and active transport meditated by carrier proteins and pumps. The rate r may depend on temperature ( The flux can then be written explicitly combining Equation (7) and Equation (9) to obtain, Taking the above observations into account, it is possible to combine Equation (4) and Equation (5), to write an expression similar to Equation (8) for the steady state balance between the forward and backward transport of all the molecules in S (see Table 1 for examples of different transport mechanisms with their energies and total charge movements).

Flux and current
After substitution of the formulas for ∆G from Equation (4) and Equation (5) where v T = kT/q and represents the net number of charges moved across the membrane. Note that v Ext should be v ATP for ATPases. The first, more complex, form of the flux in Equation (11) could be useful when working with models for which changes in the concentrations of different molecules are relevant.
If the transport is electrogenic, then the product qη (in Coulombs) represents the net charge moved across the membrane, relative to the extracellular compartment. Non electrogenic transport yields η = 0, which means the flow does not depend on the transmembrane potential, and If only ions are involved in the transport, the flux simplifies to where Ext Ext The quantity v o /η can be thought of as a reversal potential. If η < 0, then positive charge is transported inward, or negative charge is transported outward. In contrast, η > 0 means that positive charge is transported outward or negative charge transported inward. For instance, inward electrodiffusion of single Na + ions gives η = −1, which can be thought of as loosing one positive charge from the extracellular compartment in each transport event.
Transmembrane current. A flux that results in electrogenic transport (Equation (11) and Equation (14)) can be converted to current density after multiplication by qη. In short form, with qr in Amperes/m 2 or equivalent units.
Substitution of Equation (11) or Equation (14) into Equation (16) yields a general formula for the current generated by transmembrane ionic flux (Figure 2), that uses the same functional form for channels (protein or lipid) and pumps. Recall that Equation (16) can also be written explicitly in terms of the transmembrane concentrations of one or more of the ions involved using Equation (11). It is possible to derive expressions for r that take into account biophysical variables like temperature and the shape and length of the pore through which the molecules cross (Endresen et al., 2000; Pickard, 1969).

Special cases and examples
A number of nontrivial and important properties of transmembrane ionic currents, including rectification, are also described by Equation (16)  The energy ∆G S necessary can be found on the right axis if the transport is primary active, and on the left axis when transport is secondary active or passive. Extra energy from ATP hydrolysis or other sources added to ∆G S shift the total energy to the left axis and make the transport thermodynamically possible. See Table 1 for examples.

Conductance-based currents (Hodgkin & Huxley, 1952) can then be obtained as linear approximations of the general current
where g = η 2 qr/v T is in units of nS/µm 2 . For instance, the linear approximation for the current through open sodium channels around v Na in Equation (18) gives g Na = qr Na /v T , and v o /η Na = v Na , with Note that if the truncation order is larger than 1, then approximations from Equation (17) include the rectification parameter b. One already known consequence is that conductance-based currents do not capture rectification.
The Goldman constant field approximation is a particular case of the general formulation. The Goldman-Hodgkin-Katz equation describing the transmembrane current carried by ions of type x is given by where x j , j ∈ {0, 1}, represent the extra and intracellular concentrations of an ion of type x.
Multiplying the numerator and denominator of Equation (19) by , and algebraically rearranging terms yields, where is an amplitude term that can be approximated by a constant . This was first called "anomalous rectification" by Katz (1949), who noticed that K + flows through muscle membranes more easily in the inward, than in the outward direction (Adrian, 1969; Armstrong & Binstock, 1965). It was later found that some K + channels display the bias in the opposite direction (Woodbury, 1971). The former type of K + current rectification is called inward, and the latter outward ( Figure 2).
Rectification is thus a bias in one of the two directions of transport. The type of rectification (inward or outward) depends on the molecules being transported and on the structure of the proteins mediating the transport. Flux rectification is not only displayed by ions, as shown by molecules like glucose, which may cross membranes via GLUT transporters bidirectionally and asymmetrically, even if the glucose concentration is balanced across the membrane (Lowe & Walmsley, 1986).
Rectification can be described by setting b in Equation (11) to values different from 1/2 and it becomes more pronounced as b is closer to 0 or 1. These values represent biases in the transport toward the source, or the target compartment, respectively. As a consequence, rectification yields an asymmetry in the graph of α − β as a function of ∆G ( Figure 1). For electrogenic transport, rectification can be thought of as an asymmetric relationship between current flow and voltage, with respect to the reversal potential v o . The particular case b = 1/2 (non rectifying) yields a functional form for current similar to that proposed by Pickard (1969) From here on, subscripts will be used to represent different transport mechanisms. For instance, the current for a Na-Ca pump will be written as i NaCa .
Electrodiffusion of K + through channels (η = 1 and v o = v K ), is outward for v > v K , and inward for v < v K . The K + current through the open pore is therefore Current flow through inward rectifier K-channels (Riedelsberger et al., 2015) can be fit to values of b K < 1/2. For instance, describes a current generated by a K + flow that is limited in the outward direction, similar to the currents described originally by Katz (1949). Analogously, b K > 1/2 limits the inward flow. For example, the current describes outward rectification (Riedelsberger et al., 2015).
Based on the work of Riedelsberger et al. (2015) on K + channels, inward (outward) rectification arises when the S4 segment in K + channels is located in the inner (outer) portion of the membrane. These two general configurations can be thought of in terms of ranges for the parameter b K , namely, b K < 1/2 for inward, and b K > 1/2 for outward rectification (Figure 2).
In general, ion channels are typically formed by different subunits, that combined produce structural changes that may result in rectified flows. For instance, non-NMDA glutamatergic receptors that can be activated by kainic acid and α-amino-3-hydroxy-5-methyl-4-isoxazole propionic acid (AMPA), display different permeabilities to Na + , K + , and Ca 2+ depending on the subunits that form the receptor (Hollmann et al., 1991). In particular, the Ca 2+ currents recorded in oocytes injected with combinations of GluR1 and GluR3 cRNA have different steady state amplitudes and show different levels of rectification (Figure 3).
Primary active transport. The Na-K ATPase is a primary active transporter that uses the energy from the hydrolysis of one molecule of ATP for the uphill transport of 3 Na + ions outward In response to steroids like strophandin, the voltage-dependence of the Na-K ATPase current has been reported to show a plateau as v increases past the reversal potential for the current (Nakao & Gadsby, 1989). In such cases, the Na-K Table 1. Energy required for transmembrane transport mediated by different passive and active mechanisms.  .

Pump or channel
The rectification for the Na-K pump ATPase has also been reported to occur in small neurons of the dorsal root ganglion in rats (Hamada et al., 2003). The alternative expression (28) also explains qualitatively different behaviors of the Na-K current as a function of the transmembrane concentrations of Na + and K + . For instance, if either [Na] 1 or [K] 0 increase and v > v NaK , then the amplitude of i NaK would increase at a smaller rate of change in comparison to when v < v NaK , which grows exponentially in size. This is also in line with reports of non significant changes in the transport by Na-K ATPases in response to elevated intracellular Na + during heart failure (Despa et al., 2002), in which the transmembrane potential is likely to be depolarised.
Secondary active transport. An example of a pump that mediates secondary active transport is the Na-Ca exchanger, which takes 3 Na + ions from the extracellular compartment in exchange for one intracellular Ca 2+ ion in forward mode (Pitts, 1979; Reeves & Hale, 1984). The reversal potential for the current is v NaCa = 2v Ca − 3v Na , with η NaCa = −1. Assuming b NaCa = 1/2, the Na-Ca current can be rewritten as NaCa NaCa NaCa The driving force v − v NaCa could reverse in sign with large enough increases in the intracellular concentration of Ca 2+ , or in the membrane potential. As a result, the current could have a dual contribution to the change in transmembrane potential, as predicted by some theoretical models of cardiac pacemaker activity (Rasmusson et al., 1990a; Rasmusson et al., 1990b).
Electrodiffusive transport. Consider transmembrane electrodiffusive transport of a single ionic type x, with z x and v x representing the valence and the Nernst potential for ions of type x, respectively. In this case, the reversal potential satisfies which means that η x can be factorised in the argument for the exponential functions and the general expression (16) can be rewritten as .
In the absence of rectification (b x = 0.5) the formula simplifies to See Jacquez & Schultz (1974); Pickard (1969) Table 1 for other examples.
The applicability of the general formulations described above is illustrated next with models of cardiac and neuronal membrane potential.

Transmembrane potential dynamics
To illustrate applications for the formulations discussed earlier, let us build a general model of transmembrane potential dynamics. For simplification purposes, consider only one such mechanism, labelled as l ∈ {1, ..., M}, with p l N l active sites (e.g. channels), where N l is the number of membrane sites where the proteins mediating the lth transport mechanism are found, and p l is the proportion of active sites (activation might depend on voltage or a ligand). Let η l , and r l respectively represent the net number of charges moved in each transport event, and the rate of molecular transport per unit area (molecs ms -1 ⋅ µm -2 ) of a single protein mediating the lth current, l ∈ {1, ..., M}. Then the product η l r l q represents the current density of a single protein of type l (Coulombs/sec µm 2 ), and the total current mediated by the lth mechanism in a patch of membrane can be written as ā l p l l ϕ (v) with ā l = η l N l r l q (in pA/µm 2 ), and T T ( ) exp exp (  1) , where v l /η l represents the reversal potential for the lth current, l ∈ {1, ..., M}. There is experimental evidence for some ion channels that supports the replacement of η l r l q for a constant (Nonner & Eisenberg, 1998). For instance, it is reasonable to assume that the single channel current for many types of Na + channels is 1 pA (Aldrich et al. (1983)). Assuming that the change in charge density with respect to voltage (typically referred to as membrane capacitance in conductance-based models) is a constant C m (pF/µm 2 , Herrera-Valdez (2020)), the time-dependent change in transmembrane potential can then be written as 1 ( ), with v in mV and a l = ā l /C M in mV/mS (pA/pF) represents the normalized current amplitude for the lth transport mechanism, l ∈ (1, ..., N). Note that only electrogenic transport mechanisms are included.
The temporal evolution for v can then be described by where w represents the proportion of activated K + channels (and also, the proportion of inactivated where x ϕ is a difference of exponential functions as defined in Equation (16), with x ∈ {NaK, NaCa, KD, CaL}.
The steady state for the activation of voltage-dependent channels in the model is described by a generic function of the form with labels m and w used for L-type Ca 2+ and delayed rectifier K + channels, respectively. F u is increasing as a function of , where b w represents a bias in the conformational change underlying channel activation. The function R w has the shape of a hyperbolic cosine when b w is 1/2.  Figure 4A, blue line). Between the initial depolarisation and until the maximum downstroke rate, approximately, v NaCa < v, which means that J NaCa > 0. Then, Ca 2+ extrusion by the Na-Ca exchanger occurs only for a brief period of time during the downstroke and also after each action potential ( Figure 4C, blue line). Second, as previously reported in different studies involving spiking dynamics, the time course of the Ca 2+ current shows a partial inactivation with a double peak ( Figure 4B,   double peak can be much simpler. The calcium current J CaL is a negative-valued, non monotonic function for v < v Ca , which can be thought of as a product of a amplitude term that includes gating and the function φ CaL . The normalised current J CaL has a local minimum (maximum current amplitude) around -10 mV ( Figure 4B, blue line and Figure 5A, blue line), after which the current decreases, reaching a local maximum as the total current passes through zero, at the peak of the action potential around 10 mV ( Figure 4B, where ∂ t v = 0). The first peak for the Ca 2+ current occurs when v reaches the maximum depolarisation rate ( Figure 5B). As v increases (e.g. upstroke of the action potential). The second peak for the current occurs as the membrane potential decreases, and passes again through the region where the maximal current occurs (local minimum for J Ca ). The two local minima for J CaL occur at different amplitudes because of the difference in the evolution of v during the upstroke and the downstroke of the action potential ( Figure 5A, blue line, and Figure 5B, where ∂ t v = 0). It is important to remark that the dual role played by w is not the cause of the double activation. This is illustrated by analysing the behaviour of a non-inactivating J CaL without the inactivation component ( Figure 5A, gray line). The double activation can also be observed in models in which the activation of K + channels and the inactivation of Ca 2+ or Na + channels are represented by different variables (Rasmusson et al., 1990a)  The double peak in the Ca 2+ current reflects on the intracellular Ca 2+ concentration ( Figure 6, gray line), and by extension, on the Nernst potential for Ca 2+ (Figure 6, blue line), which display two increasing phases and two decreasing phases, respectively. The first and faster phase in both cases occur during the initial activation of the L-type channels. The second phase occurs during the downstroke, as second peak of the Ca 2+ current occurs. As a consequence, the reversal potential for the Na-Ca exchanger, v NaCa = 2v Ca − 3v Na (Figure 6, orange line) also has two phases, this time increasing. Increasing the intracellular Ca 2+ (Figure 6, gray line) concentration decreases the Nernst potential for Ca 2+ , and vice versa. By extension, v NaCa , becomes larger when c increases. Ca 2+ enters the cell in exchange for Na + that moves out when v > v NaCa , during most of the increasing phase and the initial depolarisation phase of the action potential (blue lines in Figure 4A and C, and Figure 6).

Fast spiking interneuron dynamics
To construct of a simple model for the dynamics of a fast spiking (FS) striatal interneuron, assume the transmembrane potential depends on three currents respectively mediated by Na-K pumps, non-inactivating K + channels, and Na + channels with transient dynamics, with voltage-dependent gating in both channels. It is also assumed that the proportion of activated K + is represented by a variable w ∈ [0, 1], which also represents the proportion of inactivated Na + channels (Av-Ron et al., 1991;Rinzel, 1985). That is, 1 − w represents the proportion of non-inactivated Na + channels. Since w models the activation of a population of channels, it makes sense assume that its dynamics follow a logistic scheme, without adding extra powers to w. This also follows the experimentally observed dynamics which have been reported repeatedly, including those of delayed-rectifier K + currents recorded in voltage clamp (see for instance Figure   Explicitly, the dynamics for the FS-interneuron membrane can be described by a system with two coupled differential equations of the form The activation rate for K + channels depends is a voltage-dependent functions R w and F w as defined for the cardiac pacemaking model. It is also assumed that the activation of sodium channels is at a quasi steady state as a function of v. Striatal FS interneurons display maximum ∂ t v between 100 and 200 V/s. In current clamp mode, most neurons are silent, and show transitions between rest and repetitive spiking at a rheobase current of approximately 90 pA, with initial firing rates between 50 and 60 Hz and a delay to first spike in the transition that decreases as the stimulus amplitude increases (Figure 7, parameters in Table 4).
To include these properties into the model, the membrane capacitance was specified first, then the maximum ∂ t v was adjusted by fitting the parameter a NaT , and then the contributions for the K + channels and the Na-K ATPase are set to obtain spiking and fit the rheobase. Even though FS-interneuron spiking is conditional to the reception of input, the dynamics of the system can be thought of as qualitatively similar to those observe in the SAN. Within the time of a single spike, and by extension during the repetitive spiking regime, if v increases, w also increases, but at a slower rate in comparison to v. This is because the activation w is always moving toward its steady state value, which increases as v increases. Once w increases, the Na + current tends to decrease and the K + current tends to increase, thereby causing a decrease in v. The slower dynamics in w relative to those in v capture the delay between the amplification caused by the Na + current and the recovery caused by the negative feedback of the K + current. The current mediated by Na/K-ATPase acts as an extra attracting force toward v NaK that increases nonlinearly as the distance between v and v NaK increases.

Discussion
A general, macroscopic model for transmembrane fluxes has been derived by directly calculating the work required to transport molecules across the membrane. The derivation is based on a general thermodynamic scheme that takes into account the rate, stoichiometry, and the direction in which the molecules are transported across the membrane. These biophysical parameters are then combined to write expressions for directional fluxes based on van't Hoff (1884) and Arrhenius (1889) formulations, weighted as in the Butler/Erdey-Gruz/Volmer equation (Butler, 1924;Erdey-Grúz & Volmer, 1930). The result is a general description of the transmembrane molecular flux as a difference of exponential functions (Equation (16)) that describes the transport dynamics of molecules in the "forward" and "backward" directions, relative to a source compartment. This general description describes transport due to electrodiffusion mediated by channels, and also translocation mediated by pumps ( Table 1). The two exponential functions depend on a common expression involving the Figure 6. Calcium dynamics during pacemaking. Time courses of the intracellular calcium concentration (gray, left axis), the Nernst potential for Ca 2+ (orange, right axis), and the reversal potential for the Na-Ca exchanger (blue, right axis). Notice the two phases of calcium increase that occur in agreement with the double peak observed in the calcium current (see Figure 4B, blue trace).  Table 4. v NaK = 3v Na -2v K + v ATP -72 mV Reversal potential for the for Na + -K + ATPase current v K -89 mV Nernst potential for K + v Na 60 mV Nernst potential for Na + v mT -17 mV Half-activation potential for the transient Na + -current v w -5 mV Half-activation potential for the transient K + -current g mT 5 -Activation slope factor for the transient Na + -current g w 4 -Activation slope factor for the K + -current r w 2 s -1 Activation rate for the neuronal K + -current k w 1 -Activation exponent for the K + -current b w 0.3 -Activation slope factor for the K + -current b NaK 0.5 -Non-rectification for the Na + -K + -current b K 0.5 -Non-rectification for the transient K + -current b Na 0.5 -Non-rectification for the transient Na + -current transmembrane concentrations of the molecules being transported, and possibly the transmembrane potential when transport is electrogenic.
Rectification, an asymmetry in the flow during transport, is typically modelled modifying the dynamics of the gating variables for the current. The general formulas for transmembrane transport include a bias term b that controls the relative contribution of inward and outward fluxes the transport. Hence, different types of rectification can be described by favouring one of the directions for transport, conceptually in line with the "anomalous rectification" originally reported by Katz (1949) for K + in muscle cells. The bias term is not part of any gating mechanism. Instead, it represents the asymmetry in bidirectional flux. For instance, in K + channels, the inward (outward, respectively) rectification occurs when the fourth transmembrane segment of the channel (S4) is located closer to the intracellular (extracellular) portion of the membrane in its open configuration (Riedelsberger et al., 2015). There are other reports that show that asymmetries in bidirectional transport occur as a consequence of changes in the three dimensional structure of the protein mediating the transport (Halliday & Resnick, 1981;Quistgaard et al., 2013). Therefore, the rectification term can be thought of as representing a structural component of the transmembrane protein through which molecules move (Figure (3)).
Outward rectification in K + channels can be explained, for instance, by biasing the flux of K + the forward (outward) direction (b K > 1/2). Instead, inward rectification can be obtained by biasing the transport in the backward (inward) direction (b K < 1/2). It is important to remark that non-rectifying currents with b = 1/2 are nonlinear functions of ∆G, which shows that the nonlinearity of the current-voltage relationships is not the defining characteristic of rectification; as argued in some textbooks (see Kew & Davies, 2010).
The formulation for transmembrane flux may be rewritten in different alternative forms that can be found throughout the literature (see Equation (11) and Equation (14) Possibly of interest to mathematicians working on bifurcation theory, a third order approximation (Equation (17)), can be used to construct models resembling the Fitz-Hugh system (FitzHugh, 1955;FitzHugh, 1961;Fitz-Hugh, 1966), which yield very close approximations to the full model, while keeping biophysical characteristics of real systems like rectification and specific ionic permeabilities, and the multiplicative interaction between the slow variable w and the fast variable v; properties that the Fitz-Hugh polynomials do not have. Third order approximations also open the possibility of expanding on the analysis of dynamical systems based on these general formulas to study normal forms and bifurcations. Depending on the ions involved in each transmembrane transport mechanism, the third order approximations for current can be shown to be very close to the full function in Equation (16) (Figure 2). Another possible use of the third order approximations could be to construct network models that take into account nonlinearities included in the general formulation, but at a reduced computational cost in comparison to the full model. This possibility is currently being tested and will be reported in the near future. A similar comparison has been made between the full model and the conductance-based approximation taking a dynamical systems perspective and also by means of computational simulations in Herrera-Valdez (2012).
One question of interest because of its possible impact on the interpretation of results from existing modelling studies is how do the excitability and the resulting dynamics in a model of membrane dynamics change when using the thermodynamic transmembrane currents or their approximations? The question has been addressed in a study in which two simple neuronal models with currents mediated by Na + and K + , each equipped with the same biophysical gating properties and the same relative contributions for the currents, but one with currents as in Equation (22), the other with conductance-based currents.
The two models display a number of qualitative and quantitative differences worth considering while making the choice of a model in theoretical studies (Herrera-Valdez, 2012). For a start, the two models are not topologically equivalent across many ratios of the relative contributions of K + and Na + channels (Herrera-Valdez, 2012); as would be expected by the fact that conductance-based formulations are only linear approximations of the general currents, around the reversal potential for each current. One of the most notable differences between the general formulation and the conductance-based formula is the contribution of the nonlinear, high order terms from Equation (16), which results in more realistic upstrokes for action potentials and an overall increased excitability; in this case characterised in terms of the minimum sustained current necessary to produce at least one action potential (Herrera-Valdez, 2012). The increased excitability of the membrane with the general formulation is due, in part, to the large, exponential contribution of the open Na + and Ca 2+ channels, but not the K + channels, to the change in the transmembrane potential near rest. The time course of the Na + current during the beginning of the action potential with the general model is much sharper than that of the conductance-based formulation, resulting in a faster upstroke of the action potential; and in better agreement with observations in cortex and other tissues (Naundorf et al., 2006). It is important to remark that the sharper increase in the change of the membrane potential shown using the general formulation is a consequence of the nonlinear driving force terms of the current in the general model (the flux term in the general formulation), and not in the activation dynamics for the transient Na + current. The nonlinearities unravelled by the general formulation could thus be part of the explanation for the observed sharpness at the beginning of the action potential upstroke observed in different experiments. However, such nonlinearities do not rule out other contributions, such as cooperation between Na channels, or the effects of spatial differences in Na-channel densities, as pointed out by Brette (2013) for the case of action potentials in cortical pyramidal cells. In the models presented here, the sharpness in the upstrokes of neuronal action potentials combine nonlinear factors including the flux given by the general formulation, channel densities, and gating that does not follow linear, but logistic-like dynamics.
The general formulation for both passive and active transmembrane transport can be thought of as a tool that facilitates the construction and analysis of models of membrane potential dynamics. The generality and versatility of the thermodynamic transmembrane transport formulations is illustrated with models of cardiac pacemaking interneuron fast spiking. One important advantage of the general formulation is that it includes the possibility of explicitly estimating the number of channels or pumps mediating each of the transport mechanisms of interest. This has proven to be useful to study the relative contributions of different currents to the excitability of neurons (see Herrera-Valdez et al., 2013) and cardiocytes (Herrera-Valdez, 2014).
Another extension of possible interest is that of modelling the transmembrane transport between organelles and the cytosolic compartment, which can be done by directly replacing the difference c s − d s with -1 or 1, in Equation (1), accounting for the direction of transmembrane motion of molecules relative to the outer compartment. This and other generalisations enable the possibility of studying the interdependence between electrical excitability across tissues and animal species (Herrera-Valdez et al., 2013), and its cross-interactions with metabolism and other processes of physiological importance, all from a general theoretical framework with common formulations.

Implications for experimentalists.
One of the main advantages of the general expressions is that fits to ionic currents can be made straight from the voltage-clamp data without much effort, and without having to calculate conductances, which amounts to imposing the assumption that the current-voltage relationship is linear. Fits to experimental currents can then be directly put into equations describing the change in the membrane potential, and model membrane dynamics of interest without having to make many extra adjustments, as it is the case for most conductance-based models restricted to data.
The model for current in Equation (22) has been used to construct simplified models for the membrane dynamics of different cell types using experimental data. Examples include fast spiking interneurons in the mice striatum (Figure 7

Conclusions
A general model that describes physiological transmembrane transport of molecules has been derived by considering basic thermodynamical principles. The model unifies descriptions of transport mediated by channels and pumps, it can model biases in either one of the directions of flow, and it can be easily converted into a model for current in the case of electrogenic transport. As it is desirable in all models, the general expressions can be thought of as extensions of some previous models. In particular, it is shown that the conductance-based model for current turns out to be a first order approximation of the general formulation.
The general formulation presented here can be used to build general models of phenomena involving transmembrane transport using a unified framework (Shou et al., 2015).

Data availability
All data underlying the results are available as part of the article and no additional source data are required to reproduce the results presented here. p9, after Eq. (39): Presumably, the author here means that the intracellular Ca2+ concentration changes (rather than the extracellular). This should be clarified. Also, the text below Eq. (40) appears to imply a distinction between free and bound Ca 2+ , but there is no buffering process described in the model, which may cause confusion.

25.
p9, Fig. 2: Please state what i m and i CaL13 are in these plots. 26.
p10, Table 2: The units of \bar{a} are inconsistent with the main text. 27.
p12, "Another possible use of the third order approximations is in the construction of network models.": Please could the author state what they mean here? 28.
p13, discussion on action potential shape: I was wondering if the author could here compare their argument surrounding the shape of the action potential near initiation with those put forward here (Brette, 2013 1 ), which puts forward a case for multi-compartment models to describe the sharpness of the spike.
SI2: Headers are required in these data files so that the values herein can be related to the physical quantities they represent.

31.
It would be useful for the author to upload any code used to generate the plots in the article.
Changed log to ln 8. p3, 2nd paragraph: It would be good to clarify in this first sentence the distinction between \Delta G and \Delta G_s (or avoid using the former entirely).
Added a phrase to explain it.
11. p3, Eq. (6): Terms in this equation need defining properly in several regards. Firstly, \Delta GATP0 should be defined. Secondly, it should be made clear that concentrations are intracellular (as opposed to in Eq. (1) (3), which correspond to both extra-and intracellular compartments). This is particularly confusing for Pi in which the subscript could easily be misread as being an intracellular label. Finally, it is not clear how the equality \Delta GATP= qvATP holds. Perhaps, as in footnote (1), it could be made clear that this relies on the definition of the reversal potential for ATP.
Added square brackets and subindices to indicate that the intracellular concentrations of ADP, ATP, and Pi are intracellular. Also corrected a typo, there was an extra q in the term for $\Delta G_0$ 12. p3, Eq. (6): I would also here like to take the opportunity to raise a point about the general framework. The models here account for the energy required to move molecules with and against their concentration gradients. In cases of active transport, as indicated by Eq. (10), ATP is then hydrolysed to meet the energy deficit required to power the pump. This then means that the concentration of ATP, relative to ADP is then changing. Indeed, at maximal rates of exercise, the ATP concentration in skeletal muscle can fall by up to 20% of its initial value. Moreover, in certain cell types, local ratios between ATP and ADP affect the opening and closing of ion channels (e.g. KATP). Eq. (10) implies that ATP and ADP in the cell are constant, in spite of the necessity that it must change. One must, therefore, conclude that an assumption is being made that ATP changes are small relative to the total amount of ATP. This assumption should be stated, or a discussion around how changes to ATP concentrations effect Eq. (10) should be made. The currents in equations 35-39 are normalized by the membrane capacitance. Therefore, to avoid abuse of notation, added a note remarking the change in the letter for the functions, explaining it is due to normalization.
25. p9, after Eq. (39): Presumably, the author here means that the intracellular Ca2+ concentration changes (rather than the extracellular). This should be clarified. Also, the text below Eq. (40) appears to imply a distinction between free and bound Ca2+, but there is no buffering process described in the model, which may cause confusion.
Added sentences to clarify these points.
26. p9, Fig. 2: Please state what im and iCaL13 are in these plots.im represents the total current.
Added a sentence to the figure caption.
27. p10, Table 2: The units of \bar{a} are inconsistent with the main text.
Corrected by deleting the bars in the table for normalized amplitude entries. \bar{a} is in pA, a = \bar{a}/Cm in pA/pF 28. p12, "Another possible use of the third order approximations is in the construction of network models.": Please could the author state what they mean here?
Expanded the sentence to clarify the point. The idea is to use cubic approximations to build network models taking advantage of the possible reduction in computational cost. However, as is explained in the paragraph, cubic approximations need to be tested first to address possible issues that may arise from using approximations.
29. p13, discussion on action potential shape: I was wondering if the author could here compare their argument surrounding the shape of the action potential near initiation with those put forward here (Brette, 20131), which puts forward a case for multi-compartment models to describe the sharpness of the spike.Added a comment about the results in the paper by Brette, 2013.
Thank you for pointing it out. Added a reference to his work too.
30. SI1, after (A3): There appears to be a reference missing here.
Edited the reference, it is for Equation 31. 31. SI2: Headers are required in these data files so that the values herein can be related to the physical quantities they represent.
Edited data files to include headers.
32. It would be useful for the author to upload any code used to generate the plots in the article.
Uploaded the code in an jupyter notebook.

Competing Interests:
The author declares no competing interests.
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias • You can publish traditional articles, null/negative results, case reports, data notes and more • The peer review process is transparent and collaborative • Your article is indexed in PubMed after passing peer review • Dedicated customer support at every stage • For pre-submission enquiries, contact research@f1000.com