*latexText page 0 contains only floats

# Supporting Information for
Biological free energy transduction is an Achilles heel of mean-field
transport theory

Kiriko Terai  Department of Chemistry, Duke University, Durham, NC, 27708
Jonathon L. Yuly  Lewis-Sigler Institute for Integrative Genomics, Princeton
University, Princeton, NJ, 08540  Corresponding authors. Email:
<EMAIL_ADDRESS>or<EMAIL_ADDRESS>Peng Zhang  Department
of Chemistry, Duke University, Durham, NC, 27708  David N. Beratan  Department
of Chemistry, Duke University, Durham, NC, 27708  Department of Physics, Duke
University, Durham, NC, 27708  Department of Biochemistry, Duke University,
Durham, NC, 27710  Corresponding authors. Email<EMAIL_ADDRESS>or<EMAIL_ADDRESS>
##  1 Full master equation and mean-field kinetics for nanoscale particle
transport

Each microstate for a transport network with $N$ sites can be described using

$\boldsymbol{\sigma}=[\sigma_{1}\hskip 2.27621pt\sigma_{2}...\hskip
2.27621pt\sigma_{i}\hskip 2.27621pt...\hskip 2.27621pt\sigma_{j}\hskip
2.27621pt...\hskip 2.27621pt\sigma_{N}]$ (1)

where $\sigma_{i}$ and $\sigma_{j}$ are non-negative integers that indicate
the number of particles on sites $i$ and $j$ in microstate
$\boldsymbol{\sigma}$. The probability $P_{\boldsymbol{\sigma}}$ that the
network is in microstate $\boldsymbol{\sigma}$ as a function of time is
described by the master equation [1]:

$\frac{dP_{\boldsymbol{\sigma}}(t)}{dt}=-\sum_{\boldsymbol{\sigma^{\prime}}}(K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}P_{\boldsymbol{\sigma}}(t)-K_{\boldsymbol{\sigma^{\prime}}\boldsymbol{\sigma}}P_{\boldsymbol{\sigma^{\prime}}}(t))$
(2)

where $K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}$ is the rate
constant for the transition from microstate $\boldsymbol{\sigma}$ to
$\boldsymbol{{\sigma^{\prime}}}$. $P_{\boldsymbol{\sigma}}$ can also be
written as an $N$-point correlation function:

$P_{\boldsymbol{\sigma}}=p_{1\sigma_{1},2\sigma_{2},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}$
(3)

which is the probability of simultaneously finding sites 1, 2, …, $i$, …, $j$,
…, and $N$ in occupation states $\sigma_{1}$, $\sigma_{2}$, …, $\sigma_{i}$,
…, $\sigma_{j}$, …, and $\sigma_{N}$, respectively.

Site occupancy constraints are enforced by the values of the rate constants
$K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}$, such that a rate
constant $K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}$ is zero if
either microstate $\boldsymbol{\sigma}$ or $\boldsymbol{\sigma^{\prime}}$
violates the occupancy constraints for any site in the network. The following
indicates a transition between microstate $\boldsymbol{\sigma}$ and
$\boldsymbol{\sigma^{\prime}}$ associated with a typical $n$-particle
transport process within the network:

$\boldsymbol{\sigma}:[\sigma_{1}\hskip 2.27621pt\sigma_{2}...\hskip
2.27621pt\sigma_{i}\hskip 2.27621pt...\hskip 2.27621pt\sigma_{j}\hskip
2.27621pt...]\rightleftharpoons\boldsymbol{\sigma^{\prime}}:[\sigma_{1}\hskip
2.27621pt\sigma_{2}...\hskip 2.27621pt(\sigma_{i}-n)\hskip 2.27621pt...\hskip
2.27621pt(\sigma_{j}+n)\hskip 2.27621pt...]$ (4)

The rate constant for an $n$-particle transfer (beginning in microstate
$\boldsymbol{\sigma}$) from site $i$, initially in occupation state
$\sigma_{i}$, to site $j$, initially in occupation state $\sigma_{j}$, is
denoted by
$k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}$.

In addition to $n$-particle transport between sites, particles are allowed to
enter or leave the system through interactions with reservoirs (Eq 5), $n$
particles at a time. The following is a transition between microstate
$\boldsymbol{\sigma}$ and $\boldsymbol{\sigma^{\prime}}$ associated with a
typical $n$-particle transfer with a reservoir at site $i$:

$\boldsymbol{\sigma}:[\sigma_{1}\hskip 2.27621pt\sigma_{2}...\hskip
2.27621pt\sigma_{i}\hskip 2.27621pt...\hskip 2.27621pt\sigma_{j}\hskip
2.27621pt...]\rightleftharpoons\boldsymbol{\sigma^{\prime}}:[\sigma_{1}\hskip
2.27621pt\sigma_{2}...\hskip 2.27621pt(\sigma_{i}+n)\hskip 2.27621pt...\hskip
2.27621pt\sigma_{j}\hskip 2.27621pt...]$ (5)

The rate constant of an $n$-particle transfer from site $i$ (in microstate
$\boldsymbol{\sigma}$) , with occupancy $\sigma_{i}$, to a reservoir (or from
a reservoir to site $i$, initially in occupation state $\sigma_{i}$) is
denoted
$k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}$
($k^{\boldsymbol{\sigma}}_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}}$).

$K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}$, written in terms of
$k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}$,
$k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}$,
and
$k^{\boldsymbol{\sigma}}_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}}$,
is:

$\begin{split}K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}&=\sum_{\begin{subarray}{c}i\\\
j\neq i\\\
n\end{subarray}}k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}\hskip
2.27621pt\delta_{\sigma_{i},\sigma^{\prime}_{i}+n}\delta_{\sigma_{j}+n,\sigma^{\prime}_{j}}\prod_{l\neq
i,j}\delta_{\sigma_{l}\sigma^{\prime}_{l}}\\\
&+\sum_{i,n}k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}\hskip
2.27621pt\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\prod_{l\neq
i}\delta_{\sigma_{l}\sigma^{\prime}_{l}}\\\
&+\sum_{i,n}k^{\boldsymbol{\sigma}}_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}}\hskip
2.27621pt\delta_{\sigma_{i}+n,\sigma^{\prime}_{i}}\prod_{l\neq
i}\delta_{\sigma_{l}\sigma^{\prime}_{l}}\end{split}$ (6)

The probability of finding site $i$ with occupancy $q$ is denoted $p_{iq}$.
This average quantity is:

$p_{iq}=\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}P_{\boldsymbol{\sigma}}$
(7)

where $\delta$ is the Kronecker delta. Therefore, the exact master equation
(Eq 2) gives the exact time derivative of $p_{iq}$ as [1]:

$\begin{split}\frac{dp_{iq}}{dt}&=\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\frac{dP_{\boldsymbol{\sigma}}(t)}{dt}\\\
&=-\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\\{\sum_{\boldsymbol{\sigma^{\prime}}}(K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}P_{\boldsymbol{\sigma}}-K_{\boldsymbol{\sigma^{\prime}}\boldsymbol{\sigma}}P_{\boldsymbol{\sigma^{\prime}}})\\}\end{split}$
(8)

###  1.1 An exact formula for the time derivative of $p_{iq}$ from the master
equation

When site-site interactions are present, the rate of $n$ particle transfer
from site $i$ to site $j$ depends on the occupancy of all the $n$-particle
transfer rate constant depends on the initial and the final microstates.

