# Thermodynamic origins of topological protection in nonequilibrium stochastic
systems

Pankaj Mehta<EMAIL_ADDRESS>Department of Physics, Boston University, Boston,
MA 02215, USA  Biological Design Center, Boston University, Boston, MA 02215,
USA  Faculty of Computing and Data Science, Boston University, Boston, MA
02215, USA  Jason W. Rocks  Department of Physics, Boston University, Boston,
MA 02215, USA

###### Abstract

Topological protection has emerged as an organizing principle for
understanding and engineering robust collective behavior in electronic and
material systems. Recent work suggests that topology may also play a role in
organizing stochastic processes relevant to biology and self-assembly. Here,
we show that topological protection in chemical networks can be understood
entirely in terms of nonequilibrium thermodynamics. We illustrate these ideas
using simple examples inspired by the literature.

Biochemical processes in living systems operate far from equilibrium. Such
nonequilibrium stochastic processes play vital roles in diverse biological
settings ranging from information processing to the self-assembly of dynamic
structures such as microtubules Beard and Qian (2008); Qian and Beard (2005);
Qian (2007); Hill (2012); Lan  _et al._ (2012); Mehta and Schwab (2012);
Murugan  _et al._ (2012); Lang  _et al._ (2014); Mehta  _et al._ (2016);
Ouldridge  _et al._ (2017); Hess and Ross (2017); McInally  _et al._ (2021);
Das  _et al._ (2021). In all of these phenomena, cells transduce free energy
by coupling chemical reactions to high-energy molecules (e.g., ATP and ADP)
that are maintained out-of-equilibrium. Outside of biology, similarly driven
stochastic chemical processes play a central role in material self-assembly
and the design of active materials Zeravcic  _et al._ (2014); Murugan  _et
al._ (2015); Nguyen and Vaikuntanathan (2016); Bisker and England (2018). In
both natural and synthetic systems, it is crucial that the behavior and
function of these systems be robust to the presence of disorder and thermal
fluctuations.

Recently, topological protection has emerged as an important design principle
for understanding and engineering robust behavior in out-of-equilibrium
stochastic systems. Inspired by the immense progress made on this topic,
several authors have translated many ideas originally developed in the context
of quantum mechanics and discrete meta-materials to stochastic networks
Murugan and Vaikuntanathan (2017); Dasbiswas  _et al._ (2018); Knebel  _et
al._ (2020, 2020); Tang  _et al._ (2021); Yoshida  _et al._ (2021). On a
technical level, many of these works exploit the formal analogy between the
Master Equation and the Schröndiger equation in which the state-to-state
transition matrix of the former plays an analogous role to the Hamiltonian of
the latter Van Kampen (1992). As in quantum systems, the existence of a
topological phase transition can be related to the closing of a gap in the
spectrum (now of the transition matrix instead of the Hamiltonian) and an
associated change in the value of a topological invariant, usually a
$\mathbb{Z}_{2}$-valued Zaks phase Zak (1989); Cohen  _et al._ (2019).

Despite the considerable insights achieved by these works, the approaches they
employ make it difficult to understand the physical origin of topological
protection. The use of transition matrices (or closely related and even more
mathematically complex objects such as tilted current matrices) seems to
suggest that topological protection is a property of the _kinetics_ of the
system under consideration Murugan and Vaikuntanathan (2017); Dasbiswas  _et
al._ (2018); Tang  _et al._ (2021). Furthermore, the relationship between
topological protection and nonequilibrium properties such as fluxes and
thermodynamic driving forces remains obscure.

Here, we revisit topological protection in chemical systems and show that this
phenomenon can be understood in terms of the theory of nonequilibrium
thermodynamics of steady-states (NESS) in open systems Qian and Beard (2005);
Qian (2007); Beard and Qian (2008). By reformulating the physical description
of these systems in terms of stoichiometric matrices and driving forces, we
show that the origin of topological protection in these systems is primarily
_thermodynamic_ rather than kinetic. Our formalism allows us to directly link
the existence of topologically protected edge modes and edge currents to
driving forces across elementary reactions that make up the system. We
demonstrate these ideas by revisiting simple examples inspired by earlier
works Knebel  _et al._ (2020); Tang  _et al._ (2021).

To make a closer connection to physical kinetics, we focus primarily on open
chemical systems (e.g., cells) where a subset of chemical species (e.g., ATP
and ADP) are maintained out-of-equilibrium and used to drive chemical
reactions. However, most of these ideas naturally transfer to a more general
setting based on nonequilibrium steady-states described by Master Equations Ge
and Qian (2010).

To begin, we consider a simple chemical network composed of chemical species
(indexed by $i$) with average concentrations $x_{i}$ that can participate in
_reversible_ reactions (indexed by $\alpha$). Reversibility of the reactions
is necessary for thermodynamic consistency. The dynamics of such system are
described by the kinetic equations

${dc_{i}\over dt}=\sum_{\alpha}S_{i\alpha}J_{\alpha}+b^{\mathrm{ext}}_{i},$
(1)

where $J_{\alpha}$ is the net flux of reaction $\alpha$,
$b^{\mathrm{ext}}_{i}$ is the rate at which species $i$ is produced and
degraded by external sources or sinks (we will generally assume all
$b^{\mathrm{ext}}_{i}=0$), and $S_{i\alpha}$ is the stoichiometric matrix
which encodes the number of each species $i$ consumed (negative entries) or
produced (positive entries) by reaction $\alpha$ Qian and Beard (2005); Qian
(2007); Beard and Qian (2008). Since all reactions are reversible, we further
divide the net flux into a forward flux $J_{\alpha}^{+}$ and a backward flux
$J_{\alpha}^{-}$ with ${J_{\alpha}=J_{\alpha}^{+}-J_{\alpha}^{-}}$. These two
fluxes can then be related to the concentrations via the law of mass action,

$J^{\pm}_{\alpha}=k^{\pm}_{\alpha}\prod_{i,S_{i\alpha}\lessgtr
0}c_{i}^{\absolutevalue{S_{i\alpha}}},$ (2)

where $k^{\pm}_{\alpha}$ are the forward and backward rate constants of
reaction $\alpha$.

As an example, consider a reaction $\gamma$ of the form

$2X_{l}\xrightleftharpoons[k_{\gamma}^{-}]{k_{\gamma}^{+}}3X_{m}+X_{n}$

that consumes two molecules of $X_{l}$ and produces three molecules of $X_{m}$
and one molecule of $X_{n}$. The associated entries in the stoichiometric
matrix are ${S_{l\gamma}=-2}$, ${S_{m\gamma}=3}$, ${S_{n\gamma}=1}$ or zero
otherwise with fluxes ${J_{\gamma}^{+}=k_{\gamma}^{+}c_{l}^{2}}$ and
${J_{\gamma}^{-}=k_{\gamma}^{-}c_{m}^{3}c_{n}}$.

We now exploit some beautiful results from the theory of NESS to reexpress
these kinetic equations in terms of thermodynamic driving forces. To begin, we
write the ratio of the forward and backward fluxes, $J_{\alpha}^{+}$ and
$J_{\alpha}^{-}$, in terms of the thermodynamic driving force
$\Delta\mu_{\alpha}$ (chemical potential difference) using the NESS identity

${J_{\alpha}^{-}\over J_{\alpha}^{+}}=e^{{\Delta\mu_{\alpha}\over k_{B}T}},$
(3)

where $k_{B}$ is Boltzmann’s constant and $T$ is the temperature (see SI and
Refs. Qian and Beard, 2005; Qian, 2007; Beard and Qian, 2008).

Next, we decompose $\Delta\mu_{\alpha}$ into three distinct terms:

$\Delta\mu_{\alpha}=\Delta\mu_{\alpha}^{0}+\Delta\mu^{\mathrm{drive}}_{\alpha}+k_{B}T\sum_{i}S_{i\alpha}\ln{c_{i}},$
(4)

with $\Delta\mu_{\alpha}^{0}=-k_{B}T\sum_{i}S_{i\alpha}\ln{c_{i}^{eq}}$ the
chemical potential difference at equilibrium (in the absence of coupling to
any high-energy molecules such as ATP/ADP, GTP/GDP, etc.) and
$\Delta\mu^{\mathrm{drive}}_{\alpha}$ the nonequilibrium free energy
associated with driving reaction $\alpha$ by coupling to high-energy molecular
processes such as phosphorylation. The last term is the usual concentration-
dependent contribution to the chemical potential difference. In kinetic
descriptions, coupling to the high-energy molecules
$\Delta\mu^{\mathrm{drive}}_{\alpha}$ is often included _implicitly_ and for
this reason, we can identify the first two terms in the driving force with the
kinetic parameters $k_{\alpha}^{\pm}$ in Eq. (2) for the law of mass action
(see SI and Ref. Qian, 2007),

$\Delta\mu_{\alpha}^{0}+\Delta\mu^{\mathrm{drive}}_{\alpha}=k_{B}T\ln{k_{\alpha}^{-}\over
k_{\alpha}^{+}}.$ (5)

Using the NESS identity of Eq. (3), we then write the net flux
${J_{\alpha}=J_{\alpha}^{+}-J_{\alpha}^{-}}$ in terms of only the backward
fluxes $J_{\alpha}^{-}$ and driving force $\Delta\mu_{\alpha}$ to get

$J_{\alpha}=(e^{-{\Delta\mu_{\alpha}\over
k_{B}T}}-1)J_{\alpha}^{-}\equiv\theta_{\alpha}J_{\alpha}^{-},$ (6)

where the second equality defines the parameter $\theta_{\alpha}$. This
parameter, which we call the _thermodynamic drive_ , will play a central role
in what follows, and for this reason, it is worth gaining some basic intuition
for its meaning. When $\theta_{\alpha}>0$ the net reaction flux is in the
forward direction. In contrast, when $\theta_{\alpha}<0$ the net flux is in
the negative direction. When $\theta_{\alpha}=0$, there is no net current
across reaction $\alpha$ and the reaction is locally at detailed balance.
Consequently, all $\theta_{\alpha}$ are identically zero for a system at
equilibrium. From these considerations, it is clear that $\theta_{\alpha}$ is
a measure of the nonequilibrium drive associated with reaction $\alpha$.

