A simple recipe for setting up the flux equations of cyclic and linear reaction schemes of ion transport with a high number of states: The arrow scheme

abstract The calculation of flux equations or current-voltage relationships in reaction kinetic models with a high number of states can be very cumbersome. Here, a recipe based on an arrow scheme is presented, which yields a straightforward access to the minimum form of the flux equations and the occupation probability of the involved states in cyclic and linear reaction schemes. This is extremely simple for cyclic schemes without branches. If branches are involved, the effort of setting up the equations is a little bit higher. However, also here a straightforward recipe making use of so-called reserve factors is provided for implementing the branches into the cyclic scheme, thus enabling also a simple treatment of such cases.


Introduction
Analyzing measured current-voltage relationships on the basis of a reaction kinetic model is a powerful, long-serving approach. [1][2][3][4][5][6] Once the equations are set up, the computational costs are very low. The benefits are the delivery of model parameters, which are related in a straightforward way to the biophysical or structural parameters of the transport protein and the kinetics of the transport process. Reaction kinetic analysis of ion transport on the basis of such a model yields the rate constants of loading of the transport site, the velocity of the electro-sensitive translocation step, the influence of substrate concentrations or of ATP on these parameters, 7-13 the occupation probability of intermediate (structural) states of the protein, 11 charge of the translocated complex, 14 binding order in the case of cotransporters and pumps. 6,15,16 Hills et al. 5 even developed a computer program to model an ensemble of transport molecules in order to describe the transport physiology of an entire guard cell.
Furthermore, the reaction kinetic analysis can also reveal molecular and structural information. For instance, Hagedorn et al. 17 employed kinetic modeling to bacteriorhodopsin and could describe the fast photoisomerization of retinal with simultaneous proton transfer to an acceptor, reprotonation of retinal from the intracellular face via a proton donor, and proton release to the extracellular space via a release complex. In the case of the calmodulin-binding domain of plasma membrane Ca 2C -ATPase isoform 4b, Penniston et al. 18 assigned the states in the reaction scheme to structural changes induced by association and dissociation of CaM and C28. Below in Figure 1, the analysis of proton transport in the viral M2 channel is given as an example in order to illustrate the reaction kinetic approach and the conclusions, which can be drawn from fitting measured IV curves. 11 Of special interest are cyclic reaction schemes. It is immediately obvious that they apply to electrogenic pumps 3,8,12,13,16,19,22 and cotransporters. 6,14,15,23,24 However, it turned out that cyclic reaction schemes also provided the best fits of IV curves from ion channels, 7,25-27 from bacteriorhodopsin, 9,17 the viral M2 channel 11 or from Polytheonamide B, a marine cytotoxic b-helical peptide. 28 Nowadays, MD simulations of ion transport proteins become more and more powerful. [29][30][31] In order to test their predictions, experimental results have to be presented in quantitative terms. The kinetic analysis provides a bunch of such parameters as mentioned above.
However, setting up the rate equations for cyclic systems with many states can become extremely cumbersome. Fortunately, there is a simple procedure for generating the equations using the so-called arrow scheme for linear and cyclic reaction models. It immediately delivers a minimum form of the desired equations without any computational effort.
We have used parts of this arrow scheme in some previous publications (e.g. ref. 15). However, we never have given a complete description, especially no proof of its validity. This is done here, hoping that it may save to lot of time and effort for researchers who want to evaluate ion transport on a quantitative basis.
Those, who are only interested in setting up the equations for evaluating their data, do not need to read the whole text. They just can use Eq. 22 for calculating gross rate constants in linear reaction schemes. In simple cyclic schemes, the occupation probabilities are obtained from Eqs. 29 or 30, adapted to the adequate number of states. The transport rate or the current are given by Eq. 33c. For branched cyclic schemes, a step by step recipe is given at the end of the article. Others, who are interested in the mathematical background, and especially in the proofs of this approach, should read the whole text.