Beginning with Eq 8, the exact time derivative of $p_{iq}$ when site-site
interactions are present is:

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\\{\sum_{\boldsymbol{\sigma^{\prime}}}(K_{\boldsymbol{\sigma}\boldsymbol{\sigma^{\prime}}}P_{\boldsymbol{\sigma}}-K_{\boldsymbol{\sigma^{\prime}}\boldsymbol{\sigma}}P_{\boldsymbol{\sigma^{\prime}}})\\}\\\
&=-\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\\\
&\indent\\{\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}[\sum_{\boldsymbol{\sigma^{\prime}}}(k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}P_{\boldsymbol{\sigma}}-k^{\boldsymbol{\sigma^{\prime}}}_{j,\sigma_{j}+n\underset{n}{\rightarrow}i,\sigma_{i}-n}P_{\boldsymbol{\sigma^{\prime}}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\delta_{\sigma_{j}+n,\sigma^{\prime}_{j}}\prod_{l\neq
i,j}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\\
&\indent+\sum_{n}[\sum_{\boldsymbol{\sigma^{\prime}}}(k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}P_{\boldsymbol{\sigma}}-k^{\boldsymbol{\sigma^{\prime}}}_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}-n}P_{\boldsymbol{\sigma^{\prime}}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\prod_{l\neq
i}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\}\\\
&=-\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\\\
&\indent\\{\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}[\sum_{\boldsymbol{\sigma^{\prime}}}(k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}-k^{\boldsymbol{\sigma^{\prime}}}_{j,\sigma_{j}+n\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{1\sigma^{\prime}_{1},...,i\sigma^{\prime}_{i},...,j\sigma^{\prime}_{j},...,N\sigma^{\prime}_{N}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\delta_{\sigma_{j}+n,\sigma^{\prime}_{j}}\prod_{l\neq
i,j}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\\
&\indent+\sum_{n}[\sum_{\boldsymbol{\sigma^{\prime}}}(k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}-k^{\boldsymbol{\sigma^{\prime}}}_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{1\sigma^{\prime}_{1},...,i\sigma^{\prime}_{i},...,j\sigma^{\prime}_{j},...,N\sigma^{\prime}_{N}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\prod_{l\neq
i}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\},\\\ \end{split}$ (9)

where the first step uses Eq 6 and last step uses Eq 3.

###  1.2 Derivation of Eq. 8 of the main text

In the absence of site-site interactions, the rate for $n$ particle transfer
from site $i$ to site $j$ does not depend on the occupancy of other sites in
the transport network. The $n$-particle transfer rate constant is independent
of the initial and the final microstates. In the absence of site interactions,
the rate constant for $n$-particle transfer from site $i$, initially in
occupation state $\sigma_{i}$, to site $j$, initially in occupation state
$\sigma_{j}$ is $k_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}$.
Beginning with Eq 9, the exact time derivative of $p_{iq}$ in the absence of
site-site interactions is:

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\\\
&\indent\\{\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}[\sum_{\boldsymbol{\sigma^{\prime}}}(k_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}-k_{j,\sigma_{j}+n\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{1\sigma^{\prime}_{1},...,i\sigma^{\prime}_{i},...,j\sigma^{\prime}_{j},...,N\sigma^{\prime}_{N}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\delta_{\sigma_{j}+n,\sigma^{\prime}_{j}}\prod_{l\neq
i,j}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\\
&\indent+\sum_{n}[\sum_{\boldsymbol{\sigma^{\prime}}}(k_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}-k_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{1\sigma^{\prime}_{1},...,i\sigma^{\prime}_{i},...,j\sigma^{\prime}_{j},...,N\sigma^{\prime}_{N}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\prod_{l\neq
i}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\}\\\ \end{split}$

Summing over the $\boldsymbol{\sigma}^{\prime}$ variables (applying the
Kronecker deltas) and using
$\sum_{\boldsymbol{\sigma}}=\sum_{\sigma_{1}}\cdot\cdot\cdot\sum_{\sigma_{i}}\cdot\cdot\cdot\sum_{\sigma_{j}}\cdot\cdot\cdot\sum_{\sigma_{N}}$
gives

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\sigma_{i}}\delta_{\sigma_{i},q}\\\
&\indent\\{\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}\sum_{\sigma_{j}}\prod_{l\neq
i,j}\sum_{\sigma_{l}}(k_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}-k_{j,\sigma_{j}+n\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{1\sigma_{1},...,i\sigma_{i}-n,...,j\sigma_{j}+n,...,N\sigma_{N}})\\\
&\indent+\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}\sum_{\sigma_{j}}\prod_{l\neq
i,j}\sum_{\sigma_{l}}(k_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}-k_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{i\sigma_{i}-n,j\sigma_{j},l\sigma_{l},...,N\sigma_{N}})\\}\\\
\end{split}$

Using the identities of Eq. 12 below gives

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\sigma_{i}}\delta_{\sigma_{i},q}\\\
&\indent\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}\sum_{\sigma_{j}}\left(k_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}p_{i\sigma_{i},j\sigma_{j}}-k_{j,\sigma_{j}+n\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{i\sigma_{i}-n,j\sigma_{j}+n}+(k_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}p_{i\sigma_{i},j\sigma_{j}}-k_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{i\sigma_{i}-n,j\sigma_{j}})\right).\end{split}$

The second term in the double sum above depends only on the probabilities
$p_{i\sigma_{i}}$ and $p_{i\sigma_{i}-n}$ when the sum over $\sigma_{j}$ is
carried out (because of the identity
$\sum_{\sigma_{j}}p_{i\sigma_{i},j\sigma_{j}}=p_{i\sigma_{i}}$). However, we
choose to keep the sum in the form of Eq 10 to obtain $w_{i,j,n,s}^{r,u}$ of
Eq 8 in the main text. Performing the sum over $\sigma_{i}$ and substituting
$\sigma_{j}$ with $s$ gives

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}\sum_{s}\left(k_{i,q\underset{n}{\rightarrow}j,s}p_{iq,js}-k_{j,s+n\underset{n}{\rightarrow}i,q-n}p_{iq-n,js+n}+(k_{i,q\underset{n}{\rightarrow}\text{out}}p_{iq,js}-k_{\text{in}\underset{n}{\rightarrow}i,q-n}p_{iq-n,js})\right)\\\
&=-\sum_{\begin{subarray}{c}j\neq i\\\
n,s\end{subarray}}\sum_{r,u}w_{i,j,n,s}^{r,u}\hskip 3.0ptp_{ir,ju}\hskip
10.0pt\end{split}$ (10)

where $w_{i,j,n,s}^{r,u}$ is defined by

$\begin{split}w_{i,j,n,s}^{r,u}&=\delta_{r,q}\delta_{u,s}k_{i,r\underset{n}{\rightarrow}j,u}-\delta_{r,q-n}\delta_{u,s+n}k_{j,u\underset{n}{\rightarrow}i,r}+\delta_{r,q}\delta_{u,s}k_{i,r\underset{n}{\rightarrow}\text{out}}-\delta_{r,q-n}\delta_{u,s}k_{\text{in}\underset{n}{\rightarrow}i,q-n}.\end{split}$
(11)

The derivation of Eq. 10 uses the identities

$\begin{split}\prod_{l\neq
i,j}\sum_{\sigma_{l}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}&=p_{i\sigma_{i},j\sigma_{j}}\\\
\sum_{\sigma_{j}}\prod_{l\neq
i,j}\sum_{\sigma_{l}}p_{1\sigma_{1},...,i\sigma_{i},...,j\sigma_{j},...,N\sigma_{N}}&=p_{i\sigma_{i}}\\\
\end{split}$ (12)

The mean-field kinetics in the absence of interactions (Eq. 9 in the main
text) follows immediately from Eq. 10 on using the mean-field approximation
$p_{ir,ju}\approx p_{ir}p_{ju}$ (Eq. 13), as shown in Section 1.4.

###  1.3 Mean-field kinetics with site-site interactions

For a transport network with $N$ sites for the particles, the mean-field
approximation neglects all statistical correlations among site occupancies.
That is, any $X$-point correlation function (Eq 3) with $X\leq N$, can be
factored into individual probabilities [1, 2, 3]:

$p_{i_{1}q_{1},i_{2}q_{2},...,i_{X}q_{X}}\approx\prod_{x=1}^{X}p_{i_{x}q_{x}}\hskip
28.45274pt(X\leq N)$ (13)

Eq 13 defines the mean-field approximation.