One subtlety worth emphasizing is that $\theta_{\alpha}$ depends not only on
the kinetic parameters via the driving potentials
$\Delta\mu^{\mathrm{drive}}_{\alpha}$, but also on the steady-state
concentrations $c_{i}$ [see Eqs. (4) and (6)]. For this reason, we cannot
directly control $\theta_{\alpha}$ in our systems by changing kinetic
parameters. However, it is still possible to choose
$\Delta\mu^{\mathrm{drive}}_{\alpha}$ to prejudice $\theta_{\alpha}$ towards a
particular direction. Physically, this corresponds to changing the amount of
free energy a reaction can transduce when it is coupled to high-energy
processes such as phosphorylation Hill (2012); Qian (2007). Finally, we note
that the choice of expressing Eq. (6) in terms of $J^{-}_{\alpha}$ or
$J^{+}_{\alpha}$ is arbitrary (i.e. a choice of gauge). Although each choice
results in a slightly different definition of $\theta_{\alpha}$, it will
follow the same sign convention in either case, so this choice does not affect
our analysis.

With this basic picture in mind, we write the steady-state condition for Eq.
(1) entirely in terms of negative fluxes,

$0=\sum_{\alpha}S_{i\alpha}^{-}J_{\alpha}^{-},$ (7)

where we have set ${b^{\mathrm{ext}}_{i}=0}$ and defined the drive-dependent
stoichiometric matrix ${S_{i\alpha}^{-}\equiv S_{i\alpha}\theta_{\alpha}}$. We
show below that $S_{i\alpha}^{-}$ encodes all the information necessary to
identify topological phase transitions. To motivate this, we note that
$S_{i\alpha}^{-}$ plays a similar role to the compatibility matrix in
topological mechanics, but for a system with directed interactions. We also
emphasize that $S_{i\alpha}^{-}$ is a purely thermodynamic quantity; it only
depends on the topology of the chemical network and the thermodynamic drives
of each reaction. Effectively, Eq. (7) is analogous to Kirchhoff’s node rule
for chemical networks expressed purely in terms of backward fluxes.

Figure 1: Topologically-protected edge modes in a network of chemical
reactions on a triangular lattice. (a) The primitive cell is composed of two
types of chemical species and three types of reversible reactions.
Forward/backward directions correspond to large/small arrows. (b) The full
chemical network is formed by repeating the basic primitive cell. (c) The
driving chemical potentials
$\Delta\mu^{0}_{\alpha}+\Delta\mu^{\mathrm{drive}}_{\alpha}$ for two different
systems with $N=10$ lattice cells. The size and direction of the arrows
indicate the magnitude and direction of chemical potential drops across each
edge. The kinetic parameters are $k_{1}^{\pm}=k_{2}^{+}=k_{3}^{+}=1$,
$k_{2}^{-}=2$, and $k_{3}^{-}=1$ ($0.25$) for the upper (lower) configuration,
corresponding to $\Delta\mu^{0}_{1}+\Delta\mu^{\mathrm{drive}}_{1}=0$,
$\Delta\mu^{0}_{2}+\Delta\mu^{\mathrm{drive}}_{2}=k_{B}T\log 2$, and
$\Delta\mu^{0}_{3}+\Delta\mu^{\mathrm{drive}}_{3}=0$ ($-2k_{B}T\log 2$). (d)
Numerical simulations of the systems in (c) showing the existence of left- and
right-localized edge modes. Magnitudes of net nonequilibrium fluxes
$J_{\alpha}$ are indicated by thickness of arrows and steady-state
concentrations by sizes of nodes. (e) Thermodynamic drives $\theta_{\alpha}$
of the two systems indicated by edge color with $\theta_{\alpha}>0$ ($<0$)
corresponding to a drive in the forward (backward) direction.

To demonstrate topological protection within this framework, we start by
analyzing a simple triangular lattice network inspired by a recent paper on
topological phase transitions in a game-theoretic model of a coupled rock-
paper-scissor system Knebel  _et al._ (2020). As shown in Fig. 1(a), each unit
cell consists of three chemical reactions (edges) that interconvert two
different chemical species (nodes). Repeating this unit $N$ times [Fig. 1(b)]
and closing the system with an additional boundary node forms a chain of
$2N+1$ chemical species coupled by $3N$ reactions. All reactions conserve
particle number with each reaction $\alpha$ between species $i$ and $j$ simply
converting one molecule of species $X_{i}$ to one molecule of species $X_{j}$
and vice versa,

$X_{i}\xrightleftharpoons[k_{\alpha}^{-}]{k_{\alpha}^{+}}X_{j}.$

Furthermore, we assume that each unit cell has the same three types of
reactions, $A=1,2,3$ with kinetic parameters $k_{A}^{\pm}$ and thermodynamic
drive $\theta_{A}$ [Eqs. (4) and (6)].

Using standard arguments, we write the full stoichiometric matrix in terms of
the $2$-by-$3$ stoichiometric matrix of the primitive lattice cell in Fourier
space (see SI),

$\hat{S}(k)=\matrixquantity[1-e^{ik}&-1&e^{ik}\\\ 0&1&-1],$ (8)

where rows correspond to chemical species and columns to fluxes. We interpret
$e^{ik}$ as a “shift operator” that translates the flux by a unit cell to the
right, allowing for a simple physical interpretation of the elements of
$\hat{S}(k)$ as subtracting and adding fluxes in adjacent unit cells. The
resulting Fourier transform of $S_{i\alpha}^{-}$ [Eq. (7)] is simply the
$2\times 3$ matrix

$\hat{S}^{-}(k)=\matrixquantity[(1-e^{ik})\theta_{1}&-\theta_{2}&e^{ik}\theta_{3}\\\
0&\theta_{2}&-\theta_{3}],$ (9)

where once again rows correspond to chemical species, but columns now
correspond to backward fluxes ($J_{\alpha}^{-}$).

To identify topological phase transitions, we analyze the spectrum of
$\hat{S}^{-}(k)$ or equivalently, the Hermitian matrix
${\hat{H}(k)=\hat{S}^{-}(k)\hat{S}^{-}(k)^{\dagger}}$ whose eigenvalues are
the squares of the singular values of $\hat{S}^{-}(k)$,

$\displaystyle\hat{H}(k)$ $\displaystyle=$
$\displaystyle\matrixquantity[2(1-\cos{k})\theta_{1}^{2}+\theta_{2}^{2}+\theta_{3}^{2}&\,-\theta_{2}^{2}-\theta_{3}^{2}e^{ik}\\\
-\theta_{2}^{2}-\theta_{3}^{2}e^{-ik}&\,\theta_{2}^{2}+\theta_{3}^{2}].$ (10)

We also further restrict our analytics to the special case $\theta_{1}=0$
where reactions of type $A=1$ are in detailed balance (attainable in an
infinite system by setting
${\Delta\mu^{0}_{1}+\Delta\mu^{\mathrm{drive}}_{1}=0}$, see SI). In this case,
we can write ${\hat{H}(k)=a_{0}I+\vec{a}\vec{\sigma}}$ where $I$ is the
identity matrix and $\vec{\sigma}$ is the set of Pauli matrices with
${a_{0}=\theta_{2}^{2}+\theta_{3}^{2}}$ and
${\vec{a}=(-\theta_{2}^{2}-\theta_{3}^{2}\cos{k},-\theta_{3}^{2}\sin{k},0)}$
(see SI). In terms of $\vec{a}$ and $a_{0}$, the spectrum is

$\lambda_{\pm}=a_{0}\pm\sqrt{\vec{a}\cdot\vec{a}}=a_{0}\pm\sqrt{\theta_{2}^{4}+\theta_{3}^{4}+2\theta_{2}^{2}\theta_{3}^{2}\cos{k}},$
(11)

which becomes degenerate at $k=\pi$ when $\theta_{2}^{2}=\theta_{3}^{2}$.

To see that this corresponds to a topological phase transition, we follow
standard arguments and plot the $x$\- and $y$-components of $\vec{a}$ as $k$
varies from $0$ to $2\pi$ (see SI Fig. A1). For all parameters, the resulting
curve is a circle of radius $\theta_{3}^{2}$ centered at the point
$(-\theta_{2}^{2},0)$ in the $x$-$y$ plane. When
$\theta_{3}^{2}<\theta_{2}^{2}$, the circle does not enclose the origin,
whereas if $\theta_{3}^{2}>\theta_{2}^{2}$, this curve always contains the
origin with the transition occurring exactly when
$\theta_{3}^{2}=\theta_{2}^{2}$ and the gap closes. Since the curve
$(a_{x}(k),a_{y}(k))$ can be viewed as a map from $S^{1}$ to $S^{1}$ (recall
$k$ is periodic with period $2\pi$), this corresponds to a change in the
winding number and hence defines a $\mathbb{Z}_{2}$ topological invariant
Nakahara (2018). Alternatively, one can show the existence of a Zak phase by
interpreting $\hat{H}(k)$ as a quantum Hamiltonian and making use of well-
established classifications of topological phases Chiu  _et al._ (2016). As in
the original game theoretic example considered in Ref. Knebel  _et al._ ,
2020, our analytics suggest that the topological transition is controlled by a
single parameter related to the strength of reactions of types 2 and 3.
However, in our system, this transition depends on the kinetic rate constants
only through the ratio of nonequilibrium thermodynamic drives.

To test these predictions, we ran numerical simulations on a triangular
lattice with closed boundary conditions. We numerically integrated the kinetic
equations of Eq. (1) assuming mass action [Eq. (2)] to calculate steady-state
concentrations $c_{i}$ and fluxes $J_{\alpha}$ (see SI for details). Fig.
1(c)-(e) shows results from simulations where we chose
$k_{1}^{\pm}=k_{2}^{+}=k_{3}^{+}=1$, $k_{2}^{-}=2$ and varied a single
parameter, $k_{3}^{-}$. Thermodynamically, this is equivalent to controlling
the thermodynamic driving force of reaction 3,
$\Delta\mu_{3}^{0}+\Delta\mu^{\mathrm{drive}}_{3}$. Using Eq. (6), we
estimated $\theta_{\alpha}$ for each edge and separately computed the mean-
squared average over all reactions of the same type,
$\langle\theta_{A}^{2}\rangle$. We find that the system robustly supports two
types of edge states localized to different ends of the lattice. Consistent
with our analytics, when
$\langle\theta_{3}^{2}\rangle\gg\langle\theta_{2}^{2}\rangle$, the systems
exhibits right localized edge states and when
$\langle\theta_{3}^{2}\rangle\ll\langle\theta_{2}^{2}\rangle$ the system has
left localized states. However, our numerics are inadequate to probe the
transition region since all reactions in our system are reversible resulting
in large finite-size effects.

