## Abstract

This paper seeks to develop a more thermodynamically sound pedagogy for students of biological transport than is currently available from either of the competing schools of linear non-equilibrium thermodynamics (LNET) or Michaelis–Menten kinetics (MMK). To this end, a minimal model of facilitated diffusion was constructed comprising four reversible steps: *cis-*substrate binding, *cis*→*trans* bound enzyme shuttling, *trans*-substrate dissociation and *trans*→*cis* free enzyme shuttling. All model parameters were subject to the second law constraint of the probability isotherm, which determined the unidirectional and net rates for each step and for the overall reaction through the law of mass action. Rapid equilibration scenarios require sensitive ‘tuning’ of the thermodynamic binding parameters to the equilibrium substrate concentration. All non-equilibrium scenarios show sigmoidal force–flux relations, with only a minority of cases having their *quasi*-linear portions close to equilibrium. Few cases fulfil the expectations of MMK relating reaction rates to enzyme saturation. This new approach illuminates and extends the concept of rate-limiting steps by focusing on the free energy dissipation associated with each reaction step and thereby deducing its respective relative chemical impedance. The crucial importance of an enzyme's being thermodynamically ‘tuned’ to its particular task, dependent on the *cis-* and *trans-*substrate concentrations with which it deals, is consistent with the occurrence of numerous isoforms for enzymes that transport a given substrate in physiologically different circumstances. This approach to kinetic modelling, being aligned with neither MMK nor LNET, is best described as intuitive non-equilibrium thermodynamics, and is recommended as a useful adjunct to the design and interpretation of experiments in biotransport.

## 1. Introduction

*Aim.* The purpose of this paper is to remedy pedagogical error in matters pertaining to the thermodynamics and kinetics of biological transport phenomena. This is done by using a minimal model of facilitated diffusion, based on the mass action law and constrained by the probability isotherm, to provide straightforward thermodynamic and kinetic insights that elude the competing approaches of Michaelis–Menten kinetics (MMK) and linear non-equilibrium thermodynamics (LNET). This new approach is aptly named intuitive non-equilibrium thermodynamics (INET) for what, it is hoped, will become obvious reasons.

It is a persistent and dreary oral tradition in the folklore of chemical energetics pedagogy that ‘thermodynamics has nothing to say about reaction rates except at equilibrium’. This quite erroneous myth persists despite occasional attempts to dispel it by establishing consistency between non-equilibrium kinetics and thermodynamics. A significant attempt was published by Boudart in the mid-1970s for the case of complex unbranched chemical reactions [1] and was extended a decade later by Wagg to include all chemical, osmotic and chemiosmotic reactions, however complex, whether branched or unbranched [2]. This line of enquiry was pursued further by Wagg and co-workers in the 1990s for theoretical analysis of membrane transport [3–5], although an important limitation affecting experimental application of such theory is the inability of radioisotopic fluxes to distinguish between the various entry and exit points of ions involved in branched transport mechanisms [6].

Meanwhile, another attempt to relate thermodynamics to non-equilibrium reaction kinetics was initiated by Porter in 1983 [7] who derived from the van't Hoff isotherm a curiously expressed formula relating a reaction's molar free energy change (Δ*µ*) to the net rate (*j*) and unidirectional forward rate (*J*) of the reaction, thus
*R* and *T* have their usual meanings. While we shall have no further use for this formulation—it will be seen that it is identical to the rate isotherm defined below as equation (2.3), though much less intuitively expressed—Porter's derivation was important because it set down the required condition for assumption of Boltzmann equilibrium in liquid-phase media that the reaction times be longer than 10^{−11} s. This condition, permitting application of the mass action law using equilibrium-determined rate constants into the non-equilibrium domain, is amply satisfied for processes occurring under biological conditions of temperature and pressure, where the second-order diffusion-limited binding constants are orders of magnitude slower than 10^{11} s^{−1}, and first-order protein conformational changes are orders of magnitude slower again (see table 1 and associated text).

Six years ago, three strands of pedagogical error in bioenergetics that have persisted since the 1970s were identified [8] and tracked historically [9] in papers that (i) identified the probability isotherm as an intuitive non-equilibrium thermodynamic framework for biochemical kinetics and (ii) stressed the importance of distinguishing between the frequently confounded concepts of *entropy creation* and *entropy exchange*. Last year, the probability isotherm was applied to offer an alternative interpretation of experimental data obtained on the F_{o}F_{1}-ATPase [10], echoing a much earlier instance of its application in 1980 to re-interpreting the kinetics of *n*, *m* and *h* in the Hodgkin–Huxley equations pertaining to electrical excitation in nerves [11]. While such application of the probability isotherm (though not under that name) is now commonplace in biological modelling of voltage gating mechanisms, it is rarer than it should be in relation to modelling of osmotic and chemiosmotic systems in biomembrane transport. It is perhaps a reflection of this rarity that the authors of the more recent application of the probability isotherm to the F_{o}F_{1}-ATPase [10] were pleasantly surprised by being prevailed upon in the review process to provide a pedagogical appendix on the distinction between entropy creation and entropy exchange and its relation to thermodynamic efficiency.

Thus, there would seem to be a need to develop a more robust pedagogy for scientists at the cutting edge of research into biological transport. The present paper uses a numerical simulation of a minimal model of facilitated diffusion to illustrate how such a pedagogy might be developed.

## 2. Technical definition of the probability isotherm and the rate isotherm