Beginning with Eq. 9 and using the mean-field approximation, the mean-field
time derivative of $p_{iq}$ including site-site interactions is:

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\\\
&\indent\\{\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}[\sum_{\boldsymbol{\sigma^{\prime}}}(k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}p_{1\sigma_{1}}...p_{i\sigma_{i}}...p_{j\sigma_{j}}...p_{N\sigma_{N}}\\\
&\indent-k^{\boldsymbol{\sigma^{\prime}}}_{j,\sigma_{j}+n\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{1\sigma^{\prime}_{1}}...p_{i\sigma_{i}^{\prime}}...p_{j\sigma^{\prime}_{j}}...p_{N\sigma^{\prime}_{N}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\delta_{\sigma_{j}+n,\sigma^{\prime}_{j}}\prod_{l\neq
i,j}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\\
&\indent+\sum_{n}[\sum_{\boldsymbol{\sigma^{\prime}}}(k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}p_{1\sigma_{1}}...p_{i\sigma_{i}}...p_{j\sigma_{j}}...p_{N\sigma_{N}}-k^{\boldsymbol{\sigma^{\prime}}}_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}-n}p_{1\sigma^{\prime}_{1}}...p_{i\sigma^{\prime}_{i}}...p_{j\sigma^{\prime}_{j}}...p_{N\sigma^{\prime}_{N}})\delta_{\sigma_{i}-n,\sigma^{\prime}_{i}}\prod_{l\neq
i}\delta_{\sigma_{l},\sigma^{\prime}_{l}}]\\}.\\\ \end{split}$

Summing over the $\boldsymbol{\sigma}^{\prime}$ variables (applying the
Kronecker deltas) and using
$\sum_{\boldsymbol{\sigma}}=\sum_{\sigma_{1}}\cdot\cdot\cdot\sum_{\sigma_{i}}\cdot\cdot\cdot\sum_{\sigma_{j}}\cdot\cdot\cdot\sum_{\sigma_{N}}$
gives

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\sigma_{i}}\delta_{\sigma_{i},q}\\\
&\indent\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}\sum_{\sigma_{j}}\\{\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}j,\sigma_{j}}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{i\sigma_{i}}p_{j\sigma_{j}}-\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}^{f}}_{j,\sigma_{j}+n\underset{n}{\rightarrow}i,\sigma_{i}-n}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{i\hskip 1.42262pt\sigma_{i}-n}p_{j\hskip
1.42262pt\sigma_{j}+n}\\}\\\ &\indent+\sum_{\begin{subarray}{c}j\neq i\\\
n\end{subarray}}\sum_{\sigma_{j}}\\{\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}}_{i,\sigma_{i}\underset{n}{\rightarrow}\text{out}}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{i\sigma_{i}}p_{j\sigma_{j}}-\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}^{f^{\prime}}}_{\text{in}\underset{n}{\rightarrow}i,\sigma_{i}-n}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{i\hskip 1.42262pt\sigma_{i}-n}p_{j\sigma_{j}}\\},\\\
\end{split}$

where $\boldsymbol{\sigma}^{f}$ and $\boldsymbol{\sigma}^{f^{\prime}}$ are
given by

$\begin{split}\boldsymbol{\sigma}^{f}&=[\sigma_{1},...,\sigma_{i}-n,...,\sigma_{j}+n,...,\sigma_{N}]\\\
\boldsymbol{\sigma}^{f^{\prime}}&=[\sigma_{1},...,\sigma_{i}-n,...,\sigma_{j},...,\sigma_{N}].\end{split}$
(14)

Performing the sum over $\sigma_{i}$ and substituting $\sigma_{j}$ with $s$
gives

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\begin{subarray}{c}j\neq i\\\
n,s\end{subarray}}\\{\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}}_{i,q\underset{n}{\rightarrow}j,s}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{iq}p_{js}-\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}^{f}}_{j,s+n\underset{n}{\rightarrow}i,q-n}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{i\hskip 1.42262ptq-n}p_{j\hskip 1.42262pts+n}\\}\\\
&\indent+\sum_{\begin{subarray}{c}j\neq i\\\ n,s\end{subarray}}\\{\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}}_{i,q\underset{n}{\rightarrow}\text{out}}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{iq}p_{js}-\prod_{l\neq
i,j}[\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}^{f^{\prime}}}_{\text{in}\underset{n}{\rightarrow}i,q-n}\prod_{u\neq
i,j}p_{u\sigma_{u}}]p_{i\hskip 1.42262ptq-n}p_{js}\\}\\\
&=-\sum_{\begin{subarray}{c}j\neq i\\\
n,s\end{subarray}}[k^{\text{MF}}_{i,q\underset{n}{\rightarrow}j,s}p_{iq}p_{js}-k^{\text{MF}}_{j,s+n\underset{n}{\rightarrow}i,q-n}p_{i\hskip
1.42262ptq-n}p_{j\hskip 1.42262pts+n}]\\\
&\indent+\sum_{\begin{subarray}{c}j\neq i\\\
n,s\end{subarray}}[k^{\text{MF}}_{i,q\underset{n}{\rightarrow}\text{out}}p_{iq}p_{js}-k^{\text{MF}}_{\text{in}\underset{n}{\rightarrow}i,q-n}p_{i\hskip
1.42262ptq-n}p_{js}]\\\ &=-\sum_{\begin{subarray}{c}j\neq i\\\
n,s\end{subarray}}\sum_{r,u}w^{r,u(\text{MF})}_{i,j,n,s}p_{ir}p_{ju}\end{split}$
(15)

where the mean-field approximated $n$-particle transfer rate constants are:

$\begin{split}k^{\text{MF}}_{i,q\underset{n}{\rightarrow}j,s}&=\prod_{l\neq
i,j}\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}}_{i,q\underset{n}{\rightarrow}j,s}\prod_{u\neq
i,j}p_{u\sigma_{u}}\\\
k^{\text{MF}}_{j,s+n\underset{n}{\rightarrow}i,q-n}&=\prod_{l\neq
i,j}\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}^{f}}_{j,s+n\underset{n}{\rightarrow}i,q-n}\prod_{u\neq
i,j}p_{u\sigma_{u}}\\\
k^{\text{MF}}_{i,q\underset{n}{\rightarrow}\text{out}}&=\prod_{l\neq
i,j}\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}}_{i,q\underset{n}{\rightarrow}\text{out}}\prod_{u\neq
i,j}p_{u\sigma_{u}}\\\
k^{\text{MF}}_{\text{in}\underset{n}{\rightarrow}i,q-n}&=\prod_{l\neq
i,j}\sum_{\sigma_{l}}k^{\boldsymbol{\sigma}^{f^{\prime}}}_{\text{in}\underset{n}{\rightarrow}i,q-n}\prod_{u\neq
i,j}p_{u\sigma_{u}}\end{split}$ (16)

and $w^{r,u(\text{MF})}_{i,j,n,s}$ is defined as

$\begin{split}w^{r,u(\text{MF})}_{i,j,n,s}=\delta_{r,q}\delta_{u,s}k^{\text{MF}}_{i,r\underset{n}{\rightarrow}j,u}-\delta_{r,q-n}\delta_{u,s+n}k^{\text{MF}}_{j,u\underset{n}{\rightarrow}i,r}+\delta_{r,q}\delta_{u,s}k^{\text{MF}}_{i,r\underset{n}{\rightarrow}\text{out}}-\delta_{r,q-n}\delta_{u,s}k^{\text{MF}}_{\text{in}\underset{n}{\rightarrow}i,q-n}\hskip
10.0pt.\end{split}$ (17)

###  1.4 The mean-field kinetics without site-site interactions

Using the mean-field approximation (Eq 13) with Eq 10 gives the mean-field
time derivative of $p_{iq}$ when site-site interactions are absent:

$\begin{split}\frac{dp_{iq}}{dt}&=-\sum_{\begin{subarray}{c}j\neq i\\\
n,s\end{subarray}}\left(k_{i,q\underset{n}{\rightarrow}j,s}p_{iq}p_{js}-k_{j,s+n\underset{n}{\rightarrow}i,q-n}p_{iq-n}p_{js+n}+(k_{i,q\underset{n}{\rightarrow}\text{out}}p_{iq}p_{js}-k_{\text{in}\underset{n}{\rightarrow}i,q-n}p_{iq-n}p_{js})\right)\\\
&=-\sum_{\begin{subarray}{c}j\neq i\\\
n,s\end{subarray}}\sum_{r,u}w^{r,u(\text{MF})}_{i,j,n,s}p_{ir}p_{ju}=-\sum_{\begin{subarray}{c}j\neq
i\\\ n,s\end{subarray}}\sum_{r,u}w^{r,u}_{i,j,n,s}p_{ir}p_{ju}\end{split}$
(18)