Figure 2: Topologically-protected edge current in a simple model of biopolymer
self-assembly where the state of the biopolymer is parameterized by the number
of monomers of two different types, along with one of four possible internal
states. (a) Schematic of eight possible elementary reactions where the
internal reactions (1-4) transition between the internal states of the
biopolymer and the external reactions (5-8) add or remove monomers. (b)
Primitive lattice cell with four chemical species corresponding to the four
internal states (nodes) and eight fluxes (edges). Forward/backward directions
correspond to large/small arrowheads. (c) A chemical reaction network with
$N_{x}=N_{y}=3$. (d) Stochastic simulations using the Gillespie algorithm
showing a typical nonequilibrium steady-state flux configurations with a
topologically protected chiral edge current ($k_{a}^{+}=9k_{a}^{-}$) and (e)
without an edge current ($k_{a}^{+}=k_{a}^{-}/3$). Magnitudes of
nonequilibrium fluxes are indicated by thickness of arrows and steady-state
concentrations by sizes of nodes. For both simulations:
$k_{a}^{-}=k_{b}^{-}=0.01$, $k_{b}^{+}=3k^{-}$.

Next, we analyze an elegant toy model for self-assembly of biopolymers shown
Fig. 2, previously introduced in Ref. Tang  _et al._ , 2021. This model was
previously shown to support topological edge currents using an approach based
on the Master Equation and techniques inspired by non-Hermitian quantum
systems Shen  _et al._ (2018); Kawabata  _et al._ (2019); Lieu (2018). In this
model, a biopolymer molecule can expand or contract by adding or losing two
different types of monomers, $X$ and $Y$. At a given time, the biopolymer is
described by an ordered pair of integers $(x,y)$ representing the number of
monomers of type $X$ and $Y$ currently incorporated into the molecule, along
with a parameter $z=A,B,C,D$ that can take on one of four values corresponding
to four different “internal states” of the biopolymer. Depending on the
internal state, the biopolymer can add or remove a monomer of $X$ and $Y$
[reactions 5-8 in Fig. 2(a)]. The biopolymer molecule can also cyclically
switch between internal states through four additional reactions [labeled 1-4
in Fig. 2(a)]. Each of these eight reactions $\alpha$ has two associated
kinetic parameters $\\{k_{\alpha}^{\pm}\\}$ and a thermodynamic driving force
$\theta_{\alpha}$. We also assume that the biopolymer can incorporate at most
$N_{x}$ of the $X$ monomers and $N_{y}$ of the $Y$ monomers so that $0\leq
x\leq N_{x}$ and $0\leq y\leq N_{y}$. As shown in Fig. 2(c), this chemical
reaction network naturally lives on a two-dimensional lattice in which each
node represents an internal state and each edge represents a chemical
reaction. In Ref. Tang  _et al._ , 2021, the analytics were largely confined
to the “chiral” case in which all reactions were strictly irreversible. Here,
we instead consider a more general case where all reactions are reversible.

As shown in Fig. 2(b), the elementary lattice cell for this chemical reaction
network is composed of four chemical species (corresponding to each of the
four possible internal states of the biopolymer at a fixed monomer
composition) and eight elementary reactions that add or remove monomers or
change the internal state. By using a procedure analogous to that used for the
triangular network, it is straightforward to calculate the Fourier transform
for the drive-dependent stoichiometric matrix $\hat{S}^{-}(\vec{k})$, where
$\vec{k}=(k_{x},k_{y})$ includes the momentum in both the horizontal and
vertical directions (see SI). To look for topological phase transitions, we
once again analyze the spectrum of
$\hat{H}(\vec{k})=\hat{S}^{-}(\vec{k})\hat{S}^{-}(\vec{k})^{\dagger}$ and ask
where an eigenvalue gap closes. For simplicity, we restrict our analysis to
the case where the thermodynamic drives associated with each of the four
internal reactions are identical with $\theta_{\alpha}=\theta_{a}$ for
$\alpha=1,2,3,4$ and similarly the external reactions are identical with
$\theta_{\alpha}=\theta_{b}$ for $\alpha=5,6,7,8$ (attainable by separately
setting $\Delta\mu^{0}_{\alpha}+\Delta\mu^{\mathrm{drive}}_{\alpha}$ equal for
the internal and external reactions, see SI). We find that the gap in the
spectrum closes when the internal and external driving forces are equal with
$\theta_{a}^{2}=\theta_{b}^{2}$. This closing of the gap is also associated
with a change in a value of a $\mathbb{Z}_{2}$ topological invariant (i.e., a
nontrivial Zaks phase) (see SI). Furthermore, in analogy with the integer Hall
effect, we expect to see a chiral edge current emerge in one of these phases
(as discussed in Ref. Tang  _et al._ , 2021).

To check this prediction numerically, we performed stochastic simulations
using the Gillespie algorithm Gillespie (1976). In our simulations, we set
$k_{a}^{-}=k_{b}^{-}=0.1$, $k_{b}^{+}=0.27$ and varied the ratio of
$\theta_{a}$ and $\theta_{b}$ by varying $k_{a}^{+}$ (or equivalently, varying
$\Delta\mu_{a}^{0}+\Delta\mu^{\mathrm{drive}}_{a}$). As predicted, we found
our system robustly supports chiral edge currents over a large range of
parameters. Typical nonequilibrium steady-state configurations observed in our
simulations with and without edge currents are shown in Figs. 2(d) and (e).

In conclusion, we have revisited the idea of topological protection in open
chemical systems and shown that the origin of topological protection can be
understood directly in terms of nonequilibrium thermodynamic quantities. In
our analysis, thermodynamic driving forces play a fundamental role: they
control the sign and magnitude of nonequilibrium fluxes, as well as the
stoichiometric matrix which encodes the topology/geometry of interactions.
This is consistent with the well-established result from the theory of
nonequilibrium thermodynamics that steady-states can be understood by
decomposing fluxes into cycles whose magnitude is controlled by chemical
potential differences/theromodynamic forces across the cycle Hill (2012).
Here, we have limited our analysis to some simple systems with just a few
types of reactions. We expect that for more complicated systems (e.g., the
case where all eight reactions in the biopolymer example are distinct) systems
will likely exhibit a much richer phase diagram with multiple topological
phase transitions. Our analysis also raises many new questions. Previous work
explicitly makes use of the non-Hermitian nature of the transition matrices
defining the Markov process Tang  _et al._ (2021). Here, we instead analyzed
the spectrum of a Hermitian matrix $\hat{H}(k)$. Given that non-Hermitian
systems exhibit unique topological properties Shen  _et al._ (2018); Kawabata
_et al._ (2019); Lieu (2018), much more work is needed to reconcile the
thermodynamic picture presented here with one based on kinetics. Doing so will
likely lead to considerable physical insight and illuminate deep mathematical
structures relating thermodynamics, stochastic systems, and topology.

## Acknowledgments

We thank Evelyn Tang for helpful discussions and Wenping Cui for comments on
the manuscript. This work was supported by NIH R35GM119461 to PM.

## References

  * Beard and Qian (2008) D.A. Beard and H. Qian,  _Chemical biophysics: quantitative analysis of cellular systems_ (Cambridge Univ Pr, 2008).
  * Qian and Beard (2005) H. Qian and D.A. Beard, “Thermodynamics of stoichiometric biochemical networks in living systems far from equilibrium,” Biophysical chemistry 114, 213–220 (2005).
  * Qian (2007) Hong Qian, “Phosphorylation energy hypothesis: open chemical systems and their biological functions,” Annu. Rev. Phys. Chem. 58, 113–142 (2007).
  * Hill (2012) Terrell Hill,  _Free energy transduction in biology: the steady-state kinetic and thermodynamic formalism_ (Elsevier, 2012).
  * Lan  _et al._ (2012) G. Lan, P. Sartori, S. Neumann, V. Sourjik, and Y. Tu, “The energy-speed-accuracy trade-off in sensory adaptation,” Nature Physics 8, 422–428 (2012).
  * Mehta and Schwab (2012) Pankaj Mehta and David J Schwab, “Energetic costs of cellular computation,” Proceedings of the National Academy of Sciences 109, 17978–17982 (2012).
  * Murugan  _et al._ (2012) Arvind Murugan, David A Huse, and Stanislas Leibler, “Speed, dissipation, and error in kinetic proofreading,” Proceedings of the National Academy of Sciences 109, 12034–12039 (2012).
  * Lang  _et al._ (2014) Alex H Lang, Charles K Fisher, Thierry Mora, and Pankaj Mehta, “Thermodynamics of statistical inference by cells,” Physical review letters 113, 148103 (2014).
  * Mehta  _et al._ (2016) Pankaj Mehta, Alex H Lang, and David J Schwab, “Landauer in the age of synthetic biology: energy consumption and information processing in biochemical networks,” Journal of Statistical Physics 162, 1153–1166 (2016).
  * Ouldridge  _et al._ (2017) Thomas E Ouldridge, Christopher C Govern, and Pieter Rein ten Wolde, “Thermodynamics of computational copying in biochemical systems,” Physical Review X 7, 021004 (2017).
  * Hess and Ross (2017) H Hess and Jennifer L Ross, “Non-equilibrium assembly of microtubules: from molecules to autonomous chemical robots,” Chemical Society Reviews 46, 5570–5587 (2017).
  * McInally  _et al._ (2021) Shane G McInally, Jane Kondev, and Bruce L Goode, “Scaling of subcellular actin structures with cell length through decelerated growth,” Elife 10, e68424 (2021).
  * Das  _et al._ (2021) Krishnendu Das, Luca Gabrielli, and Leonard J Prins, “Chemically fueled self-assembly in biology and chemistry,” Angewandte Chemie International Edition 60, 20120–20143 (2021).
  * Zeravcic  _et al._ (2014) Zorana Zeravcic, Vinothan N Manoharan, and Michael P Brenner, “Size limits of self-assembled colloidal structures made using specific interactions,” Proceedings of the National Academy of Sciences 111, 15918–15923 (2014).
  * Murugan  _et al._ (2015) Arvind Murugan, James Zou, and Michael P Brenner, “Undesired usage and the robust self-assembly of heterogeneous structures,” Nature communications 6, 1–10 (2015).
  * Nguyen and Vaikuntanathan (2016) Michael Nguyen and Suriyanarayanan Vaikuntanathan, “Design principles for nonequilibrium self-assembly,” Proceedings of the National Academy of Sciences 113, 14231–14236 (2016).
  * Bisker and England (2018) Gili Bisker and Jeremy L England, “Nonequilibrium associative retrieval of multiple stored self-assembly targets,” Proceedings of the National Academy of Sciences 115, E10531–E10538 (2018).
  * Murugan and Vaikuntanathan (2017) Arvind Murugan and Suriyanarayanan Vaikuntanathan, “Topologically protected modes in non-equilibrium stochastic systems,” Nature communications 8, 1–6 (2017).
  * Dasbiswas  _et al._ (2018) Kinjal Dasbiswas, Kranthi K Mandadapu, and Suriyanarayanan Vaikuntanathan, “Topological localization in out-of-equilibrium dissipative systems,” Proceedings of the National Academy of Sciences 115, E9031–E9040 (2018).
  * Knebel  _et al._ (2020) Johannes Knebel, Philipp M Geiger, and Erwin Frey, “Topological phase transition in coupled rock-paper-scissors cycles,” Physical Review Letters 125, 258301 (2020).
  * Tang  _et al._ (2021) Evelyn Tang, Jaime Agudo-Canalejo, and Ramin Golestanian, “Topology protects chiral edge currents in stochastic systems,” Physical Review X 11, 031015 (2021).
  * Yoshida  _et al._ (2021) Tsuneya Yoshida, Tomonari Mizoguchi, and Yasuhiro Hatsugai, “Chiral edge modes in evolutionary game theory: A kagome network of rock-paper-scissors cycles,” Physical Review E 104, 025003 (2021).
  * Van Kampen (1992) Nicolaas Godfried Van Kampen,  _Stochastic processes in physics and chemistry_ , Vol. 1 (Elsevier, 1992).
  * Zak (1989) J Zak, “Berry’s phase for energy bands in solids,” Physical review letters 62, 2747 (1989).
  * Cohen  _et al._ (2019) Eliahu Cohen, Hugo Larocque, Frédéric Bouchard, Farshad Nejadsattari, Yuval Gefen, and Ebrahim Karimi, “Geometric phase from aharonov–bohm to pancharatnam–berry and beyond,” Nature Reviews Physics 1, 437–449 (2019).
  * Ge and Qian (2010) Hao Ge and Hong Qian, “Physical origins of entropy production, free energy dissipation, and their mathematical representations,” Physical Review E 81, 051133 (2010).
  * Nakahara (2018) Mikio Nakahara,  _Geometry, topology and physics_ (CRC press, 2018).
  * Chiu  _et al._ (2016) Ching-Kai Chiu, Jeffrey CY Teo, Andreas P Schnyder, and Shinsei Ryu, “Classification of topological quantum matter with symmetries,” Reviews of Modern Physics 88, 035005 (2016).
  * Shen  _et al._ (2018) Huitao Shen, Bo Zhen, and Liang Fu, “Topological band theory for non-hermitian hamiltonians,” Physical review letters 120, 146402 (2018).
  * Kawabata  _et al._ (2019) Kohei Kawabata, Ken Shiozaki, Masahito Ueda, and Masatoshi Sato, “Symmetry and topology in non-hermitian physics,” Physical Review X 9, 041015 (2019).
  * Lieu (2018) Simon Lieu, “Topological phases in the non-hermitian su-schrieffer-heeger model,” Physical Review B 97, 045106 (2018).
  * Gillespie (1976) Daniel T Gillespie, “A general method for numerically simulating the stochastic time evolution of coupled chemical reactions,” Journal of computational physics 22, 403–434 (1976).
  * Virtanen  _et al._ (2020) Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, _et al._ , “Scipy 1.0: fundamental algorithms for scientific computing in python,” Nature methods 17, 261–272 (2020).

