# Phase boundaries promote chemical reactions through localized fluxes

Alexandra Shelest  Institute of Physics, École Polytechnique Fédérale de
Lausanne (EPFL), 1015 Lausanne, Switzerland  Hugo Le Roy  Institute of
Physics, École Polytechnique Fédérale de Lausanne (EPFL), 1015 Lausanne,
Switzerland  Daniel M. Busiello  Max Planck Institute for the Physics of
Complex Systems, 01187 Dresden, Germany  Paolo De Los Rios  Institute of
Physics, École Polytechnique Fédérale de Lausanne (EPFL), 1015 Lausanne,
Switzerland  Institute of Bioengineering, École Polytechnique Fédérale de
Lausanne (EPFL), 1015 Lausanne, Switzerland

###### Abstract

One of the hypothesized functions of biomolecular condensates is to act as
chemical reactors, where chemical reactions can be modulated, i.e. accelerated
or slowed down, while substrate molecules enter and products exit from the
condensate. Likewise, the components themselves that take part in the
architectural integrity of condensates might be modified by active (energy
consuming, non-equilibrium) processes, e.g. by ATPase chaperones or by kinases
and phosphatases. In this work, we study how the presence of spatial
inhomogeneities, such as in the case of liquid-liquid phase separation,
affects active chemical reactions and results in the presence of directional
flows of matter, which are one of the hallmarks of non-equilibirum processes.
We establish the minimal conditions for the existence of such spatial
currents, and we furthermore find that these fluxes are maximal at the
condensate interface. These results propose that some condensates might be
most efficient as chemical factories due to their interfaces rather than their
volumes, and could suggest a possible biological reason for the the observed
abundance of small non-fusing condensates inside the cell, thus maximizing
their surface and the associated fluxes.

##  I Introduction

Biomolecular condensates are membraneless cellular organelles that can
specifically control their own composition, being enriched in some components
and depleted in others. This barrierless partitioning offers, besides other
properties, the possibility to concentrate specific enzymes and their
substrates, and to potentially accelerate the rates of the reactions they
catalyze, possibly resulting in fluxes of molecules both to and from
condensates. Furthermore, the molecules that are structurally responsible for
the condensate formation might themselves undergo transformations, driven by
chemical reactions, that tune their conformational states and their
interactions, consequently modulating the condensate properties. For example,
Narayanan  _et al._ [1] found that ruvBL, an AAA+ ATPase, is involved in the
size control of Synphilin1 condensates and Guilhas  _et al._ [2] showed that
the ATPase ParA is actively involved in preventing the fusion of parABS
condensates. Another typical hallmark of active processes, the presence of
fluxes, has been observed in the nucleolus, which is instrumental to the
assembly of ribosomes from rRNA [3]. Directional mRNA fluxes have also been
measured in nuclear speckles [4, 5]. There are several further systems where
the presence of fluxes is hypothesized because mRNA must be sequentially
modified while transiting from one condensate to another [6], e.g. between
P-bodies and stress granules [7]. The presence of fluxes dramatically
underscores the necessity of non-equilibrium driving for the molecular
transformations taking place in condensates.

Physically, condensates are often modelled in the framework of liquid-liquid
phase separation, according to the Flory-Huggins theory ([8], [9], [10]). This
description is only valid at equilibrium and must be extended to account for
active effects. In most works (e.g., [10, 11]), the kinetics is derived close
to equilibrium from linearized fluxes which depend on the chemical potentials
of the species, calculated from the free energy, and on the chemical potential
of the fuel.

A number of studies also investigates the interplay between active condensate
maintenence and the presence of fluxes [10, 12, 11, 13, 14], yet a clear
analysis of the conditions that are necessary for fluxes to be present, and of
their spatial structure, is still lacking. This is precisely what we address
in this work.

In Section II we describe in detail our mathematical framework, where we
derive and extend the customary free-energy approach from a microscopic
perspective. In Section III, we show that the simple presence of activity is
not sufficient for the presence of fluxes, and we thus establish the minimal
set of conditions that must be satisfied in a cell (or chemical system) to
maintain steady-state spatial fluxes in the presence of an interface between
two phases. In Section IV, we study how these fluxes depend on space and on
some of the salient parameters of the system. We find that diffusive fluxes
are always maximal close to the interface, decreasing exponentially away from
it, with the decay governed by the reaction-diffusion length. As a
consequence, depending on this length and on the condensate size, the entire
condensate, or just a region close to its surface, contributes to the fluxes,
with relevant consequences on the current view of these membraneless
organelles as chemical factories.

##  II Derivation of the formalism

We model a system with a fixed number of particles, $N$, where each particle
can be in any of $K$ chemical states. The different chemical states can
correspond to, e.g., different protein conformations that change the way the
particles interact, and hence their preferences for either the condensate or
the dilute phase. We assume that the diffusion constant $D$ does not depend on
the chemical state of the particle (this condition can be relaxed without loss
of generality of our results).
The total energy of the system comprises a contribution from external
potentials, $V(\vec{x}_{i},\sigma_{i})$, and a part due to inter-particle
interactions, $U(\vec{x}_{i},\vec{x}_{j},\sigma_{i},\sigma_{j})$, both
depending on the positions and states of the particles:

$E(\vec{x},\vec{\sigma})=\sum_{i=1}^{N}V(\vec{x}_{i},\sigma_{i})+\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}U(\vec{x}_{i},\vec{x}_{j},\sigma_{i},\sigma_{j})$
(1)

where $\vec{x}=(\vec{x}_{1},...,\vec{x}_{N})$ and
$\vec{\sigma}=(\sigma_{1},...,\sigma_{N})$ denote the ensemble of particle
positions and states, respectively, and for the sake of simplicity, we are
restricting the analysis to one- and two-body interactions.

We describe the system through the probability distribution
$P(\vec{x},\vec{\sigma},t)$, starting from the full $N$-particle Fokker-Planck
equation [15], which is a known approach for single state systems undergoing
phase separation and here we endow particles with multiple internal states
governed by transition rates between each other:

$\begin{split}\partial_{t}P(\vec{x},\vec{\sigma},t)=&\sum_{i=1}^{N}\Big{[}D~{}\nabla_{i}\big{[}\nabla_{i}P(\vec{x},\vec{\sigma},t)\\\
&+\beta P(\vec{x},\vec{\sigma},t)\nabla_{i}E(\vec{x},\vec{\sigma})\big{]}\\\
&+\sum_{\\{{\sigma_{i}^{\prime}}\\}}^{K}\Big{(}k_{\sigma_{i}^{\prime}\sigma_{i}}(\vec{x})P(\vec{x},\vec{\sigma}_{/i},\sigma_{i}^{\prime},t)\\\
&-k_{\sigma_{i}\sigma_{i}^{\prime}}(\vec{x})P(\vec{x},\vec{\sigma}_{/i},\sigma_{i},t)\Big{)}\Big{]}\end{split}$
(2)

where
$\vec{\sigma}_{/i}=(\sigma_{1},...,\sigma_{i-1},\sigma_{i+1},...,\sigma_{N})$
indicates the set of all chemical states apart from $\sigma_{i}$, and
$\beta=1/(k_{B}T)$, with $k_{B}$ the Boltzmann constant.

Under the assumption that all reactions can be seen as activation processes
over energy barriers, the rate of the chemical transition of particle $i$ from
state $\sigma_{i}$ to state $\sigma_{i}^{\prime}$ via a transition state
$[\sigma_{i}\sigma^{\prime}_{i}]$ is

$k_{\sigma_{i}\sigma_{i}^{\prime}}(\vec{x})=k^{0}_{\sigma_{i}\sigma_{i}^{\prime}}e^{-\beta[E(\vec{x},\vec{\sigma}_{/i},[\sigma_{i}\sigma_{i}^{\prime}])-E(\vec{x},\vec{\sigma}_{/i},\sigma_{i})]}e^{\beta\mu^{f}_{\sigma_{i}\rightarrow\sigma_{i}^{\prime}}}\quad,$
(3)

where, for clarity, we separated the state of the $i$th particle from all the
others. This rate is composed of three factors. The first one,
$k^{0}_{\sigma_{i}\sigma_{i}^{\prime}}$, is the intrinsic rate of the
transition as it would happen if the particle were not affected by external
interactions, and it depends on the internal energies of each state and of the
barrier. Yet, in the presence of external potentials and interactions, as
detailed in Eq. 1, both the energies of the states and of the barriers change
(Fig. 1A), resulting in the modifications captured by the second factor in
(3). Here, the rates take the customary Kramers form with an energy barrier
between $\sigma_{i}$ and $\sigma_{i}^{\prime}$, two local minima of the
modified energy landscape (see App. A for details). Finally, the third factor
describes the further modification of the rate in the presence of an active
process that affects the reaction in selected directions, associated to an
effective chemical potential
$\mu^{f}_{\sigma_{i}\rightarrow\sigma_{i}^{\prime}}$ (where $f$ stands for
fuel).

In general, we are interested in the evolution of the concentrations of
particles, which are just proportional to their probabilities obtained by
marginalizing the $N$-body Fokker-Planck equation (2) over the coordinates and
states of all the particles that are not of interest:

$\begin{split}&P(\vec{x}_{i},\sigma_{i},t)=\sum_{\\{\sigma_{j}=1,...K\\}_{j\neq
i}}\int\prod_{j\neq i}d\vec{x}_{j}P(\vec{x},\vec{\sigma},t)\end{split}$ (4)

In the following, for brevity we define the notation

