Optimization of nonequilibrium free energy harvesting illustrated on bacteriorhodopsin

Many molecular systems, both biological and artificial, harvest free energy from their environments. Biological organisms require free energy to grow and replicate [1, 2], and many species undergo selection for increased harvesting [3, 4, 5, 6]. Artificial harvesting systems have also been constructed and optimized in the field of synthetic biology [7, 8, 9, 10, 11, 12, 13, 14]. The optimization of free energy harvesting is thus a central problem both in biology and engineering.

As an example, consider a harvesting system such as a biological metabolic network that converts glucose to ATP [15]. Suppose that the kinetic and thermodynamic parameters of one or more reactions can be optimized, either by natural selection or artificial design. What is the maximum rate of free energy harvesting that the network can achieve, and what are the kinetic and thermodynamic parameters that achieve it? These questions are relevant both for design of optimal harvesting devices and for quantifying the efficiency of existing systems.

In this paper, we use techniques from stochastic thermodynamics to derive bounds on maximum rate of free energy harvesting. We consider a harvesting system in nonequilibrium steady state which is coupled to an external source of free energy, an internal free energy reservoir, and a heat bath. The setup is well-suited for studying the kinds of small-scale systems usually considered in stochastic thermodynamics [16], where assumptions of local detailed balance and steady state are justified. The steady-state assumption is reasonable in many molecular systems, where there is a separation of timescales between internal relaxation and environmental change.

We suppose that the system’s dynamics can be separate into two kinds of processes, termed baseline and control. The baseline processes, which are held fixed, mediate the coupling to the external source of free energy. Control refers to all other processes which can be optimized for increasing the harvesting rate at which free energy flows into the internal reservoir. The particular separation of baseline/control generally depends on domain knowledge about the system and the scientific question at hand. For example, to study the efficiency of a particular reaction in a metabolic network, that reaction may be treated as control while the other reactions are baseline. The baseline/control separation is similar to the distinction in control theory between “plant” (a given system with fixed dynamical laws) and “controller” (the part that undergoes optimization) [17].

In our first set of results, we demonstrate that the optimization of the harvesting rate can be expressed as a convex optimization problem. The optimal solution of this problem determines both the maximum harvesting rate and the specific control processes that achieve that maximum. Importantly, the optimization can also account for various types of constraints on the possible network topology, kinetic timescales, and thermodynamic forces of the control processes.

We illustrate our results on bacteriorhodopsin (Fig. 1), a proton-pumping membrane protein. Bacteriorhodopsin is a photosynthetic system found in Archaea, with close relatives in bacteria [20, 21]. It is also used in many artificial light-harvesting systems [7, 9, 8, 11]. Using published thermodynamic and kinetic data, we develop a thermodynamically consistent stochastic model of bacteriorhodopsin. We demonstrate that, under normal operating conditions, the bacteriorhodopsin system is remarkably efficient.

Our main result is formulated as a convex optimization problem which must be solved numerically in the general case. In the second part of this paper, we derive closed-form solutions of this problem for three physically meaningful regimes: the weakly-driven linear response regime, the irreversible deterministic regime, and the intermediate near-deterministic regime. These solutions illustrate how optimal harvesting reflects the “alignment” between free energy input and relaxation dynamics. We illustrate these closed-form solutions on two unicyclic systems, which may be interpreted as idealized models of two types of nonequilibrium harvesting devices.

We finish our paper with a brief Discussion. There we relate our approach to previous work, including maximization of power output in steady-state engines and flux balance analysis. We also propose directions for future research.

We consider a system with $n$ mesostates described by the distribution $\bm{p}=(p_{1},\dots,p_{n})\in\mathbb{R}_{+}^{n}$. The distribution evolves according to the master equation $\dot{p}_{i}=\sum_{j}R_{ij}p_{j}$, where $R_{ji}$ is the transition rate $i\to j$ ($R_{ii}=-\sum_{j}R_{ji}$). Usually $\bm{p}$ represents a probability distribution of a stochastic system with Markovian dynamics [22, 23]. However, under an appropriate choice of units, it may also represent chemical concentrations in a deterministic chemical reaction network with pseudo-first-order reactions, such as an enzymatic cycle [24, 25].

The system is coupled to a heat bath at inverse temperature $\beta=1/k_{B}T$, an internal free energy reservoir, and a nonequilibrium environment that serves as an external source of free energy. For example, in a metabolic network, one may consider an internal reservoir of ATP and an external source of glucose. The system has nonequilibrium free energy