##  Appendix A Review of basic nonequilibrium thermodynamic steady-states

We start by reviewing some basic ideas on how to describe the nonequilibrium
steady-states (NESS) of open chemical systems. In doing so, we largely follow
the formalism developed by H. Qian and D. Beard in Beard and Qian (2008); Qian
and Beard (2005); Qian (2007).

The fundamental objects in our formalism for describing a chemical reaction
network are the concentrations of chemical species $c_{i}$ that can undergo
different _reversible_ chemical reactions, labeled by $\alpha$. Reversibility
is necessary in order to have a consistent thermodynamic description of the
system.

###  A.1 Fluxes in Open Systems

Each reaction $\alpha$ is characterized by a net flux $J_{\alpha}$, that can
be further decomposed into a forward flux $J_{\alpha}^{+}$ and a backward flux
$J_{\alpha}^{-}$ with

$J_{\alpha}=J_{\alpha}^{+}-J_{\alpha}^{-}.$ (12)

In equilibrium, we have detailed balance so that
$J_{\alpha}^{+}=J_{\alpha}^{-}=J_{\alpha}^{\pm eq}$ and there is no net
current $J_{\alpha}=0$. However, first emphasized by Hill Hill (2012), and
then Qian and colleagues Beard and Qian (2008); Qian and Beard (2005); Qian
(2007), in open chemical systems like cells where the concentration of some
chemical species are held fixed (e.g., ATP and ADP concentrations inside a
cell), reactions can be driven by consuming energy. In this case, we know in
general that $J_{\alpha}\neq 0$ and we can define a thermodynamic driving
force/chemical potential $\Delta\mu_{\alpha}$ associated with each reaction,

$\Delta\mu_{\alpha}=k_{B}T\ln{J_{\alpha}^{-}\over
J_{\alpha}^{+}}=k_{B}T\ln{\left(1-{J_{\alpha}\over J_{\alpha}^{+}}\right)}.$
(13)

We can combine Eqs. (12) and (13) to write

$J_{\alpha}=J_{\alpha}^{+}-J_{\alpha}^{-}=(e^{-{\Delta\mu_{\alpha}\over
k_{B}T}}-1)J_{\alpha}^{-}.$ (14)

Near equilibrium, $J_{\alpha}\ll J_{\alpha}^{\pm eq}$ so that Eq. (13) reduces
to the usual linear relationship between drives and fluxes for linear
response,

$\Delta\mu_{\alpha}\approx-{k_{B}T\over J_{\alpha}^{\pm eq}}J.$ (15)

Finally, we note that we can define the entropy production rate $epr$ as

$epr=-\sum_{\alpha}\Delta\mu_{\alpha}J_{\alpha}\geq 0,$ (16)

where the inequality follows from substituting Eqs. (12) and (13) into the
expression above.

Eq. (13) relating $\Delta\mu_{\alpha}$ to the ratio of the forward and
backward flux will play a central role in what follows.

###  A.2 Stoichiometric Matrix

The second important ingredient needed to describe NESS in our formalism is
the stoichiometric matrix $S_{i\alpha}$, where $i$ runs over chemical species
and $\alpha$ runs possible reactions. As the name suggests, $S_{i\alpha}$
encodes the stoichiometry of the reactions. In particular, $S_{i\alpha}$
encodes how many molecules $i$ are consumed (negative entries) or produced by
reaction $\alpha$ (positive entries). If a reaction $\alpha$ consumes two
molecules of species $l$ to produces a molecule of species $m$ and a molecule
of species $n$ then $S_{l\alpha}=-2$, $S_{m\alpha}=1$ and $S_{n\alpha}=1$ with
$S_{i\alpha}=0$ if $i\neq l,m,n$.

By definition, the dynamics of species $i$ are described by the differential
equation

${dc_{i}\over dt}=\sum_{i\alpha}S_{i\alpha}J_{\alpha}+b^{\mathrm{ext}}_{i},$
(17)

where $b^{\mathrm{ext}}_{i}$ is the rate at which $c_{i}$ is produced or
depleted by external sources or sinks (i.e., produced/depleted by reactions
not included in the network or through exchange with the environment). For a
NESS, we can set the left side to zero and we have

$0=\sum_{i\alpha}S_{i\alpha}J_{\alpha}+b^{\mathrm{ext}}_{i}$ (18)

For chemical species that are are not externally replenished or depleted, this
reduces to the familiar Flux Balance equation,

$0=\sum_{\alpha}S_{i\alpha}J_{\alpha}.$ (19)

###  A.3 Chemical potential differences for reactions

We can also associate each reaction $\alpha$ with a chemical potential
difference $\Delta\mu_{\alpha}$ which can be decomposed into three parts:

$\Delta\mu_{\alpha}=\Delta\mu_{\alpha}^{0}+\Delta\mu^{\mathrm{drive}}_{\alpha}+k_{B}T\sum_{i}S_{i\alpha}\ln{c_{i}},$
(20)

where $\Delta\mu_{\alpha}^{0}$ is just the equilibrium chemical potential
difference characterizing the reactions in the _absence_ of nonequilibrium
driving (i.e., the difference in the _absence_ of any coupling to high-energy
molecules maintained out of equilibrium by cells such as ATP/ADP, GTP/GDP,
etc.), the term $\Delta\mu^{\mathrm{drive}}_{\alpha}$ represents the
transduction of free energy by coupling the reaction to a to high-energy
molecular processes such as phosphorylation Hill (2012), and the last term
encodes the usual dependence of the free energy differences on the
concentrations of reactants and products.

It is worth better understanding this decomposition in greater detail. To do
so, first consider the two following simple chemical reactions:

$\displaystyle A$ $\displaystyle\xrightleftharpoons[k_{0}^{-}]{k_{0}^{+}}$
$\displaystyle B$ (21) $\displaystyle A+\mathrm{ATP}$
$\displaystyle\xrightleftharpoons[k^{-}]{k^{+}}$ $\displaystyle
B+\mathrm{ADP}+\mathrm{P}$ (22)

For Eq. (21) without a driving force, we know that
$\Delta\mu^{\mathrm{drive}}_{\alpha}=0$ and
$\Delta\mu=\Delta\mu^{0}+k_{B}T\ln{[B]\over[A]}$. Since at equilibrium, we
know that $\Delta\mu=0$ and detailed balance is satisfied,
$k_{0}^{-}[B]_{eq}=k_{0}^{+}[A]_{eq}$. For this reason, we conclude that we
must have the usual relation between kinetic constants and thermodynamic
chemical potentials,

$\Delta\mu^{0}=k_{B}T\ln{k_{0}^{-}\over k_{0}^{+}}.$ (23)

Now consider the second driven reaction in Eq. (22). In this case, we have by
definition that Qian (2007)

$\Delta\mu=\Delta\mu^{0}-k_{B}T\ln{[ATP]\over[ADP][P]}+k_{B}T\ln{[B]\over[A]}.$
(24)

Thus, we can identify

$\Delta\mu^{\mathrm{drive}}=-k_{B}T\ln{[ATP]\over[ADP][P]}.$ (25)

As emphasized by Qian, Beard, and colleagues, we can actually just implicitly
treat the high-energy phosphorylation reactions be defining effective forward
and backward rates $k_{+}=k_{0}^{+}[ATP]$ and $k_{-}=k_{0}^{-}[ADP][P]$. In
terms, of these effective rates for Eq. (22), we have