An example for the kinetic analysis of membrane transport
In order to illustrate the rationale of the kinetic analysis of ion transport, we give a short outline of the application of the cyclic 3-state model to the evaluation of measured IV curves of proton flux through the M2 channel of the influenza A virus (Fig. 1A). 11 The selection of the number of states has to be done along the following rules: Experiments reveal only those states, which are affected by the experimental conditions. In the experiments on the viral M2 channel, internal pH H C influences the binding of protons to His37, and membrane potential V m controls the exchange with the external medium (Fig. 1B). These two reactions require three states. Figure 1C shows these three states: C D His37 ready to bind an internal proton, B D His37 having bound a proton (actually 3, but only one is exchanged) and O the external site having accepted a proton.
The model in Figure 1C leads to the following rate equations for the occupation probabilities P i of the states with i D C, B and O. . This is legitimate if the individual reactions cannot be distinguished by the experimental procedure. 4 The binding and release of protons from and to the external side is thus included into k CO and k OC , respectively.
Macroscopic IV curves were measured under steady-state conditions, as indicated by the zeros in Eqs. 1a to c. Solving the above set of equations leads to the occupation probabilities of the three states In those experiments, internal proton concentration, H C , and membrane voltage, V m , were changed leading to the following dependence of the related rate constants on those experimental conditions with s being the location of the Eyring barrier and u T the characteristic voltage.
The single-channel current is calculated from the difference between the unidirectional currents between states B and O With the charge ze D 1e for protons, insertion of P B and P O leads the following macroscopic current with N being the number of M2 proteins and P Open the open-probability in a macro-patch. In those experiments or in whole-cell experiments, the product NP Open cannot be determined and influences as a common scaling factor the values of the rate constants. Thus, the absolute values have to be revealed by single-channel experiments. In the case of M2, the number N and thus the single-channel conductivity can for example be determined from the amount of protein added to a suspension of liposomes. 32 Nevertheless, the occupation probabilities and the ratios between the rate constants are independent on this scaling factor, and thus provide valuable information also without the knowledge of N.
In the case of the M2 channel, IV curves measured at different internal pH were subject to a global fit., i.e., u T , s, and the rate constants k BC , k CB,1 , k OB,0 , k BO,0 , k OC , k CO , were free fit parameters, but had to be the same for all curves, and H C and V m were fixed parameters. The key findings obtained from the fit are as follows: very strong binding of the cytosolic protons to His37, only a minority of them are finally transferred to the outside. Most of them return to the cytosol via k BC . The binding rate constant k CB is proportional to the internal proton concentration. The voltage-sensitive reactions k BO and k OB are very slow and asymmetric. A relevant finding was that good fits of the IV curves could not be obtained without the assumption of k OC and k CO . These reactions imply a cyclic reaction scheme and are related to the resetting of the protonation-induced tilting of the helices. 33,34 For further details see. 11 The above example serves to discuss the crucial questions of such kind of analysis, namely, what is the gain of the analysis, and how reliable are the results? The first question is already addressed above by listing the insights, which are obtained from the numerical values of the rate constants determined by the fit.
The second question has to deal with two aspects: the legitimacy and the implications of using models with a minimum number of states and the power of the analysis to reveal many details of the transport process. These two questions are strongly related, thus we discuss them together. An individual single-channel IV curve can be described by a cyclic 2-state model. 3,4 The states O and B are connected by voltage-sensitive reactions k OB D k BO,0 exp((1-s)V m /u T ) and k BO D k BO,0 exp(-sV m /u T ) and the neutral recycling reactions k OB and k BO . s is the location of an Eyring barrier between O and B, and u T the characteristic voltage (mostly u T D 25 mV for ions with one charge). The 2-state equation of the singlechannel current I is If the range of voltage V m is wide enough (e.g., C/¡ 200 mV, depending on the transport protein) all unknown parameters of Eq. 5 can be obtained from a fit of the IV curve. Often, two parameters can already be determined from a visual inspection. 3,35 At very high positive or negative potentials, the current saturates, and takes the values I sat,neg D -zek OB or I sat,pos D zek BO . Furthermore, the shape of the IV curve provides enough features to determine also the other unknown parameters. s and u T determine the dependence of current on membrane potential, which may be different for inward and outward currents if s is different from 0.5. In the neighborhood of the reversal potential, the IV curve increases linearly with membrane potential if k BO,0 and k OB,0 are much greater than k OB and k BO , and exponentially when they are much smaller. Thus, the bending of the IV curve provides the information about the values of k BO,0 and k OB,0 . Furthermore, the numerator of Eq. 5 has to be zero at the reversal potential, thus providing another constraint on the values of the rate constants. If the regions of saturation are not reached, the determination of the parameters may be still correct as the bending of the IV curve toward the saturation region can still provide enough information. Whether this information is enough can be tested by repeating the fit, with that parameter, which is doubted to be unreliable, set to a different fixed value. As mentioned above, the absolute values of the rate constants are only obtained from single-channel IV curves.
Data from macro-patches or whole-cell experiments include a normally unknown scaling factor N P OPEN with N D number of transport proteins and P OPEN the average open probability. However, also in this case valuable information can be obtained from the relative magnitude of the rate constants showing, which reactions are fast or slow or whether they are asymmetric.
For fitting IV curves with a 3-state model (as in Fig. 1C), IV curves have to be measured at different concentrations of the transported ion (e.g. internal pH in the experiments of Fig. 1). Then, a set of IV curves is available, providing additional information. There are two different approaches. The analysis can be done using a 2-state model (O-B, not shown) or using the 3-state model shown in Figure 1C. In the case of the 2-state model, the evaluated rate constants depend on the concentrations, because now k OB and k BO are gross rate constants of the O-C-B branch. The relationship between the individual reactions k OC , k CB , k BC , k CO and the gross rate constants can be calculated by means of the arrow scheme described below (Eq. 12). This calculation would result in a unique determination of the individual rate constants because the relationships between the ion concentration, [ion] i , in the binding reaction k CB D k CB,1 [ion] I and k OB and k BO are non-linear (see Eq. 12, below).
However, the conversion of the 2-state results to the 3-state parameters is a little bit more complicated, because the law of mass conservation applied to the occupation probabilities P C C P B C P O D 1 (Eq. 1d, the protein can be only in one state) has to be taken into account. P C changes with the ion concentration and so do P O and P B . This effect can be handled by the usage of reserve factors, 4,35 which account for hidden states. C would be a hidden state in the 2-state analysis, but not in the 3-state analysis.
However, the states WcVi3, WcVo2Pn and WcVc2 in Figure 1B, which are suggested by structural analysis as discussed by ref. 11, remain hidden also in the 3state analysis. An access to one of these hidden reactions could be obtained by measuring at different external pH as already indicated in Figure 1B at the rate constants between O D WcVo2Pm and WcVo2Pn (Tryptophan gate closed, Valine gate open with m or n protons in the cavity near Val27). Then, the same considerations hold as described for the effect of internal pH. The 3-state model including the effect of internal pH would be extended to a 4-state model including the effects of internal and external pH. Thus, it always has to be kept in mind that also the 3-state parameters include unknown reserve factors. However, this holds for any kind of kinetic analysis as a protein has a plethora of different states. The theory of reserve factors implies that an experimentally determined state mostly presents a pool of states, which are connected by interactions not influenced by the actual experimental procedure. Details of the philosophy of reserve factors are provided in the last part of this article.
The message of the above considerations is: 1. The introduction of a new experimental approach provides enough information to reveal more details of the transport process. In the example of changing internal pH, it is the non-linear relationship between k CB D k CB,1 [ion] i and the gross rate constants k OB and k BO . 2. The calculation of the 3-state parameters from a 2-state analysis is mathematically complex because of the reserve factors, but could be solved with the tools provided below. 3. The fit results are more stable and more directly connected to the involved mechanism when a model is used, which includes all transitions between states which are influenced by the experimental procedure.
Setting up the equation for the current (Eq. 4) by solving the rate equations (Eq. 1) manually is still simple for a 3-state model like that in Figure 1C. However, doing so we did not realize at that time that k OB Ck BC in Eq. 4b could be factored out also in the denominator, which would have reduced the number of terms by a factor of 2. Setting up the equations becomes still more and more cumbersome, when the number of states is increased. For instance, when we tried to apply the cyclic 5-state model of single-file transport through KcsA 31 for the analysis of IV curves from a viral potassium channel (unpublished data), the step-wise elimination of the intermediate states from equations similar to Eq. 1 led to terms with 125 products of 10 rate constants in the numerator and 600 in the denominator.
Fortunately, there exists the arrow scheme, which enables an immediate setting-up of the equations for the occupation probability P j of the involved states S j and for the fluxes between the states. This arrow scheme is described below, and the proof of its validity is given. It provides the correct solutions with the minimum number of terms.