where $S(\bm{p}):=-\sum_{i}p_{i}\ln p_{i}$ is the Shannon entropy and $f_{i}$ is the internal free energy of mesostate $i$ [26].

As mentioned in the Introduction, we suppose that the dynamics of the system can be separated into baseline and control processes. We make one important assumption in our analysis: the control processes are only coupled to the heat bath and internal free energy reservoir, but not directly to the external source of free energy. This means that control can only increase harvesting by interacting with the baseline, rather than directly increasing the inflow of free energy from the external source. For example, in a metabolic network, control processes cannot directly increase the import of glucose, but they can affect the rate at which glucose is converted into ATP by optimizing intermediate reactions and mechanisms. Control processes may be driven by the internal reservoir (e.g., coupled to ATP hydrolysis) or undriven (e.g., enzymes).

To formalize the baseline/control distinction, we write the rate matrix as $R=R^{b}+R^{c}$, where $R^{b}_{ji}$ and $R^{c}_{ji}$ represent the transition rate of $i\to j$ due to baseline and control. Given distribution $\bm{p}$, the increase of system free energy due to baseline processes is

The increase due to control processes is defined analogously but using rate matrix $R^{c}$,

The harvesting rate is the rate at which free energy flows to the internal reservoir. The harvesting rate due to baseline processes is

where ${g}^{b}_{ji}$ is the free energy increase in the internal reservoir due to a single baseline transition $i\to j$ and $\dot{g}^{b}_{i}$ is the rate of free energy flow to the internal reservoir due to internal transitions within $i$ (assuming $i$ is a coarse-grained mesostate). Similarly, the harvesting rate due to control processes is

where ${g}^{c}_{ji}$ is the free energy increase in the internal reservoir due to control transition $i\to j$. For simplicity, we assume that control cannot exchange free energy with the internal reservoir due to internal transitions within $i$. Negative values of ${g}^{b}_{ji},\dot{g}^{b}_{i},{g}^{c}_{ji}$ indicate driving done on the system by the internal reservoir.

For a concrete example of how $(R^{b},R^{c},\bm{f},\bm{g}^{b},\dot{\bm{g}}^{b},\bm{g}^{c})$ are defined for a real biomolecular system, see the bacteriorhodopsin example below and SM-II [27].

As standard in stochastic thermodynamics, we assume that control processes obey local detailed balance (LDB) [16, 28],

Eq. (6) guarantees that the irreversibility of each control transition is determined by the amount of free energy dissipated by that transition [29]. Observe that the right side accounts for free energy changes of the system ($f_{i}-f_{j}$) and the internal reservoir (${g}^{c}_{ji}$), but not the external source. This formalizes the assumption that control processes do not exchange free energy with the external source.

We do not require that the baseline rate matrix obeys LDB, although it will often do so for reasons of thermodynamic consistency.

Suppose that the combined rate matrix $R=R^{b}+R^{c}$ has the steady-state distribution $\bm{\pi}$, which satisfies $R^{b}\bm{\pi}+R^{c}\bm{\pi}=0$. The total steady-state harvesting rate due to baseline and control is

We seek to maximize this harvesting rate by varying the parameters of the control processes $(R^{c},\bm{g}^{c})$ while holding the baseline parameters ($\bm{f},R^{b},\bm{g}^{b},\dot{\bm{g}}^{b}$) fixed. Finding this maximum would allow us to determine fundamental bounds on harvesting and to evaluate the efficiency of existing harvesting systems.

However, $\dot{\mathcal{G}}^{\mathrm{tot}}$ is not a concave function of the parameters $(R^{c},\bm{g}^{c})$, therefore maximization of (7) is not a convex optimization problem and is not generally intractable. In the following, we reformulate this maximization as a convex optimization with a physically interpretable objective. This allows us to solve the optimization numerically and, at least for some special cases, also in closed form.

To begin, we rewrite (7) as

where we introduced the Schnackenberg formula for the entropy production rate (EPR) [22],

where ${J^{c}_{ji}}=\pi_{i}R^{c}_{ji}\geq 0$ is the one-way probability flux due to control transition $i\to j$.