$\Delta\mu=k_{B}T\ln{k_{-}\over k_{+}}+k_{B}T\ln{[B]\over[A]}.$ (26)

Notice that this is just the special case of our general formula in Eq. (13)
that

$\Delta\mu_{\alpha}=k_{B}T\ln{J^{-}\over J^{+}}.$ (27)

More generally, we see that a very similar argument gives

$\Delta\mu_{\alpha}=k_{B}T\ln{k_{\alpha}^{-}\over
k_{\alpha}^{+}}+\sum_{i}S_{i\alpha}\ln{c_{i}},$ (28)

where we have defined “effective rate constants” $k_{\alpha}^{-}$ and
$k_{\alpha}^{+}$ that implicitly include any high-energy molecules involved in
the reactions. Furthermore, we know that we can identify the driving force
with the chemical potential difference,

$\Delta\mu^{\mathrm{drive}}_{\alpha}=k_{B}T\ln{k_{\alpha}^{-}\over
k_{\alpha}^{+}}-k_{B}T\ln{k_{\alpha 0}^{-}\over k_{\alpha 0}^{+}},$ (29)

or equivalently,

$\Delta\mu^{\mathrm{drive}}_{\alpha}+\Delta\mu_{\alpha}^{0}=k_{B}T\ln{k_{\alpha}^{-}\over
k_{\alpha}^{+}}.$ (30)

##  Appendix B Analytics: Edge modes in one-dimensional triangular lattice

We now give a detailed analysis of the triangular lattice discussed in the
main text.

###  B.1 Defining generalized stoichiometric matrix

Consider a periodic triangle network of the type shown in Fig. 1. The
primitive lattice for this consists of three fluxes and two nodes which are
repeated $N$ times. For a such a lattice, we can consider the stoichiometric
matrix in Fourier space. This matrix is given by

$\hat{S}(k)=\matrixquantity[1-e^{ik}&-1&e^{ik}\\\ 0&1&-1],$ (31)

where without loss of generality, we have set the lattice spacing to be one so
that $e^{ik}$ represents a shift by a single lattice spacing. In Fourier
space, the steady-state condition

$\sum_{\alpha}S_{i\alpha}J_{\alpha}=0$ (32)

becomes

$\hat{S}(k)\hat{J}(k)=0,$ (33)

where $\hat{J}(k)$ is just the usual Fourier transform of the flux
distributions,

$\hat{J}_{A}(k)=\sum_{\alpha=1}^{N}e^{ikj}J_{A,\alpha},$ (34)

where $A$ indicates the reaction type ($A=1,2,3$) and $\alpha$ ranges over the
unit cells.