The probability isotherm affords a thermodynamic definition of a *probability ratio* for forward and backward reaction, *p*_{f}/*p*_{b}, thus
*G*_{Diss} is the molar free energy dissipation of the reaction. It should be noted that the unidirectional probabilities are individually undefined; only their *ratio* is defined according to the second law of thermodynamics as expressed in equation (2.1). As shown elsewhere [8], the probability ratio given in equation (2.1) is the product of an *intrinsic* probability ratio (equal to the equilibrium constant, *K*_{eq}) determined by the *nature* of the reactants and products, and an *extrinsic* probability ratio determined by the *composition* (concentrations of reactants and products). As with all thermodynamic relationships, equation (2.1) is totally aloof from mechanism, so much so that it remains true regardless of whether or not any reaction mechanism actually exists.

For example, the oxidation of glucose according to
*in vivo* conditions in respiring animals and plants, the value is somewhat less at around 2.85 MJ mol^{−1}. The probability isotherm determines a probability ratio for forward to reverse directions of reaction (2.2) at 37 K *in vivo* equal to *e*^{1108}, a number with a magnitude beyond the RAM capacity of today's personal computers; its value might be made more accessible to current computing technology and human appreciation by representing it numerically as the alternative expression

However, this does not mean that glucose (or a glucose-containing organism) forms an explosive mix with air. Fortunately, there is no known mechanism for the direct oxidation of glucose by molecular oxygen; while glucose may be heated in air, only residual carbon is oxidized following thermal dehydration. Nonetheless, the probability ratio for glucose oxidation is thermodynamically defined according to equation (2.1) and sets a benchmark into which all thermodynamically determined probability ratios for the known partial reactions of intermediary metabolism of glucose must ‘fit’. While this ‘fitting’ requirement often goes under the quaint name of ‘detailed balance’, the present author prefers the ‘second law’ designation of what is essentially a thermodynamic requirement.

In cases where there are known mechanisms for reaction, the molar free energy dissipation also determines a rate isotherm giving a *rate ratio* for forward and backward reaction, *r*_{f}/*r*_{b}, identical to the probability ratio determined by the probability isotherm, thus

In this case, the unidirectional rates, unlike the unidirectional probabilities in equation (2.1), are actually defined and are quantitatively determined according to whatever level of catalysis might be present for any given composition of reactants and products. However, while the level of catalysis may influence the absolute values of the unidirectional rates, it cannot possibly influence their *ratio* which is always and everywhere determined thermodynamically. Even though equation (2.3) is a kinetic expression of the more fundamental equation (2.1), its thermodynamic message is that the level of catalysis cannot influence the rate ratio any more than it can influence the position of thermodynamic equilibrium. In this sense, therefore, thermodynamics has a great deal to say about non-equilibrium kinetics, and equation (2.3) affords a quantitative constraint that can be brought into play to enhance the design and interpretation of studies of membrane transport in which the thermodynamic boundary conditions are under experimental control.

Equations (2.1) and (2.3) apply at all times to all chemical, osmotic and chemiosmotic reactions; in the case of complex reactions, they apply equally to overall reactions and to individual reaction steps. As mentioned above, the rate isotherm expressed in equation (2.3) is equivalent to Porter's formulation shown in equation (1.1), where Porter's Δ*μ* is equal to the rate isotherm's –Δ*G*_{Diss} and his expression in parenthesis is equal to *r*_{b}/*r*_{f}, i.e. the reciprocal of the corresponding expression in the rate isotherm. However, the rate isotherm is more intuitively expressed in that it relates an increase in molar free energy dissipation directly to an increasing ratio of forward to backward rates of reaction.

## 3. A minimal model of facilitated diffusion (uniport)

Consider the osmotic reaction
*cis-*) and 2 (*trans-*), respectively. The translocation mechanism involves a transmembrane uniporter enzyme, E, having two conformational states, E_{1} and E_{2}, respectively, and catalysing the reaction according to the following mechanism comprising four elementary steps:

where *k _{i}* and

*k*

_{−i}are the respective forward and backward rate coefficients of the

*i*th elementary step.

For a passive process of facilitated diffusion, the standard molar free energy change of the overall reaction,

In order to work with a model that can be related to the actual performance of real transport processes reported in the biomedical literature, we shall aim to ‘design’ a uniporter enzyme that can achieve a membrane transport density of the order of 1 pmol cm^{−2} s^{−1} with a maximum molecular turnover number of 1000 s^{−1} and requiring an enzyme site density of around 1 fmol cm^{−2}. In accordance with the simpler enzymatic function of a uniporter, this turnover number is roughly ten times that reported for the more complex chemiosmotic Na^{+},K^{+}-ATPase [12] while being comparable to that found for the glucose transporter 1 (GLUT1) that mediates facilitated diffusion [13].

### 3.1. Assumptions about the rate coefficients

#### 3.1.1. Enzyme substrate binding

We assume that the second-order association constants *k*_{1} and *k*_{−3} are both 10^{7} M^{−1} s^{−1}, i.e. a ‘middle-order’ value in the commonly encountered range of 10^{6}–10^{8} M^{−1} s^{−1} for diffusion-limited ligand–receptor binding [14]. This means that the first-order dissociation constants *k*_{−1} and *k*_{3} will be thermodynamically determined by the standard free energies,

#### 3.1.2. Conformational changes (translocation)

We assume that the first-order turnover numbers *k*_{2}, *k*_{−2}, *k*_{4} and *k*_{−4} are limited to being no greater than 10^{4} s^{−1} (corresponding to the upper limit of reported protein conformational turnover numbers). When exploring the effects of variation in the standard free energies,

Table 1 thus shows the kinetic and thermodynamic relationships of the parameters for each of the reaction steps. It may be seen that, within the stated constraints on the upper limit for the rate coefficients of either association or translocation, specification of a numerical value for the standard free energy, *i*th step will determine the absolute values of the forward and reverse rate coefficients, *k _{i}* and

*k*, for that step.

_{−i}## 4. Simulation methods

### 4.1. Steady-state equations

The following four equations for the rates of change in concentration of the respective four forms of the enzyme apply in the steady state:

Of these four equations in four unknowns, only three are independent. The fourth independent equation required for steady-state solution is the conservation equation for the total concentration of all forms of the enzyme, thus

### 4.2. Mathematical methods

These steady-state equations were solved numerically using the algorithms built into Microsoft Excel, the resulting concentrations of the four forms of the enzyme then being used to calculate the unidirectional rates of each of the four reaction steps as given in table 2. The formulae for calculating the overall unidirectional rates were derived according the method given by Boudart [1], assuming a value of unity for the stoichiometric number (Boudart's *χ _{i}*) of each of the reaction steps 1–4, as given in the ‘overall’ row of table 2. The concentration of total enzyme was set at 1 fmol cm

^{−2}of membrane between compartments 1 and 2, comparable to the density of Na

^{+},K

^{+}-ATPase in plasmalemmal membranes reported in the literature [15]. The Excel spreadsheet (INET.xlsx) used for the calculations reported in this paper is supplied as the electronic supplementary material .

### 4.3. Choosing and sweeping the independent variable(s); display of data

The studies described here examine the influence of (i) the affinity of the enzyme for the substrate, (ii) the absolute concentration of the substrate and (iii) the concentration gradient of the substrate. This was done by using ToolBook to set the respective boundary conditions and sweep the independent variable(s) by passing them as input variables to Excel via dynamic data exchange. The primary solutions of the four linear equations and the secondary calculations of rates, enzyme saturation and free energy dissipation were all performed within Excel, and selected results were extracted by ToolBook via dynamic data exchange and plotted graphically within ToolBook. ToolBook's graphical displays were saved as screen captures and converted to PDFs for printing.

Using ToolBook to pass parameters to Excel and to retrieve solutions therefrom through dynamic data exchange was much faster than performing the primary and secondary operations within ToolBook, provided that the open Excel file was kept minimized during any given simulation run. As recorded in the Acknowledgements, the accuracy of these unorthodox numerical methods was initially cross-checked using Matlab by colleagues at the University of Auckland. Once this accuracy was thus confirmed, all further virtual experimental scenarios were simulated using this ToolBook–Excel regime. The heuristic and pedagogical advantages of constructing the numerical simulation in Excel should become obvious on perusal of the INET.xlsx file provided.

## 5. Numerical simulation of model behaviour

### 5.1. Rapid equilibration across a barrier—effects of binding affinity

Rapid equilibration across biomembranes is not commonly found for transport of ions or metabolites, but it is common for water transport subserved by several different aquaporin proteins. Nonetheless, our heuristic purposes are well served by examining the factors that optimize the ability of a uniporter to facilitate rapid equilibration. The design specification for such a uniporter might reasonably be supposed to involve (i) perfect symmetry for the *cis-* and *trans-*binding affinities (i.e.

We begin by sweeping the standard free energies of binding for steps 2 and 4 identically over the domain of –3*RT*ln(10) (high affinity) to +3*RT*ln(10) (low affinity) with unit substrate concentration (1 M) in both compartments. It is important to note that the thermodynamic constraint on the rate coefficients for the binding steps was realized by fixing the second-order diffusion-limited association constants *k*_{1} and *k*_{−3} at 10^{7} M^{−1} s^{−1} while allowing the enzyme-determined dissociation constants *k*_{−1} and *k*_{3} to take on various values constrained by the formulae given in table 1. This is equivalent to simulating the behaviour of a continuum of enzyme isoforms, each isoform having a different enzyme-specific dissociation constant.

The plots displayed in figure 1 show that, with equal *cis-* and *trans*-concentrations of substrate, the percentages of the two bound forms of the enzyme superimpose at all levels of binding affinity, as do the percentages of the two free (unbound) forms. Because all the conditions are at equilibrium, the forward and backward rate curves also superimpose, while the net rate is zero throughout. In this specific case, where the substrate concentration is unity (1 M) in both compartments, the enzyme is exactly half saturated at zero standard free energy of binding, with each of the four enzyme forms being equally represented at 25%. Curiously, the maximum unidirectional rates, *not* achieved at exactly this point of 50% saturation with unit substrate concentration each side of the membrane, but at a saturation level of 49.9750%, with ^{−4}*RT*ln(10) or +2.5776 J mol^{−1}, such that *k*_{−1} = *k*_{3} = 1.0010 × 10^{7} s^{−1}. This discrepancy between the affinity level corresponding to exactly 50% enzyme saturation and that corresponding to maximum unidirectional rates becomes progressively larger as the equilibrium substrate concentration is lowered and is owing to the fact that, while 50% saturation exactly follows the reciprocal changes in binding affinity and concentration, the denominator used for calculating the overall unidirectional rates (table 2) does not show such constancy in the face of reciprocal variation in substrate concentration and binding affinity.

This effect is illustrated in figure 2 which shows the influence of binding affinity on equilibrium saturation levels (figure 2*a*) and unidirectional forward rates (figure 2*b*) at fixed levels of equilibrium substrate concentration ranging from 10^{−6} to 10 M. The conditions determining the maximum equilibrium unidirectional rate (

Note that the maximum unidirectional rate achieved by different binding affinities is relatively constant over a substrate concentration domain of 10 M to 10 mM and is achieved at saturation levels of around 50% or slightly less; however, over the substrate concentration domain of 1 mM to 1 µM, *RT*ln(10) unit of change in

To assist with a perspective on the simulations recorded in figures 1 and 2, the continuous variation of

### 5.2. Rapid equilibration across a barrier—effects of substrate concentration

The results shown in figure 2 and table 3 confirm that, for any given equilibrium substrate concentration, there exists a unique binding affinity at which the equilibrium unidirectional rates are maximal. The complementary expectation—that, for any given binding affinity, there exists a unique equilibrium substrate concentration at which the equilibrium unidirectional rates are maximal—is confirmed by the results displayed in figure 3 and table 4, showing the effect of sweeping the equilibrium substrate concentrations over the domain of 1 µM to 10 M with the free energy of binding fixed at exact integer multiples, *m*, of *RT*ln(10), where –6 ≤ *m* ≤ 1. The curves shown in figure 3 thus pertain to eight different virtual isoforms of the uniporter.

A particularly interesting observation arising from these ‘complementary’ experiments is that, at each *RT*ln(10) unit of variation in the binding affinity or 10-fold variation of the substrate concentration, the equilibrium value of *R*_{max} is identical (within the limits of numerical precision) even though the enzyme saturation level at which it is achieved is different, with the differences in saturation becoming very large at low substrate concentrations. Moreover, there is a clear ‘symmetry’ between the optimal substrate concentration and the optimal binding affinity to achieve *R*_{max} at each decade of variation, as reflected by the essentially identical numerical value of the logarithm of optimal substrate concentration (table 4, column 3) and the number of *RT*ln(10) units of binding affinity (table 3, column 2). This ‘symmetry’ also extends to the deviation from 50% enzyme saturation for the *R*_{max} achieved at each decade of variation: the optimal saturation levels recorded for each decade in tables 3 and 4 deviate by identical amounts from 50% such that the average enzyme saturation recorded in the two tables for *R*_{max} is exactly 50% at every decade, regardless of how far the individual optimal saturation deviates from 50%.

### 5.3. Rapid equilibration across a barrier—disclosure of isokines

The symmetry observed between the sets of data presented in figures 2 and 3 and tables 3 and 4 suggests that there exist pathways in the concentration–affinity space upon which the equilibrium unidirectional rate remains constant; we shall call such pathways ‘isokines’. The isokine pertaining to data rows 4 of tables 3 and 4 (i.e. for log[A]_{1} = log[A]_{2} = –2(M) in table 3 and for *a* over the domain −1.8,−1.8 to −2.0,−2.0 (figure 4*a*) and over the domain 0,0 to −3.0,−3.0 (figure 4*b*, red curve). Figure 4*b* also shows isokines pertaining to rows 3 (blue curve) and rows 5 (green curve) of tables 3 and 4. The data for figure 4*a* were obtained by stepping log[A]_{1} = log[A]_{2} over the required domain and, at each step, stepping *R*_{max} = 1.1387 pmol cm^{−2} s^{−1} recorded in data rows 4 of tables 3 and 4. While the methods of generating the isokines in figure 4*a,b* were identical, the tolerance used for the curve in figure 4*a* was 10^{−4}% while that used to generate the curves in figure 4*b* varied between 10^{−4}% and 0.05% according to the required grid resolution.

The isokine shown in figure 4*a* has the appearance of a parabola tilted at 45° to the perpendicular, although it is clearly not parabolic (figure 4*b*, red curve). As one proceeds clockwise round the left-hand red isokine, the enzyme saturation increases monotonically from 42.52% to 57.48%. For the right-hand red isokine, the saturation increases monotonically clockwise from 35.25% to 64.75%. For each of the red isokines, the point corresponding to exactly 50.00% saturation occurs at the point of intersection of each curve with its axis of symmetry (the 45° line given by *R*_{max} for the virtual isoform for which *R*_{max} = 1.1410 pmol cm^{−2} s^{−1} when the equilibrium substrate concentration is 11.19 mM at 52.23% saturation.

Similar behaviour is seen for the blue and green isokines of figure 4*b*. For the blue isokine (Rate = 1.2377 pmol cm^{−2} s^{−1}), as one proceeds clockwise the enzyme saturation increases monotonically from 45.41% to 54.59%, with exact 50% saturation occurring for a virtual isoform having *R*_{max} = 1.2379 pmol cm^{−2} s^{−1} when the equilibrium substrate concentration is 103.32 mM at 50.24% saturation.

For the green isokine (Rate = 0.6699 pmol cm^{−2} s^{−1}), the enzyme saturation increases monotonically clockwise from 16.00% to 84.00%, with exact 50% saturation occurring for a virtual isoform having *R*_{max} = 0.7105 pmol cm^{−2} s^{−1} when the equilibrium substrate concentration is 1.91 mM at 62.30% saturation.

The concept of theoretical isokines developed here demonstrates that different uniporter isoforms are capable of operating at similar absolute rates in the face of different substrate concentrations found in different environments. This kind of concept could be explored experimentally by comparing and contrasting different uniporter enzymes inserted in lipid vesicles, as is becoming increasingly possible through the evolving techniques for micro-measurement of fluxes and conformational transitions of enzymes.

It is worth noting that the concentration/affinity isokine and its associated continuum of variation of enzyme saturation are theoretical concepts of potentially significant experimental consequence that find no resonance in pedagogies deriving from LNET or MMK.

### 5.4. Rapid equilibration across a barrier—effects of translocation affinity

The data produced so far have all been obtained on the assumption that there is no intrinsic translocational preference between *cis-* and *trans-*conformations for either the bound or unbound forms of the enzyme (i.e.

Figure 5 shows the variation of equilibrium unidirectional rate, enzyme saturation and enzyme distribution when *RT*ln(10) to +3 *RT*ln(10) under three different sets of conditions. The data of figure 5*a* were obtained under the same conditions that yielded *R*_{max} = 1.1387 pmol cm^{−2} s^{−1} in figure 2 and row 4 of table 3, i.e. [A]_{1} = [A]_{2} = 10 mM and _{2} forms of the enzyme while positive values for _{1} forms of the enzyme. These variations are reciprocal in that the level of enzyme saturation remains constant throughout, despite the large variation in unidirectional rate. Figure 5*b* shows the ‘complementary’ results obtained under conditions yielding the same *R*_{max} = 1.1387 pmol cm^{−2} s^{−1} in figure 3 and row 4 of table 4, i.e. [A]_{1} = [A]_{2} = 10.9544 mM and *c* shows the results obtained for the suboptimal intermediate condition in which [A]_{1} = [A]_{2} = 10 mM and *R*_{max} = 1.1364 pmol cm^{−2} s^{−1}; under these conditions the two conformational states of the enzyme are exactly half distributed between the bound and free forms.

The ‘symmetry’ already noted in relation to the data of tables 3 and 4 and figure 4 is reinforced by the graphical data shown in figure 5*a*,*b* where there is exact numerical equality between figure 5*a* solutions for the percentages of the enzyme forms E_{1}, E_{1}A, E_{2}A and E_{2} and figure 5*b* solutions for the percentages of E_{1}A, E_{1}, E_{2} and E_{2}A, respectively.

Within the functional context of rapid equilibration between compartments, the results shown in figure 5 invite the teleological supposition that there is a significant disadvantage in a uniporter's having a preferred conformation; this is because the rate-limiting constraint inherent in protein conformational changes becomes amplified on one of the respective unidirectional rate coefficients as the standard free energy of conformational change departs from zero.

We now turn our attention to non-equilibrium simulations involving facilitated diffusion down a concentration gradient.

### 5.5. Force–flux relations

As noted earlier, it is a central expectation of LNET [16] that, in the absence of complicating factors, the force–flux relation of a process should be linear, provided that the process is occurring not too far from equilibrium (a proviso that is commonly very elastic in its application). While it is a trivial expectation that the force–flux relation of a saturable system such as a uniporter enzyme will be sigmoidal, with a quasi-linear inflection point occurring between the saturating asymptotes in either direction, it is not clear as to whether the inflection point will always be found at or near equilibrium. The next simulations address this question by examining the model's kinetic behaviour under non-equilibrium conditions in which log[A]_{1} is swept over the molar domain of −6 to +1 at eight fixed integer values of log[A]_{2} over the same domain for three different fixed values of binding affinity.

Figures 6*a*–*f* and 7*a,c,e* show plots of enzyme saturation (figures 6*a,b* and *7a*), unidirectional forward rate, A_{1→2} (figures 6*c,d* and *7c*) and net reaction rate (A_{1→2} – A_{2→1}, figures 6*e,f* and *7e*) as functions of *cis-*substrate concentration, [A]_{1}, for eight different fixed values of *trans-*substrate concentration, [A]_{2}, with *a,c,e*), –2 (figure 7*a,c,e*) and –3 (figure 6*b,d,f*) times *RT*ln(10).

Each graph of net rate versus log[A]_{1} shown in figures 6*e,f* and 7*e* are thermodynamic force–flux relations that all fulfil the intuitive expectation of being sigmoidal and showing saturating behaviour in both directions. These results thus confirm and extend the earlier demonstration by Chapman & Loiselle [10] that linear force–flux relations are not generally to be expected. Of all the non-equilibrium force–flux relations simulated here, only those obtained with *RT*ln(10) times the respective equilibrium value of log[A]_{1} = log[A]_{2} have their *quasi*-linear inflection points approximately centred at equilibrium; all other force–flux relations have their inflection points centred both away from equilibrium and increasingly displaced to the right of the value of log[A]_{1} = log[A]_{2} equal to the multiple of *RT*ln(10) for the respective

The data plotted in figures 6 and 7*a,c,e* show the results obtained from three different virtual isoforms (three columns), differing from each other by their respective binding affinities. Nonetheless, the bottom panel of each column shows data pertaining to a single virtual isoform, demonstrating that the *cis-*force–flux relation for a particular enzyme is extremely variable, depending on the *trans-*substrate concentration. This, of course, is a trivial observation, but it raises doubt about the interpretation to be placed upon the shape of a force–flux relationship and its continuously variable ‘phenomenological coefficient’ expressing the continuously changing ratio of the flux to its respective force. By contrast, consider the interpretation to be placed on the eight different equilibrium slopes of the lines shown in figures 6*e,f* and 7*e*, and the fact that these slopes are apparently maximal when the binding affinity of the enzyme is closely ‘tuned’ to its equilibrium substrate concentration (cf. figures 2 and 3). Moreover, given that linear force–flux relations do not occur for this simplest of models of facilitated diffusion, it does not seem that linear force–flux relations are likely to be generally encountered elsewhere in experimental studies of biotransport mechanisms.^{1} Indeed, the application of LNET to chemistry, biochemistry and molecular physiology should never have been admissible. While linear force–flux relations have clear application to the work of electricians and plumbers, they have no place in any domain governed by the law of mass action. But, of course, it is far worse than that; the saturability of enzyme-catalysed reactions is yet another source of nonlinearity in biological force–flux relations, compounding the already inescapable nonlinear consequences of the law of mass action. Proponents of LNET might claim that the present results vindicate their expected near-equilibrium linearity for enzymes that are teleologically tuned to their respective substrate concentrations; but readers of LNET-inspired speculations are not generally informed about any kind of clearly defined or openly admitted domain of near-equilibrium validity.

Two features of some of the net rate curves shown in figures 6*e,f* and 7*e,f* warrant some cautionary remarks, particularly in relation to the uppermost clusters of curves arising from close to zero net rate. If such curves were actual experimental data, they might be taken as evidence of asymmetric diffusion on the one hand, or enzyme activation on the other.

Asymmetric diffusion is a difficult concept that has achieved some traction in the ‘origin-of-life literature’ [17] but is of no relevance to the present studies which are predicated on the applicability of the probability isotherm to any overall transport reaction and to its individual steps. In all the studies reported here, the overall equilibrium constant is unity for every condition simulated, meaning that the overall ratio of forward to backward transport rates is always equal to the *cis *: *trans* concentration ratio.

This leads immediately to the apparent depiction of enzyme activation in these same curves as they seem to rise *quasi*-exponentially from near-zero to reach their maximum slopes at their points of inflection. The portions of such curves lying to the left of their inflection points are similar in shape to those reported elsewhere from studies of the ΔpH-dependence of the rate of ATP synthesis by the F_{o}F_{1}-ATP synthase [18–20]. The authors of these studies were influenced by the LNET school, leading them to require an explanation for the fact that the curves were nonlinear; the explanation offered was that the pH gradient itself might serve to activate the enzyme. In the context of the present simulations depicted in figures 6 and 7, there is no enzyme activation because all eight rate coefficients are constant *within* any given curve, and they are also constant *between* the eight curves present in any given panel of graphs. The various shapes of all the curves shown in figures 6 and 7 are due solely to the playing out of the law of mass action as determining the unidirectional rates of each step and of the overall reaction in a saturable system. Upward concavity of force–flux relations in saturable systems, therefore, cannot be taken as evidence for enzyme activation in the absence of other kinetic and thermodynamic information.

### 5.6. Saturation kinetics and the Michaelis–Menten *K*_{m}

Michaelis–Menten kinetic modelling has long afforded a useful basis on which to develop some simple principles of catalysis and to allow analytical solutions of various enzymatic scenarios to be obtained. The main problem with such models is that they are unable to afford any useful reconciliation between kinetics and thermodynamics, except at equilibrium. Indeed, the only so-called ‘thermodynamic’ equation that has appeared in these contexts with any regularity has been the van't Hoff isotherm, an expression more appropriately described as a ‘thermostatic’ equation. However, its thermodynamic equivalent form, the probability isotherm [8,10], leads directly to the powerful thermodynamically constrained kinetic equation, appropriately called the rate isotherm [8] and which is immediately available for use in thermodynamically constrained computer simulations of enzyme kinetics such as that developed here.

Figure 7 shows logarithmic (figure 7*a,c,e*) and linear (figure 7*b,d,f*) plots of the dependence on *cis*-substrate concentration, [A]_{1}, of enzyme saturation (figure 7*a*,*b*), unidirectional forward rate (A_{1→2}, figure 7*c*,*d*) and net reaction rate (A_{1→2}–A_{2→1}, figure 7*e*,*f*) for eight different values of *trans-*substrate concentration, [A]_{2}, with *RT*ln(10). In these graphs, the only curves that can be related to MMK are those for the unidirectional rate shown in figure 7*c*,*d*. In these instances, the logarithmic plots (figure 7*c*) show that the ‘*K*_{m} value’ for half-maximal unidirectional rate is relatively constant at approximately 10 mM [A]_{1} for all fixed values of [A]_{2}. However, this does not correspond to any kind of association with 50% saturation of the enzyme except for the case of [A]_{2} = 10 mM, as is evident in the logarithmic saturation plots. Moreover, the concept of half-maximal rate finds resonance in only the *net* rate curves for the cases of −6 ≤ log[A]_{2} (M) ≤ −3.

While the data plotted in figure 7*c* suggest at least a near concordance with MMK for the *cis*-substrate concentration at which half-maximal unidirectional rates occur, this apparent ‘*K*_{m}’ also varies somewhat with the fixed *trans*-substrate concentration at which the value is determined, and these situations do not correspond to 50% enzyme saturation. This is demonstrated by the data shown in table 5 for the values of [A]_{1} at which the half-maximal unidirectional rates occur. These values were determined with *RT*ln(10) for two different values of *k*_{4} = *k*_{−4}. The ‘*K*_{m}’ was determined by forcing *R*_{max} on the enzyme by setting [*A*]_{1} = 10^{5} M and then finding the respective value of [A]_{1} to yield 0.5 *R*_{max}. Thus, the apparent *cis*-‘*K*_{m}’ varies with both *trans-*substrate concentration and with the relative rate-limitation among the steps.

This kinetic behaviour generally lacks any useful correspondence with classical MMK. Although the saturation behaviour is roughly comparable, it finds no resonance in the corresponding rate behaviour as shown in figure 7. This suggests that Michaelis–Menten analysis is of limited value for understanding the physiological function of transporter enzymes either *in situ* or in the increasingly sophisticated nano-experimental protocols being used today. In particular, it seems unlikely that the classical Michaelis–Menten framework will be any more useful than the expectations of LNET in assisting with the design and interpretation of experiments in this expanding field.

### 5.7. Effect of *trans*-substrate concentration on *cis*-substrate transport

It is also shown in figure 7 that increasing the concentration of *trans*-substrate is inhibitory to both the unidirectional rate of *cis*-substrate transport and the net rate. However, this is not necessarily in conflict with the well-known observation that *trans*-substrate can be stimulatory towards unidirectional isotopic flux of *cis*-substrate, a phenomenon often referred to as *trans-*acceleration [13]. This is because the experimental measurement of unidirectional flux of *cis*-substrate is made, not across all of reaction steps 1–4, but only across steps 1–3. This distinction may be confirmed in the electronic supplementary material (the Excel file ‘INET.xlsx’) by contrasting the formula of cell F7 (overall unidirectional forward rate) with that of cell F8 (labelled ‘Isotopic A Influx across Steps 1 to 3’). Figure 8 shows the virtual experimental unidirectional flux rate, A_{1→2}, as a logarithmic function of *cis-*substrate concentration, [A]_{1}, for different values of *trans-*substrate concentration, [A]_{2}, with *RT*ln(10) and *k*_{4} = *k*_{−4} set to 10^{4} s^{−1} (figure 8*a*) and 10^{3} s^{−1} (figure 8*b*). These graphs show the effects of two different degrees of rate limitation of the conformational change in step 4 relative to that in step 2. In figure 8*a*, where *k*_{4} = *k*_{−4} is set to 10^{4} s^{−1}, increasing the *trans*-substrate concentration has a small *negative* effect on unidirectional flux, while it has a marked *positive* effect if *k*_{4} = *k*_{−4} is set to 10^{3} s^{−1} (figure 8*b*). Thus, when there is no relative rate limitation of step 4, increasing the *trans*-substrate concentration is slightly inhibitory to the flux of *cis-*substrate, as shown by the small reductions in *R*_{max} as the *trans*-substrate concentration is raised from 1 µM to 10 M (figure 8*a*). On the other hand, when the rate-limitation of step 4 is present, this flux is greatly stimulated as *trans-*substrate is raised over the same domain (figure 8*b*). Note that the curves for the three highest *trans*-substrate concentrations form the lower (almost superimposed) cluster in figure 8*a*, while they form the upper cluster in figure 8*b*.

These results confirm that, in principle, *trans*-substrate can be either inhibitory or stimulatory on *cis*-substrate flux, depending on the presence or absence of rate limitation in the translocation of free enzyme relative to translocation of bound enzyme. When there is no relative rate limitation on the translocation of free enzyme, step 4, then *trans*-substrate competes for free enzyme on an equal footing with *cis*-substrate through the rapidly reversible step 3; under these conditions, *trans*-substrate can reduce the unidirectional flux of *cis*-substrate by reducing the availability of unbound E_{2} to participate in step 4 as the *trans*-substrate concentration increases. However, when step 4 is rate limiting, the quickest way to regenerate E_{1} from E_{2} to bind *cis*-substrate in step 1 is via reverse unidirectional fluxes of steps 3, 2 and 1 rather than from forward unidirectional flux of step 4; under these conditions, *trans*-substrate will increase the unidirectional flux of *cis*-substrate by increasing the availability of unbound E_{1} to participate in forward unidirectional flux through step 1.

It is worth noting that the *trans*-acceleration demonstrated for the present model arises purely from the relative rate-limitation of step 4 relative to step 2. Thus, while elsewhere it has been suggested that tetramer formation by the glucose transporter GLUT1 might be ‘important for understanding the *trans* acceleration observed for GLUT1-mediated transport’ [13], the present studies show that this is not a necessary condition for *trans*-acceleration.

### 5.8. Dissipation of free energy determined by the rate isotherm

The rate isotherm determines the kinetics of the actual reaction rates according to the second law and may be applied as such not only to the overall uniport process comprising steps 1–4, but also to the individual steps themselves. Thus, the free energy dissipation of the *i*th individual step, *i*th step. For the facilitated diffusion represented by steps 1–4, we have, at all times

Figure 9 shows the results of simulations in which *RT*ln(10) and the *trans-*substrate concentration, [A]_{2}, fixed at 10 mM while the *cis*-substrate concentration, [A]_{1}, was swept over the domain from 1 µM to 10 M. Figure 9*a,c,e* and *b,d,f* differs in that the rate coefficients for step 4 were set at 10^{4} s^{−1} (figure 9*a,c,e*) and 10^{3} s^{−1} (figure 9*b,d,f*). Figure 9*a,b* shows the variation in proportions of the four forms of the enzyme, including the total saturation, while figure 9*c,d* shows the overall unidirectional and net rates of transport. Figure 9*e,f* shows solutions to equation (5.6) as straight 45° lines indicating the total free energy dissipation, while solutions to equations (5.2)–(5.5) are shown as differently shaded areas indicating their respective contributions to the total free energy dissipation. Note that, in this context, negative free energy dissipation for overall transport from *cis*-compartment 1 to *trans-*compartment 2 simply means that the net transport is proceeding exergonically in the *reverse* direction whenever [A]_{1} < [A]_{2} as occurs for all values of [A]_{1} < 10 mM for the conditions represented in figure 9.

The dissipation plots shown in figure 9*e,f* demonstrate that there is no such thing as a unique rate-limiting step, even in this irreducibly simplified model of uniport. All that can be claimed with any consistency is that, for the conditions simulated, the dissociation of E_{2}A in step 3 is the least rate-limiting of the four steps (indicated by very little free energy dissipation, *e,f*). Although the rate-limitation of step 4 forced in figure 9*b,d,f* resulted in much reduced reaction rates (figure 9*c,d*), the overall *molar rate* of free energy dissipation (figure 9*e,f*, 45° lines) is unaltered.

Nonetheless, the concept of ‘rate-limiting’ steps can be usefully quantified if the four reaction steps are regarded as chemical impedances in series, across which the chemical potential drops in steps, by analogy with the stepwise drop in electrical potential that occurs when electric current flows through a set of resistors in series. This idea is given graphical representation in figure 10, where the absolute free energy dissipations shown in figure 9*e,f* are re-drawn as *proportions* (i.e. percentages) of the total free energy dissipation determined by the *cis *: *trans* concentration ratio for substrate A. Figure 10*a* corresponds to figure 9*e*, i.e. that for which *k*_{4} = *k*_{−4} = 10^{4} s^{−1}, while figure 10*b* corresponds to figure 9*f*, i.e. that for which *k*_{4} = *k*_{−4} = 10^{3} s^{−1}.

As free energy dissipation at each step is an indicator of that step's chemical impedance, the proportional plots of figure 10 also provide a useful indication of the relative presence of rate-limiting behaviour for each of the four steps. As is to be expected from the rate coefficients shown in table 1, the translocation steps 2 and 3 are the most rate-limiting (highest impedance), except for step 1 at low *cis*-substrate concentrations. The main conclusion to draw from the free energy dissipation data shown in figures 9 and 10 is that the relative chemical impedance (rate-limitation) of any given step in an unbranched reaction sequence is a continuously variable function of substrate concentration, as are the unidirectional and net reaction rates.

Application of the probability isotherm affords a more realistic framework for developing the concept of rate-limiting steps for transport enzymes operating *in situ*. As demonstrated by the data shown in figure 10, the degree of rate limitation occurring in any given step can be quantified in terms of its instantaneous molar free energy dissipation as a proportion of that for the overall reaction.

### 5.9. Designing an optimal uniporter for a particular concentration gradient

The thinking behind the simulations in this section is unashamedly teleological, stimulated by the findings reported in the preceding sections. Taking overall unidirectional rate as a proxy for catalytic activity, it has been shown that the activity of a uniporter is highly sensitive to binding affinity, to translocational affinity and to substrate concentration, and that there is no simple correlation between the maximal unidirectional rate, *R*_{max}, and enzyme saturation under a variety of scenarios. There is a functional advantage for cellular economy in optimizing *R*_{max} for a uniporter through the ‘design’ of its standard free energies of binding and translocation rather than through the membrane density of the enzyme. The data shown in figures 2 and 3 indicate that different concentrations of substrate require different thermodynamic and kinetic parameters for the enzyme to be optimally ‘tuned’ to any given situation, while the data shown in figure 5 indicate that any significant translocational affinity would be kinetically disadvantageous. Because any significant asymmetry of *cis*- and *trans-*binding affinity would require a compensatory asymmetry in the translocational affinity, resulting in kinetic disadvantage, we shall proceed on the assumption of no asymmetry in binding affinity and simply search for the optimal symmetric binding affinity for transport down four different concentration gradients.

The results of these simulations are shown in figure 11 for enzyme saturation (left axis and dashed lines) and net transport rate (right axis and solid lines) plotted against the standard free energy of binding. The results are shown in four different colours, corresponding to four different *cis*-[A]:*trans*-[A] concentration gradients as given in the legend. The coordinate points for optimal rate (*R*_{max}) and the levels of binding affinity and enzyme saturation at which they are achieved, are all marked on figure 11 and recorded in table 6. These data demonstrate that, as for rapid equilibration between compartments, transport down a concentration gradient requires fine tuning of the binding affinity for optimum performance. Unsurprisingly, the optimal binding affinity for a particular concentration gradient is found to lie between the optimal affinities for rapid equilibration of the respective concentrations (figure 2 and table 3), and is generally less than the mean of the two affinities.

### 5.10. Factors contributing to the kinetic performance of a uniporter

The simulations described above demonstrate that optimal ‘tuning’ of a simple four-step uniporter enzyme for a particular task, e.g. uniport down a relatively constant concentration gradient, involves the possible adjustment of a complex set of thermodynamic and kinetic parameters, thus

(a) Thermodynamic parameters:

— Affinity of

*cis*-binding at step 1— Intrinsic probability ratio for the two bound states of the enzyme at step 2

— Affinity of

*trans*-binding at step 3— Intrinsic probability ratio for the two free states of the enzyme at step 4

— Binding imbalance between steps 1 and 3

— Translocation imbalance between steps 2 and 4

(b) Kinetic parameters:

— Assumed upper limit for the second-order diffusion-limited association coefficients (

*k*_{1}and*k*_{−3})— Assumed upper limit for the first-order translocation constants (

*k*_{2},*k*_{−2},*k*_{4}and*k*_{−4})— Relative rate limitation of steps 2 and 4, i.e.

*k*_{4}relative to*k*_{2}.

Added to these relatively simple considerations is the possibility, even the probability, that some of these parameters may be voltage-dependent as the three-dimensional structure of the protein adapts to physiological differences in transmembrane potential (e.g. allowing glucose influx to increase in response to depolarization arising from metabolic or circulatory insufficiency).

In considering the thermodynamics and kinetics of transporters operating *in situ*, it is important to account for anisotropy in every aspect, whether in theoretical concepts, experimental protocols or physiological boundary conditions involving transmembrane electrical and chemical gradients. Very few of these aspects are comprehended within the traditional approaches of LNET or MMK; however, they are all comprehended within the approach of INET, based as it is on the law of mass action constrained by the probability isotherm.

From one point of view, the virtual experimental results reported here pertain to idealized abstractions, with only the broadest attempt to relate the absolute values of any of the kinetic constants or substrate concentrations to actual experimentally determined values; only the thermodynamic constraints on such kinetic constants have been strictly honoured. Nonetheless, a consistent pattern of results has emerged, indicating that unidirectional and net rates of a simple uniport reaction are exquisitely sensitive to the above-listed thermodynamic and kinetic factors. Put another way, the thermodynamic and kinetic parameters of an enzyme may be expected to be ‘tuned’ to its specific physiological function *in situ*. Moreover, slow changes in membrane potential may influence the manifestation of both anisotropy in binding affinity and anisotropy in the shuttling of bound and unbound forms of the enzyme, leading to voltage-dependent inhibition or activation of transport. Different situations will call for different thermodynamic and kinetic ‘tuning’, such as might occur among the many isoforms that are known to exist for enzymes that transport the same molecule in different anatomical locations with different physiological boundary conditions. In this light, the abundance of information on the numerous isoforms of glucose transporters [13,21,22] would seem to provide much fertile ground for application of the INET methods demonstrated here, constrained by the probability isotherm.

## 6. Conclusion

The purpose of this study was to explore what kinds of insights might be gained by applying classical kinetic concepts, such as the law of mass action, to explicit molecular models of membrane transport constrained by the probability isotherm. This had not been possible in the world that existed before the advent of high-speed digital computers. In those days, it was necessary to resort to simplifying assumptions of one kind or another to deal with kinetic models for which there were no ready-made analytical solutions available from mathematics. Therefore, it was that investigators came to rely on the Michaelis–Menten approach to enzyme kinetics on the one hand [22], and on the near-equilibrium linear force–flux assumptions of LNET, on the other [16].

Against this pre-digital backdrop, the present paper provides an extremely belated unveiling of the possibilities that may arise if the design and interpretation of the increasingly ingenious nano-experiments being performed around the world were to be informed by what may reasonably be called the emergent discipline of INET. INET may be distinguished in the following ways:

— INET is intuitive because it is based on familiar classical concepts in chemical kinetics, such as the law of mass action;

— INET is non-equilibrium in two senses, because:

it deals with unidirectional and net rates both at the level of the elementary step and at the level of the overall multi-step reaction, and

its application is valid under all conditions, whether they be very far from equilibrium, close to equilibrium or even precisely at thermodynamic equilibrium;

— INET is thermodynamic in that it recognizes the second law constraint on all forms of kinetic modelling that inheres in the probability isotherm and its kinetic equivalent, the rate isotherm. This effectively neutralises the unhelpful myth that thermodynamics has nothing to say about kinetics except at equilibrium.

The kinetic model studied in this paper is irreducibly minimalist as a model of facilitated diffusion, yet it has already indicated essential complexities of real enzyme behaviour that are not adequately accounted for by MMK in the pedagogical literature [22] or LNET principles in the research literature [16]. And these complexities are formidable, even without bringing enzyme mobilization, multimeric cooperativity or gene expression to account [21,22]. However, the list of adjustable parameters summarized above is already sufficiently large to make the model of very limited specific predictive value on its own. All that the present study can indicate is the breadth of the range of integrated kinetic and thermodynamic considerations that can, and should, be brought to the design and interpretation of experiments in biotransport enzymology. Further development of such models should ideally proceed closely in tandem with experimental work, either re-interpreted retrospectively [10,11] or designed prospectively, and be informed by the principles of INET, constrained everywhere and always by the probability isotherm.

It is also possible that, informed by this INET approach, we might begin to gain insight into the factors that determine the molecular behaviour of transport enzymes operating across the very strong transmembrane electric fields that are generally present. What actually *constrains* uniporters to do what they do (binding reactions, conformational flips, etc.) while obeying the second law exactly so as to ensure compliance with the constraint that the overall Δ*G*^{o} remains zero over the entire cycle of reaction steps? The second law is an (apparently) empirical fact, deriving from no known ‘laws’ more fundamental than itself, but its iron rule constrains even the most Heath Robinsonian collections of transmembrane amino acid sequences as they undergo their conformational changes and reveal their respective *cis*- and *trans*-binding affinities in accordance with their physiological functions. The second law informs us that this must be so and, though we accept this by inductive faith, we do not yet apprehend the mechanism by which this is determined.

## Data accessibility

All data were created using the Excel file (INET.xlsx) supplied as the electronic supplementary material. The numerical solutions were not stored but rather plotted graphically and screen captured for reproduction as text figures. The text figures are supplied as separate image files.

## Competing interests

I declare I have no competing interests.

## Funding

This work was carried out by the author *pro bono* with no salary or any other kind of funding. The author is supported by the School of Applied and Biomedical Science, Federation University Australia (Gippsland Campus), with the provision of an office and library privileges.

## Acknowledgements

This work has benefited from the kind assistance of Associate Professor Jennifer Mosse, Head of School of Applied and Biomedical Science at Federation University Australia, and Associate Professor Denis Loiselle and Dr Kenneth Tran of the University of Auckland. Associate Professors Mosse and Loiselle provided critical reviews of earlier drafts of this paper. Dr Tran performed independent cross-checking of the simulation methods using both the analytical solution to the steady-state linear algebraic equations and Matlab.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3792307.

↵1 It is also straightforward to demonstrate the continuous variability of net rate (flux) with composition at constant thermodynamic force by sweeping the

*cis-*and*trans*-substrate concentrations while holding their ratio constant, thereby demonstrating, once again, the misplaced expectation of LNET in a domain ruled by the law of mass action.

- Received April 28, 2017.
- Accepted May 16, 2017.

- © 2017 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.