Definitions
The following symbols are used in the text below: S i : state of the model P i : occupation probability of state S i k i,j : rate constant of the transitions between neighboring states S i and S j , represented by ordinary arrows in the arrow schemes, below g i,j and h i,j : gross rate constants between non-adjacent states (see Eqs. 17a, b below), represented by slightly longer arrows s i,j : gross rate constants of single loops in branched systems Several arrow schemes are used as a shorthand for sums of products of rate constants. In linear reaction models (Eq. 16, below) ranging from state S m to state S n , the arrow scheme D Ln m is used. The superscript Ln means "linear scheme with the last considered state being n." n is always at the right-hand side. The subscript m is the index of the state S m at the left-hand side where the chain starts. For the most frequent case m D 1, the matrix is (6) Each arrow represents a rate constant k i,j between the state S i and its neighbor S j with j D i C 1 or j D i¡1. The related states are given on top of the scheme over the gaps between the arrows with the direction of the flux defined by the arrows. A typical feature is the empty diagonal with arrows on either side pointing away from the diagonal. D Ln m is quadratic and its dimension is n-m (i.e., n-1 in the case of D Ln 1 ). In a cyclic reaction scheme (see Eq. 23, below), the arrow scheme D Cn m includes one step more (from S n to S 1 closing the loop, and thus has the dimension n and not n-1). The superscript Cn means "cyclic scheme with n states" (in contrast to D Ln m in Eq. 6, n is not the index of the last state). The subscript denotes the index of the state S m where the chain starts and ends. For example, for m D 1 this is (7) In the definition of D Cn m (Eqs. 7,8), all indices i greater than n are equivalent to i-n, due to the cyclic nature of the model, as illustrated by the equations for The gross rate constants of a linear reaction scheme with n states Necessary for the treatment of models with an arbitrary number of states is the calculation of gross rate constants g i,j , lumping the reactions between non-adjacent states S i and S j . The proof of the final equation for n states is done by the procedure of induction. Thus, we start with a simple chain of 3 states k i,j are individual rate constants between two adjacent states, g i,j are gross rate constants including intermediate states. The rate equations for S 1 and S 2 in steadystate are with P j being the probability of being in state S j . Considering steady-state fluxes (which determine measured IV curves), it follows: The comparison of Eqs. 10a and 11c yields the gross rate constants For the extension of these equations to chains with n states, we present the gross rate constants by arrow schemes (see Eq. 6). An arrow between the numbers j and jC1 in the headline stands for k j,jC1 if it points to the right, and k jC1,j if it points to the left. Arrows in the same row mean that the related rate constants are multiplied. A vertical stack of rows of arrows represents the sum of these products.
With a set of equations similar to Eqs. 13a, b (or by the calculations below) it can be shown that the gross rate constants between 4 states are The denominator for a chain between state S m and state S n is called D Ln m (Eq. 6). With the rule that the arrows in one row are multiplied, and the products are added, the scheme in Eq. 14a leads to the following equation for the 3-state model Now, the above procedure is extended to linear chains with nC1 states by induction, i.e., we assume that the equations hold for n states. Based on this assumption, we show that it holds for nC1 states. Having shown the validity for 2 states, it follows that it hold for 3 states. If it holds for 3 states, it holds for 4 states, etc.
In order to show that the arrow scheme applies also to a linear reaction scheme with nC1 states, we start from a chain with n states, and assume that the following equations are correct for this linear reaction scheme with n states (see also Eqs. 13 and 14) Next, it is shown that a similar equation holds for a linear scheme with nC1 states. For the extension to a linear reaction scheme with nC1 states, we replace the two reaction steps between S n-1 and S nC1 by the gross rate constants g n-1,nC1 and g nC1,n-1 .
S n ¡ 1 ÀÀÀÀÀ ! ÀÀÀÀ k n ¡ 1;n k n;n ¡ 1 S n ÀÀÀÀÀ ! ÀÀÀÀ k n;n C 1 k n C 1;n S n C 1 g n C 1;n ¡ 1 ÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀ and thus we still can use the equation 17a for n states. However, there is a problem with the nomenclature. The chain in Eq. 18 with the lumped reaction has n states, but S n is hidden and the last one is S nC1 . We make the ad hoc decision to call the respective matrix Dg Ln 1 analog to D Ln 1 (Eq. 6), because state S nC1 takes over the role of S n in the lumped chain.
j D 1 k j;j C 1 k n ¡ 1;n k n;n C 1 k n;n ¡ 1 C k n;n C 1 Dg Ln 1 (19a) with Dg Ln 1 being the denominator for n states including gross rate constants (long arrows) between S n-1 and S nC1 : Introducing Eq. 13a for the lumped reaction g n-1,nC1 (n corresponds to 2 in Eq. 13a) leads to (20) The sum in the last row of Dg Ln 1 leads to an additional row with respect to D Ln 1 , thus proving that Dg Ln 1 has n rows as required for D Ln C 1 1 : (21a) The first factor is the standard quadratic symmetrical form D Ln Inserting Dg Ln 1 (Eq. 21c) into Eq. 19b leads to the desired form of the gross rate constant for nC1 states, which is of the same form as that for n states. Thus, the induction verifies the validity of Eq. 17 (for n) or Eq. 22 (for nC1): In the case of g nC1,1 , the arrows in the numerator point into the opposite direction.
Probability P m of being in state S m in a cyclic model with n states The trick of the following calculations is the expression of P m in terms of P 1 . For this issue, we use the gross rate constants g i,j for the first part of the cycle and h i,j for the second part.
The symbols S 1 at the right-hand side and at the left-hand side denote the same state S 1 . This pseudolinear presentation renders the generation of the equations simpler than the cyclic representation.
Again, fluxes under steady-state conditions (as apply to measured IV curves) are considered. To calculate the occupation probability P m of an intermediate state, we read Eq. 23 as a pseudo 3-state model. Thus, according to Eq. 11a, P m is determined by the sum of fluxes into S m from source states via individual and via gross rate constants depending on how far the two source states are away from P m . In a cyclic scheme, the two source states are identical, e.g., S 1 . Thus, the steady-state probability P m of being in state S m is The relationships between the gross rate constants g i,j or h i,j and the individual rate constants k i,j are given above by Eqs. 17 or 22.
(Remember the definition: index nC1 D index 1 in cyclic schemes). DEN is the product of the two denominators in Eq. 25a.
One transport protein is considered. Thus, the sum over all P m (Eq. 26) is 1.
Resulting in Eq. 28 is a quite simple equation for the occupation probability of state S 1 . It is obvious that this form also holds for any other state S m of the cyclic reaction scheme, because each state can be called S 1 by rotating the indices in the circle. The The arrow scheme in the denominator is the same for all P m . Only that one of the numerator has to be adapted to the actual state P m , i.e., D Cn m .