Eq. (8) has an intuitive physical interpretation: the total steady-state harvesting rate is the rate of free energy increase in the system and internal reservoir due to baseline, minus the rate of dissipation (EPR) due to the control fluxes. The derivation of this result proceeds in two steps (see SM-I.1 [27] for details). The first step is to show that $\dot{\Sigma}({\boldsymbol{J}}^{c})=-\beta[\dot{\mathcal{F}}^{c}(\bm{\pi})+\beta\dot{\mathcal{G}}^{c}(\bm{\pi})]$, which follows by combining (9) with (3) and (6). This states that the EPR due to control is proportional to the free energy loss in the system and internal reservoir due to control. The second step is to show that $\dot{\mathcal{F}}^{b}(\bm{\pi})+\dot{\mathcal{F}}^{c}(\bm{\pi})=0$, which follows because the left side is the overall derivative of the nonequilibrium free energy $\mathcal{F}$, as defined in (1), therefore it must vanish in steady state. The result (8) then follows by combining with (7) and rearranging.

Importantly, when expressed in the form (8), the harvesting rate is a concave function of the steady-state distribution $\bm{\pi}$ and the control fluxes ${\boldsymbol{J}}^{c}$ (see SM-I.2 [27]). To find the maximum harvesting rate, we optimize (8) with respect to $\bm{\pi}$ and ${\boldsymbol{J}}^{c}$. Note that varying $\bm{\pi}$ and ${\boldsymbol{J}}^{c}$ is equivalent to varying the control rate matrix via $R^{c}_{ji}={J^{c}_{ji}}/\pi_{i}$ and control driving ${g}^{c}_{ji}$ via (6). However, when performing this optimization, we must also ensure that $\bm{\pi}$ is the steady-state distribution induced by the fluxes ${\boldsymbol{J}}^{c}$. This condition can be expressed as a linear constraint on $\bm{\pi}$ and ${\boldsymbol{J}}^{c}$ via $R^{b}\bm{\pi}+\mathbb{B}{\boldsymbol{J}}^{c}=0$. Here ${\boldsymbol{J}}^{c}$ is treated as a vector in $\mathbb{R}^{n^{2}}$ and $\mathbb{B}\in\mathbb{R}^{n\times n^{2}}$ is the incidence matrix with entries $\mathbb{B}_{k,ij}=\delta_{ki}-\delta_{kj}$, which guarantees $R^{c}\bm{\pi}=\mathbb{B}{\boldsymbol{J}}^{c}$.

Combining, we arrive at the bound $\dot{\mathcal{G}}^{\mathrm{tot}}\leq\mathscr{G}$, where

In this expression, $\Lambda$ is the feasible set of distributions $\bm{p}$ and control fluxes ${\boldsymbol{J}}$. At a minimum, $\Lambda$ ensures the validity of the distribution $\bm{p}$ and the fluxes ${\boldsymbol{J}}$ via the linear constraints $\sum p_{i}=1$, $p_{i}\geq 0$, and $J_{ji}\geq 0$. We write $\sup$ instead of $\max$ because the set of allowed fluxes is potentially unbounded. Eq. (10) implies a tradeoff between the total gain of free energy in the system and internal reservoir due to baseline (which depends only on $\bm{p}$) and the dissipation due to control fluxes (which depends only on ${\boldsymbol{J}}$).

Importantly, the feasible set $\Lambda$ can include additional convex constraints, which may introduce topological, kinetic, thermodynamic, etc. restrictions on the control processes. Topological constraints restrict which transitions can be controlled; e.g., $J_{ji}=0$ ensures that control does not use transition $i\to j$). Kinetic constraints restrict timescales of control processes, as might reflect underlying physical processes like diffusion; e.g., an upper bound on control transition rate $R^{c}_{ji}=J_{ji}/p_{i}\leq\kappa$ can be enforced via $J_{ji}\leq p_{i}\kappa$. Thermodynamic constraints bound the affinity [22] of control transitions; e.g., $J_{ji}e^{-\mathcal{A}}\leq J_{ij}\leq J_{ji}e^{\mathcal{A}}$ ensures that $|\ln(J_{ij}/J_{ji})|\leq\mathcal{A}$. The above examples all involve linear constraints. An example of a nonlinear, but still convex, constraint is an upper bound on the EPR incurred by control, $\dot{\Sigma}({\boldsymbol{J}})\leq\dot{\Sigma}^{c}_{\max}$.