Thus, in the absence of site-site interactions,

$\begin{split}k^{\text{MF}}_{i,q\underset{n}{\rightarrow}j,s}&\rightarrow
k_{i,q\underset{n}{\rightarrow}j,s}\\\
k^{\text{MF}}_{j,s+n\underset{n}{\rightarrow}i,q-n}&\rightarrow
k_{j,s+n\underset{n}{\rightarrow}i,q-n}\\\
k^{\text{MF}}_{i,q\underset{n}{\rightarrow}\text{out}}&\rightarrow
k_{i,q\underset{n}{\rightarrow}\text{out}}\\\
k^{\text{MF}}_{\text{in}\underset{n}{\rightarrow}i,q-n}&\rightarrow
k_{\text{in}\underset{n}{\rightarrow}i,q-n}\end{split}$ (19)

so $w^{r,u(\text{MF})}_{i,j,n,s}\rightarrow w^{r,u}_{i,j,n,s}$, and one
obtains Eq 9 in the main text, that is the mean-field kinetics in the absence
of site-site interactions.

##  2 The mean-field approximation is exact at equilibrium without site-site
interactions

The mean-field approximation for transport kinetics is exact at thermal
equilibrium, provided there are no site-site interactions. This can be
understood by noting that the microstates $\boldsymbol{\sigma}$, and their
total energies, map directly onto those of a non-interacting spin chain. In
particular, the state occupancies $s_{j}$ may be mapped onto the magnetic
quantum numbers $m_{s}^{j}$ of spins $j$. Non-interacting spins are known to
be uncorrelated at thermal equilibrium [4]. We provide a mathematical proof of
this well-known result [5, 6].

Each microstate of a transport network is described using

$\boldsymbol{\sigma}=[\sigma_{1}\hskip 2.27621pt\sigma_{2}...\hskip
2.27621pt\sigma_{i}\hskip 2.27621pt...\hskip 2.27621pt\sigma_{j}\hskip
2.27621pt...]$ (20)

where $\sigma_{i}$ and $\sigma_{j}$ are non-negative integers that indicate
the number of particles on sites $i$ and $j$. Therefore, microstate
$\boldsymbol{\sigma}$ thus tracks the occupation states of all sites. At
thermodynamic equilibrium in the grand canonical ensemble, the probability to
occupy microstate $\boldsymbol{\sigma}$ ($P_{\boldsymbol{\sigma}}$) and the
partition function ($Z$) of a many-site transport chain are:

$P_{\boldsymbol{\sigma}}=\frac{\exp({-\beta
G_{\boldsymbol{\sigma}})}}{Z}=\frac{\exp{(-\beta\sum_{k}G^{\boldsymbol{\sigma}}_{k})}}{Z}=\frac{\prod_{k}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}}{Z}$ (21)
$Z=\sum_{\boldsymbol{\sigma}}\exp{(-\beta\sum_{k}G^{\boldsymbol{\sigma}}_{k})}=\sum_{\boldsymbol{\sigma}}\prod_{k}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}=\prod_{k}\sum_{\boldsymbol{\sigma}}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}=\prod_{k}Z_{k}$ (22)

where

$Z_{k}=\sum_{\boldsymbol{\sigma}}\exp{(-\beta G^{\boldsymbol{\sigma}}_{k})}$
(23)

is the partition function associated with site $k$. The quantity
$G_{\boldsymbol{\sigma}}$ is the free energy associated with microstate
$\boldsymbol{\sigma}$ and $G^{\boldsymbol{\sigma}}_{k}$ is the free energy
associated with site $k$ in any microstate $\boldsymbol{\sigma}$. Thus, at
equilibrium, the probability of finding site $i$ with occupancy $q$ ($p_{iq}$)
is

$\begin{split}p_{iq}&=\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}P_{\boldsymbol{\sigma}}\\\
&=\frac{1}{Z}\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\prod_{k}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}\\\
&=\frac{1}{\sum_{\boldsymbol{\sigma}}\prod_{k}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}}\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\prod_{k}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}\\\
&=\frac{\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{i})}\prod_{k\neq i}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}}{\sum_{\boldsymbol{\sigma}}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{i})}\prod_{k\neq i}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}}\\\
&=\frac{\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i},q}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{i})}}{Z_{i}}\end{split}$ (24)

The joint probability to find sites $i_{1},...,i_{X}$ with occupancies
$q_{1},...,q_{X}$ (i.e., the $X$-point correlation function
$p_{i_{1}q_{1},...,i_{X}q_{X}}$) at equilibrium is

$\begin{split}p_{i_{1}q_{1},...,i_{X}q_{X}}&=\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i_{1}},q_{1}}\cdot\cdot\cdot\delta_{\sigma_{i_{X}},q_{X}}P_{\boldsymbol{\sigma}}\\\
&=\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i_{1}},q_{1}}\cdot\cdot\cdot\delta_{\sigma_{i_{X}},q_{X}}\frac{\prod_{k}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}}{Z}\\\
&=\frac{1}{\prod_{k}Z_{k}}\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i_{1}},q_{1}}\cdot\cdot\cdot\delta_{\sigma_{i_{X}},q_{X}}\prod_{k}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{k})}\\\
&=\frac{1}{\prod_{k}Z_{k}}\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i_{1}},q_{1}}\cdot\cdot\cdot\delta_{\sigma_{i_{X}},q_{X}}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{i_{1}})}\cdot\cdot\cdot\exp{(-\beta
G^{\boldsymbol{\sigma}}_{i_{X}})}\prod_{k\neq\\{i_{1},...,i_{X}\\}}Z_{k}\\\
&=\left(\frac{1}{Z_{i_{1}}}\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i_{1}},q_{1}}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{i_{1}})}\right)\cdot\cdot\cdot\left(\frac{1}{Z_{i_{X}}}\sum_{\boldsymbol{\sigma}}\delta_{\sigma_{i_{X}},q_{X}}\exp{(-\beta
G^{\boldsymbol{\sigma}}_{i_{X}})}\right)\\\ &=p_{i_{1}q_{1}}\cdot\cdot\cdot
p_{i_{X}q_{X}}\\\ &=\prod_{x=1}^{X}p_{i_{x}q_{x}}.\end{split}$ (25)

Thus, the mean-field approximation (Eq. 13) is exact at equilibrium for a
classical transport chain without site-site interactions.

##  3 Hill’s model for biological energy transduction

Free energy transduction occurs when one or more exergonic reactions ($\Delta
G<0$ or “downhill”) are leveraged to drive one or more endergonic reactions
($\Delta G>0$ or “uphill”) [7]. Reaction cycles that only perform exergonic
processes without driving endergonic cycles are known as slippage events,
dissipating otherwise available free energy as heat [7]. Slippage reduces the
turnover rates of the endergonic processes, so efficient free energy
transduction occurs when slippage is suppressed.

A useful model to understand biological free energy transduction was described
by T.L. Hill [7], and Fig. S1 summarizes his model. This model maps onto an
equivalent model for many-particle transport, by assigning the microstates of
the enzyme and binding sites to a particle transport network with isomorphic
kinetics. We denote the states of the sites in this mapping by ($\sigma_{L}$,
$\sigma_{E}$, $\sigma_{M}$). The occupancy of the L and M binding sites,
$\sigma_{L}$ and $\sigma_{M}$ respectively, indicate whether L and M are bound
to E (0 if no ligand is bound, 1 if a ligand is bound). An additional site is
used to denote the state of the enzyme $E$ in its two available conformations
(E:$\sigma_{E}=0$ and E*:$\sigma_{E}=1$). Importantly, particles cannot move
between the sites in this mapping, as such processes are unphysical (that is,
M is eliminated and converted into bound species L). However, the particle
sites in the mapping do influence each other through interactions. The L and M
concentration gradients across the membrane, and the energy of the
E$\rightarrow$E* conformational change, are reflected in the energies of the
reservoirs that are coupled to the three sites. There are eight microstates in
this model. The redox-coupled proton pump described in Figure 3 of the main
text is analogous to Hill’s model (Fig S1), albeit with additional transitions
between the states (diagonals of the cube faces).