To understand the topological phase transitions we must work with slightly
different objects. Instead of writing our dynamics in terms of the _net flux_
$J_{\alpha}$, it is useful to write the flux in terms of the backward flux
$J_{\alpha}^{-}$ (or alternatively, in terms of the forward flux
$J_{\alpha}^{+}$. Substituting Eq. (14) into Eq. (32) yields the modified
steady state equation

$\sum_{\alpha}S_{i\alpha}(e^{-{\Delta\mu_{\alpha}\over
k_{B}T}}-1)J_{\alpha}^{-}=0.$ (35)

Motivated by this equation, we define the parameter

$\theta_{\alpha}\equiv e^{-{\Delta\mu_{\alpha}\over k_{B}T}}-1,$ (36)

which measures the direction and magnitude of the flux across reaction
$\alpha$ and the chemical potential-dependent stochiometric matrix for the
negative fluxes

$S_{i\alpha}^{-}\equiv
S_{i\alpha}\theta_{\alpha}=S_{i\alpha}(e^{-{\Delta\mu_{\alpha}\over
k_{B}T}}-1).$ (37)

For the the periodic triangular lattice [see Eq. (31)], we can take the
Fourier transform of $S_{i\alpha}^{-}$ to get

$\hat{S}^{-}(k)=\matrixquantity[(1-e^{ik})\theta_{1}&-\theta_{2}&e^{ik}\theta_{3}\\\
0&\theta_{2}&-\theta_{3}],$ (38)

where $\theta_{A}=e^{-{\Delta\mu_{\alpha}\over k_{B}T}}-1$ is the drive
associated with the three reactions $A=1,2,3$ associated with the three fluxes
in the elementary lattices shown in Fig. 1.

###  B.2 Spectrum and Topology

We now show that we can associate a $Z_{2}$ topological invariant with the
triangular lattice by carefully analyzing the spectrum associated with
$\hat{S}^{-}(k)$ as a function of $\theta_{1}$, $\theta_{2}$, and $\theta_{3}$
(or equivalently, the driving forces $\Delta\mu_{1}$, $\Delta\mu_{2}$, and
$\Delta\mu_{3}$ associated with the three reactions).

To do so, we can look at the singular values associated with $\hat{S}^{-}(k)$
or equivalently, the eigenvalues associated with the symmetric $2\times 2$
matrix

$\hat{H}(k)=\hat{S}^{-}(k)\hat{S}^{-}(k)^{\dagger}=\matrixquantity[2(1-\cos{k})\theta_{1}^{2}+\theta_{2}^{2}+\theta_{3}^{2}&\,-\theta_{2}^{2}-\theta_{3}^{2}e^{ik}\\\
-\theta_{2}^{2}-\theta_{3}^{2}e^{-ik}&\,\theta_{2}^{2}+\theta_{3}^{2}]$ (39)

It is also helpful to decompose this Hamiltonian into a diagonal
$k$-independent term,

$H_{0}=\matrixquantity[\theta_{2}^{2}+\theta_{3}^{2}&0\\\
0&\,\theta_{2}^{2}+\theta_{3}^{2}],$ (40)

and a non-trivial $k$-dependent part ,

$\hat{H}_{1}(k)=\matrixquantity[2(1-\cos{k})\theta_{1}^{2}&\,-\theta_{2}^{2}-\theta_{3}^{2}e^{ik}\\\
-\theta_{2}^{2}-\theta_{3}^{2}e^{-ik}&\,0]$ (41)

so that

$\hat{H}(k)=H_{0}+\hat{H}_{1}(k).$ (42)

To identify the topological phase transition,we rewrite $\hat{H}_{1}(k)$ in
terms of the $2\times 2$ identity matrix $I$ as well as the three Pauli
matrices $\sigma_{x}$, $\sigma_{y}$, and $\sigma_{z}$, allowing us to connect
to standard arguments Nakahara (2018),

$\hat{H}_{1}(k)=\theta_{1}^{2}(1-\cos{k})I-(\theta_{2}^{2}+\theta_{3}^{2}\cos{k})\sigma_{x}-\theta_{3}^{2}\sin{k}\sigma_{y}+\theta_{1}^{2}(1-\cos{k})\sigma_{z}.$
(43)

In particular, for a Hamiltonian of the form
$H=a_{0}+\vec{a}\cdot\vec{\sigma}$, the two eigenvalues are given by

$\lambda_{\pm}=a_{0}\pm\sqrt{\vec{a}\cdot\vec{a}}.$ (44)

A $Z_{2}$ topological transition can occur when the gap between these two
eigenvalues closes for any $k$.

We now focus on the special case where $\Delta\mu_{1}=0$, or equivalently
$\theta_{1}=0$ for which the topological transition occurs. This corresponds
to the physical situation where reaction 1 is in detailed balance. This
scenario can be attained by setting
$\Delta\mu^{0}_{1}+\Delta\mu^{\mathrm{drive}}_{1}=0$ and noting that the
concentrations of the chemical species connected by type 1 reactions must be
uniform due to translational symmetry. As a result, $\Delta\mu_{1}$ and
$\theta_{1}$ are guaranteed to be zero. For this choice of parameters,
$\hat{H}_{1}(k)=\vec{a}\cdot\vec{\sigma}$ with

$\vec{a}=(-(\theta_{2}^{2}+\theta_{3}^{2}\cos{k}),-\theta_{3}^{2}\sin{k},0).$
(45)

Notice that

$\vec{a}\cdot\vec{a}=\theta_{2}^{4}+2\theta_{2}^{2}\theta_{3}^{2}\cos{k}+\theta_{3}^{4}$
(46)

so that for $k=\pi$ and $\theta_{2}^{2}=\theta_{3}^{2}$,
$\vec{a}\cdot\vec{a}=0$ and hence the gap in the spectrum closes.

To show that this is associated with a $\mathbb{Z}_{2}$ topological invariant,
we can plot the two-dimensional closed curve defined by the coefficients of
$\sigma_{x}$ and $\sigma_{y}$, namely $(a_{x}(k),a_{y}(k))$ where $k$ varies
from $0$ to $2\pi$. A $\mathbb{Z}_{2}$ topological transition occurs at the
values of $\theta_{2}$ and $\theta_{3}$ where this closed curve transitions
from enclosing the origin to not enclosing the origin when the gap also
closes. For $\hat{H}_{1}(k)$ in Eq. (41), this closed curve is given by

$(a_{x}(k),a_{y}(k))=(-\theta_{2}^{2}-\theta_{3}^{2}\cos{k},-\theta_{3}^{2}\sin{k}).$
(47)

These are just the coordinates for a circle with radius $\theta_{3}^{2}$
centered at the point $(0,-\theta_{2}^{2})$. Thus, we see that there is a
topological phase transition when $\theta_{2}^{2}=\theta_{3}^{2}$. For
$\theta_{2}^{2}<\theta_{3}^{2}$, the curve does not enclose the origin,
whereas for $\theta_{2}^{2}>\theta_{3}^{2}$ it does enclose the origin. This
choice of parameters, $\theta_{2}^{2}=\theta_{3}^{2}$ is also where the
spectrum in the gap closes for $k=\pi$ as expected [see Fig. A1].

##  Appendix C Analytics: Chiral currents in two-dimensional rectangular
lattice

In this section, we analyze a simple model introduced in Ref. Tang  _et al._
(2021) for biopolymer formation. In this model, biopolymers can contain
different numbers of two monomers $X$ and $Y$ and be in one of four different
internal states $z=A,B,C,D$. As shown in Fig. 2, these internal states tag a
monomer and prime the biopolymer for the addition or removal of a subunit.

The state of the system is determined by the ordered pair $(x,y)$ that
represents the number of each monomer incorporated into the biopolymer and an
internal variable $w=A,B,C,D$ that encodes the internal state. We assume that
the biopolymer can have at most $N_{x}$ of the $X$ monomers and $N_{y}$ of the
$Y$ monomers so that $0\leq x\leq N_{x}$ and $0\leq y\leq N_{y}$. This basic
dynamics can be visualized as a lattice where edges correspond to reactions
and nodes to different states of the biopolymer.

Let us denote by $c_{i}$ the number of polymers in state $i$ (where the index
$i=(x,y,z)$ refers to an ordered triplet specifying the state of the
biopolymer). We know that the dynamics of this system can also be specified by
a stoichiometric matrix $S_{i\alpha}$ where $i$ runs over biopolymer states
(i.e., nodes in the lattice) and $\alpha$ runs over possible reactions (i.e.,
edges in the lattice). In writing $S_{i\alpha}$, we must make a gauge choice
for which direction in each reaction we call forward and which reaction we
call backward. We choose the convention shown in Fig. 2. With this choice, we
have

$S_{i\alpha}=\begin{cases}1\text{ if reaction }\alpha\text{ ``produces'' a
biopolymer in state }i\\\ -1\text{ if reaction }\alpha\text{ ``consumes'' a
biopolymer in state }i\\\ 0\text{ otherwise}\end{cases}$ (48)

and the full dynamics are given by

${dc_{i}\over dt}=\sum_{\alpha}S_{i\alpha}J_{\alpha},$ (49)

where $J_{\alpha}$ is the net flux across reaction $\alpha$.

It is clear from Fig. 2 that $S_{i\alpha}$ is periodic in both the $x$\- and
$y$-directions. The elementary lattice cell consists of four chemical species
corresponding to the four internal states $z=A,B,C,D$ at a fixed $(x,y)$ and
eight elementary reactions corresponding to the four internal reactions
(modifying the state $z$) and four external reactions (adding or removing a
monomer of $X$ or $Y$). Translations in $X$($Y$) correspond to adding/removing
a monomer of type $X$($Y$). For these reasons it is natural to consider the
Fourier transform of this stoichiometric matrix,

$\hat{S}(k)=\matrixquantity[-1&1&0&0&e^{ik_{x}\over 2}&0&-e^{ik_{y}\over
2}&0\\\ 0&-1&1&0&0&-e^{ik_{x}\over 2}&e^{-{ik_{y}\over 2}}&0\\\
0&0&-1&1&0&e^{-{ik_{x}\over 2}}&0&-e^{-{ik_{y}\over 2}}\\\
1&0&0&-1&-e^{-{ik_{x}\over 2}}&0&0&e^{ik_{y}\over 2}],$ (50)

where $e^{ik_{x}}$ encodes shifts in the $x$-direction and $e^{ik_{y}}$ encode
shifts in the $y$-direction.

As before, instead of writing our dynamics in terms of the _net flux_
$J_{\alpha}$, it is useful to write the flux in terms of the backwards flux
$J_{\alpha}^{-}$. Substituting Eq. (14) into Eq. (32) yields the modified
steady-state equation

$\sum_{\alpha}S_{i\alpha}(e^{-{\Delta\mu_{\alpha}\over
k_{B}T}}-1)J_{\alpha}^{-}=0.$ (51)

Defining the parameter

$\theta_{\alpha}\equiv e^{-{\Delta\mu_{\alpha}\over k_{B}T}}-1,$ (52)

yields the chemical potential-dependent stoichiometric matrix

$S_{i\alpha}^{-}\equiv
S_{i\alpha}\theta_{\alpha}=S_{i\alpha}(e^{-{\Delta\mu_{\alpha}\over
k_{B}T}}-1).$ (53)

In Fourier space, this becomes

$\hat{S}^{-}(k)=\matrixquantity[-\theta_{1}&\theta_{2}&0&0&\theta_{5}e^{ik_{x}\over
2}&0&-\theta_{7}e^{ik_{y}\over 2}&0\\\
0&-\theta_{2}&\theta_{3}&0&0&-\theta_{6}e^{ik_{x}\over
2}&\theta_{7}e^{-{ik_{y}\over 2}}&0\\\
0&0&-\theta_{3}&\theta_{4}&0&\theta_{6}e^{-{ik_{x}\over
2}}&0&-\theta_{8}e^{-{ik_{y}\over 2}}\\\
\theta_{1}&0&0&-\theta_{4}&-\theta_{5}e^{-{ik_{x}\over
2}}&0&0&\theta_{8}e^{ik_{y}\over 2}].$ (54)

In order to identify a topological invariant, we examine the spectrum of

$\hat{H}(k)=\hat{S}^{-}(k)\hat{S}^{-}(k)^{\dagger}=\matrixquantity[\theta_{1}^{2}+\theta_{2}^{2}+\theta_{5}^{2}+\theta_{7}^{2}&-\theta_{2}^{2}-\theta_{7}^{2}e^{ik_{y}}&0&-\theta_{1}^{2}-\theta_{5}^{2}e^{ik_{x}}\\\
-\theta_{2}^{2}-\theta_{7}^{2}e^{-ik_{y}}&\theta_{2}^{2}+\theta_{3}^{2}+\theta_{5}^{2}+\theta_{7}^{2}&-\theta_{3}^{2}-\theta_{6}^{2}e^{ik_{x}}&0\\\
0&-\theta_{3}^{2}-\theta_{6}^{2}e^{-ik_{x}}&\theta_{3}^{2}+\theta_{4}^{2}+\theta_{6}^{2}+\theta_{8}^{2}&-\theta_{4}^{2}-\theta_{8}^{2}e^{-ik_{y}}\\\
-\theta_{1}^{2}-\theta_{5}^{2}e^{-ik_{x}}&0&-\theta_{4}^{2}-\theta_{8}^{2}e^{ik_{y}}&\theta_{1}^{2}+\theta_{4}^{2}+\theta_{5}^{2}+\theta_{8}^{2}].$
(55)

We now focus on the case when the “internal” reactions are identical so that
$\Delta\mu^{0}_{\alpha}+\Delta\mu^{\mathrm{drive}}_{\alpha}=a$ for
$\alpha=1,2,3,4$ and the “external” reactions are identically driven so that
$\Delta\mu^{0}_{\alpha}+\Delta\mu^{\mathrm{drive}}_{\alpha}=b$ for
$\alpha=5,6,7,8$. In this case, translational and 4-fold rotational symmetry
dictate that the concentrations for a periodic lattice must all be uniform so
that the thermodynamic drives have the same symmetry as the driving forces
with two corresponding values, $\theta_{a}$ and $\theta_{b}$. This gives us

$\hat{H}(k)=\hat{S}^{-}(k)\hat{S}^{-}(k)^{\dagger}=\matrixquantity[2\theta_{a}^{2}+2\theta_{b}^{2}&-\theta_{a}^{2}-\theta_{b}^{2}e^{ik_{y}}&0&-\theta_{a}^{2}-\theta_{b}^{2}e^{ik_{x}}\\\
-\theta_{a}^{2}-\theta_{b}^{2}e^{-ik_{y}}&2\theta_{a}^{2}+2\theta_{b}^{2}&-\theta_{a}^{2}-\theta_{b}^{2}e^{ik_{x}}&0\\\
0&-\theta_{a}^{2}-\theta_{b}^{2}e^{-ik_{x}}&2\theta_{a}^{2}+2\theta_{b}^{2}&-\theta_{a}^{2}-\theta_{b}^{2}e^{-ik_{y}}\\\
-\theta_{a}^{2}-\theta_{b}^{2}e^{-ik_{x}}&0&-\theta_{a}^{2}-\theta_{b}^{2}e^{ik_{y}}&2\theta_{a}^{2}+2\theta_{b}^{2}],$
(56)

where $\theta_{a}=e^{-{a\over k_{B}T}}-1$ and $\theta_{b}=e^{-{b\over
k_{B}T}}-1$. As before, it is helpful to separate out the $k$-dependent and
independent parts,

$\hat{H}(k)=H_{0}-\hat{H}_{1}(k),$ (57)

with

$H_{0}=\matrixquantity[2\theta_{a}^{2}+2\theta_{b}^{2}&0&0&0\\\
0&2\theta_{a}^{2}+2\theta_{b}^{2}&0&0\\\
0&0&2\theta_{a}^{2}+2\theta_{b}^{2}&0\\\
0&0&0&2\theta_{a}^{2}+2\theta_{b}^{2}],$ (58)

and

$\hat{H}_{1}(k)=\matrixquantity[0&\theta_{a}^{2}+\theta_{b}^{2}e^{ik_{y}}&0&\theta_{a}^{2}+\theta_{b}^{2}e^{ik_{x}}\\\
\theta_{a}^{2}+\theta_{b}^{2}e^{-ik_{y}}&0&\theta_{a}^{2}+\theta_{b}^{2}e^{ik_{x}}&0\\\
0&\theta_{a}^{2}+\theta_{b}^{2}e^{-ik_{x}}&0&\theta_{a}^{2}+\theta_{b}^{2}e^{-ik_{y}}\\\
\theta_{a}^{2}+\theta_{b}^{2}e^{-ik_{x}}&0&\theta_{a}^{2}+\theta_{b}^{2}e^{ik_{y}}&0].$
(59)

We can then rewrite this as

$\hat{H}_{1}(k)=I\otimes[\vec{a}(k_{y})\cdot\vec{\sigma}]+\sigma_{x}\otimes[\vec{a}(k_{x})\cdot\vec{\sigma}],$
(60)

where $\vec{a}=(\theta_{a}^{2}+\theta_{b}^{2}\cos{x},\theta_{b}^{2}\sin{x},0)$
and $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of the
Pauli matrices. This suggests that in analogy with the triangular lattice,
there should be a topological phase transition when
$\theta_{a}^{2}=\theta_{b}^{2}$, the gap closes, and $\vec{a}(k_{x})$ and
$\vec{a}(k_{y})$ go from enclosing the origin to not enclosing the origin.

This basic intuition can be verified by directly diagonalizing
$\hat{H}_{1}(k)$ and looking at the corresponding eigenvalues and
eigenvectors. To do so, it is useful to define

$A(k)=\theta_{a}^{2}+\theta_{b}^{2}e^{ik}=|A(k)|e^{i\phi(k)},$ (61)

where
$|A(k)|=\sqrt{\theta_{a}^{4}+\theta_{b}^{4}+2\theta_{a}^{2}\theta_{b}^{2}\cos{k}}$
is just the magnitude and $\phi_{k}$ is the corresponding phase. In terms of
these, we can write the four eigenvalues as

$\displaystyle\lambda=\left\\{\pm\mathrm{abs}\left(|A(k_{y})|^{2}-|A(k_{x})|^{2}\right),\pm\left(|A(k_{y})|^{2}+|A(k_{x})|^{2}\right)\right\\}.$
(62)

We see that the spectrum becomes degenerate when $|A(k_{x})|^{2}=0$ or
$|A(k_{y})|^{2}=0$ which occur when $\theta_{a}^{2}=\theta_{b}^{2}$ and
$k_{x}=\pi$ or $k_{y}=\pi$, respectively. We can also look at the
corresponding eigenvectors. To do so, we define a variable
$d=\mathrm{sgn}(\left(|A(k_{y})|^{2}-|A(k_{x})|^{2}\right)$ that takes on
values $\pm 1$. In terms of these quantities, one can show that the four
eigenvectors are

$\displaystyle|e_{1}\rangle$ $\displaystyle=$
$\displaystyle\quantity[de^{i\phi(k_{x})},-e^{i[\phi(k_{x})-\phi(k_{y})]},-de^{-i\phi(k_{y})},1]$
$\displaystyle|e_{2}\rangle$ $\displaystyle=$
$\displaystyle\quantity[-de^{i\phi(k_{x})},-e^{i[\phi(k_{x})-\phi(k_{y})]},de^{-i\phi(k_{y})},1]$
$\displaystyle|e_{3}\rangle$ $\displaystyle=$
$\displaystyle\quantity[-e^{i\phi(k_{x})},e^{i[\phi(k_{x})-\phi(k_{y})]},-e^{-i\phi(k_{y})},1]$
$\displaystyle|e_{4}\rangle$ $\displaystyle=$
$\displaystyle\quantity[e^{i\phi(k_{x})},e^{i[\phi(k_{x})-\phi(k_{y})]},e^{-i\phi(k_{y})},1].$
(63)

The corresponding gauge potential (i.e., Berry phase) for the first eigenstate
is just

$(A_{kx},A_{ky})=(\langle e_{1}|\partial_{k_{x}}e_{1}\rangle,\langle
e_{1}|\partial_{k_{y}}e_{1}\rangle)=2i\quantity({\partial\phi(k_{x})\over\partial
k_{x}},{\partial\phi(k_{y})\over\partial k_{y}}).$ (64)

We can then look at the corresponding charge over a Wilson loop (contour
integral in $k_{x}$ and $k_{y}$ space) to define the $\mathbb{Z}_{2}$
topological invariant. We see that when $\theta_{a}^{2}\geq\theta_{b}^{2}$,
this phase is trivial in both cases (the origin is not enclosed), whereas for
$\theta_{a}^{2}\leq\theta_{b}^{2}$ it is nontrivial and encloses origin. This
confirms that there is a topological phase transition when
$\theta_{a}^{2}=\theta_{b}^{2}$.

The exact nature of the phase where topological edge currents are supported
can be gleaned by numerical simulations. We find that edge currents exist in
one of these phases, but not the other [see Fig. 1].

##  Appendix D Numerical Simulations

In this section, we provide details on numerical simulations. All code used to
perform simulations and generate figures can be found in the corresponding
Jupyter Notebook on our GitHub repository at https://github.com/Emergent-
Behaviors-in-Biology/NESS-Topology.

###  D.1 One-dimensional triangular lattice

In order to simulate the triangular lattice, we directly simulate the ODE
system defined by Eqs. (1) and (2) using ODEint from SciPY package Virtanen
_et al._ (2020). We specify the six kinetic constants corresponding to the
forward and backward rates for the three types of reactions in the lattice
cells as $k_{1}^{\pm}$, $k_{2}^{\pm}$, and $k_{3}^{\pm}$. For all three
reactions, we choose backward reactions so that
$k_{1}^{-}=k_{2}^{-}=k_{3}^{-}=1$. Furthermore, we choose forward rates
$k_{j}^{+}$ to be thermodynamically consistent using Eq. (30) so that

$k_{A}^{+}=k_{A}^{-}e^{-{\Delta\mu_{A}^{0}+\Delta\mu^{\mathrm{drive}}_{A}\over
k_{B}T}},$ (65)

for $A=1,2,3$. As discussed in main text, we choose the first reaction to be
at equilibrium so that

$\Delta\mu_{1}^{0}+\Delta\mu^{\mathrm{drive}}_{1}=0$ (66)

and vary $\Delta\mu_{A}^{0}+\Delta\mu^{\mathrm{drive}}_{A}$ for reactions 2
and 3 as indicated in the main text. We choose initial conditions by sampling
$c_{i}$ uniformly in the interval $[0,5]$ (though any reasonable initial
conditions for which the integrator converges give identical results). We
compute fluxes, chemical potentials, and thermodynamic drives directly from
the steady-state concentrations. See Jupyter Notebook for detailed code.

###  D.2 Two-dimensional rectangular lattice

In order to simulate the two-dimensional rectangular lattice, we make use of
the Gillespie algorithm for simulating chemical reaction networks Gillespie
(1976). Direct numerical simulation of the corresponding ODEs proved to be
numerically unstable due to positivity constraints. In addition to the
stoichiometric matrix, simulations require us to specify the 16 kinetic
constants corresponding to the 8 types of reactions that can occur. Following
Tang et. al. Tang  _et al._ (2021), we divide this 8 reactions into two
classes: 4 “external reactions” that change the number of monomers [reactions
5-8 in Fig. 2] and 4 “internal reactions” [reactions 1-4 in Fig. 2] that
change the internal states of the biopolymer.

We choose $k_{j}^{-}=0.01$ ($j=1,\ldots,8$) for all 8 reactions and use Eq.
(30) to set

$k_{j}^{+}=e^{-{\Delta\mu_{a}^{0}+\Delta\mu^{\mathrm{drive}}_{a}\over
k_{B}T}}k_{j}^{-}$ (67)

for reactions $j=1,2,3,4$ and and

$k_{j}^{+}=e^{-{\Delta\mu_{b}^{0}+\Delta\mu^{\mathrm{drive}}_{b}\over
k_{B}T}}k_{j}^{-}$ (68)

for reactions $j=5,6,7,8$. The two parameters
$\Delta\mu_{a}^{0}+\Delta\mu^{\mathrm{drive}}_{a}$ and
$\Delta\mu_{b}^{0}+\Delta\mu^{\mathrm{drive}}_{b}$ are chosen such that
$k_{a}^{-}=k_{b}^{-}=0.01$, $k_{b}^{+}=0.03$ and $k_{a}^{+}$ is varied.

To calculate the average abundance of a chemical species $c_{i}$ we average
the stochastic trajectory corresponding to the species. The average flux is
then calculated directly from the simulated average abundances using the
formula for the currents and mass action kinetics. See Jupyter Notebook for
detailed code.

##  Appendix E Supplementary Figures

Figure A1: Analytics for triangular lattice. (Top) Singular value spectrum of
$\hat{S}^{-}(k)$. (Middle) Eigenvalue spectrum of
${\hat{H}(k)=\hat{S}^{-}(k)\hat{S}^{-}(k)^{\dagger}}$. (Bottom) Topological
classification using map from $\hat{H}_{1}(k)=a_{0}+\vec{a}\cdot\vec{\sigma}$
for $k=[0,2\pi]$ to $(\sigma_{x},\sigma_{y})$. Curve is
$(a_{x}(\vec{k}),a_{y}(\vec{k}))$. Topological invariant measures whether
origin is enclosed.


nt\rfloor+j}\left(-\frac{1}{2}w_{\mu}(n\eta,j)\cdot\hat{\mu}^{(d)}(Z)\right)$
$\displaystyle\qquad~{}+2\left(-\bar{Z}^{(d)}\right)\left(-\frac{1}{2}w_{\mu}(n\eta,j)\cdot\hat{\mu}^{(d)}(Z)\right)\bigg{)}$
$\displaystyle=\sum_{d=1,2}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\sum_{j=1}^{\lfloor n\eta\rfloor}\left(Z^{(d)^{2}}_{\lfloor
nt\rfloor+j}+\bar{Z}^{(d)^{2}}+\frac{1}{4}w_{\mu}(n\eta,j)^{2}\hat{\mu}^{(d)}(Z)^{2}\right.$
$\displaystyle\qquad~{}-2Z^{(d)}_{\lfloor
nt\rfloor+j}\bar{Z}^{(d)}-w_{\mu}(n\eta,j)Z^{(d)}_{\lfloor
nt\rfloor+j}\hat{\mu}^{(d)}(Z)+w_{\mu}(n\eta,j)\bar{Z}^{(d)}\hat{\mu}^{(d)}(Z)\bigg{)}$
$\displaystyle=\sum_{d=1,2}\Bigg{(}\underbrace{\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor n\eta\rfloor-4}Z^{(d)^{2}}_{\lfloor
nt\rfloor+j}}_{=:S_{1}}+\underbrace{\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\bar{Z}^{(d)^{2}}}_{=:S_{2}}$
$\displaystyle\qquad~{}+\underbrace{\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\frac{1}{4}w_{\mu}(n\eta,j)^{2}\hat{\mu}^{(d)}(Z)^{2}}_{=:S_{3}}-\underbrace{\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor n\eta\rfloor-4}2Z^{(d)}_{\lfloor
nt\rfloor+j}\bar{Z}^{(d)}}_{=:S4}$
$\displaystyle\qquad~{}-\underbrace{\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}w_{\mu}(n\eta,j)Z^{(d)}_{\lfloor
nt\rfloor+j}\hat{\mu}^{(d)}(Z)}_{=:S_{5}}$
$\displaystyle\qquad~{}+\underbrace{\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}w_{\mu}(n\eta,j)\bar{Z}^{(d)}\hat{\mu}^{(d)}(Z)}_{=:S_{6}}\Bigg{)}$
(61)