Eq. (10) is our first main result. Importantly, $\mathscr{G}$ is defined purely in terms of the thermodynamic and kinetic properties of the baseline process, along with desired constraints encoded in $\Lambda$. Thus, $\mathscr{G}$ is the maximum steady-state harvesting rate that can be achieved by any allowed control processes, given a fixed baseline. In addition, Eq. (10) involves the maximization of a concave objective given convex constraints. This is equivalent to the minimization of a convex objective, thus Eq. (10) is a convex optimization problem that can be efficiently solved using standard numerical techniques [30]. The optimization also identifies an optimal steady-state distribution $\bm{p}^{*}$ and control fluxes ${{\boldsymbol{J}}}^{*}$ that achieve the maximum harvesting rate $\mathscr{G}$ (or come arbitrarily close to achieving it). These fix the optimal control rate matrix via $R^{c*}_{ji}={J_{ij}^{*}}/p^{*}_{i}$. Thus, our optimization specifies an upper bound on harvesting as well as the optimal control strategy that achieves this bound.

There is an important special case in which our optimization problem is simplified. Suppose that $\Lambda$ does not enforce additional constraints on $\bm{p}$ and ${\boldsymbol{J}}$ (more generally, we permit topological constraints if the graph of allowed transitions is connected and contains all $n$ states). Then, the objective is maximized in limit of fast control, ${\boldsymbol{J}}\to\infty$ and $\dot{\Sigma}({\boldsymbol{J}})\to 0$. We can then write (10) as an optimization over steady-state distributions:

The optimal $\bm{p}^{*}$ is unique as long as the baseline rate matrix is irreducible. The optimal control rate matrix is very fast ($R^{c*}\to\infty$) and obeys detailed balance for $\bm{p}^{*}$, $R^{c*}_{ji}{p}^{*}_{i}=R^{c*}_{ij}{p}^{*}_{j}$. For details, see SM-I.3 and SM-I.4 [27].

We illustrate our results using bacteriorhodopsin, a light-driven proton pump from Archaea [18].

We define a thermodynamically consistent model of the bacteriorhodopsin cycle using published thermodynamic [31] and kinetic [32] data (see SM-II [27]). The system operates in a cyclical manner, absorbing a photon and pumping a proton during each turn of the cycle (Fig. 1, right). Specifically, the transition $M_{1}\to M_{2}$ pumps a proton out of the cell. This stores free energy in the internal reservoir (the membrane electrochemical potential),

where $e$ is the elementary charge constant, $\Delta\psi$ is the membrane electrical potential, and $\Delta\text{pH}$ is the membrane pH difference. The other transitions in the cycle do not affect the free energy of the internal reservoir ($g_{ij}=0$ and $\dot{g}_{i}=0$).

During the transition $\mathrm{bR}\to K$, the system leaves the ground state by absorbing a photon at 580nm, thereby harvesting free energy from the external source and dissipating some heat to the environment at $T=293^{\circ}$ K. This transition is much faster (picoseconds) than the other transitions in the photocycle (micro- to milliseconds). As commonly done in photochemistry [33], we coarse-grain transitions $O\to\mathrm{bR}$ and $\mathrm{bR}\to K$ into a single effective transition $O\to K$.

To explore the performance of bacteriorhodopsin under different conditions, we vary the membrane electrical potential $\Delta\psi$ between $-75$ and $350$ mV, while using a realistic fixed $\Delta\text{pH}=-0.6$ [34]. We show the actual harvesting rate ($\dot{\mathcal{G}}^{\mathrm{tot}}$ in units of $k_{B}T/$sec) at different potentials as a black line in Fig. 2 (a). At a plausible in vivo $\Delta\psi=120$ mV [34], the model exhibits a steady-state current of 11.5 protons/sec, each proton carrying 6.1 $k_{B}T$ of free energy, corresponding to ${\dot{\mathcal{G}}^{\mathrm{tot}}}\approx 70\,k_{B}T$/sec. The largest output is achieved near the in vivo potential: at lower potentials, the cycle current saturates while the free energy per proton drops, and at higher potentials the pump stalls.

Next, we quantify the maximum harvesting rate that can be achieved by optimizing the parameters of individual transitions. This analysis is relevant for understanding limits on increasing bacteriorhodopsin output, whether via natural selection or biosynthetic methods [35, 36, 37, 38]. Interestingly, such transition-level optimization may be feasible in bacteriorhodopsin, as the properties of several transitions are known to be individually controlled by particular amino acid residues in the bacteriorhodopsin protein [35, 39, 40, 41].

For each reversible transition in the cycle, for instance $N\leftrightarrow O$, we define the baseline as the rest of the cycle without that transition. We then optimize control under the topological constraint that only the relevant transition (e.g., $N\leftrightarrow O$) is allowed. The analysis is repeated for all transitions, except for the (coarse-grained) photon-absorbing transition $O\leftrightarrow K$, which is in accordance with our assumption that control cannot directly exchange free energy directly with the external source.