##  4 Electron transfer through redox cofactor chains

Electron transport through proteins is often mediated by a chain of multiple
redox cofactors between two or more substrates that interact with the cofactor
chain. We model the cofactors as hopping sites and the substrates as electron
reservoirs. In studies of electron transport, the electron hopping rate from
site $i$ to site $j$ may be approximated by [1, 8, 9, 10, 11]

$k_{ij}=\frac{2\pi}{\hbar}|H_{ij}|^{2}\frac{1}{\sqrt{4\pi\lambda_{ij}k_{B}T}}\exp(-\frac{((\epsilon_{j}-\epsilon_{i})+\lambda_{ij})^{2}}{4\lambda_{ij}k_{B}T})$
(26)

where $H_{ij}$ (Å) and $\lambda_{ij}$ (eV) are the distance dependent
electronic coupling and reorganization energy for electron transfer between
sites $i$ and $j$. $k_{B}T$ is the thermal energy. The electronic coupling is
approximated as $|H_{ij}|^{2}\approx V^{2}\exp(-\beta R)$ ($\beta$ is the
tunneling decay factor) and $V$ is the electronic coupling at contact between
cofactors $i$ and $j$. The parameters were approximated as $\beta$ = 1 Å-1,
$V$ = 0.1 eV, $\lambda$ = 1.0 eV, and $\nicefrac{{1}}{{k_{B}T}}$ = 39.06
$\text{eV}^{-1}$, which are typical for protein electron transfer [11, 12,
13]. The distances between cofactors in the electron bifurcation model were
set to typical values for electron bifurcating enzymes. We set $R$ = 10 Å for
nearest-neighbor cofactors, $R$ = 20 Å for second nearest-neighbors, and $R$ =
30 Å for third nearest-neighbors, as in a previous model [14].

In the electron bifurcation model (Fig. 5 of the main text), the high- and
low-potential substrates $A_{H}$ and $A_{L}$, respectively, were modeled as
one electron exchange reservoirs, and the two-electron substrate $D$ was
modeled as a two-electron exchange reservoir ($n=2$ in Eq 6). The rate
constants for electron transfer from the system to the reservoirs were set to
be fast ($k_{i,q\underset{n}{\rightarrow}\text{out}}=10^{7}s^{-1}$) so they
are not not rate-limiting, and the rate constants for electron transfer into
the system from these reservoirs were determined using Eq 27. The rate
constant for $n$ electrons to transfer into a reservoir from site $i$ in redox
state $q$, $k_{\text{in}\underset{n}{\rightarrow}i,q-n}$, is [14]:

$k_{\text{in}\underset{n}{\rightarrow}i,q-n}=k_{i,q\underset{n}{\rightarrow}\text{out}}\exp\left[\frac{n\times(\mu_{\text{res}}+|e^{-}|E_{i}^{\circ\text{mid}})}{k_{B}T}\right]$
(27)

where $\mu_{\text{res}}$ is the chemical potential of electrons in the
reservoir, and $E_{i}^{\circ\text{mid}}$ is the average midpoint redox
potential of cofactor $i$ (our model with rate constants given by Eq 26 and Eq
27 satisfy local detailed balance [15]) over all $n$-electron transfers with
the reservoir.

The exact kinetics within this theoretical framework is described by the
master equation (Eq 2 in main text), using rate constants defined by Eq 26 and
27. The mean-field approximated kinetics is described by the approximate rate
equation (Eq 18) with rate constants defined by equations 26 and 27,
particularly noting that no interactions between redox cofactors are included
in this model. To simulate the linear transport chain of Fig. 1 of the main
text, the parameters of Eq 26 above were taken from calculations [8] and
experiment [16], summarized in Table S1.

##  5 Summary of the kinetics model of Jiang et al.

This section describes the kinetic approach used by Jiang et al. [8] and its
differences from the mean-field approximation in the presence of interactions
between particle binding sites. Jiang et al. included interactions between
sites by making the driving force for electron transfer depende on site
occupancies. Jiang et al. made $\epsilon_{i}$ (and $\epsilon_{j}$) in Eq. 26
dependent on the site occupations using:

$\epsilon^{\circ}_{i}=\epsilon_{i}^{r}+\sum_{j\neq
i}(\Delta\epsilon_{ij})(1-p_{j1})$ (28)

where $\epsilon_{i}^{r}$ is the redox potential of site $i$ when all sites
(including site $i$) are occupied, $\Delta\epsilon_{ij}$ is the change in
redox potential of site $i$ as site $j$ changes from occupied to unoccupied,
and $p_{j1}$ is the probability of site $j$ being occupied. Table S1 shows
$\Delta\epsilon_{ij}$ values taken from experiments [16] (and used for both
the mean-field model, and the model of Jiang et al. [8]) for the STC protein.
The approach of Jiang et al. uses the equations

$\begin{split}\frac{dp_{i1}}{dt}&=-\sum_{j\neq
i}[k^{\circ}_{i,1\underset{1}{\rightarrow}j,0}p_{i1}p_{j0}-k^{\circ}_{j,1\underset{1}{\rightarrow}i,0}p_{i0}p_{j1}]\\\
&\indent+k^{\circ}_{i,1\underset{1}{\rightarrow}\text{out}}p_{i1}-k^{\circ}_{\text{in}\underset{1}{\rightarrow}i,0}p_{i0}\end{split}$
(29)

where $k^{\circ}_{i,q\underset{n}{\rightarrow}j,s}$ is the electron transfer
rate (Eq. 26) based on the modified reduction potentials
$\epsilon^{\circ}_{i}$ and $\epsilon^{\circ}_{j}$ values of Eq. 28 (the
reservoir rates $k^{\circ}_{i,q\underset{1}{\rightarrow}\text{out}}$ and
$k^{\circ}_{\text{in}\underset{1}{\rightarrow}i,q-1}$ were modeled as
irreversible rates that are independent of the site occupancies, which we
reproduce in the limit of strong driving forces, infinite $\Delta G$ in Fig.
S2). The kinetics of Eq 29 is clearly distinct from the mean-field kinetics
(Eq. 15 and 16), except when there are no interactions (Eq. 18). This case is
relevant for our electron bifurcation model, where the mean-field
approximation and the approach of Jiang et al. are identical. The results for
the linear transport chain that are derived from the exact master equation and
the mean-field approximation are shown in Fig. 2 of the main text, and the
results using the approach of Jiang et al. are shown in Fig. S2.

We extended the approach of Jiang et al. to study the redox-coupled proton
pump. In this model (Fig. 3 of the main text), the rate constant for a proton
to transfer from site $i$ to site $j$ is [17, 18]

$k_{i,1\underset{1}{\rightarrow}j,0}=\kappa_{ij}\exp{(-\frac{E_{j}-E_{i}}{2k_{B}T})}$
(30)

The reverse rate constant is

$k_{j,1\underset{1}{\rightarrow}i,0}=\kappa_{ij}\exp{(-\frac{E_{i}-E_{j}}{2k_{B}T})},$
(31)

and the rate constant for a single proton or electron transfer from a
reservoir to site $i$ is

$k_{\text{in}\underset{1}{\rightarrow}i,0}=\kappa_{ii}\exp{(\frac{E_{i}}{2k_{B}T})}$
(32)

$\kappa_{ij}$ is the rate constant between sites $i$ and $j$ when all sites,
other than site $i$, are unoccupied. $\kappa_{ii}$ is the rate constant for a
single proton or electron transfer from a reservoir to site $i$ when all sites
are unoccupied, and $E_{i}$($E_{j}$) is the energy of site $i$($j$). As in the
approach of Jiang et al., the electrostatic interactions between sites are
included by making the driving force for a proton or electron transfer
dependent on the site occupancies (i.e. making $E_{i}$ (and $E_{j}$) in Eq. 30
and 32 dependent on site occupancies). Studies by Kim et al. [17, 18] modified
$E_{i}$ to be

$E^{\circ}_{i}=E_{i}^{*}+\sum_{j\neq
i}\epsilon_{ij}p_{j1}+\frac{q_{i}Z_{i}V_{m}}{L}$ (33)