$\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}\equiv\sum_{\\{\sigma_{j}=1,...K\\}_{j\neq
i}}\int\prod_{j\neq i}d\vec{x}_{j}\quad.$ (5)

and we drop the vector sign for single-particle positions. The consequent
marginalization of Eq. (2) leads to the following equation:

$\begin{split}\partial_{t}&P(x_{i},\sigma_{i},t)=D~{}\nabla_{i}\Big{[}\nabla_{i}P(x_{i},\sigma_{i},t)\\\
&+\beta
P(x_{i},\sigma_{i},t)\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}\Big{[}P(\vec{x}_{/i},\vec{\sigma}_{/i}|{x_{i},\sigma_{i}})\nabla_{i}E(\vec{x}_{/i},x_{i},\vec{\sigma}_{/i},\sigma_{i})\Big{]}\\\
&+\sum_{\sigma_{i}^{\prime}\neq\sigma_{i}}^{K}\Big{[}-k^{eff}_{\sigma_{i}\sigma_{i}^{\prime}}(x_{i})P(x_{i},\sigma_{i},t)+k^{eff}_{\sigma_{i}^{\prime}\sigma_{i}}(x_{i})P(x_{i},\sigma_{i}^{\prime},t)\Big{]}\end{split}$
(6)

where we have defined

$\begin{split}k^{eff}_{\sigma_{i}\sigma_{i}^{\prime}}(x_{i})=&k_{\sigma_{i}\sigma_{i}^{\prime}}^{0}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}\Big{(}P(\vec{x}_{/i},\vec{\sigma}_{/i}|{x_{i},\sigma_{i}})\\\
&e^{-\beta
E(\vec{x}_{/i},x_{i},\vec{\sigma}_{/i},[\sigma_{i}\sigma_{i}^{\prime}])+\beta
E(\vec{x}_{/i},x_{i},\vec{\sigma}_{/i},\sigma_{i})}\Big{)}\end{split}$ (7)

and we have used
$P(\vec{x},\vec{\sigma})=P(\vec{x}_{/i},\vec{\sigma}_{/i}|{x_{i},\sigma_{i}})P(x_{i},\sigma_{i})$.

The presence of the conditional probability
$P(\vec{x}_{/i},\vec{\sigma}_{/i}|{x_{i},\sigma_{i}})$ maintains the same
complexity as the full $N$-particle Fokker-Planck equation, and we must thus
use approximations to reduce it to a more manageable form. While it is for
sure not the only possibility, here we are postulating that the $N-1$
particles that we are tracing over adapt to the evolved state of particle $i$
with a much faster rate, eventually settling in their equilibrium state given
$x_{i}$ and $\sigma_{i}$.
As a consequence (see App. A),

$P(\vec{x}_{/i},\vec{\sigma}_{/i}|{x_{i},\sigma_{i}})=\frac{e^{-\beta
E(\vec{x}_{/i},x_{i},\vec{\sigma}_{/i},\sigma_{i})}}{Z(x_{i},\sigma_{i})}$ (8)

where

$Z(x_{i},\sigma_{i})=\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}{e^{-\beta
E(\vec{x}_{/i},x_{i},\vec{\sigma}_{/i},\sigma_{i})}}.$ (9)

From this partition function of $N-1$ particles given the state of particle
$i$, we can define a corresponding conditional free energy:

$F_{\text{cond}}(x_{i},\sigma_{i})=-k_{B}T\ln Z(x_{i},\sigma_{i}).$ (10)

We similarly define the barrier partition function

$\displaystyle Z(x_{i},[\sigma_{i}\sigma_{i}^{\prime}])$ $\displaystyle=$
$\displaystyle\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}{e^{-\beta
E(\vec{x}_{/i},x_{i},\vec{\sigma}_{/i},[\sigma_{i}\sigma_{i}^{\prime}])}}$
(19) $\displaystyle=$ $\displaystyle e^{-\beta
F_{\text{cond}}(x_{i},[\sigma_{i}\sigma_{i}^{\prime}])}\quad.$ (20)

Using these expressions, after some trivial algebra the effective transition
rate Eq. (7) becomes

$k^{eff}_{\sigma_{i},\sigma_{i}^{\prime}}(x_{i})=k^{0}_{\sigma_{i},\sigma_{i}^{\prime}}e^{-\beta\left[F_{\text{cond}}(x_{i},[\sigma_{i}\sigma_{i}^{\prime}])-F_{\text{cond}}(x_{i},\sigma_{i})\right]}$
(21)

Taking into account that all particles in the same state $\sigma$ have the
same probability distribution, and that their concentration is
$c(x,\sigma)=NP(x,\sigma)$, Eq. (6) eventually reduces to

$\begin{split}\partial_{t}&c(x,\sigma,t)=D~{}\nabla\Big{[}\nabla
c(x,\sigma,t)\\\ &+\beta c(x,\sigma,t)\nabla
F_{\text{cond}}(x,\sigma))\Big{]}\\\
&+\sum_{\sigma^{\prime}\neq\sigma}^{K}\Big{[}-k^{eff}_{\sigma\sigma^{\prime}}(x)c(x,\sigma,t)+k^{eff}_{\sigma^{\prime}\sigma}(x)c(x,\sigma^{\prime},t)\Big{]}\end{split}$
(22)

as described in more detail in Appendix A. The diffusive part of Eq.22 is the
usual generalized Cahn-Hilliard equation, where $F_{\text{cond}}$ is the
chemical potential, since it is the derivative of the free energy with respect
to the concentration.

Note that the Kramers form of the rates with free energies instead of energies
emerges as a consequence of choosing Kramers rates for the microscopic
reactions in the full, $N$-body description of the system. In this respect, it
is worth highlighting that while the ratio between the forward and backward
effective reaction rates in Eq. (21) expectedly corresponds to the exponential
of the free-energy difference, the form of the rates themselves cannot be
chosen arbitrarily with the sole constraint of obeying detailed balance.
Instead, a careful treatment of the microscopic reaction is in principle
necessary to know which functional form the rates take. In particular, we
emphasise the relevance of the barrier, as it will be shown in the remainder
of the work to play a crucial role. Furthermore, it is also worth noting that
the same free energy, $F_{\text{cond}}(x_{i},\sigma_{i})$, shaped by the local
environment, dictates both the chemical rates and the drift force,
underscoring the intimate connection between the different parts of Eq. 22.

In this work, we are interested in the relation between spatial
inhomogeneities, non-equilibrium conditions and the presence of spatially
extended fluxes at steady state ($\partial_{t}c(x,\sigma)=0$). In what
follows, without loss of generality, we will describe one-dimensional systems
in a finite box of size $2L$.

Figure 1: The model. Left: Potential landscape of a chemical reaction changing
state i to state j and vice versa. Dashed line shows the potential landscape
without the influence of the background potential. Middle: Active chemical
cycle with three chemical species. Right: Coarse-grained chemical cycle where
one of the species is getting replaced by a second pathway with different
reaction rates than on the first one.

##  III Necessary conditions for spatial fluxes

While the full steady-state solution of the system must satisfy Eq. (22) with
$\partial_{t}c(x,\sigma)=0$, our strategy for this Section is to solve first
the chemical steady state locally, and then check whether the ensuing solution
also separately satisfies the purely diffusion-drift part. If that is the
case, then there are no spatially extended fluxes, otherwise they are present
by necessity.

Chemical reaction network without cycles: Let us first consider a system in
which particles can be in two states, $\sigma=\\{A,B\\}$, such that the
conversion of A into B is fuelled by a chemical potential $\mu_{A\rightarrow
B}$. The local purely chemical steady state is

$\begin{split}c_{A}(x)&=c_{B}(x)\frac{k_{BA}(x)}{k_{AB}(x)}\\\
&=c_{B}(x)\frac{k^{0}_{BA}}{k^{0}_{AB}}e^{-\beta(F_{A}(x)-F_{B}(x))}e^{-\beta\mu_{A\rightarrow
B}}\quad.\end{split}$ (23)

where we have dropped ”cond” from the free energy and written the chemical
state in the subscript for the sake of simplicity. It is important to stress
that, in the absence of cycles, there can be no fluxes on the network at
steady-state, and the solution must respect the mathematical detailed balance,
even if there is an underlying energy-consuming process ($\mu_{A\to B}$)
driving some of the transitions.

Substituting this relation in the steady-state Eq. 22 for $A$, the chemical
part vanishes by construction, and the spatial part is solved by the same
spatial no-flux solution that solves the equation for B (see Appendix B),
ultimately resulting in

$\displaystyle c_{B}(x)$ $\displaystyle=$ $\displaystyle c_{0}e^{-\beta
F_{B}(x)}$ $\displaystyle c_{A}(x)$ $\displaystyle=$ $\displaystyle
c_{0}e^{-\beta F_{A}(x)}e^{-\beta\mu_{A\rightarrow B}}$ (24)

where $c_{0}$ is fixed by the overall normalisation. The presence of a non-
equilibrium driving that does not depend on space is thus not sufficient to
generate fluxes because it simply changes some apparent energies of the system
by a constant multiplicative term.

Similar arguments can be used in the case of more complex chemical reaction
networks, as long as they lack cycles: the concentrations of the states on the
leaves of the network can be expressed in terms of the concentrations of the
nodes they are connected to, in analogy to (23), restoring to the effective
local detailed balance condition used above The iteration of this procedure
reduces all equations to the equation of the root of the tree-like reaction
network (see Appendix C), which can then be solved for spatial no-flux
conditions. We have thus highlighted the first requirement that is necessary
for the presence of spatially extended fluxes: the chemical reaction network
must have cycles. Yet, the presence of cycles is not sufficient in itself for
the presence of fluxes.