The current in a cyclic reaction scheme
We consider the net current starting from state S 1 . This is not a restriction of general applicability, because in a cyclic scheme the indices can be rotated in such a way that the state of interest can always be called S 1 . The current can be calculated from the difference of the overall forward and backward currents of the whole circular reaction scheme or from the difference of the fluxes between two adjacent states. The result is the same as shown below. We start with the circular currents for one transport protein. If macroscopic currents are considered, the expressions have to be multiplied by the number of active transport proteins. From the circular forward and backward gross rate constants g 1 ! 1 and g 1 1 , respectively, (Eqs. 17 and 22) (spanning over the whole loop from S 1 to S 1 ) and the probability P 1 (Eq. 28) the current is calculated with z being the charge of the electro-sensitive translocation step, e the elementary charge, and index nC1 D index 1, because of the cyclic scheme (Eq. 23).
If we calculate the net current from the flow between two adjacent states (For the sake of clarity we use n D 5 however, the algorithms can easily be transferred to any other number of states.): The difference of fluxes between two states gives the same value of current as that calculated from the circular gross rate constants starting from one state. Including k 1,2 and k 2,1 in the (now asymmetric) matrices leads to (33a) In the numerator of Eq. 33a, only the first row in the left-hand matrix and the last row in the right-hand matrix are not canceled out.
In a product, the sequence of factors can be rearranged, thus the arrow of k 1,2 is moved from the end to the beginning of the first term in the numerator.

5
(33c) The denominator in Eq. 33 is represented by the same arrow scheme as the denominator of Eq. 30. The comparison between Eq. 31 and Eq. 33c shows that the net fluxes are the same regardless of whether they are calculated from the whole circular gross reactions or from the fluxes between two states. However, the unidirectional fluxes are quite different, because the four rows in the matrices of Eq. 33a cancel each other only if net fluxes are considered.