where $E_{i}^{*}$ is the free energy to place a particle on site $i$ while all
other sites are unoccupied and the membrane potential $V_{m}$ is zero,
$\epsilon_{ij}$ is the electrostatic interaction between sites $i$ and $j$,
$p_{j1}$ is the probability to find site $j$ occupied, $q_{i}$ is the charge
of site $i$, $Z_{i}$ is the distance from the N-side to site $i$, and $L$ is
the membrane thickness. Table S2 shows the relevant parameters. The kinetics
are computed using Eq. 29 with the modified site energies $E^{\circ}_{i}$ and
$E^{\circ}_{j}$ of Eq. 33. The results for the analysis of the redox-coupled
proton pump based on the exact master equation and the mean-field
approximation are shown in Fig. 4 of the main text, and the results using the
approach of Jiang et al. are shown in Fig. S3.

##  6 Effective activation energies for short-circuiting in the mean-field
approximation

This section explores the short-circuit flux in the electron bifurcation
network of Fig. 5 when $\Delta G_{bifurc}=0$ (i.e., “idle short-circuiting”
[14]). We find that mean-field kinetic analysis fails to produce two separate
short-circuiting regimes as a function of temperature (see Fig. S5 below).
These regimes are determined by two separate short-circuit channels
(illustrated in Fig. S6a and S6b) that have different activation energies and
hence dominate the kinetics at different temperatures. The electron transfer
steps that occur in three possible short-circuit channels are shown in Fig.
S6. In contrast, the mean-field approximation to the kinetics predicts one
regime with an effective activation energy that is the average of the
activation energies of the regimes from the exact master equation (purple
curve in Fig. S5).

Earlier studies found that the short-circuiting flux predicted by the exact
master equation agrees well with the estimated flux through the dominant
short-circuit channel at a given temperature, given by [14, 19]

$\text{flux}\sim k_{rl}\exp(-\frac{\Delta G^{\ddagger}}{k_{B}T}),$ (34)

where $\Delta G^{\ddagger}=\Delta G^{e^{-}}+\Delta G^{h^{+}}$ is the free
energy that moves an electron and hole to the cofactors required for short-
circuiting to occur, and $k_{rl}$ is the rate-limiting electron transfer rate
constant of the short-circuiting channel(Fig. S6). Thus, short-circuits are
thermally-activated processes [14]. We define $J1$ and $J2$ as the short-
circuiting electron fluxes through the orange and green channels indicated in
Fig. 5 of the main text, respectively. The factor of $2$ in the equation for
$J1$ arises from the fact that there are two-short circuits $B^{-}\rightarrow
H_{1}$ and $L_{1}^{-}\rightarrow B^{-}$ that both contribute equally to $J1$.
Using Eq. 34, the flux $J1$ is approximated as

$\begin{split}J1&=\text{flux}_{B^{-}\rightarrow
H1}+\text{flux}_{L1^{-}\rightarrow B^{-}}\\\ &\approx k_{L1,1\rightarrow
B,0}\exp(-\frac{\Delta G_{B^{-}\rightarrow
H1}^{e^{-}}}{k_{B}T})\exp(-\frac{\Delta G_{B^{-}\rightarrow
H1}^{h^{+}}}{k_{B}T})\\\ &\hskip 15.0pt+k_{B,2\rightarrow
H1,0}\exp(-\frac{\Delta G_{L1^{-}\rightarrow
B^{-}}^{e^{-}}}{k_{B}T})\exp(-\frac{\Delta G_{L1^{-}\rightarrow
B^{-}}^{h^{+}}}{k_{B}T})\\\ &=k_{L1,1\rightarrow B,0}\exp(-\frac{\Delta
G_{B^{-}\rightarrow H1}^{\ddagger}}{k_{B}T})+k_{B,2\rightarrow
H1,0}\exp(-\frac{\Delta G_{L1^{-}\rightarrow B^{-}}^{\ddagger}}{k_{B}T})\\\
&=2k^{J1}_{rl}\exp(-\frac{\Delta G^{\ddagger}_{J1}}{k_{B}T})\end{split}$ (35)

where $\Delta G^{\ddagger}_{SC}=\Delta G^{e^{-}}_{SC}+\Delta G^{h^{+}}_{SC}$
$(SC=B^{-}\rightarrow H1,L1^{-}\rightarrow B^{-})$ is the energy required to
move a hole from $A_{H}$ and an electron from $A_{L}$ to realize the precursor
state from which the rate-limiting process with rate constant $k_{rl}$
initiates the short-circuiting (see Fig S6). $\Delta G^{\ddagger}_{SC}$ is a
product of Boltzmann weights for particle occupancies in each branch. There
are two short-circuits that contribute equally to $J1$ (see Fig. 5 of the main
text), and their contributions are identical for the parameters of the
electron bifurcation model used (see Fig. 5 of the main text). Thus,
$k^{J1}_{rl}=k_{L1,1\rightarrow B,0}=k_{B,2\rightarrow H1,0}$ and $\Delta
G^{\ddagger}_{J1}=\Delta G_{B^{-}\rightarrow H1}^{\ddagger}=\Delta
G_{L1^{-}\rightarrow B^{-}}^{\ddagger}$ Similarly, the $J2$ flux is
approximated as

$\begin{split}J2&=\text{flux}_{L1^{-}\rightarrow H1}\\\ &\approx
k_{L1,1\rightarrow H1,0}\exp(-\frac{\Delta G_{L1^{-}\rightarrow
H1}^{e^{-}}}{k_{B}T})\exp(-\frac{\Delta G_{L1^{-}\rightarrow
H1}^{h^{+}}}{k_{B}T})\\\ &=k_{L1,1\rightarrow H1,0}\exp(-\frac{\Delta
G_{L1^{-}\rightarrow H1}^{\ddagger}}{k_{B}T})\\\ &=k_{L1,1\rightarrow
H1,0}\exp(-\frac{\Delta G^{\ddagger}_{J2}}{k_{B}T})\end{split}$ (36)

These expressions for the $J1$ and $J2$ fluxes are plotted as orange and green
dashed lines, respectively, in Fig. S5. The activation energies $\Delta
G^{\ddagger}_{J1}$ and $\Delta G^{\ddagger}_{J2}$ are estimated by analyzing
the energy landscape for the electrons in the system, again as illustrated in
Fig. S6. Since the activation energy for $J2$ is lower than that for $J1$, the
$J2$ flux dominates at low temperature (large $1/k_{B}T$).

We were intrigued to find that the mean-field approximated short-circuiting
flux has an effective activation energy (slope of the purple line in Fig. S5)
that is approximately the mean of the slopes of $J1$ and $J2$ on a log scale
(supplementary Fig. S5). Defining $J3$ as the geometric mean of $J1$ and $J2$,

$J3=(J1J2)^{1/2}=\sqrt{2k_{L1,1\rightarrow D,0}k_{L1,1\rightarrow
H1,0}}\exp(-\frac{\Delta G^{\ddagger}_{J1}+\Delta
G^{\ddagger}_{J2}}{2k_{B}T})$ (37)

we find that $J3$ has the activation energy $\frac{\Delta
G^{\ddagger}_{J1}+\Delta G^{\ddagger}_{J2}}{2k_{B}T}$ that is the average of
the activation energies of $J1$ and $J2$ ($J3$ is indicated with the pink line
in Fig. S5) Since the purple curve (mean-field flux) is nearly parallel to
$J3$ (supplementary Fig. S5), the effective activation energy of the flux
calculated with the mean-field approximation is approximately the geometric
mean of $J1$ and $J2$, but the mean-field flux is not equal to $J3$, having a
different $Y$-intercept in Fig. S5.

The underlying cause of the mean-field approximation to fail to capture the
different short-circuiting regimes as a function of temperature remains
unclear. As mentioned in the main text, mean-field failure may involve
extending the near-equilibrium statistics within each branch unnaturally
across the two branches. Understanding why the effective activation energy of
the mean-field flux is approximately the average of the activation energies
for the orange and green short-circuiting channels in Fig. S5 and is of
interest for future analysis.

##  7 Supplementary tables and figures