General chemical reaction networks: For a general chemical reaction network,
the local chemical steady-state for particles in state $i$ can be written
using the spanning-tree form of the solution [16]:

$c_{\sigma}(x)=\frac{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}k_{l}^{\sigma}(x)}{\mathcal{N}(x)}$
(25)

where $\\{\mathcal{T}\\}$ is the set of spanning trees over the network,
$k_{l}^{\sigma}$ is the rate of the reaction on the network edge $l$ belonging
to $\mathcal{T}$, directed toward $\sigma$, and $\mathcal{N}(x)$ is a general
normalisation that does not depend on the state $\sigma$. Using (21), after
some algebra it is possible to write this expression relative to a reference
state concentration (say, the concentration of state $1$) as

$\displaystyle c_{\sigma}(x)$ $\displaystyle=$ $\displaystyle
c_{1}(x)e^{-\beta\left(F_{\sigma}(x)-F_{1}(x)\right)}\frac{\sum_{\\{\mathcal{T}\\}}e^{-\beta\mu_{1\sigma}}{\prod_{l\in\mathcal{T}}}k_{l}^{1}(x)}{\sum_{\\{\mathcal{T}\\}}\prod_{l\in\mathcal{T}}k_{l}^{1}(x)}$
(26) $\displaystyle=$ $\displaystyle
c_{1}(x)e^{-\beta\left(F_{\sigma}(x)-F_{1}(x)\right)}S_{1\sigma}(x)$

where $\mu_{1\sigma}$ is the total chemical potential due to the fuel
accumulated over the chemical reaction pathway from states $1$ and $\sigma$
(for the detailed derivation see Appendix D).

According to Eq. (21), the product of the rates over the edges of each
spanning tree depends on the exponential of the sum of all (free) energies,
which is the same for each tree and thus cancels in Eq. (26). Consequently
$S_{1\sigma}(x)$ depends only on the fuel chemical potentials associated to
the reactions, and on the barriers of the transitions. Substituting this
expression in the steady-state equation for $c_{i}(x)$, the chemical reactions
part vanishes by construction, while the spatial part acts solely on the
reference concentration, i.e., $c_{1}(x)$ in this case. The resulting equation
coincides with the stationary condition of $c_{1}(x)$, and thus can be
consistently solved, only if $\nabla{S_{1\sigma}(x)}=0$, that is,
$S_{1\sigma}(x)$ is constant in space (Appendix D). This condition is
satisfied in three cases: i) the system is at equilibrium
($\mu_{1\sigma}=0,\forall\sigma$) and thus $S_{1\sigma}(x)=1$, or ii) the
barriers do not depend on space, or iii) the barrier energies all depend on
space in the same way, thus leading to no spatial dependence when computing
the rate (see Eq. (26)). The two-state system and the cycle-less cases are
trivially recovered because there is only one spanning tree and, as such, the
system satisfies an effective equilibrium, as detailed above.

We can thus summarise the necessary requirements for diffusive fluxes:

  * •

Activity: At least one of the chemical reactions has to be out of equilibrium.

  * •

Chemical cycles: Without a chemical cycle, the systems adjusts into a pseudo-
equilibrium state.

The presence of activity and chemical cycles ensures that there will be
chemical fluxes, i.e., the chemical steady state does not locally satisfy
detailed balance. However, these fluxes only govern the local conversion of
species and do not necessarily turn into spatial fluxes, unless a third
condition is satisfied:

  * •

Spatial dependency of energies. The energy landscape of the chemical reaction,
and in particular the energy barriers of the transitions, must depend on
space.

The minimal system with a cycle contains three species, $A$, $B$ and $C$. This
cycle can describe, for example, the conversion of $A$ into $B$ via an active
and a passive pathway, so that the third species $C$ represents a composite
state of $A$ interacting with an enzyme. This three-state system can be
studied either as just sketched or coarse-grained to a two-state system, with
a passive pathway and an active pathway that is absorbing all interaction with
the enzyme and the $C$ state (see Appendix E and Fig. 1). It is important to
note that, in the presence of multiple reaction pathways, a two-state system
shows cycles, and can thus exhibit fluxes.

##  IV Spatial structure of the fluxes

Having established our formalism, we can now address the problem of the
structure of fluxes in the presence of a spatial partition of the system
induced by non-homogeneous free energies $F(x,\sigma)$ and
$F(x,[\sigma^{\prime}\sigma])$. Here, we assume that the functional space-
dependent form of the free energies is given as an external potential:
$F(x,\sigma)=V(x,\sigma)$ and
$F(x,[\sigma^{\prime}\sigma])=V(x,[\sigma^{\prime}\sigma])$. From a biological
perspective, such a scenario might capture, for example, the physics of a
chemical reaction system where the different molecules (e.g. substrates,
enzymes and products) are sensitive to the presence of phase-separated
(condensate) regions that can selectively attract or repel them.
Alternatively, the active maintenance of liquid-liquid phase separation might
be captured by our framework, with the caveat that the inhomogeneous free
energies would depend on some of the concentrations themselves, thus leading
in principle to self-consistent equations. At steady-state, and within a
certain degree of approximation, we are confident that our conclusions are
valid also in this setting.

In what follows, we study a three-state system with the states
$\sigma=\\{A,B,C\\}$. To illustrate our results, we restrict our analysis to
the specific case when only one reaction, from A to C, is actively driven.

###  IV.1 Step potentials

As the simplest model for an interface we use a step potential of the form

$V_{\sigma}(x)=s_{\sigma}V_{\sigma}^{0}\Theta(x)$ (27)

where $V_{\sigma}^{0}$ determines the amplitude of the potential, $s_{\sigma}$
determines its sign and $\Theta(x)$ is the Heaviside Theta function. We write
the energy of the barrier as

$V_{[\sigma\sigma^{\prime}]}(x)=V_{[\sigma\sigma^{\prime}]}^{0}(C_{[\sigma\sigma^{\prime}]}+\Theta(x))$
(28)

with the offset $C_{[\sigma\sigma^{\prime}]}$ that guarantees that the energy
barriers are always higher than the initial and final energies in the
reactions, a condition necessary for the Kramers approximation to be
legitimate.
The analytical solutions for Eq. (22) are obtained by solving it for $x<0$ and
$x>0$ separately and connecting the solutions at $x=0$ (for more details, see
App. G). The solution on the left of the interface reads:

$\vec{c}^{L}(x)=A^{L}_{0}\vec{v}^{L}_{0}+\sum_{i=1,2}A^{L}_{i}\cosh{\Big{(}\sqrt{\lambda^{L}_{i}}(x+L)\Big{)}}\vec{v}^{L}_{i}$
(29)

and similarly on the right of the interface, where the argument of the
hyperbolic cosine changes to $x-L$ and $A^{L}_{i}\to A^{R}_{i}$ are constants
that are fixed by the boundary conditions at the interface. The parameters
$\lambda_{i}^{L/R}$ for $i=1,2$ are the two non-zero eigenvalues of the
chemical transition matrix on each side of the interface (the same approach
can be used to find spherically symmetric solutions in three dimensions, which
is more pertinent to real condensates, see App. H).

Representative concentration profiles are shown in Fig. 2A (the values of the
used parameters are reported in the figure caption). Because of the structure
of the solution (29), the largest deviations from a constant solution are
found at the interface (inset in Fig. 2A), and as a consequence the spatial
fluxes, that for the step potential are just the first derivatives of the
concentrations, are maximal at the interface and decay away from it (Fig. 2B).
It is important to stress here that the absence of spatial fluxes far away
from the interface does not imply the absence of fluxes over the chemical
reaction network (inset in Fig. 2B). Rather, away from the interface reactions
are chemically compensated and there is no ”spill-over” generating fluxes in
space.