Since the summands in the two dimensions $d=1,2$ are independent, we focus on
one dimension and suppress the index $d$. As $\hat{\sigma}^{2}(nt,n\eta)$ is
unbiased we can split the expectation $\sigma^{2}$ into the expectations of
each of the summands. We show that each summand in the brackets in (61)
centred with its expectation converges to $0$ almost surely.
For summand 1 we have

$\displaystyle S_{1}-\mathbb{E}[S_{1}]$ $\displaystyle=\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\left(Z^{(d)^{2}}_{\lfloor
nt\rfloor+j}-\mathbb{E}\bigg{[}Z^{(d)^{2}}_{\lfloor
nt\rfloor+j}\bigg{]}\right)$ $\displaystyle=\sum_{j=1}^{\lfloor
n\eta\rfloor}v(n,j)\left(\underbrace{Z^{(d)^{2}}_{\lfloor
nt\rfloor+j}-\mathbb{E}\bigg{[}Z^{(d)^{2}}_{\lfloor
nt\rfloor+j}\bigg{]}}_{\mathbb{E}[*]=0}\right),$

where we here define $v(n,j):=\frac{\sigma^{2}}{2\lfloor n\eta\rfloor-4}$.
According to Proposition  B.1, we need to verify
$\sum_{j=1}^{m_{n}}v(n,j)^{2}\in\mathrm{O}(n^{-1})$. This holds because