Table S1:  Parameters for the STC linear transport chain taken from previous
studies [8, 20]. $\lambda_{ij}$ is the reorganization energy (meV), and
$H_{ij}$ is the electronic coupling (meV) between cofactors $i$ and $j$.
$\epsilon_{i}^{r}$ is the redox potential (meV) of cofactor $i$ when all other
sites are occupied, and $\Delta\epsilon_{ij}$ is the change in redox potential
(meV) of cofactor $i(j)$ when $j(i)$ is unoccupied.

$\lambda_{12}$ | 1080 meV
---|---
$\lambda_{23}$ | 760 meV
$\lambda_{34}$ | 880 meV
$H_{12}$ | 2.17 meV
$H_{23}$ | 3.08 meV
$H_{34}$ | 2.08 meV

$\epsilon_{1}^{r}$ | -243 meV
---|---
$\epsilon_{2}^{r}$ | -222 meV
$\epsilon_{3}^{r}$ | -189 meV
$\epsilon_{4}^{r}$ | -171 meV
$\Delta\epsilon_{12}$ | 28 meV
$\Delta\epsilon_{13}$ | 21 meV
$\Delta\epsilon_{14}$ | 11 meV
$\Delta\epsilon_{23}$ | 72 meV
$\Delta\epsilon_{24}$ | 11 meV
$\Delta\epsilon_{34}$ | 29 meV

Table S2:  Parameters for the redox-coupled proton pump model used in earlier
studies [18, 17, 21]. $E_{i}^{*}$ is the free energy associated with site $i$
when site $i$ is occupied while all other sites are unoccupied and the
membrane potential $V_{m}$ is zero. $\epsilon_{ij}$ is the electrostatic
interaction between sites $i$ and $j$. $\kappa_{ii}$ is the rate constant for
a single proton or electron transfer from a reservoir to site $i$ when all
other sites are unoccupied. $V_{m}$ is the membrane potential. $q_{i}$ is the
elementary charge of site $i$. $Z_{i}$ is the distance from the N-side to site
$i$. $L$ is the membrane thickness.

$E_{1}^{*}$ |  $98.4$ meV
---|---
$E_{2}^{*}$ |  $229$ meV
$E_{3}^{*}$ |  $386$ meV
$\epsilon_{12}$ |  $319$ meV
$\epsilon_{13}$ |  $-386$ meV
$\epsilon_{23}$ |  $-578$ meV
$\kappa_{11}$ |  $7.58\times 10^{6}$ sec-1
$\kappa_{22}$ |  $2.84\times 10^{4}$ sec-1
$\kappa_{33}$ |  $2.77\times 10^{4}$ sec-1

$V_{m}$ | 100
---|---
$q_{1}$ |  1 $|e^{-}|$
$q_{2}$ |  1 $|e^{-}|$
$q_{3}$ |  -1 $|e^{-}|$
$Z_{1}$ | 15 Å
$Z_{2}$ | 25 Å
$Z_{3}$ | 20 Å
$L$ | 30 Å

Figure S1: Hill’s model for free energy transduction. (a) Schematic
representation of the model. Two substrates, L and M, have different
concentrations (chemical potentials) on the two sides of a membrane. The
concentration of molecule M inside ($c_{M_{i}}$) is greater than outside
($c_{M_{o}}$), while molecule L has greater concentration outside [7]. An
energy transducing enzyme, E, with binding sites for L and M, is bound to the
membrane. The enzyme pumps L from the inside to the outside against its
concentration gradient (“L-cycle”: uphill transport of L), using the energy
released by moving M from inside to outside down its concentration gradient
(“M-cycle”: downhill transport of M) [7]. E may exist in two conformations (E
and E*) [7]. In conformation E (E*), the binding sites are accessible only for
L and M that are inside (outside) [7]. The thicker arrows show pumping of L
against a concentration gradient and transport of M down the concentration
gradient. (b) Hill’s model mapped onto particle occupancy microstates. The
microstates are ($\sigma_{L}$, $\sigma_{E}$, $\sigma_{M}$). $\sigma_{L}$ and
$\sigma_{M}$ indicate the site occupancies of the L and M binding sites.
Different enzyme conformations are denoted by the fictitious E site, using
$\sigma_{E}=0$ and E*:$\sigma_{E}=1$. The schematics adjacent to the
microstates correspond to the binding sites and enzyme shown in panel (a) and
illustrate occupied sites and conformations of E for each microstate.
$\text{M}_{i}$ ($\text{L}_{i}$) and $\text{M}_{o}$ ($\text{L}_{o}$) represent
reservoirs of molecule M (L) that are inside and outside, respectively. The
red and orange arrows represent the free energy transduction cycle, which
transports one M and one L species [7]. The green dashed arrow represents
slippage, namely the transport of M from the inside to the outside without
coupling to L transport, and the blue dashed arrow represents slippage
associated with the flow of L from the outside to the inside [7].

Figure S2: Kinetics of a linear transport chain inspired by the STC protein.
The plot shows the steady-state net electron fluxes into reservoir IV as a
function of $\Delta G$ varied from -300 to 300 meV, computed using the mean-
field approximation (black solid line), the approach of Jiang $et$ $al.$ (blue
dashed line), and exact master equation (red dashed line). Panel (a) shows the
case with site-site interactions and panel (b) shows the case without site-
site interactions (i.e., $\Delta\epsilon_{ij}$ = 0, where $i$ and $j$ are
labels of the redox cofactors shown in main text Fig. 1). The mean-field
approximation, the approach of Jiang et al. and the exact master equation all
agree (in both regimes) within an order of magnitude. In (b), the mean-field
limit approximates the kinetics without site-site interactions (error $<1\%$)
very well. Without site-site interactions, the mean-field approximation and
the approach of Jiang et al. are equivalent.

Figure S3: Steady-state flux of protons pumped to the P-side
($J_{pump}(s^{-1})$), calculated using (a) Jiang $et$ $al.$’s approach
including site-site interactions, (b) the exact master equation without site-
site interactions, (c) the mean-field model without site-site interactions,
and (d) Jiang $et$ $al.$’s approach without site-site interactions. Positive
$J_{pump}$ values indicate energy transduction since $J_{pump}>0$ indicates
spontaneous proton flow uphill from the N- to the P-side of the membrane
(against the the membrane potential or proton motive force $V_{m}=$ 100 mV).
Negative $J_{pump}$ indicates that slippage reactions dominate, causing proton
leakage from the P-side to the N-side. Comparing panel (a) here with panel (a)
of Fig. 4 in the main text, we find that in regimes ($k_{p}\gtrsim 10^{8}$
sec-1) where the exact master equation predicts no free energy transduction,
the approach of Jiang $et$ $al.$ also predicts no free energy transduction.
However, in regimes with $k_{p}\lesssim 10^{8}$ sec-1, the exact master
equation predicts free energy transduction, and the approach of Jiang $et$
$al.$ does not predict free energy transduction. Thus, the approach of Jiang
$et$ $al.$ is not able to describe free energy transduction for this redox-
coupled proton pump model, similar to the failure of the mean-field kinetics
to describe the energy transduction. We observe that site-site interactions
are required for free energy transduction in this proton-pump model,
consistent with conclusions from earlier studies [17, 18]. Panels (c) and (d)
are exactly the same, as expected because the mean-field kinetics and the
approach of Jiang et al. [8] are equivalent in the absence of site-site
interactions.

Figure S4: Kinetic simulations of the electron bifurcation network in regimes
not explored in Fig. 6 of the main text. The solid curves indicate the
electron fluxes computed with the mean-field approximation, and the dashed
curves are the fluxes computed using the exact master equation. The sum of the
values represented by the blue and red curves matches the values indicated by
the purple curve, consistent with the conservation of electrons at steady-
state. The efficiency of electron bifurcation is indicated by the gap between
the blue and red curves, which arises from short-circuiting. Panel (c) shows
the case when the high and low-potential reservoirs have different chemical
potentials ($\Delta G_{slope}$=250 meV for the low-potential branch and
$\Delta G_{slope}$=100 meV for the high-potential branch), but have lower
chemical potentials than the two-electron reservoir ($\Delta G_{HP}<0,\Delta
G_{LP}<0$). Therefore, case (c) does not produce energy transduction because
energy transduction requires an uphill and downhill process, or $\Delta
G_{HP}<0$ and $\Delta G_{LP}>0$. A significant discrepancy remains between the
blue curves, indicating that the mean-field approximation of the kinetics may
fail in systems that do not transduce energy.