Since the potentials are constant on both sides of the interface, there is no
drift term, and as a consequence the fluxes are purely diffusive and decay in
space according to the reaction-diffusion length $l^{rd}_{L/R}$, which is
asymptotically dominated by the smallest (non-zero) eigenvalue of the rate
matrix of the chemical reactions, $l^{rd}_{L/R}=(\lambda_{min}^{L/R})^{-1/2}$
(for simplicity, here we consider all diffusion constants to be identical, and
the reaction-diffusion length is actually related to the diffusion-normalized
chemical rates matrix; the same treatment can be generalized to systems with
state-dependent diffusion constants. In regions where chemical rates are
sufficiently slow compared to the diffusive rate, the diffusive fluxes show an
almost linear behaviour from the system boundaries ($x=\pm L$) to the
interface $x=0$ (Fig. 2C). This result stems from Eq. (29): the fluxes are
proportional to the derivative of the concentrations, so they will be
proportional to $\sinh((x\pm L)/l_{L/R}^{rd})$, which is roughly linear over
the full spatial range if $L/l_{L/R}^{rd}\ll 1$. In such a situation, fluxes
can “penetrate” fully in at least one of the two regions.

The behavior of the fluxes with respect to the chemical driving is shown in
Fig. 2D, where we report their value at the interface. Fluxes vanish at
equilibrium ($\mu_{A\rightarrow C}=0$), as it must be, and increase (in
absolute value) proportionally to $\mu_{A\rightarrow C}$. While this relation
is expected for small values of $\mu_{A\to C}$ as from a first-order
perturbative approach, it breaks down for higher values of $\mu_{A\rightarrow
C}$, becoming non-monotonic: one of the fluxes ($J_{A}$) vanishes for large
non-equilibrium driving, and one ($J_{C}$) changes signs, namely, it changes
direction. Concomitantly, one of the chemical species approaches extinction
(Fig. 2E). While this behavior can be understood because as $\mu_{A\to C}$
grows state A is depleted and thus it cannot give rise to significant spatial
fluxes, it highlights how, in more general scenarios and for more complex
chemical reaction networks, the behavior for intermediate and large driving
can be highly non-trivial.

Since in most of the literature, spatial inhomogeneities and chemical
reactions are studied for two-state systems. In Appendix E we show how to
coarse-grain a three-state chemical reaction network with cycles to a two-
state ones, with cycles due to the presence of two different reaction pathways
between the two remaining states, and we further re-derive the results of this
section for that setup.

Figure 2: A: Concentration profiles of the three species for a step potential.
Unless mentioned otherwise, we set the parameters for all calculations to be
$E_{A}=-0.5$, $E_{B}=2.$, $E_{C}=2.5$, $E_{AB}=2.$, $E_{BC}=2.5$, $E_{AC}=3.$,
$s_{B}=s_{C}=-s_{A}=1$, $\mu_{A\rightarrow C}=2$ and
$C_{[AB]}=C_{[BC]}=C_{[AC]}=0.5$. Inset: Profiles of ‘excess’ concentration,
calculated as the difference of the concentration value at x and the
concentration value at the respective boundary of the system. B: Spatial
fluxes of the three species. Inset: Chemical fluxes of the three species. C:
Fluxes for $E_{AB}=4.$, $E_{BC}=4.5$, $E_{AC}=5.$, resulting in a linear
profile for $x>0$. D: Diffusive fluxes at $x=0$ as a function of chemical
potential $\mu_{A\rightarrow C}\equiv\mu$. E: Concentrations at $x=L$ (solid)
and $x=-L$ (dotted) as a function of $\mu_{A\rightarrow C}$, normalized to the
respective equilibrium concentration values ($\mu=0$). As seen in A, $\rho(x)$
varies little on each side of the step, so it suffices to study its behaviour
at the extremities of the box.

###  IV.2 Sigmoidal potentials

To study the effects of a finite interface width, and to create a more
realistic model of a biological interface, we turn to a sigmoidal potential
$V(x,\sigma)$ with the functional form here of a hyperbolic tangent. This
choice, that introduces a smooth interface of controlled width, is inspired by
the profile of the interface of phase-separated regions [10], which would
provide the relevant background energy for the molecules involved in the
chemical reaction network.

We thus choose potentials of the form

$V_{\sigma}(x)=0.5~{}V^{0}_{\sigma}(1+s_{\sigma}\tanh(n\cdot x))$ (30)

where $V^{0}_{\sigma}$ is controlling the potential amplitude and $n$ is its
slope. The sign of the hyperbolic tangent is set by $s_{\sigma}\in\\{-1,1\\}$.
For the energy barriers, we choose

$V_{[\sigma^{\prime}\sigma]}(x)=0.5~{}V^{0}_{[\sigma^{\prime}\sigma]}(C_{[\sigma^{\prime}\sigma]}+\tanh(n\cdot
x))$ (31)

where $C_{[\sigma^{\prime}\sigma]}$ is an additional offset which guarantees
that the energy barriers are always higher than the initial and final energies
in the reactions, a condition necessary for the Kramers approximation to be
legitimate.

In the absence of analytical solutions of Eq. (22) with these potentials, we
solve it numerically (see Appendix F for a description of the method). The
resulting concentrations have, as expected, a sigmoidal shape, centred at
$x=0$ (Fig. 3A), different from the equilibrium ones ($\mu_{AC}=0$, Fig. 3A
dashed lines). The behavior of the concentrations as a function of the
sharpness of the interface is shown in Fig. 3B (only the results for the
concentration of species A are shown) and expectedly, they approach a limit
profile for $n\to\infty$, corresponding to step potentials (here we have
chosen the same $n$ for all species, while of course it could be different for
each one of them). As the interface becomes sharper the fluxes also approach
the step-function limit behaviour.

We again observe linear fluxes for the same parameter regime of energies and
energy barriers as for the step function (Fig. 3D). This parameter range
corresponds to reaction-diffusion lengths much bigger than the system size,
highlighting the consistency with the solution with the step potential.

In order to assess whether a smoother interface is increasing or decreasing
the net flux with respect to the step potential case, we study the excess
average total flux magnitude. It is calculated as an integral of the absolute
flux value over the system length for a slope $n$ minus the same quantity for
$n\rightarrow\infty$ (step function). We find that it is positive and
proportional to $n^{-t}$ with exponents $t$ very close to 1 (for Fig. 3E,
$t_{A}=0.95$, $t_{B}=0.98$, $t_{C}=0.87$). This shows that a smoother
interface corresponds to higher spatial fluxes, a result that can possibly be
assessed in experimental observations of known active condensates.

Figure 3: Sigmoidal potentials for the same parameters as the step potentials,
with $C_{[AB]}=C_{[BC]}=C_{[AC]}=3$ (a change to ensure the barrier potentials
have the same values at the system boundaries). A: Comparison of the
concentration profiles out of equilibrium with the equilibrium case (dashed
lines) for the three species A, B and C for $n_{A}=n_{B}=n_{C}=n=1$. B:
Concentration profiles of species A with different tanh slopes $n$. C: Flux
profiles of species A with different tanh slopes $n$. D: Concentration
profiles of species A for $E_{AB}=4.$, $E_{BC}=4.5$, $E_{AC}=5.$, resulting in
an approximately linear profile for $x>0$. E: Integral of the absolute value
of the flux excess over the flux found for a step potential, averaged over the
system length, as a function of slope $n$. The dashed line shows the fit to
the last 6 points of the function $y(n)=n^{a}10^{b}$. The retrieved fitting
parameters a are very close to -1.

##  V Discussion

For condensates to act as chemical factories, there must be an influx of
substrates and outflux of products, implying that fluxes are necessarily
present across their interfaces. Likewise, when structural components of
condensates are actively modified, so to change their state and condensation
propensity, outgoing and ingoing fluxes of molecules in different states have
to be expected.

Here we focused on a spatial description of the condensates to explore the
conditions that are necessary for the presence of such fluxes, and their
spatial structure. We have found that two non-trivial conditions for the
presence of fluxes are the presence of cycles in the chemical reaction
networks, and the space and state dependence of the energy barriers governing
the transition rates. Of course, and more trivially, when these conditions are
met, fluxes are present only if reaction are driven away from thermodynamic
equilibrium by energy consumption (some $\mu_{\sigma\sigma^{\prime}}\neq 0$).

If this set of conditions is satisfied, then fluxes are present, and are
mostly localized at the condensate interface, extending on both sides by,
essentially, the reaction-diffusion length $l^{rd}$. As a consequence, regions
of the system that are far away (relative to $l^{rd}$) from the interface do
not contribute to the transport of substrates and products, because the non-
equilibrium fluxes over the chemical reaction network are fully chemically
compensated, and do not leak in space. Therefore, since only a limited region
close to the interface can act as a chemical factory, it is tempting to
speculate that maintaining the size of condensates limited, instead of letting
them all coalesce as they would if they where just following the laws of
equilibrium thermodynamics, might be useful to increase the overall condensate
surface and the associated effective reactive region. Furthermore, our results
show that the net flux is decreasing with increasing interface slope, allowing
to speculate that condensates functioning as chemical factories will be
characterized by a smooth, fuzzy interface and negligible surface tension.

In this respect, renewed attention should be devoted to studying the
properties of condensates close to their interfaces in more detail to better
characterize their molecular arrangements. Indeed, numerous works, both
numerical [17, 18] and experimental [5, 4, 19, 20, 21], have shown that the
surfaces of condensates have distinct physical properties from the bulk and
the surroundings of the condensate. Furthermore, there is an increasing number
of observations that aggregation and fibril formation rates are enhanced at
the condensate interface [22, 23, and others]. Simulations at different levels
of detail provide the same results ([24, 18]).

From a modelling perspective, our work shows the relevance of accounting for
the spatial dependence of the chemical rates, as well as the explicit
treatment of the energy barriers. Some authors (e.g., [11, 25]) choose
$k^{0}_{ij}$ and $\mu_{i\rightarrow j}$ to be space-dependent, absorbing
$\epsilon_{ij}(x)$ and $V_{i}(x)$, respectively, into their definitions.
Comparing to such works, it should be kept in mind how the constraints we
derived here are translated with differing definitions of the variables.
Furthermore, it is important to stress that there is an intimate connection
between the spatial and chemical parts of the governing equation of the
system, Eq.22: space-dependent chemical transition rates should be accompanied
by corresponding drift terms.

While this work has been devoted to elucidating the physical rules that are
necessary for the appearance of spatial fluxes, and we have started to
highlight some of their properties, it ultimately calls for a better
experimental characterization both of the biochemical networks underlying
active processes in phase-separated systems, and of the interactions that
determine the phase separation itself and, as a consequence, the perturbation
of the chemical rates. The consequent implementation of the present findings
in more detailed models, possibly closer to biochemical/biological instances,
will provide avenues to understand the functional advantages of certain
condensates, of their compositions and of their size distributions.

## Acknowledgements

We thank Alessandro Barducci and David Zwicker for fruitful discussions, as
well as Francesco Piazza for a critical reading of the manuscript. This work
was supported by the Swiss National Science Foundation under grant
CRSII5_193740.

## References

  * Narayanan  _et al._ [2019] A. Narayanan, A. Meriin, J. O. Andrews, J.-H. Spille, M. Y. Sherman, and I. I. Cisse, A first order phase transition mechanism underlies protein aggregation in mammalian cells, eLife 8, e39695 (2019).
  * Guilhas  _et al._ [2020] B. Guilhas, J.-C. Walter, J. Rech, G. David, N. O. Walliser, J. Palmeri, C. Mathieu-Demaziere, A. Parmeggiani, J.-Y. Bouet, A. Le Gall, and M. Nollmann, Atp-driven separation of liquid phase condensates in bacteria, Molecular Cell 79, 293 (2020).
  * Correll  _et al._ [2019] C. C. Correll, J. Bartek, and M. Dundr, The nucleolus: A multiphase condensate balancing ribosome synthesis and translational capacity in health, aging and ribosomopathies, Cells 8, 10.3390/cells8080869 (2019).
  * Wang  _et al._ [2018] K. Wang, L. Wang, J. Wang, S. Chen, M. Shi, and H. Cheng, Intronless mrnas transit through nuclear speckles to gain export competence, Journal of Cell Biology 217, 3912 (2018).
  * Ishov  _et al._ [2020] A. M. Ishov, A. Gurumurthy, and J. Bungert, Coordination of transcription, processing, and export of highly expressed rnas by distinct biomolecular condensates, Emerging Topics in Life Sciences 4, 281 (2020).
  * Hondele  _et al._ [2020] M. Hondele, S. Heinrich, P. De Los Rios, and K. Weis, Membraneless organelles: phasing out of equilibrium, Emerging topics in life sciences 4, 343 (2020).
  * Kedersha  _et al._ [2005] N. Kedersha, G. Stoecklin, M. Ayodele, P. Yacono, J. Lykke-Andersen, M. J. Fritzler, D. Scheuner, R. J. Kaufman, D. E. Golan, and P. Anderson, Stress granules and processing bodies are dynamically linked sites of mrnp remodeling, The Journal of cell biology 169, 871 (2005).
  * Flory [1942] P. J. Flory, Thermodynamics of High Polymer Solutions, J. Chem. Phys. 10, 51 (1942).
  * Huggins [1942] M. L. Huggins, Some properties of solutions of long-chain compounds., The Journal of Physical Chemistry 46, 151 (1942).
  * Weber  _et al._ [2019] C. A. Weber, D. Zwicker, F. Jülicher, and C. F. Lee, Physics of active emulsions, Reports on Progress in Physics 82, 064601 (2019).
  * Kirschbaum and Zwicker [2021] J. Kirschbaum and D. Zwicker, Controlling biomolecular condensates via chemical reactions, Journal of The Royal Society Interface 18, 20210255 (2021).
  * Bauermann  _et al._ [2022] J. Bauermann, S. Laha, P. M. McCall, F. Jülicher, and C. A. Weber, Chemical kinetics and mass action in coexisting phases, Journal of the American Chemical Society 144, 19294 (2022), pMID: 36241174.
  * Osmanovic and Franco [2022] D. Osmanovic and E. Franco, Chemical reaction motifs driving non-equilibrium behaviors in phase separating materials (2022), arXiv:2207.10135 [cond-mat.soft] .
  * Zwicker [2022] D. Zwicker, The intertwined physics of active chemical reactions and phase separation, Current Opinion in Colloid Interface Science 61, 101606 (2022).
  * Archer and Evans [2004] A. J. Archer and R. Evans, Dynamical density functional theory and its application to spinodal decomposition, The Journal of chemical physics 121, 4246 (2004).
  * Schnakenberg [1976] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Reviews of Modern physics 48, 571 (1976).
  * Cho and Jacobs [2023] Y. Cho and W. M. Jacobs, Tuning nucleation kinetics via nonequilibrium chemical reactions, Physical Review Letters 130, 10.1103/physrevlett.130.128203 (2023).
  * Farag  _et al._ [2022] M. Farag, S. R. Cohen, W. M. Borcherds, A. Bremer, T. Mittag, and R. V. Pappu, Condensates of disordered proteins have small-world network structures and interfaces defined by expanded conformations, bioRxiv , 2022 (2022).
  * Böddeker  _et al._ [2022] T. J. Böddeker, K. A. Rosowski, D. Berchtold, L. Emmanouilidis, Y. Han, F. H. T. Allain, R. W. Style, L. Pelkmans, and E. R. Dufresne, Non-specific adhesive forces between filaments and membraneless organelles, Nature Physics 18, 571 (2022).
  * Folkmann  _et al._ [2021] A. W. Folkmann, A. Putnam, C. F. Lee, and G. Seydoux, Regulation of biomolecular condensates by interfacial protein clusters, Science 373, 1218 (2021).
  * Wu  _et al._ [2023] T. Wu, M. King, M. Farag, R. Pappu, and M. Lew, Single fluorogen imaging reveals spatial inhomogeneities within biomolecular condensates., Biorxiv: the Preprint Server for Biology (2023).
  * Lipiński  _et al._ [2022] W. P. Lipiński, B. S. Visser, I. Robu, M. A. A. Fakhree, S. Lindhoud, M. M. A. E. Claessens, and E. Spruijt, Biomolecular condensates can both accelerate and suppress aggregation of alpha-synuclein, Science Advances 8, eabq6495 (2022).
  * Shen  _et al._ [2022] Y. Shen, A. Chen, W. Wang, Y. Shen, F. S. Ruggeri, S. Aime, Z. Wang, S. Qamar, J. R. Espinosa, A. Garaizar, P. S. George-Hyslop, R. Collepardo-Guevara, D. A. Weitz, D. Vigolo, and T. P. J. Knowles, Solid/liquid coexistence during aging of fus condensates, bioRxiv 10.1101/2022.08.15.503964 (2022).
  * Toprakcioglu  _et al._ [2022] Z. Toprakcioglu, A. Kamada, T. C. T. Michaels, M. Xie, J. Krausser, J. Wei, A. Saric, M. Vendruscolo, and T. P. J. Knowles, Adsorption free energy predicts amyloid protein nucleation rates, Proceedings of the National Academy of Sciences 119, e2109718119 (2022).
  * Laha  _et al._ [2024] S. Laha, J. Bauermann, F. Jülicher, T. C. Michaels, and C. A. Weber, Chemical reactions regulated by phase-separated condensates, arXiv preprint arXiv:2403.05228 (2024).
  * Uhl  _et al._ [2021] M. Uhl, V. Weissmann, and U. Seifert, Propagator for a driven Brownian particle in step potentials, Journal of Physics A Mathematical General 54, 065002 (2021), arXiv:2009.03201 [cond-mat.stat-mech] .
  * Piazza [2022] F. Piazza, The physics of boundary conditions in reaction–diffusion problems, The Journal of Chemical Physics 157, 234110 (2022).

##  Appendix A Connecting a Fokker-Planck description to a free energy
description

In order to determine the conditional probability
$P(\vec{x}_{/i},\vec{\sigma}_{/i}|{x_{i},\sigma_{i}})$, we study the system
close to equilibrium, so that all particles but the one in consideration are
equilibrated. Then, we can connect the Fokker-Planck description of the
problem to a statistical physics description [15]. For better readability, in
this section we denote $P(x_{i},\sigma_{i})=P_{i}$,
$P_{/i}=P(\vec{x}_{/i},\vec{\sigma}_{/i}|x_{i},\sigma_{i})$ and
$E_{i,\sigma_{i}}\equiv E(\vec{x}_{/i},x_{i},\vec{\sigma}_{/i},\sigma_{i})$.
We write a general free energy as the sum of internal energy, entropy and a
normalization factor for the conditional probability.

$\begin{split}F[P]=&\sum_{\sigma_{i}}\int
dx_{i}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}\Big{[}E_{i,\sigma_{i}}\cdot
P_{/i}P_{i}+k_{B}TP_{/i}P_{i}\cdot\ln[P_{/i}P_{i}]\Big{]}\\\
&+\lambda(x)(1-\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}P_{/i})\end{split}$
(32)