Fig. 2 (a) shows $\mathscr{G}$, the maximum $\dot{\mathcal{G}}^{\mathrm{tot}}$ achievable by optimizing each reversible transition. In Fig. 2 (b), we also show the efficiency $\dot{\mathcal{G}}^{\mathrm{tot}}/\mathscr{G}\leq 1$ for each transition, that is the ratio of the actual and maximum harvesting rate.

Several transitions, such as $K\leftrightarrow L$,$L\leftrightarrow M_{1}$,$M_{1}\leftrightarrow M_{2}$, are remarkably efficient ($\geq 85\%$) near in vivo membrane potentials. The reprotonation step $N\leftrightarrow O$ is the least efficient ($\sim 40\%$) and also has the slowest kinetics of the 5 transitions studied in Fig. 2. This suggests that $N\leftrightarrow O$ is a bottleneck whose optimization can have a big impact on the harvesting rate, while optimization of other non-bottleneck transitions has a more limited effect.

Observe that $\mathscr{G}$ for $M_{1}\leftrightarrow M_{2}$ does not depend on $\Delta\psi$. This is because $\mathscr{G}$ is a function of baseline properties, which do not depend on the membrane potential when $M_{1}\leftrightarrow M_{2}$ is treated as control. Conversely, $M_{1}\leftrightarrow M_{2}$ as control transition can be optimized by varying the membrane potential and/or scaling up the forward/backward rates. Our results show that this transition is very close to optimal at in vivo membrane potentials and kinetic timescales.

Optimal distributions $\bm{p}^{*}$ are also obtained, with a typical one shown in Fig. 2 (c). We find a consistent shift toward state $O$, which accelerates the reset of the cycle and increases the flux across the photon-absorbing transition $O\to K$.

In SM-II [27], we illustrate how the efficiency of bacteriorhodopsin transitions can be evaluated under other types of constraints, including constraints on thermodynamic affinity, dynamical activity, and kinetics.

As a final analysis, instead of optimizing individual existing transitions in the bacteriorhodopsin cycle, we ask to what extent the harvesting rate can be increased by any additional control processes. For example, this could involve an additional enzyme that shifts the cycle’s steady state by permitting a new transition between distant states (e.g. $L\leftrightarrow N$), possibly yielding an increase of the proton pumping rate.

In this case, we treat the entire bacteriorhodopsin system as the baseline, and we do not introduce any additional constraints on the steady-state distribution or the control fluxes. Then, as shown in Sec. SM-I.4 [27], the objective in (10) is achieved in the limit of fast control, and the maximum harvesting rate can be found by solving the simplified optimization problem (11).

For this setup, Fig. 3 (a) shows the baseline (actual) harvesting rate $\dot{\mathcal{G}}^{\mathrm{tot}}$ and the maximum harvesting rate $\mathscr{G}$ at varying $\Delta\psi$. Interestingly, both peak at around the in vivo values of the membrane potential. In Fig. 3 (b), we show that the actual bacteriorhodopsin cycle harvests approximately 50% of the fundamental bound given by $\mathscr{G}$ (at in vivo values of $\Delta\psi$). This suggests that bacteriorhodopsin is remarkably close to optimal, relative to improvements that could be achieved by introducing any additional control processes.

We also show the actual steady-state distribution and the optimal distribution $\bm{p}^{*}$ in Fig. 3 (c). The optimal distribution increases the probability of state $O$, similar to the optimal distribution found by optimizing the $N\leftrightarrow O$ transition, shown in Fig. 2 (c). However, unlike Fig. 2 (c), where most of the extra probability is taken from state $N$, in Fig. 3 (c) the probability is drawn more uniformly from other states in the cycle, indicating the presence of distributed control.

Our results are stated via an optimization problem that generally does not have a closed-form solution. In our second set of results, we identify closed-form expressions in three physically meaningful regimes. For simplicity, here we focus on the simplified objective (11). See SM-III [27] for detailed derivations, including analysis of the conditions under which each of these three approximation are be valid.

For convenience, we first rewrite (11) as