Figure S5: Temperature-dependent short-circuiting fluxes at zero driving force
$\Delta G_{bifurc}=0$ meV plotted on a log scale. The exact master equation
calculated flux (black), mean-field approximated flux (purple), $J1$
calculated using Eq. 35 (orange), $J2$ calculated using Eq. 36 (green), and
$J3$ calculated using Eq. 37 (pink). These data were computed for the network
with ($\Delta G_{slope},\Delta G_{HL}$) = (300 meV, 200 meV), as shown in Fig.
5 of the main text. The slope of the purple line that represents the mean-
field calculated flux is approximately the same as the slope of the line
representing $J3$. Since the plot is on a log scale, this indicates that the
mean-field calculated flux has an effective activation energy that is the
average of the activation energies of the short circuiting channels $J1$ and
$J2$, as described in Section 6 of the SI appendix.

Figure S6: Illustration of electron transfer steps associated with three
different short-circuiting channels: (a) $B^{-}\rightarrow H1$ (orange arrow),
(b) $L1^{-}\rightarrow B^{-}$ (orange arrow) and (c) $L1^{-}\rightarrow H1$
(green arrow) [14, 19]. When $\Delta G_{bifurc}=0$, there is no energy cost to
refill the two-electron cofactor $B$. $k_{rl}$ labels the rate-limiting
electron-transfer rate constant in each short-circuiting channel. $\Delta
G^{e^{-}}_{SC}$ ($\Delta G^{h^{+}}_{SC}$) is the free energy to bring an
electron (hole) to the site required for the short-circuiting channel shown.
(a) An electron is transferred from the low-potential reservoir with energy
cost $\Delta G_{B^{-}\rightarrow H1}^{e^{-}}$, and a hole is transferred from
the high-potential reservoir to $H_{1}$ with energy $\Delta
G_{B^{-}\rightarrow H1}^{h^{+}}$, followed by rate-limiting electron transfer
from $L_{1}^{-}$ to $B$. (b) A hole is transferred from the high-potential
reservoir to $B^{=}$, forming $B^{-}$, costing energy $\Delta
G_{L1^{-}\rightarrow B^{-}}^{h^{+}}$. Next, an electron is transferred from
the low-potential reservoir to $L_{1}$ with energy cost $\Delta
G_{L1^{-}\rightarrow B^{-}}^{e^{-}}$, followed by a rate-limiting electron
transfer to $B^{-}$ from $L_{1}^{-}$. (c) An electron (hole) is transferred
from the low-potential (high-potential) reservoir with $\Delta
G_{L1^{-}\rightarrow H1}^{e^{-}}$ ($\Delta G_{L1^{-}\rightarrow H1}^{h^{+}}$)
to $L1$ ($H1$) followed by the short-circuit electron transfer from
$L_{1}^{-}$ to $H_{1}$ that is the rate-limiting electron transfer for the
short-circuit channel shown in (c).

## References

  * [1] L. Bányai and A. Aldea “Master Equation Approach to the Hopping Transport Theory” In _Fortschritte der Phys._ 27.9, 1979, pp. 435–462 DOI: 10.1002/prop.19790270903
  * [2] B. Derrida, E. Domany and D. Mukamel “An exact solution of a one-dimensional asymmetric exclusion model with open boundaries” In _J. Stat. Phys._ 69.3, 1992, pp. 667–687 DOI: 10.1007/BF01050430
  * [3] A.. Shumilin and Y.. Beltukov “Effect of Nonequilibrium Correlations on the Effective Resistances between Sites in the Theory of Hopping Conductivity” In _Phys. Solid State_ 61.11, 2019, pp. 2090–2095 DOI: 10.1134/S1063783419110349
  * [4] David Chandler “Introduction to modern statistical mechanics” New York: Oxford University Press, 1987
  * [5] Peter M Richards “Theory of one-dimensional hopping conductivity and diffusion” In _Phys. Rev. B_ 16.4 APS, 1977, pp. 1393
  * [6] Vinay Ambegaokar, BI Halperin and JS Langer “Hopping conductivity in disordered systems” In _Phys. Rev. B_ 4.8 APS, 1971, pp. 2612
  * [7] Terrell L Hill “Free energy transduction and biochemical cycle kinetics” New York: Springer, 2013
  * [8] Xiuyun Jiang et al. “Cysteine Linkages Accelerate Electron Flow through Tetra-Heme Protein STC” In _J. Am. Chem. Soc._ 139.48, 2017, pp. 17237–17240 DOI: 10.1021/jacs.7b08831
  * [9] JJ Hopfield “Electron transfer between biological molecules by thermally activated tunneling” In _Proc. Natl. Acad. Sci. U.S.A._ 71.9, 1974, pp. 3640–3644 URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC433831/pdf/pnas00072-0326.pdf
  * [10] R.. Marcus and Norman Sutin “Electron transfers in chemistry and biology” In _Biochim. Biophys. Acta - Bioenerg._ 811.3, 1985, pp. 265–322 DOI: https://doi.org/10.1016/0304-4173(85)90014-X
  * [11] Harry B. Gray and Jay R. Winkler “Electron tunneling through proteins” In _Q. Rev. Biophys._ 36.3, 2003, pp. 341–372 DOI: 10.1017/S0033583503003913
  * [12] David N Beratan et al. “Charge Transfer in Dynamical Biosystems, or The Treachery of (Static) Images” In _Acc. Chem. Res._ 48.2, 2015, pp. 474–481 DOI: 10.1021/ar500271d
  * [13] Harry B. Gray and Jay R. Winkler “Long-range electron transfer” In _Proc. Natl. Acad. Sci. U.S.A._ 102.10, 2005, pp. 3534–3539 DOI: 10.1073/pnas.0408029102
  * [14] J.. Yuly et al. “Universal free-energy landscape produces efficient and reversible electron bifurcation” In _Proc. Natl. Acad. Sci. U.S.A._ 117.35, 2020, pp. 21045 DOI: 10.1073/pnas.2010815117
  * [15] Michel Bauer and Françoise Cornu “Local detailed balance: a microscopic derivation” In _J. Phys. A Math. Theor._ 48.1 IOP Publishing, 2014, pp. 015008
  * [16] Marian Breuer, Piotr Zarzycki, Jochen Blumberger and Kevin M. Rosso “Thermodynamics of Electron Flow in the Bacterial Deca-heme Cytochrome MtrF” In _J. Am. Chem. Soc._ 134.24, 2012, pp. 9868–9871 DOI: 10.1021/ja3027696
  * [17] Young C Kim, Mårten Wikström and Gerhard Hummer “Kinetic models of redox-coupled proton pumping” In _Proc. Natl. Acad. Sci. U.S.A._ 104.7 National Acad Sciences, 2007, pp. 2169–2174
  * [18] Young C Kim and Gerhard Hummer “Proton-pumping mechanism of cytochrome _c_ oxidase: A kinetic master-equation approach” In _Biochim. Biophys. Acta - Bioenerg._ 1817.4 Elsevier, 2012, pp. 526–536
  * [19] Jonathon L. Yuly, Peng Zhang and David N. Beratan “Energy transduction by reversible electron bifurcation” In _Curr. Opin. Electrochem._ 29, 2021, pp. 100767 DOI: https://doi.org/10.1016/j.coelec.2021.100767
  * [20] Bruno M Fonseca et al. “The tetraheme cytochrome from _Shewanella oneidensis_ MR-1 shows thermodynamic bias for functional specificity of the hemes” In _J. Biol. Inorg. Chem._ 14.3 Springer, 2009, pp. 375–385
  * [21] Young C Kim, Mårten Wikström and Gerhard Hummer “Kinetic gating of the proton pump in cytochrome _c_ oxidase” In _Proc. Natl. Acad. Sci. U.S.A._ 106.33 National Acad Sciences, 2009, pp. 13707–13712