Taking the variation of this term with respect to $P_{/i}$, we find:

$\begin{split}&P_{/i}=\frac{e^{-\beta
E_{i,\sigma_{i}}-1-\beta\frac{\lambda(x)}{P_{i}}}}{P_{i}}\end{split}$ (33)

which we can re-write using the normalization condition:

$\begin{split}P_{/i}=\frac{\exp\big{[}-\beta
E_{i,\sigma_{i}}\big{]}}{Z(x_{i},\sigma_{i})}\end{split}$ (34)

with the conditional partition function

$\begin{split}Z(x_{i},\sigma_{i})=\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}e^{-\beta
E_{i,\sigma_{i}}}\end{split}$ (35)

This is the full partition function for the remaining $N-1$ particles. The
corresponding conditional free energy is
$F_{\text{cond}}(x_{i},\sigma_{i})=-k_{B}T\ln Z(x_{i},\sigma_{i})$.

Plugging this expression in Eq. 32, the free energy becomes

$\displaystyle F[P_{i}]$ $\displaystyle=$ $\displaystyle\sum_{\sigma_{i}}\int
dx_{i}\Big{[}P_{i}F_{\text{cond}}(x_{i},\sigma_{i})+k_{B}TP_{i}\ln
P_{i}\Big{]}$ (36) $\displaystyle=$ $\displaystyle U[P_{i}]-TS[P_{i}]$

where the last expression highlights the natural split of the free energy into
the internal energy and entropy for particle $i$. We can use these results to
re-express the diffusive part of Eq. 6:

$\begin{split}\partial_{t}P_{i}&=D~{}\nabla_{i}\big{[}\nabla_{i}P_{i}+\beta
P_{i}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}\frac{e^{-\beta
E_{i,\sigma_{i}}}}{Z(x_{i},\sigma_{i})}\nabla_{i}E_{i,\sigma_{i}}\big{]}\\\
&=D~{}\nabla_{i}\big{[}P_{i}\nabla_{i}\big{(}\ln(P_{i})+\ln(Z(x_{i},\sigma_{i}))\big{)}\big{]}\\\
&=D~{}\nabla_{i}\big{[}P_{i}\nabla_{i}\frac{\delta F[P_{i}]}{\delta
P_{i}}\big{]}=D~{}\nabla_{i}\big{[}\nabla_{i}P_{i}+\nabla_{i}\frac{\delta
U[P_{i}]}{\delta P_{i}}\big{]}\end{split}$ (37)

A full description of the system also requires us to specify
$k^{eff}_{\sigma_{i}\sigma_{i}^{\prime}}(x_{i})$ as a function of
thermodynamic quantities. As outlined in Section II, we apply Kramers
approximation, assuming that both the initial and the final points of the
chemical reaction form local minima in the energy landscape along the reaction
coordinate, separated by the energy barrier which we denote as
$[\sigma_{i}\sigma_{i}^{\prime}]$. Then:

$\begin{split}&k^{eff}_{\sigma_{i}\sigma_{i}^{\prime}}(x_{i})=k^{0}_{\sigma_{i}\sigma_{i}^{\prime}}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}e^{-\beta
E_{i,[\sigma_{i}\sigma_{i}^{\prime}]}+\beta E_{i,\sigma_{i}}}P_{/i}\\\
&=k^{0}_{\sigma_{i}\sigma_{i}^{\prime}}\exp\Big{[}-\beta\Big{\\{}-\frac{1}{\beta}\ln\Big{[}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}e^{-\beta
E_{i,[\sigma_{i}\sigma_{i}^{\prime}]}+\beta
E_{i,\sigma_{i}}}P_{/i}\Big{]}\Big{\\}}\Big{]}\end{split}$ (38)

Consistently with the previous approximation, we use again that all particles
but the $i$th are in local equilibrium, and therefore we can express the
conditional probability by Eq. 34. Inserting it in Eq. 38, the denominator of
Eq. 34 cancels with one of the two energies from Kramer’s rates.

$\begin{split}k^{eff}_{\sigma_{i}\sigma_{i}^{\prime}}(x_{i})&=k^{0}_{\sigma_{i}\sigma_{i}^{\prime}}\exp\Big{[}-\beta\Big{\\{}-\frac{1}{\beta}\ln\Big{[}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{/i}\frac{e^{-\beta
E_{i,[\sigma_{i}\sigma_{i}^{\prime}]}}}{Z(x_{i},\sigma_{i})}]\Big{]}\Big{\\}}\Big{]}\end{split}$
(39)

By analogy with Eq. 34, we can define the conditional partition function of
the energy barrier and the corresponding conditional barrier free energy
$F_{\text{cond}}(x_{i},[\sigma_{i}\sigma_{i}^{\prime}])=-k_{B}T\ln
Z(x_{i},[\sigma_{i}\sigma_{i}^{\prime}])$. With this,

$k^{eff}_{\sigma_{i}\sigma_{i}^{\prime}}(x_{i})=k^{0}_{\sigma_{i}\sigma_{i}^{\prime}}\exp\left\\{-\beta\left[F_{\text{cond}}(x_{i},[\sigma_{i}\sigma_{i}^{\prime}])-F_{\text{cond}}(x_{i},\sigma_{i})\right]\right\\}$
(40)

Remarkably, the conditional free energies that modify the transition rates in
Eq.40 can also be written as functional derivatives of the internal energy,
just as the drift term, stressing the intimate connection between drift and
chemical reactions.

Since the probability $P_{i}=P(x_{i},\sigma_{i})$ is the same for all
particles in state $\sigma=\sigma_{i}$, we drop the index $i$ and, using the
previous results, we can rewrite the full reaction-diffusion equation for the
probability $P(x,\sigma)$ as:

$\displaystyle\partial_{t}P(x,\sigma)=D~{}\nabla\left[P(x,\sigma)\nabla\frac{\delta
F[P]}{\delta P(x,\sigma)}\right]+$ (41) $\displaystyle+$
$\displaystyle\sum_{\sigma^{\prime}}\left[k^{eff}_{\sigma^{\prime}\sigma}(x)P(x,\sigma^{\prime})-k^{eff}_{\sigma\sigma^{\prime}}(x)P(x,\sigma)\right]$

where the usual form of the generalized Cahn-Hilliard equation for the
diffusive part.

To further emphasize the connection with the generalized Cahn-Hilliard
equation, it is possible to express Eq.41 in terms of concentrations, by
multiplying it by $N$ and using $c(x,\sigma)=NP(x,\sigma)$. We obtain

$\displaystyle\partial_{t}c(x,\sigma)=D~{}\nabla\left[c(x,\sigma)\nabla\frac{\delta
F[c]}{\delta c(x,\sigma)}\right]+$ (42) $\displaystyle+$
$\displaystyle\sum_{\sigma^{\prime}}\left[k^{eff}_{\sigma^{\prime}\sigma}(x)c(x,\sigma^{\prime})-k^{eff}_{\sigma\sigma^{\prime}}(x)c(x,\sigma)\right]$

where we have used that, apart from an irrelevant additive constant, Eq.36
becomes

$F[c]=\frac{1}{N}\sum_{\sigma}\int
dx\Big{[}c(x,\sigma)F_{\text{cond}}(x,\sigma)+k_{B}Tc(x,\sigma)\ln
c(x,\sigma)\Big{]}\\\ $ (43)

and that taking its functional derivative with respect to $P$ becomes, due to
the $1/N$ term, the functional derivative with respect to $c$. Furthermore,
also the conditional free energies in the exponent in Eq.40 can also be
written as functional derivatives of the internal energy with respect to
concentrations, making Eq.42 only dependent on $c$.

With this result, we can rewrite Eq. 6 as Eq. 22, as is done in Section II.

##  Appendix B Derivation of no-flux-condition for two species out of
equilibrium

For two species A and B, we can write the time evolution of their
concentrations as:

$\begin{split}\partial_{t}c_{A}(x)=&D\partial_{x}[\partial_{x}c_{A}(x)+\beta
c_{A}(x)\partial_{x}V_{A}(x)]\\\ &-k_{AB}c_{A}(x)+k_{BA}c_{B}(x)\end{split}$
$\begin{split}\partial_{t}c_{B}(x)=&D\partial_{x}[\partial_{x}c_{B}(x)+\beta
c_{B}(x)\partial_{x}V_{B}(x)]\\\ &+k_{AB}c_{A}(x)-k_{BA}c_{B}(x)\end{split}$
(44)

Solving the chemical part of the equation at steady state out of equilibrium,
we get:

$\frac{c_{A}(x)}{c_{B}(x)}=\frac{k_{BA}(x)}{k_{AB}(x)}=\frac{k^{0}_{BA}}{k^{0}_{AB}}e^{-\beta(V_{A}(x)-V_{B}(x)+\mu_{A\rightarrow
B})}$ (45)

We can express $c_{B}(x)$ as a function of $c_{A}(x)$:

$c_{B}(x)=\frac{k^{0}_{AB}}{k^{0}_{BA}}e^{\beta(V_{A}(x)-V_{B}(x)+\mu_{A\rightarrow
B})}c_{A}(x)$ (46)

We can plug this into Eq. 44 at steady state. The chemical part of the
equation disappears by construction, so we are left only with the diffusive
flux:

$\begin{split}0&=\frac{k^{0}_{AB}}{k^{0}_{BA}}e^{\beta(V_{A}(x)-V_{B}(x)+\mu_{A\rightarrow
B})}\\\ &\cdot\partial_{x}[\partial_{x}c_{A}(x)+\beta
c_{A}(x)\partial_{x}V_{A}(x)]\end{split}$ (47)

Canceling
$\frac{k^{0}_{AB}}{k^{0}_{BA}}e^{\beta(V_{A}(x)-V_{B}(x)+\mu_{A\rightarrow
B})}$ on both sides, we recover the equation for $c_{A}(x)$. In this case, we
can find a steady-state solution without fluxes for
$c_{A}(x)=\hat{c}_{A}e^{-\beta V_{A}(x)}$, which determines also the solution
for $c_{B}(x)=\hat{c}_{A}\frac{\hat{k}_{AB}}{\hat{k}_{BA}}e^{-\beta
V_{B}(x)}$, where $\hat{c}_{A}$ can be found through the normalization. We
find that the chemical and the diffusive part of Eq. 44 are solved by the same
probabilities and hence there are no diffusive fluxes for two species even
when the system is out of equilibrium.

##  Appendix C Calculating the flux-less case for a tree system

For a chemical system without cycles, every two concentrations of neighbouring
nodes $\sigma$ and $\sigma^{\prime}$ are connected via the detailed balance:

$c_{\sigma}^{\prime}(x)=\frac{k_{\sigma\sigma^{\prime}}(x)}{k_{\sigma^{\prime}\sigma}(x)}c_{\sigma}(x)$
(48)

Therefore, any two concentrations for states $\sigma$ and $\sigma^{\prime}$
will be connected as

$c_{\sigma}^{\prime}(x)=K_{\sigma^{\prime}\sigma}(x)c_{\sigma}(x)$ (49)