where $\dot{S}^{b}(\bm{p})=-\sum_{i,j}R^{b}_{ij}p_{j}\ln p_{i}$ is the increase of the Shannon entropy of $\bm{p}$ under $R^{b}$ and for convenience we defined $\phi_{i}:=\dot{g}^{b}_{i}+\sum_{j}R^{b}_{ji}(f_{j}-f_{i}+{g}^{b}_{ji})$. The objective (13) contains a nonlinear term $-\dot{S}^{b}(\bm{p})/\beta$ quantifying the decrease of information-theoretic entropy and a linear term $\sum_{i}p_{i}\phi_{i}$ quantifying the flow of thermodynamic free energy.

Next, we consider three regimes.

Linear response (LR) applies when the optimal distribution $\bm{p}^{*}$ is close to the steady-state distribution of the baseline rate matrix $R^{b}$. Suppose that $R^{b}$ is irreducible and has a unique steady state ${\bm{\pi}^{b}}$ with full support. We introduce the “additive reversibilization” of $R^{b}$,

The rate matrix $A$ obeys detailed balance (DB) for the steady-state distribution ${\bm{\pi}^{b}}$ and has the same dynamical activity [42] on all edges as $R^{b}$. $A$ may be considered as a DB version of $R^{b}$, and it is equal to $R^{b}$ when the latter obeys DB [43, 44]. Let $\boldsymbol{u}^{\alpha}$ indicate the $\alpha$-th right eigenvector of $A$ normalized as $\sum_{i}({u^{\alpha}_{i}})^{2}/\pi^{b}_{i}=1$, and $\lambda_{\alpha}$ the corresponding real-valued eigenvalue ($\lambda_{1}=0$). The LR solution for the maximum harvesting rate and the optimal distribution is

where $\Omega_{\alpha}=(\boldsymbol{\phi}+{\beta}^{-1}{R^{b}}^{T}\ln{\bm{\pi}^{b}})^{T}\boldsymbol{u}^{\alpha}/2$ quantifies the harvesting “amplitude” for mode $\alpha$.

Eq. (15) has a simple interpretation. In addition to the baseline harvesting rate $\dot{\mathcal{G}}^{b}({\bm{\pi}^{b}})$, $\mathscr{G}$ contains contributions from the relaxation modes of $A$, with each mode weighed by its squared amplitude and relaxation timescale $-1/\lambda_{\alpha}$. All else being equal, $\mathscr{G}$ is large when slow modes have large harvesting amplitudes. The optimal $\bm{p}^{*}$ shifts the baseline steady state ${\bm{\pi}^{b}}$ toward mode $\alpha$ in proportion to that mode’s harvesting amplitude and relaxation timescale, thereby optimally balancing the tradeoff between harvesting and dissipation.

The Deterministic (D) regime applies when the nonlinear information-theoretic term in (13) is much smaller than the linear thermodynamic term. We can then ignore the former, turning (13) into a simple linear optimization. This gives the approximate solution

where ${i}^{*}=\arg\max_{i}\phi_{i}$ is the optimal mesostate. This solution concentrates probability on a single mesostate, effectively ignoring the cost of maintaining this low-entropy distribution.

The Near-Deterministic (ND) regime lies between Linear Response and Deterministic ones. By perturbing $\bm{p}^{*}$ around $\delta_{{i}^{*}i}$, we can decouple the values of $p_{i}$ in the objective function (11). The maximal harvesting rate and optimal distribution in this regime are then given by

The ND solution also has a simple interpretation. It perturbs the Deterministic solution by shifting probability towards states with high transition rates away from the optimal state (large $R^{b}_{i{i}^{*}}$) and small decreases in harvesting ($\phi_{{i}^{*}}-\phi_{i}$). This balances the benefit of harvesting against the cost of pumping probability against $R^{b}_{i{i}^{*}}$.

We illustrate our closed-form solutions using two simple models, both based on a unicyclic system with $n$ states. The baseline dynamics involve diffusion across a 1-dimensional ring, with left and right jump rates set to 1. The baseline steady state is a uniform distribution, $\pi^{b}_{i}=1/n$, with no cyclic current. We assume a uniform free energy function, $f_{i}=0$ for all $i$.

We consider two different scenarios. In the first scenario, shown schematically in Fig. 4 (a), $\Theta$ of free energy is harvested each time the system carries out the transition $1\to 2$, so

and ${g}^{b}_{ij}=\dot{g}^{b}_{i}=0$ otherwise. This scenario may be interpreted as an idealized model of a biomolecular harvesting cycle, such as a transporter. In the second scenario, shown schematically in Fig. 4 (b), free energy is harvested at a rate of $\theta$ per unit time when the system is located in one particular mesostate ${i}^{*}=1$, so