Unidirectional fluxes and tracer fluxes
The knowledge of the rate constants and occupation probabilities as obtained from fitting an appropriate set of IV curves by means of Eq. 33c also enables the determination of the unidirectional currents. Overall outward current is given by the first (positive) term in Eq. 31 and inward current by the second (negative) term in Eq. 31. Caution: They cannot be calculated from the first and second term in Eq. 32. Even though the net current is the same in Eq. 31 and Eq. 32, the fluxes related to the exchange between adjacent or (more distant) states may be higher than unidirectional current obtained from the whole cycle. For instance in the example of Figure 1C, the unidirectional fluxes between C and B are extremely high. 11 This inward current is carried by k BC . However, since most of the protons return to B, only those few protons continuing to O contribute to the measured outward current.
In a general tracer experiment, there is (for instance in a cotransporter, see Figure 2, below) a cycle where the tracer is bound at state S a and released from state S b on the other side of the membrane. In this case, Eq. 32 becomes important. Using Eq. 17a for g a,b and Eq. 29 for P a , the unidirectional flux from S a to S b is Thus, the reaction kinetic evaluation of IV curves also provides a means of comparing electrophysiological experiments with tracer measurements. However, the number of transport molecules again remains an unknown scaling factor of the rate constants.

Branched reaction schemes
Some transport proteins require a model, which is a little bit more complicated than the cyclic scheme, as a third path of transmembrane transport has to be incorporated. One example is the Na C /glucose cotransporter 36 with the following transmembrane movements: In the model of Figure 2A, the reaction S 1 -S 6 is the translocation of the unloaded "carrier." In cotransporters, the "movement of the carrier" means conformational change of the protein from "open to the outside" to "open to the inside." 24 S 2 -S 5 is the translocation of the "carrier" with 2 Na C bound, S 3 -S 4 is the real cotransport with 2 Na C and one glucose bound. Another example is the association and dissociation of the Calmodulin-binding domain of the plasma membrane Ca 2C -ATPase. 18 Branched reactions schemes mostly apply to cotransport systems. Single-channel experiments are not available for cotransport systems because of the low transfer rates, and a reliable determination of the number of channels in macroscopic experiments (as e.g. done in ref. 32) might not always be possible. Nevertheless, the kinetic analysis provides valuable information about the (a)symmetry of forward and backward reactions, the relative affinities of the binding reactions (S 1 -S 2 , S 2 -S 3 , S 6 -S 5 , S 5 -S 4 in Fig. 2A), the rate limiting steps, and the ratio of fluxes through different branches like S 2 -S 5 and S 3 -S 4 in Figure 2A.
A prerequisite for extracting all this information again is the usage of many experimental approaches, e.g., changing membrane potential to obtain IV curves and studying the influence of the cytosolic and external concentrations (activities) of all transported species for the global fit of many IV curves. As discussed above for the M2 system in Figure 1, this provides enough information for determining all (relative) rate constants of a 6-state model in Figure 2A and the occupation probabilities, which are not affected by the incomplete states, as, for instance applies to a 2Na C /glucose cotransporter described by. 36 Binding reactions glucose (G) and Na C (2Na) are marked. (B) Illustration of the 2 possible paths between the states S 2 and S 5 , which are presented by the gross rate constants g 5,2 and g 2,5 . The longer single path via S 3 and S 4 is described by the gross rate constants s 5,2 and s 2,5 . (C) Reduced 4-state model, in which the reactions between states S 2 and S 5 are lumped as indicated by the box in (A). Double arrows indicate pseudo-4-state rate constants, which are different from the corresponding ones of the 6-state model in (A). Apparent states S Ã i and rate constants k i,j Ã different from those of the full 6-state model as changed by lumping the transitions between state S 2 and S 5 are marked by stars, g i,j Ã are gross rate constants comprising both paths in (B).
knowledge of the number of transporters N. The occupation of the states may become valuable when occupation probabilities will be predicted by MD simulations, thermodynamic considerations or new techniques (e.g., voltage-clamp fluorometry as used in gating studies 37,38 ) are able to show different structures for different binding states. Forrest et al. 39 complain that despite the wealth of structural information, not enough information about dynamics and dwell-times of the intermediate states is available. Dwell-times correspond to occupation probabilities. They can be provided by a kinetic analysis of IV curves or tracer experiments. However, the derivation of the flux equations of branched systems is even more cumbersome than in the case of the schemes with an unbranched cycle. Here, we make use of the concept of reserve factors and gross rate constants in order to derive the equations of the full system in Figure 2A. First, the reserve factors help us to derive the equations of a branched system ( Fig. 2A) from the results obtained above for unbranched systems. For this purpose, the calculation starts with merging the 6-state model in Figure 2A into the pseudo 4-state model of Figure 2C. Second, the calculations below also give an introduction into the philosophy of reserve factors.
From the equations of the pseudo 4-state model (Fig. 2C), the equations of the full 6-state model can be calculated. Because of the great variety of such models, it is not possible to give a general scheme. However, from the calculations below some rules should become obvious, which can considerably ease the creation of the equations for the current or the occupation probability of the involved states. We will first derive the equations for the example in Figure 2, and then give a recipe for the generation of the equations for a general scheme at the end.
The reduced reaction scheme "Lumped" reactions are used to describe the transport in the "full" model in Figure 2A in terms of a simpler "pseudo-4-state" cyclic model (Fig. 2C). For this purpose, the reactions between state S 2 and state S 5 are merged in gross rate constants g 2,5 and g 5,2 as indicated by the box in Figure 2A. However, in contrast to the matrix of the linear scheme in Eq. 19c, we have to account for the states, which are hidden in the reduced model (S 3 and S 4 in Figure 2B). In order to remember that we deal with a pseudo 4-state model, the symbols D Ã g, double arrows g Ã i,j (longer double arrows in Eq. 35) and states S i Ã with an asterisk are used for the arrow scheme of the reduced model (see "Extended matrix definitions," below). Nevertheless, also for this reduced model, the arrow scheme of a 4-state model applies.
Next, the equations, which give the relationships between the k i,j (ordinary arrows) of the 6-state model and the k i,j Ã and g i,j Ã (double arrows, starting from states S 2 Ã or S 5 Ã ) of the pseudo four-state model are derived. The double arrows k 5,6 Ã and k 2,1 Ã are calculated in Eq. 42a, g 2,5 Ã and g 5,2 Ã in Eq. 42b, and the occupation probabilities of the states S 2 Ã and S 5 Ã in Eqs. 43 and 44a.