$\displaystyle\sum_{j=1}^{\lfloor n\eta\rfloor}v(n,j)^{2}=\sum_{j=1}^{\lfloor
n\eta\rfloor}\left(\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\right)^{2}=\sigma^{4}\frac{\lfloor n\eta\rfloor}{(2\lfloor
n\eta\rfloor-4)^{2}}\in\mathrm{O}(n^{-1}).$

For summand 2, we use $\mathbb{E}\bar{Z}^{(d)^{2}}=1/\lfloor n\eta\rfloor$ and
get

$\displaystyle S_{2}-{\mathbb{E}[S_{2}]}$ $\displaystyle=\frac{\lfloor
n\eta\rfloor\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\left(\bar{Z}^{(d)^{2}}-\mathbb{E}\bar{Z}^{(d)^{2}}\right)$
$\displaystyle=\frac{\lfloor n\eta\rfloor\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\bar{Z}^{(d)^{2}}-\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\xrightarrow[n\rightarrow\infty]{}0\text{ a.s.}$

For summand 3 we first rewrite $S_{3}$ as

$\displaystyle S_{3}$ $\displaystyle=\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}\frac{1}{4}w_{\mu}(n\eta,j)^{2}\left(\frac{6}{\lfloor
n\eta\rfloor^{3}-\lfloor n\eta\rfloor}\sum_{k=1}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,k)Z^{(d)}_{\lfloor nt\rfloor+k}\right)^{2}$
$\displaystyle=\sum_{k,l=1}^{\lfloor
n\eta\rfloor}\underbrace{v(n,k,l)}_{\text{ in part 2 of Proposition
\ref{proposition:SLLNnew}}}Z^{(d)}_{\lfloor nt\rfloor+k}Z^{(d)}_{\lfloor
nt\rfloor+l},$

with

$\displaystyle v(n,k,l):=\frac{6^{2}\sigma^{2}}{4(2\lfloor
n\eta\rfloor-4)(\lfloor n\eta\rfloor^{3}-\lfloor
n\eta\rfloor)^{2}}\left(\sum_{j=1}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,j)^{2}\right)w_{\mu}(n\eta,k)w_{\mu}(n\eta,l)$

such that we have

$\displaystyle S_{3}-\mathbb{E}[S_{3}]$ $\displaystyle=\sum_{k,l=1}^{\lfloor
n\eta\rfloor}v(n,k,l)\left(Z^{(d)}_{\lfloor nt\rfloor+k}Z^{(d)}_{\lfloor
nt\rfloor+l}-\mathbb{E}\left[Z^{(d)}_{\lfloor nt\rfloor+k}Z^{(d)}_{\lfloor
nt\rfloor+l}\right]\right).$

To show almost sure convergence to $0$ we need to verify that the term in (57)
is in $\mathrm{O}(n^{-2})$. We thus note that for the present definition of
$v(n,k,l)$, the term in (57) corresponds to the following term

$\displaystyle\underbrace{\frac{6^{8}\sigma^{8}}{4^{4}(2\lfloor
n\eta\rfloor-4)^{4}(\lfloor n\eta\rfloor^{3}-\lfloor
n\eta\rfloor)^{8}}\left(\frac{\lfloor n\eta\rfloor^{3}-\lfloor
n\eta\rfloor)}{3}\right)^{4}}_{\in\mathrm{O}(n^{-28}n^{12})=\mathrm{O}(n^{-16})}$
$\displaystyle~{}\times\bigg{[}\underbrace{\sum_{j,k=1}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,j)^{4}w_{\mu}(n\eta,k)^{4}}_{\lfloor
n\eta\rfloor^{2}\cdot\mathrm{O}(n^{8})=\mathrm{O}(n^{10})}$
$\displaystyle\qquad~{}+\underbrace{\sum_{\begin{subarray}{c}i,j,k,l=1\\\
(i,j)\neq(k,l)\end{subarray}}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,i)^{2}w_{\mu}(n\eta,j)^{2}w_{\mu}(n\eta,k)^{2}w_{\mu}(n\eta,l)^{2}}_{\in\lfloor
n\eta\rfloor^{4}\cdot\mathrm{O}(n^{8})=\mathrm{O}(n^{12})}$
$\displaystyle\qquad~{}+\underbrace{\sum_{\begin{subarray}{c}i,j,k,l=1\\\
(i,j)\neq(k,l)\end{subarray}}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,i)^{3}w_{\mu}(n\eta,j)^{3}w_{\mu}(n\eta,k)w_{\mu}(n\eta,l)}_{\in\lfloor
n\eta\rfloor^{4}\cdot\mathrm{O}(n^{8})=\mathrm{O}(n^{12})}$
$\displaystyle\qquad~{}+\underbrace{\sum_{\begin{subarray}{c}i,j,k,l,o,p=1\\\
(i,j)\neq(k,l)\\\ (i,j)\neq(o,p)\\\ (o,p)\neq(k,l)\end{subarray}}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,i)w_{\mu}(n\eta,j)w_{\mu}(n\eta,k)w_{\mu}(n\eta,l)w_{\mu}(n\eta,o)^{2}w_{\mu}(n\eta,p)^{2}}_{\in\lfloor
n\eta\rfloor^{6}\cdot\mathrm{O}(n^{8})=\mathrm{O}(n^{14})}$
$\displaystyle\qquad~{}+\sum_{\begin{subarray}{c}i,j,k,l,o,p,r,s=1;\\\
(i,j),(k,l),\\\ (o,p),(r,s)\\\ \text{ pairwise
different}\end{subarray}}^{\lfloor
n\eta\rfloor}\prod_{q=i,j,k,l,o,p,r,s}w_{\mu}(n\eta,q)\bigg{]}.$ (62)

All sums in the latter display except for the last one are easily upper
bounded and seen to be in
$\mathrm{O}(n^{14}/n^{16})=\mathrm{O}\left(n^{-2}\right)$. The last sum
requires an exact, tedious calculation. We used the Python module Sympy and
obtained $\mathrm{O}(n^{12}/n^{16})$.

For summand 4 we use $\mathbb{E}[Z^{(d)}_{\lfloor
nt\rfloor+i}\bar{Z}]=1/\lfloor n\eta\rfloor$ and get

$\displaystyle S_{4}-\mathbb{E}[S_{4}]$ $\displaystyle=\sum_{j=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{\lfloor n\eta\rfloor-2}\left(Z^{(d)}_{\lfloor
nt\rfloor+j}\bar{Z}^{(d)}-\frac{1}{\lfloor
n\eta\rfloor}\right)\xrightarrow[n\rightarrow\infty]{}0$ (63)

by the conventional SLLN.

For summand 5 we again first rewrite $S_{5}$ as

$\displaystyle S_{5}$ $\displaystyle=\sum_{i=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}w_{\mu}(n\eta,i)Z^{(d)}_{\lfloor
nt\rfloor+i}\left(\frac{6}{\lfloor n\eta\rfloor^{3}-\lfloor
n\eta\rfloor}\sum_{j=1}^{\lfloor n\eta\rfloor}w_{\mu}(n\eta,j)Z^{(d)}_{\lfloor
nt\rfloor+j}\right),$

which yields

$\displaystyle S_{5}-\mathbb{E}[S_{5}]$ $\displaystyle=\sum_{i,j=1}^{\lfloor
n\eta\rfloor}v(n,i,j)\left(Z^{(d)}_{\lfloor nt\rfloor+i}Z^{(d)}_{\lfloor
nt\rfloor+j}-\mathbb{E}\left[Z^{(d)}_{\lfloor nt\rfloor+i}Z^{(d)}_{\lfloor
nt\rfloor+j}\right]\right)$

where we define

$\displaystyle v(n,i,j):=\frac{6\sigma^{2}}{(2\lfloor n\eta\rfloor-4)(\lfloor
n\eta\rfloor^{3}-\lfloor n\eta\rfloor)}w_{\mu}(n\eta,i)w_{\mu}(n\eta,j).$

Now we are in the setting of the SLLN as in Proposition C.4, part 2 because
the summands are not necessarily independent. Therefore we need to verify that
the term in (57) is in $\mathrm{O}(n^{-2})$. We define the expression in the
square brackets in (62) as $L$, for which we showed $L\in\mathrm{O}(n^{14})$
and find that the term in (57) in the current setting can be written as

$\displaystyle\underbrace{\frac{6^{4}\sigma^{8}}{(2\lfloor
n\eta\rfloor-4)^{4}(\lfloor n\eta\rfloor^{3}-\lfloor
n\eta\rfloor)^{4}}}_{\in\mathrm{O}(n^{-16})}\cdot\underbrace{L}_{\in\mathrm{O}(n^{14})}\in\mathrm{O}(n^{-2}).$

For summand 6 we first observe that

$\displaystyle\mathbb{E}\bigg{[}\bar{Z}^{(d)}\hat{\mu}^{(d)}(Z)\bigg{]}$
$\displaystyle=\frac{6}{\lfloor n\eta\rfloor(\lfloor n\eta\rfloor^{3}-\lfloor
n\eta\rfloor)}\sum_{k=1}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,k)\underbrace{\mathbb{E}\bigg{[}Z^{(d)^{2}}_{\lfloor
nt\rfloor+k}\bigg{]}}_{=1}$ $\displaystyle=\frac{6}{\lfloor
n\eta\rfloor(\lfloor n\eta\rfloor^{3}-\lfloor
n\eta\rfloor)}\underbrace{\sum_{k=1}^{\lfloor
n\eta\rfloor}w_{\mu}(n\eta,k)}_{=0}=0,$ (64)

such that

$\displaystyle{S_{6}-\mathbb{E}[S_{6}]}=\sum_{i=1}^{\lfloor
n\eta\rfloor}\frac{\sigma^{2}}{2\lfloor
n\eta\rfloor-4}w_{\mu}(n\eta,i)\bar{Z}^{(d)}\hat{\mu}^{(d)}(Z)$

This term vanishes because

$\displaystyle\sum_{i=1}^{\lfloor n\eta\rfloor}w_{\mu}(n\eta,i)=0.$

In total, each summand of $\hat{\sigma}^{2}(nt,n\eta)-\sigma^{2}$ converges to
$0$ almost surely, which implies

$\displaystyle\hat{\sigma}^{2}(nt,n\eta)\xrightarrow[n\rightarrow\infty]{}\sigma^{2}\text{
a.s.}$