where $K_{\sigma^{\prime}\sigma}(x)$ is the product of the rate ratios for all
nodes along the path between $\sigma$ and $\sigma^{\prime}$.

Substituting Eq. 49 into Eq. 22 at steady state, we find that the chemical
part is always satisfied as detailed balance was obeyed all along the path.
Therefore, we are solving

$0=\nabla\Big{[}\nabla c_{\sigma}^{\prime}(x)+\beta
c_{\sigma}^{\prime}(x)\nabla F_{\sigma^{\prime}}(x))\Big{]}$ (50)

When we substitute Eq. 49, the first term becomes two:

$\nabla K_{\sigma^{\prime}\sigma}(x)\cdot
c_{\sigma}(x)+K_{\sigma^{\prime}\sigma}(x)\cdot\nabla c_{\sigma}(x)$ (51)

According to Eq. 21, the structure of $K_{\sigma^{\prime}\sigma}(x)$ is the
product of constant ratios of the relevant intrinsic transition rates and of
the exponentials of the corresponding free energies. If the same state appears
both as a final state for one reaction and an initial state for another, its
energy will enter in the exponentials with opposite signs. Therefore, only
$e^{-\beta(F_{\sigma^{\prime}}(x)-F_{\sigma}(x))}$ will survive. Substituting
this into the first term of Eq. 51, we recover $\beta
K_{\sigma^{\prime}\sigma}(x)c_{\sigma}(x)\nabla[F_{\sigma}(x,\sigma)-F_{\sigma^{\prime}}(x)]$.
Substituting this back into Eq. 50, we find that the resulting equation is
identical to the equation for $c_{\sigma}(x)$. Therefore, for a tree-like
reaction structure, the solution never exhibits fluxes as a consequence of the
detailed balance for every chemical reaction.

##  Appendix D Calculating the flux-less case for a system with cycles

For a reaction network structure with cycles, we can use the Kirchhoff
approach to express the concentration $c_{\sigma}(x)$ through the products of
all rates along all edges $l$ leading to this state, summed over all possible
spanning trees including it ($\mathcal{T}$) [16].

$c_{\sigma}(x)=\frac{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}k_{l}^{\sigma}(x)}{\mathcal{N}(x)}$

We can relate this concentration to another concentration, for example
$c_{1}(x)$ via $c_{\sigma}(x)=\frac{c_{\sigma}(x)}{c_{1}(x)}c_{1}(x)$. Then,

$\begin{split}c_{\sigma}(x)&=c_{1}(x)\frac{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}k_{l}^{\sigma}(x)}{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}k_{l}^{1}(x)}\\\
&=c_{1}(x)\frac{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}\frac{k_{l}^{\sigma}(x)}{k_{l}^{\sigma}(x)}\prod_{l\in\mathcal{T}}k_{l}^{1}(x)}{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}k_{l}^{1}(x)}\end{split}$
(52)

Note that every path can be divided into the path
$\gamma^{\mathcal{T}}_{\sigma 1}$ that connects $\sigma$ and $1$ and
sidebranches that are the same whether one is looking at all edges leading to
$\sigma$ or to $1$. Therefore, for the ratio of rates we can replace the
product over $l\in\mathcal{T}$ by one over $l\in\gamma^{\mathcal{T}}_{\sigma
1}$. When substituting the Kramer form of $k_{l}^{\sigma}(x)$ we find
analogously to before, $\prod_{l\in\gamma^{\mathcal{T}}_{\sigma
1}}\frac{k_{l}^{\sigma}(x)}{k_{l}^{1}(x)}=e^{-\beta(F_{\sigma}(x)-F_{1}(x))}e^{-\beta\mu^{l}_{\sigma
1}}$, where the free energy difference is independent of the taken path and
$\mu^{l}_{\sigma 1}$ is the total fuel chemical energy along the path.
Therefore, we obtain

$\begin{split}c_{\sigma}(x)=c_{1}(x)e^{-\beta(F_{\sigma}(x)-F_{1}(x))}S_{1\sigma}(x)\end{split}$
(53)

where we identified

$S_{1\sigma}(x)=\frac{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\gamma^{\mathcal{T}}_{\sigma
1}}}e^{-\beta\mu^{l}_{\sigma
1}}\prod_{l\in\mathcal{T}}k_{l}^{1}(x)}{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}k_{l}^{1}(x)}$
(54)

Substituting the Kramer form of $k_{l}^{1}(x)$, we note how every node except
$1$ will appear as an initial state in $\prod_{l\in\mathcal{T}}k_{l}^{1}(x)$.
Therefore, if we write
$F_{\sigma^{\prime}}(x)=F^{b}_{\sigma^{\prime}}(x)+F^{init}_{\sigma^{\prime}}(x)$,
where $F^{b}_{\sigma^{\prime}}(x)$ is the energy of the barrier in this
direction along the edge and $F^{init}_{\sigma^{\prime}}(x)$ is the energy of
the initial state, then the product of all $F^{init}_{\sigma^{\prime}}(x)$
will cancel and only the barrier energies will remain in Eq. 54:

$S_{1\sigma}(x)=\frac{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\gamma^{\mathcal{T}}_{\sigma
1}}}e^{-\beta\mu^{l}_{\sigma 1}}\prod_{l\in\mathcal{T}}e^{-\beta
F^{b}_{l}(x)}e^{-\beta\mu^{l}_{1}}}{\sum_{\\{\mathcal{T}\\}}{\prod_{l\in\mathcal{T}}}e^{-\beta
F^{b}_{l}(x)}e^{-\beta\mu^{l}_{1}}}$ (55)

If the system is in equilibrium and all $\mu=0$, $S_{1\sigma}(x)$ reduces to
one and detailed balance is obeyed. Note that if the energies of the barriers,
$F^{b}_{l}(x)$ are constant, even if offset from one another by a constant
factor, $S_{1\sigma}(x)$ will be a constant, affecting the distribution, but
resulting in a quasi-equilibrium state.
Now, plugging Eq. 53 into Eq. 22, we obtain

$\begin{split}\dot{c}_{\sigma}(x)=&\nabla[c_{1}(x)e^{-\beta(F_{\sigma}(x)-F_{1}(x))}S_{1\sigma}(x)]\\\
&+c_{1}e^{-\beta(F_{\sigma}(x)-F_{1}(x))}S_{1\sigma}(x)\nabla
F_{\sigma}(x)\end{split}$ (56)

which transforms into

$\begin{split}\dot{c}_{\sigma}(x)=&S_{1\sigma}(x)e^{-\beta(F_{\sigma}(x)-F_{1}(x))}\Big{[}\nabla
c_{1}(x)+\beta c_{1}(x)\nabla F_{1}(x)\Big{]}\\\
&+c_{1}(x)e^{-\beta(F_{\sigma}(x)-F_{1}(x))}\nabla S_{1\sigma}(x)\end{split}$
(57)

The term in the square brackets is $\dot{c}_{1}(x)$, therefore
$\dot{c}_{1}(x)=0$ implies $\dot{c}_{\sigma}(x)=0$ only if $\nabla
S_{1\sigma}(x)=0$.

##  Appendix E Coarse-graining three states into two states with two pathways

Consider the chemical cycle in Fig. 1 in the case when C is a state that is so
short-lived, so that its chemical reactions are much faster than its diffusion
and it is always in steady state with respect to A and B. In that case, we can
expess its concentration through the concentrations of the other two species,

$c_{C}(x)=\frac{k_{AC}}{k_{CA}+k_{CB}}c_{A}(x)+\frac{k_{BC}}{k_{CA}+k_{CB}}c_{B}(x)$
(58)

We can substitute this solution in the equations for $c_{A}(x)$ and
$c_{B}(x)$, which are not individually equilibrated:

$\begin{split}\partial_{t}c_{A}&=-\partial_{t}c_{B}\\\
&=-(k_{AB}+\frac{k_{AC}k_{CB}}{k_{CA}+k_{CB}})c_{A}+(k_{BA}+\frac{k_{BC}k_{CA}}{k_{CA}+k_{CB}})c_{B}\end{split}$
(59)

We can identify two pathways between A and B, where one corresponds to the
rates $k_{AB}$ and $k_{BA}$, and the other – to the rates

$\begin{split}&\tilde{k}_{AB}=\frac{k_{AC}k_{CB}}{k_{CA}+k_{CB}}\\\
&\tilde{k}_{BA}=\frac{k_{BC}k_{CA}}{k_{CA}+k_{CB}}\end{split}$ (60)

Substituting the definitions of the rates, we find

$\begin{split}\tilde{k}_{AB}&=\frac{e^{-\beta(V_{AC}-V_{A})}e^{-\beta(V_{BC}-V_{C})}}{e^{-\beta(V_{AC}-V_{C})}}+e^{-\beta(V_{BC}-V_{C})}e^{\beta\mu_{A\rightarrow
C}}\\\ &=e^{\beta
V_{A}}\frac{e^{-\beta(V_{AC}+V_{CB}-V_{C})}}{e^{-\beta(V_{AC}-V_{C})}}+e^{-\beta(V_{BC}-V_{C})}e^{\beta\mu_{A\rightarrow
C}}\end{split}$ (61)

where we can define the fraction in the last line to be
$\exp(-\beta\tilde{V}_{AB})$, the energy barrier of the second pathway. As
this energy in general is not equal to $V_{AB}$, we indeed have two pathways
that depend differently on the parameters of the system.
Note that this coarse-graining does not alter the affinity of the cycle:

$\frac{k_{AB}\tilde{k}_{BA}}{k_{BA}\tilde{k}_{AB}}=\frac{k_{AB}k_{BC}k_{CA}}{k_{BA}k_{AC}k_{CB}}$
(62)