Incorporation of gross rate constants
First, the gross rate constants g 2,5 and g 5,2 between S 2 and S 5 (of the 6 state model) are calculated. They are the sum of the rate constants of the two paths between S 2 and S 5 (Fig. 2A). These paths are (1) the direct path via k 2,5 and k 5,2 and (2) the indirect path via S 3 and S 4 (Fig. 2B). For the latter reactions from S 2 to S 5 via S 3 and S 4 , the two gross rate constants s 2,5 and s 5,2 have to be calculated by means of Eq. 22. (36a) However, these gross rate constants are not yet sufficient to describe the transitions between S 2 Ã and S 5 Ã of the reduced 4-state model in Figure 2C. In order to obtain g 2,5 Ã and g 5,2 Ã , Eqs. 36a, b still have to be divided by so-called reserve factors. Reserve factors account for the fact that the occupation probabilities of S 2 Ã and S 5 Ã are higher than those of S 2 and S 5 , respectively, because they also have to represent the hidden states S 3 and S 4 .

Mass conservation and reserve factors
In the arrow scheme of Eq. 35, asterisks are attached to the states S 2 Ã and S 5 Ã , because S 2 and S 5 have a different meaning depending on whether we use them in the reduced model (Fig. 2C) or in the model of Figure 2A, which includes the knowledge of the 2-3-4-5 branch. This follows from the equation of mass conservation. In the case of a transport protein, the sum of all occupation probabilities has to be 1. The protein can be in only one of the states of the reaction scheme. The sum of occupation probabilities is in the reduced model : with P Ã 1 D P 1 and P Ã 6 D P 6 , because S 1 and S 6 are not connected to hidden states. in the full model : 1 D P 1 C P 2 C P 3 C P 4 C P 5 C P 6 (37b) The comparison of Eq. 37a, b shows that P 2 Ã C P 5 Ã D P 2 C P 3 C P 4 C P 5 D r 2 P 2 C r 5 P 5 (38) In Eq. 38, the states S 3 and S 4 in Figure 2A are merged with the states S 2 and S 5 via the so-called reserve factors r 2 and r 5 , respectively (ref. 4, do not use the inverse definition of ref. 3). The values of the reserve factors are calculated from equations similar to Eq. 11 (which gives the occupation probability of an intermediate state), but now with gross rate constants s i,j (of single branches) between the representative states S 2 or S 5 and the intermediate states S 3 and S 4 . P 3 D k 2;3 P 2 C s 5;3 P 5 k 3;2 C s 3;5 D r 3;2 P 2 C r 3;5 P 5 (39a) P 4 D s 2;4 P 2 C k 5;4 P 5 s 4;2 C k 4;5 D r 4;2 P 2 C r 4;5 P 5 (39b) The comparison between the second and third term of Eq. 39a, b yields the reserve factors e.g., r 3,2 r 3;2 D k 2;3 k 3;2 C s 3;5 (39c) We have to distinguish between reserve factors with one index and those with two indices. Those with two indices describe that part of state S 3 or S 4 which is attributed to state S 2 or S 5 (Eq. 39).
The full reserve factor r j of a state S j (j D 2 or 5), which is connected to a pool of merged states, (S 3 and S 4 in Fig. 2A and B) is obtained from the following identities according to Eqs. 39: leading to Correction of the rate constants of the reduced model (Fig. 2C), which point away from a representative state as indicated by an asterisk (e.g. S Ã 2 and S Ã 5 ) The "real" unidirectional flux between two states in the full model in Fig. 2A is k i,j P i with j D iC1 or j D i-1. However, in the reduced model (Fig. 2C), P Ã i is calculated, and P i is not explicitly shown. Nevertheless, both models have to yield the same unidirectional fluxes for each pair of states, because those are the parameters, which can potentially be obtained from experiments, e.g., from fitting IV curves or from tracer experiments.
In order to yield the same unidirectional flux as in the full model a different value of k i,j , namely k Ã i;j (as indicated by the double arrows in Eq. 35) has to be used. The value of k Ã i;j is obtained from the condition that the fluxes in both models have to be the same, i.e.
Equation 41 implies that the "real" rate constants k i,j of the extended model are obtained from those of the reduced model by multiplying them with the reserve factor of the states from where the double arrow starts. The same argument can be made for the gross rate constant g i,j (Eq. 36) In Figure 2C, the k Ã i;j are those reactions, which start from i D 2 or 5, as indicated by the double arrows.
Since the reserve factors are greater than 1 (Eq. 40), P Ã i is greater than P i . Consequently, k Ã i;j has to be smaller than k i;j in order to result in the same flux.