and ${g}^{b}_{ij}=\dot{g}^{b}_{i}=0$ otherwise. This scenario may be interpreted as an idealized model of error correction or self-assembly, where free energy can only be harvested when the system is in some particular functional mesostate.

For both scenarios, we evaluate the maximum harvesting rate $\mathscr{G}$ and the optimal distributions achievable by adding any possible control to the system, without constraints. We report exact values found by numerical optimization of Eq. (13), as well as the LR, ND, and D approximations described above. To calculate the LR values, we exploit the fact that the baseline unicyclic rate matrix is a circulant matrix with a simple eigendecomposition [45]. Full details of the derivations for the two scenarios are provided in SM-IV [27] and SM-IV.3 [27], respectively.

We first report results for the first scenario from Fig. 4 (a), where free energy is harvested during the transition $2\to 1$. Observe that baseline harvesting rate vanishes, $\dot{\mathcal{G}}^{b}({\bm{\pi}^{b}})=0$, since harvesting free energy requires a cyclic current. In Fig. 5 (a), we plot the maximum harvesting rate $\mathscr{G}$ and its approximations as a function of the supplied free energy $\Theta$. For small $\Theta$, LR applies and the maximum harvesting rate is

The optimal distribution in the LR regime, shown in Fig. 5 (c), builds up in equal increments starting from $i={i}^{*}+1$ until the optimal state ${i}^{*}=1$, after which it drops sharply. For large $\Theta$, the D regime is relevant and the optimal distribution concentrates on the optimal state ${i}^{*}=1$, so

At intermediate $\Theta$, the ND regime applies, which gives

The optimal distribution in the ND regime, shown in Fig. 5 (b), allocates ${p}^{*}_{{i}^{*}-1}=1/(\beta\Theta-2)$, ${p}^{*}_{{i}^{*}+1}=1/(2\beta\Theta-2)$ and the rest to the optimal state ${p}^{*}_{{i}^{*}}$.

Next, we consider the second scenario from Fig. 4 (b), where free energy is harvested when the system is in the optimal mesostate ${i}^{*}=1$. Observe that the uniform baseline steady state assigns $1/n$ probability to the optimal state, thus in this scenario the baseline harvesting rate is $\dot{\mathcal{G}}^{b}({\bm{\pi}^{b}})=\theta/n$. To facilitate comparison with the first scenario, we focus on the increase of the maximum harvesting rate relative to baseline,

In Fig. 6 (a), we show $\Delta\mathscr{G}$ and its approximations as a function of the free energy supply rate $\theta$. For small $\theta$, LR applies and the maximum harvesting rate is

The optimal distribution in the LR regime, shown in Fig. 6 (c), is symmetric about the optimal state ${i}^{*}=1$. For large $\theta$, the D regime is relevant and the optimal distribution concentrates on the optimal state ${i}^{*}=1$, so

At intermediate $\theta$, the ND regime applies, giving

The optimal distribution in the ND regime, shown in Fig. 6 (b), allocates ${p}^{*}_{{i}^{*}-1}={p}^{*}_{{i}^{*}+1}=1/(\beta\theta-2)$ and the rest to ${p}^{*}_{{i}^{*}}$.

There are some similarities among the two harvesting scenarios. For both scenarios, in the LR regime, the increase of the harvesting rate scales quadratically in the supplied free energy ($\Theta$ or $\theta$) and linearly in inverse temperature $\beta$. This scaling reflects the fact that the optimal strategy has to balance harvesting ($\Theta$ or $\theta$ contributions) with the thermodynamic cost of maintaining a low entropy $\bm{p}^{*}$ ($\beta$ contributions). In the ND and D regimes, $\mathscr{G}$ scales linearly in the supplied free energy and loses its linear dependence on $\beta$. Thus, outside of LR, the thermodynamic cost of maintaining a low entropy distribution has a minor effect on the optimal strategy.

There are also important differences between the two scenarios. For the first scenario, the optimal strategy maintains an asymmetric $\bm{p}^{*}$, thereby generating a net flux across the transition $2\to 1$. In the LR regime, the cost of maintaining this asymmetric distribution grows with the system size $n$, therefore the maximum harvesting rate in Eq. (19) scales as $\sim O(n^{-1})$. This is shown in Fig. 5 (d), where we plot $\mathscr{G}$ and its LR approximation at various $\Theta$ and $n$. For the second scenario, the optimal strategy maintains a peaked but symmetric $\bm{p}^{*}$. Remarkably, the cost of maintaining this distribution does not depend on system size $n$. This is shown in Fig. 6 (d), where we plot $\Delta\mathscr{G}$ at various $\theta$ and $n$.