Using this mapping, equation 22 can be solved analytically for a step
potential (analogously to App. G):

$\begin{split}&J_{A}=-J_{B}\\\ &=\begin{cases}\begin{aligned}
\frac{e^{e_{R}+\epsilon_{L}}-e^{e_{L}+\epsilon_{R}}}{\mathcal{N}}(e^{\beta\tilde{\mu}}-1)\sinh(q_{R}&L)\sinh(q_{L}(L+x))\\\
&\mathrm{if}~{}-L\leq x\leq 0\end{aligned}\\\ \begin{aligned}
\frac{-e^{e_{R}+\epsilon_{L}}+e^{e_{L}+\epsilon_{R}}}{\mathcal{N}}(e^{\beta\tilde{\mu}}-1)\sinh(q&{}_{L}L)\sinh(q_{R}(x-L))\\\
&\mathrm{if}~{}0\leq x\leq L\end{aligned}\end{cases}\end{split}$ (63)

where $e_{L/R}$ and $\epsilon_{L/R}$ are the energy barriers along the two
different chemical pathways and $\mathcal{N}$ is the overall normalisation
factor. The factors $q_{L}$ and $q_{R}$ are the square roots of the non-zero
eigenvalues of the chemical transition matrix on the left and right side of
the interface, respectively. The fluxes disappear of course in the absence of
non-equilibrium driving ($\tilde{\mu}=0$), but also when the energy barriers
of the two pathways satisfy the condition
$e_{R}-e_{L}=\epsilon_{R}-\epsilon_{L}$, which corresponds to barriers that do
not depend on space or that do, but in the same way for both pathways. This
setting could actually be re-mapped into a single pathway between A and B,
which, as shown in Section III, does not create diffusive fluxes even away
from equilibrium (see Appendix B).

The spatial decay of the fluxes on each side of the interface is modulated
only by the corresponding non-zero eigenvalue of the chemical transition
matrix, which is (where $i\in\\{L,R\\}$):

$\begin{split}q_{i}^{2}=~{}&k_{AB}^{0}e^{-\beta e_{i}+\beta
V_{A,i}+\beta\tilde{\mu}}+k_{BA}^{0}e^{-\beta e_{i}+\beta V_{B,i}}+\\\
&k_{AB}^{0}e^{-\beta\epsilon_{i}+\beta
V_{A,i}}+k_{BA}^{0}e^{-\beta\epsilon_{i}+\beta V_{B,i}}\end{split}$ (64)

This result provides the reaction-diffusion length

$l^{rd}_{L/R}=\sqrt{\frac{D}{q_{L/R}}}$ (65)

For eigenvalues small enough that the argument of the hyperbolic sinus is
small over the whole intervals $x\in[-L,0]$ and $x\in[0,L]$, the fluxes can be
Taylor expanded to the first order, consistent with the observed linear
behaviour.

##  Appendix F Numerical solution method

In order to solve Eq. 22 numerically, we convert discretize space and solve
the equation in every point. To do so, we construct a local stochastic matrix,
the entries of which include both terms coming from the discretized diffusion
and terms accounting for the chemical transformations. We then solve
$\vec{0}=\mathbf{M}(x)\vec{\rho}(x)$, where
$\vec{\rho}(x)=(\rho_{A}(1),\rho_{A}(2)...\rho_{A}(n),\rho_{B}(1),...\rho_{C}(n))$
for a system discretized to $n$ spatial points. The steady-state solution of
this equation is the non-zero eigenvector for the zero eigenvalue of the
resulting stochastic matrix. We obtain it by dividing $\mathbf{M}$ into a
reduced matrix $\mathbf{M}_{red}$ which lacks the first row and first column
of $\mathbf{M}$ and into the vector $\vec{b}$, which is the first column of
$\mathbf{M}$ without the first entry. We can choose the first entry of $\rho$,
corresponding to the first (redundant) row of $\mathbf{M}$, freely and we
choose it to be 1, so that $\vec{\rho}=(1,0,...,0)^{T}+\vec{\rho^{\prime}}$.
Then, we find that
$0=\mathbf{M}\rho=\vec{b}(1,0,...,0)^{T}+\mathbf{M}_{red}\vec{\rho^{\prime}}$
and therefore $\vec{\rho^{\prime}}=-\mathbf{M}_{red}^{-1}\vec{b}$. The final
result is renormalized so that
$\sum_{\sigma\in\\{A,B,C\\}}\sum_{x}\rho_{\sigma}(x)=1$.
In the numerical solution, we can choose the form of the background potential
to have an arbitrary functional form as long as it does not vary too much
within one space unit.

##  Appendix G Boundary conditions for step potential

On each side, the system of equations for the three species can be represented
by a stochastic matrix and written in the basis given by its eigenvectors
$\vec{v}_{i}$. The concentration then can be written in the eigenvector basis,
denoted as $\tilde{c}(x)$. By construction, the matrix has one zero
eigenvalue, $\lambda_{0}$. For it, we are solving
$\partial_{x}\tilde{c}_{0}(x)=0$, while for the other two eigenvalues we solve
$\partial_{x}^{2}\tilde{c}_{i}(x)=\lambda_{i}(x)\tilde{c}_{i}(x)$, where
$\lambda_{i}$ are the remaining two eigenvalues. In total, there are twelve
free parameters and twelve boundary conditions. Of those, nine are coming
from:

  * •

the conservation of the total particle number (1 condition)

  * •

no fluxes through the boundaries (at $x=-L$ and $x=L$) (6 conditions)

  * •

the continuity of the fluxes at the interface (2 conditions)

The last three boundary conditions come from the behaviour of the particles
immediately at the interface [see also 26, 27]. For a band around the
interface with a width $\epsilon\rightarrow 0$, the contribution of chemical
reactions is negligible and the dynamics of the particles is dominated by
diffusion. The probability of a particle crossing the interface is therefore
proportional to the ratio of the Boltzmann factors on either side.
Intuitively, it can be seen as the rule to be obeyed by particles moving
across the interface without undergoing chemical transitions.
The exact condition can be derived by integrating Eq. 22 around the interface.
One finds:

$\begin{split}0&=D\int_{-\epsilon}^{\epsilon}\Big{(}\partial_{x}c_{\alpha}(x)+\beta
c_{\alpha}(x)\partial_{x}V_{\alpha}(x)\Big{)}dx\\\
&=D\int_{-\epsilon}^{\epsilon}e^{-\beta
V_{\alpha}(x)}\partial_{x}(c_{\alpha}(x)e^{\beta V_{\alpha}(x)})dx\end{split}$
(66)

To fulfil this equation, $c_{\alpha}(x)$ has to have the form
$c_{\alpha}(x)=C~{}\exp(-\beta V_{\alpha}(x))$ with a constant $C$. It follows
that

$\frac{c_{\alpha}(0+)}{c_{\alpha}(0-)}=\frac{\exp(-\beta
V_{\alpha}(0+))}{\exp(-\beta V_{\alpha}(0-))}$ (67)

which provides the last three boundary conditions.

##  Appendix H Spherical coordinates

A condensate inside the cell is better described as a sphere in 3D. The fact
that we are describing it in 1D does not alter the results qualitatively, as
we show here for the solution of the problem for step potentials. We consider
a spherical condensate centered around the origin with its interface at the
radius R. The system is invariant under rotation, so that
$c(r,\theta,\phi,\sigma)=c(r,\sigma)$. Then, we can solve the vectorial form
of Eq. 22

$\partial_{t}\vec{c}(r)=0=D\vec{\nabla}(\vec{\nabla}\vec{c}(r)+\beta\vec{\nabla}\mathbf{V}(r)~{}\vec{c}(r))+\mathbf{K}(r)\vec{c}(r)$
(68)

where now the divergence is in spherical coordinates. We proceed as in Section
IV, changing the basis of Eq. 68 to the one spanned by the eigenvectors of the
chemical transition matrix $\mathbf{K}(r)$ and solving the equation on the
left and on the right of the interface, where the drift term is equal to zero.
For the non-zero eigenvalues, the left-hand side becomes
$\frac{1}{r^{2}}\partial_{r}(r^{2}\partial_{r}\tilde{c}(r))$, so we look for a
solution of the form of $\tilde{c}(r)=f(r)/r$, for which we retrieve the same
solutions as in the 1D case. Collecting all terms, the concentration profile
inside the interface is described by

$\begin{split}c^{in}(r)&=A^{in}_{0}\vec{v}^{in}_{0}+\left(\frac{A^{in}_{1+}}{r}e^{\lambda^{in}_{1}r}+\frac{A^{in}_{1-}}{r}e^{-\lambda^{in}_{1}r}\right)\vec{v}^{in}_{1}\\\
&+\left(\frac{A^{in}_{2+}}{r}e^{\lambda^{in}_{2}r}+\frac{A^{in}_{2-}}{r}e^{-\lambda^{in}_{2}r}\right)\vec{v}^{in}_{2}\end{split}$
(69)

and analogously outside of the interface. The boundary conditions will be
analogous to the ones for 1D, with the difference that the no-flux conditions
on the boundaries now transform into a no-flux condition at the outer boundary
and a no-sources condition at the origin. Taking this into account, the
solution can be determined as before. The result can be seen in Fig. S1. The
fluxes still are the largest at the interface and the 1D results remain
qualitatively unaltered.

Figure S1: Fluxes calculated for the step potential in 3D with the same
parameters as in 1D. The system length is 10, the potential changes at $R=5$.