Extended matrix definitions for branched models
For the sake of clarity, here the definitions of the matrices of the branched model are summarized. For the model in Figure 2A is n D 6 and for the reduced one in Figure 2C is n Ã D 4. D Cn j : arrow scheme starting with S j of the 6-state model (n D 6) in Figure 2A with all rate constants k i,j explicitly given, equal to Eq. 8. Dg CnÃ j : arrow scheme starting with S j of the 6-state model in Figure 2A, but with the rate constants between S 2 and S 5 lumped into the gross rate constants g i,j according to Eq. 36a, b (n Ã D 4). This matrix must not be used for a 4-state model in equations like Eqs. 31 or 32, because then the law of mass conservation would be violated. D Ã g CnÃ j : arrow scheme starting with S j of the reduced 4-state model (n Ã D 4) in Figure 2C with the rate constants between S 2 and S 5 lumped into the gross rate constants g i,j Ã (Eq. 42b). All rate constants starting from a representative state (S 2 Ã and S 5 Ã in Fig. 2C) are presented by a double arrow in Eq. 35 and include the reserve factor r i of that state according to Eq. 42. This matrix can be used for a 4-state model in equations like Eqs. 31 or 32.

Calculation of the occupation probabilities and of the current
Occupation probabilities, net current and unidirectional fluxes can be calculated in terms of the pseudo 4-state model or in terms of the full 6-state model. The results of IV curve fits are more stable when the calculations are done on the basis of the full model, but for the derivation of equations for the full model, we start from the pseudo-4 state model and then convert it to the full model. First, we show that the occupation of states S 1 , S 2 , S 5 , and S 6 of the 6-state model ( Fig. 2A) can be expressed by the matrices of the 4-state model (Fig. 2C): Adapting Eq. 29 to gross rate constants, we obtain (Remember, the sum of the denominator is 1 because of the mass conservation.). In the denominator of the third term, the number of states is that of the reduced model in Figure 2C: n Ã D 4. For the sake of a comprehensive description, we use r j D 1 for states not related to a pool of lumped reactions, e.g. r 1 D r 6 D 1 in the example of the 6-state model in Figure 2A. In the third term of Eq. 43, the indices of the independent states are called "ind" (S 1 and S 6 in Fig. 2) and the those of the representative states "rep" (S Ã 2 and S Ã 5 in Fig. 2C).
Calculating the occupation probabilities in terms of the reduced model in Figure 2CF ormelabschnitt (n€ achster) leads to a modified version of Eq. 29: Next, we show that Eqs. 43 and 44a are equivalent for the 6-state model in Figure 2A. To this end, we first need to express D Ã g C4 1 in terms of Dg C4 1 : (44b) In each row, double arrows start from state S 2 Ã and S 5 Ã , thus each row contains the factors 1/r 2 and 1/r 5 .
(Remember: ordinary arrows contain no reserve factors).
In contrast, the rows in the next matrix (starting from state S 2 Ã ) have only double arrows starting from S 5 Ã : (44c) Note: The long ordinary arrows between S 2 Ã and S 5 Ã in the right-hand matrix denote lumped reactions (g 2,5 and g 5,2 ) without reserve factors. In a way similar to Eqs. 44b, c we obtain Introducing Eq. 44b to e into Eq. 44a leads to The factors r m and r j serve to introduce the correct factors according to Eqs. 44b to e.
The last equality follows from Eq. 43. For the sake of a comprehensive description, we use r m D 1 and P m Ã D P m for states S 1 and S 6 , which are not connected to lumped states. Equation 45 shows that it does not matter whether we use Eq. 43 (model in Fig. 2A) or Eq. 44a (model in Fig. 2C) for the calculation of the occupation probabilities. Also for the current, the same results are obtained regardless of whether the matrices r j Dg C4 j or D Ã g C4 j are used. In terms of the 6-state model, the current is after adapting Eq. 33c to a pseudo 4-state model Note: The long single arrows between S 2 Ã and S 5 Ã denote lumped reactions (g 2,5 and g 5,2 ) without reserve factors.
In terms of the reduced 4-state model (Fig. 2C), the current is  Figure 2C is the same as Eq. 46 (or Eq. 47c) derived in terms of the model in Figure 2A.
Fitting IV curves using Eq. 46 is actually a fit in terms of the 6-state model, because the rate constants k i,j of the branch between S 2 and S 5 via S 3 and S 4 have to be free fit parameters. They are used to calculate the reserve factors r i and the lumped reactions g i,j Ã , which then are introduced into Eq. 46. In contrast, employing Eq. 47a for a 4-state fit yields the gross rate constants g i,j Ã including reserve factors as fit parameters. Then, parameters of the 6-state model have to be extracted by means of Eqs. 22, 39c, 40c, 42a,b). This is not recommended, because of the mathematical effort of these conversions, the influence of the individual 6 states is less obvious (e.g., when investigating the significance of the fitted values by introducing deviating fixed values of these parameters), and the accuracy of the calculated values of the 6-state parameters suffers.
The recipe for setting up the equations of a branched system with n states Cotransport systems (or others) may have more states than the 6 states of the cotransporter in Figure 2A. An example is the 9-state model of the Na,K-ATPase. 16 Thus, we consider schemes with additional states in the translocations steps between states S 2 and S 5 and states S 3 and S 4 . It may be difficult to find an experimental access, which influences the reactions related to these intermediate states. However, that is a problem of the particular experiment. In the case of the Na,K-ATPase it was the ATP concentration, which provided an access to the intermediate states.
Here, we provide the tool for a simple generation of the adequate equations. This tool may also be helpful if the existence of such an intermediate state is suspected from crystallographic, spectroscopic, computational or other approaches.
Step 1, gross rate constants The representative states are S a and S b (S 2 and S 5 in Fig. 2). The gross rate constants s i,j have to be calculated by means of Eq. 22 (similar to Eq. 36) (48a) (48b) Eqs. 48a, b generate many equations: One set with i D a (the representative state) and j D intermediate state (can be more than one, e.g., S 3 and S 4 in Fig. 2), one set with i D b (the representative state at the other end of the chain, and j D intermediate state, and two equations with i D a and j D b (Eq. 36).
Step 2: The gross rate constants s i,j of the individual branches between S a and S b have to be added (compare Eq. 36). In general, both (all) branches between S a and S b may require gross rate constants The superscript distinguishes between the different branches. In the model of Figure 2, Figure 2. Using gross rate constants also for the second branch (and additional ones) may become necessary, for instance, for a 9-state scheme of the Na C /K C pump. 16 Step 3: Knowing the gross rate constants g a,b and g b,a the matrices Dg CnÃ j can be generated, as shown for (50) The long arrows between S a and S b are the gross rate constants g a,b and g b,a without reserve factors.
Step 4: The gross rate constants of Eq. 49 can be used to calculate the reserve factors according to Eqs. 39 The indices j and i denote the intermediate states in branch 1 or branch 2. In the model of Figure 2A, the second term is omitted because branch 2 has no intermediate states. Again, Eq. 49 can easily be extended to more than three branches between S a and S b .
Step 5: Knowing the reserve factors, the matrices of the reduced model can be generated (52) The double arrows include reserve factors according to Eq. 42a,b. The ordinary arrows are k i,j of the full n-state model. The double arrows between S Ã a and S Ã b are g a,b /r a or g b,a /r b (Eqs. 49 and 51). The double arrows starting from S Ã a to an individual state S a-1 are k a,a-1 /r a and those from S Ã b to S bC1 are k b,bC1 /r b .
Step 7: The current can be calculated by means of Eq. 46 for n Ã states.
with c i,iC1 D k i,iC1 and c iC1,i D k iC1,i if i not equal a (and iC1 not equal b), c a,b D g a,b , c b,a D g b,a , and index (n Ã C1) D index 1 (because of the cyclic scheme).
After incorporating Eqs. 48 to 56 in a fitting routine, the rate constants can be obtained from measured IV curves, and the occupation probabilities can by calculated. The intermediate states in branches have to be obtained from equations similar to Eq. 39. Fortunately, in most cases, it is not necessary to merge these equations into one comprehensive equation. Many fitting programs allow the model to be implemented as a set of shorter equations. In a fitting routine, the search algorithms suggests a set of rate constants k i,j for the full model like that in Figure 2A. Then, this set of rate constants is used to calculate the IV curve (or the tracer fluxes) from the set of equations.

Calculation of the reversal potential
In the case of the transport of one ion species, the kinetic approach does not yield any improvement over the calculation of the Nernst potential from cytosolic and external activity of the ion. This follows from the fact that the numerator of the flux equation (Eq. 56) has to be zero at the reversal potential. However, the situation may be different in the case of transporting different ions simultaneously as in the case of cotransporters. In this case, the intersection of the IV curve with the voltage axis has to be calculated from the fit parameters of an adequate branched model (e.g., done by ref. 40 for channelrhodopsin). Of course, this value has to be the same as that one of the measured curve, however, the theoretical background is different to that of the Goldman equation. 41 Furthermore, from this equation for the branched model, reversal potentials can be calculated for different ion concentrations than those used for determining the rate constants by fitting.

Conclusion
The reaction kinetic evaluation of IV curves or tracer experiments can reveal a plethora of information if global fitting is applied to sets of IV curves measured under different conditions. This approach may become very valuable when computational approaches like MD simulations will be able to predict occupation probabilities of intermediate states and rate constants of the transitions between those states. Then, they would build the bridge between physiology and structural approaches. With increasing number of experimental modifications of the system, the complexity of the reaction scheme increases. The arrow scheme provides the simplest approach to calculate gross rate constants, state occupancies and fluxes in cyclic reaction schemes with an arbitrary number of states. The scheme can also be applied to branched reaction schemes. However, in contrast to the purely cyclic scheme the branched schemes require some more steps when the equations are to be set-up and fitting results are evaluated.

Disclosure of potential conflicts of interest
No potential conflicts of interest were disclosed.