In this paper, we consider the problem of optimizing free energy harvesting in a nonequilibrium steady-state system. We demonstrate that this problem can be formulated as a constrained convex optimization problem, and we use this formulation to study optimal harvesting and efficiency in the bacteriorhodopsin proton pump. We also solve the convex optimization problem in closed-form for three limiting regimes, as illustrated on two unicyclic models discussed above.

A key step in our analysis is to separate the dynamics of the system into separate contributions from fixed baseline processes and optimizable control processes. We note that, in stochastic thermodynamics, the baseline/control separation has been previously used to study autonomous Maxwellian demons [46, 47], counterdiabatic driving [48], and the cost of maintaining a nonequilibrium steady state [49, 50].

We derive a simplified bound on the maximum harvesting rate in (13), which is achieved in the limit of fast dissipation-less control. Interestingly, this expression involves a tradeoff between two terms, one information-theoretic and one thermodynamic. At first glance, this resembles information/free-energy tradeoffs characteristic of Maxwellian demons and other “information engines” [51, 52, 53, 54, 55, 56, 57, 58]. However, there are important differences. In a typical information engine, there is no external source of driving and information serves as fuel, which can be converted into ${\beta}^{-1}\ln 2$ of thermodynamic free energy per bit. In our case, there is an external source of free energy that in some cases can be harvested more effectively by reducing the system’s statistical entropy, e.g. by concentrating it on favorable states. Here, a bit of information can increase the harvesting rate by a very large amount (much larger than ${\beta}^{-1}\ln 2/\text{bit}$), and information acts more like a catalyst than a fuel [59, 60]. Loosely speaking, this is similar to how information encoded in the sequence of a metabolic enzyme is not itself fuel, but rather allows metabolism to harvest a large amount of fuel from elsewhere.

We finish by mentioning some connections to previous work and future directions. First, our approach may be related to prior work on optimizing power output and free energy transduction in steady-state molecular machines [61, 62, 63, 64, 65, 29, 66]. Here we consider the general problem of optimizing a set of control processes, given a fixed baseline and possible additional constraints on kinetics, topology, and thermodynamics. Previous work does not make the baseline/control distinction; instead, it is mostly concerned with the problem of optimizing system performance with respect to a small set of specific parameters or observables of interest, such as the location of free energy barriers [66, 62, 63, 65], efficiency [66], and the size of fluctuations [61, 64].

There is also an interesting relation between our work and flux balance analysis (FBA) [67, 68, 69]. The goal of FBA is to identify deterministic fluxes in biological metabolic systems that optimize biomass production, or other similar metrics of performance. This can be formulated as a linear program, which may include linear constraints that enforce thermodynamically favored reaction directions [69] (interestingly, in Ref. [70], the authors propose a version of FBA that also accounts for the entropy production rate). Our setting and optimization are different from FBA and its variants. We seek to optimize free energy harvesting at the stochastic level, and our objective involves nonlinear information-theoretic contributions to free energy. In addition, our optimization involves both the steady-state distribution $\bm{p}$ and fluxes ${\boldsymbol{J}}$, which allows us to optimize harvesting due to to internal transitions within coarse-grained mesostates, as in Fig. 4 (b). Nonetheless, investigating the relationship between our approach and FBA is an interesting direction for future work.

Another interesting direction for future work is to consider stochastic fluctuations of free energy harvesting. In particular, the thermodynamic uncertainty relation may be used to study tradeoffs between the entropy production rate, the average harvesting rate (the quantity $\dot{\mathcal{G}}^{\mathrm{tot}}$ considered here), and the fluctuations in the amount of harvested free energy [71, 64]. For biomolecular systems, large fluctuations in harvesting can lead to starvation, suggesting that minimizing fluctuations may be of significant biological importance.

Finally, an interesting direction for future work is to consider free energy harvesting in a system embedded in a fluctuating environment. For example, one may imagine a harvesting system in an environment with fluctuating sugar sources, or with a fluctuating amount of available light. In this setting, it is natural to optimize the harvesting rate under the topological constraint that control fluxes cannot directly change the state of the environment, for instance using bipartite models of Markovian dynamics [72]. It would be interesting to investigate how, under the optimal harvesting strategy, the information flow from the environment to the system varies with the abundance of free energy and complexity of the environment.