\section{Physical example of KL divergence as energy cost} \label{app:physical_systems}

\begin{figure}[!t]
    \centering

    \includegraphics[width=3in]{figures/KL_div_work_v4.pdf}

    \caption{(top) A gas molecule is found underneath the $i$-th piston with
    probability proportional to volume $V_i$. (bottom) Inserting impermeable
    partitions between the pistons, isothermally compressing to new volumes
    $V'_i$, and then removing the partitions incurs expected work on the system
    proportional to the KL divergence between the initial and final probability
    distributions of the molecule's position.}
    \label{fig:KL_div_work}
\end{figure}

We offer an example of the KL divergence as the cost of sampling from a
target distribution $p$ given a ``base'' distribution $q$ (see discussion in
section \ref{sec:KL_div_cost}).

Figure \ref{fig:KL_div_work} presents a slight generalization of the Szilard
information engine \cite{Szilard29}: molecules of an ideal gas inhabit the
space formed by movable pistons indexed by $i \in \mathcal{X}$.  Let $V_i$ denote
the volume beneath the $i$-th piston and $V \triangleq \sum_{i \in \mathcal{X}} V_i$ be
the total volume.  A molecule of gas is equally likely to be found anywhere
within the space beneath the pistons, corresponding to probability distribution
$p \triangleq (V_i/V)_{i \in \mathcal{X}}$ on the pistons.  Now imagine we add
impermeable partitions between the pistons (vertical dashed lines) and move the
pistons to new positions (green lines) at constant temperature (perhaps the
bottom of the box is in thermal contact with a heat reservoir).  The partitions
prevent mixing between different pistons during compression; we remove them
afterwards.  Let $V_i'$ be the new volume beneath the $i$-th piston, $V' \triangleq
\sum_{i \in \mathcal{X}} V'_i$ be the new total volume, and $q \triangleq (V'_i/V)_{i
\in \mathcal{X}}$ be the new piston probability distribution after this deformation.

What is the work used to perform this deformation?  The work to isothermally
compress an ideal gas is $k_\text{B} \tau \log(V_i/V'_i)$, where $k_\text{B}$
is Boltzmann's constant and $\tau$ is the temperature.  A gas molecule
occupies the space beneath the $i$-th piston with probability $V_i/V$ before
compression, so the expected work to move the pistons is
\begin{equation}\label{gaspistons}
    \sum_{i \in \mathcal{X}} \frac{V_i}{V} k_\text{B} \tau
    \log\left(\frac{V_i}{V'_i}\right) = k_\text{B} \tau (D(p||q) + \log(V/V'))
\end{equation}
per molecule of gas.  If the new pistons positions are such that the total
volume is unchanged, then $V' = V$ and the work is proportional to the KL
divergence $D(p||q)$. We can imagine a sequence of such gas boxes and
deformations, where the pre-deformation volumes $(V_{t,i})_{i \in \mathcal{X}}$ at
time $t$ are determined by drawing a single molecule from volumes
$(V_{t-1,i})_{i \in \mathcal{X}}$ at time $t-1$, forming a Markov chain with KL
divergence control cost.

\section{Proof of Theorem 1} \label{app:proof}

\subsection{Part 1}
We wish to solve the following problem:
\begin{equation}
    P^* \triangleq \argmin_{P: \pi(P)= \pi^*} D(P || Q) \label{eq:def_optimization_problem}
\end{equation}

For ease of manipulation, we work with the empirical joint transition probability distribution $N = (N[i,j])_{i,j\in\mathcal{X}}$
\begin{equation}
	N[i,j](x_0^n) \triangleq \frac{1}{n}|\{t : (x_{t-1},x_t) = (i,j)\}|.
\end{equation}

We can solve (\ref{eq:def_optimization_problem}) by setting up the Lagrangian:
\begin{align}
    \Lambda \triangleq & D(P||Q) + \alpha \left(\sum_{i,j \in \mathcal{X}} N_P[i,j] - 1\right) \nonumber \\
    &+ \sum_{i \in \mathcal{X}} \lambda_i \left(\sum_{j \in \mathcal{X}} N_P[i,j] - \pi^*[i]\right) \\
    & + \sum_{j \in \mathcal{X}} \eta_j \left(\sum_{i \in \mathcal{X}} N_P[i,j] - \pi^*[j]\right) \label{eq:def_Lagrangian}
\end{align}
where $\alpha$ is a Lagrange multiplier enforcing normalization of the joint transition probability distribution under $P$, $N_P$, and $(\lambda_i)_{i \in \mathcal{X}}$, $(\eta_i)_{i \in \mathcal{X}}$ are Lagrange multipliers enforcing the stationary distribution condition $\pi(P) = \pi^*$.  Our solution is a stationary point of the Lagrangian with respect to $N_P$:
\begin{align}
&\left.\frac{\partial \Lambda}{\partial N_P[i,j]} \right|_{N_P = N^*} = \left. \frac{\partial D(P||Q)}{\partial N_P[i,j]} + \alpha + \lambda_i + \eta_j \right|_{N_P = N^*} = 0 \nonumber \\
&\forall i,j \label{eq:Lagrangian_partial_derivative}
\end{align}
Since $\pi_P = \pi^*$, then $N^*$, $A^*$ satisfy:
\begin{align}
    \left.\frac{\partial D(P||Q)}{\partial N_P[i,j]}\right|_{N_P = N^*} &= \left.\frac{\partial}{\partial N_P[i,j]}\left(D(N_P||N_Q) - D(\pi^*||\pi)\right)\right|_{N_P = N^*} \\
    &= \left.\frac{\partial}{\partial N_P[i,j]}D(N_P||N_Q)\right|_{N_P = N^*} \\
    &= \left.1 + \log\left(\frac{N_P[i,j]}{N_Q[i,j]}\right)\right|_{N_P = N^*} \label{eq:Lagrangian_partial_derivative_expression}
\end{align}
Now using (\ref{eq:def_Lagrangian}), (\ref{eq:Lagrangian_partial_derivative}), and (\ref{eq:Lagrangian_partial_derivative_expression}) we find:
\begin{align}
    N^*[i,j] = N_Q[i,j]e^{-1-\alpha-\lambda_i-\eta_j}
\end{align}
Equivalently
\begin{equation}
    N^* = \text{diag}((e^{-\lambda_i})_{i \in \mathcal{X}}) \ N_Q \ \text{diag}((e^{-\eta_i})_{i \in \mathcal{X}}) e^{-1-\alpha}
\end{equation}
Now using the condition $\pi(P) = \pi^*$:
\begin{align}
    \pi^*[i] &= \sum_{j \in \mathcal{X}} N^*[i,j] \\
    &= e^{-1-\alpha-\lambda_i}\sum_{j \in \mathcal{X}} N_Q[i,j] e^{-\eta_j}
\end{align}
Therefore
\begin{equation}
    \lambda_i = \log\left(\sum_{j \in \mathcal{X}} N_Q[i,j] e^{-\eta_j}\right) - \log(\pi^*[i]) - 1 - \alpha \label{eq:recursive_lambda}
\end{equation}
Analogously
\begin{equation}
    \eta_i = \log\left(\sum_{j \in \mathcal{X}} N_Q[j,i] e^{-\lambda_j}\right) - \log(\pi^*[i]) - 1 - \alpha \label{eq:recursive_eta}
\end{equation}
Now using the normalization condition
\begin{equation}
    1 = \sum_{i,j} N^*[i,j]
\end{equation}
we find
\begin{equation}
    \alpha = \log\left(\sum_{i,j \in \mathcal{X}} N_Q[i,j]e^{-\lambda_i - \eta_j}\right) - 1 \label{eq:recursive_alpha}
\end{equation}

\subsection{Part 2}
In the continuous case, we work with rate matrices instead of probability transition matrices. We wish to solve the following problem:

\begin{equation}
\begin{aligned}
& \underset{\B{P}: \pi{(\B{P})} = \pi^*}{\text{min}}
& & D(\B{P}||\B{Q}) \\
\end{aligned}
\end{equation}
The Lagrangian is $D(\B{P}||\B{Q})+\sum_i \lambda_i \sum_j \B{P}_{ij} + \sum_j \eta_j \sum_i \pi_i \B{P}_{ij}.$ Differentiating it we get the conditions
\begin{equation}
      \B{P}_{ij}^* = \begin{cases} \label{eq:rec}
	\B{Q}_{ij} e^{\eta_i-\eta_j} & i \neq j \\
	-\sum_{j:j\neq i} \B{Q}_{ij} e^{\eta_i-\eta_j} & i = j \\
   \end{cases}
\end{equation}

\section{Discussion}\label{sec:discussion}

This work derives the minimum power required to maintain a target stationary
distribution given uncontrolled Markov chain dynamics in both discrete and
continuous time.  We relate KL divergence-like penalties from control theory
\cite{Todorov06,Todorov09} to the power used to control a Markov process, using
muscular molecular motors as a guiding example.  The problem of minimizing a KL
divergence subject to a constrained stationary distribution is familiar from
large deviations theory
\cite{DemboZeitouniBook,BaldiPiccioni98,CsiszarCoverChoi87}; the novelty of our
work is in relating these large deviations results to the thermodynamics of
``holding'' a distribution, and in computing the minimum-cost chains in some
important examples: the birth-and-death process in continuous time (a toy model
for a muscular molecular motor) and two-state chains in discrete time.

To the best of our knowledge, this is the first time a lower bound on average
power consumption has been studied in detail for control of the stationary
distribution.  \cite{HorowitzEngland17} study a related quantity, the minimum
entropy production rate associated with adding edges (allowing control to
increase but not decrease transition rates) to a continuous time Markov chain,
but their notion of cost has the interpretation of energy efficiency, as
opposed to ours, which is to be interpreted as total energy usage.
Unsurprisingly, different notions of cost lead to different optimal controlled
chains: the optimal controlled chain of \cite{HorowitzEngland17} depends on the
underlying uncontrolled chain only in the requirement that it be much faster,
while our minimum-cost controlled chain is a function of the uncontrolled
chain; this function is easy to compute (\ref{eq:reversible_cont_time}) in the
case of a continuous time, reversible uncontrolled chain, an important case in
modeling biological processes.

\section{Introduction} \label{sec:introduction}

Let $(X_1,X_2,\ldots)$ be a sequence of random variables forming a first-order
Markov chain on a finite set $\mathcal{X}$ with
transition probability
matrix $Q$, where $Q_{ij} = \mathbb{P}_Q(X_{t+1}=j|X_t=i)$
for all $t$ and $i,j\in\mathcal{X}$
We think of $Q$ as the ``uncontrolled'' or ``base'' chain.
If $Q$ is irreducible and aperiodic, then there exists a unique stationary
distribution $\pi$, viewed as a row vector:
\begin{equation}

    \pi Q = \pi

\end{equation}
Let the initial state $X_1$ have distribution $\mu$.  Then, denoting $X_1^T
\triangleq (X_1, \ldots X_T)$, the distribution of $X_1^T$ is:
\begin{align}
    p_{Q,\mu}(x_1^T) &\triangleq \mathbb{P}_{Q,\mu}(X_1^T = x_1^T) = \mu_{x_1}
    \prod_{t=2}^T Q_{x_{t-1},x_t} \label{eq:def_MC_distribution}
\end{align}

Let $\pi^*$ be some probability distribution on $\mathcal{X}$ viewed as a row
vector.  We study the nearest Markov chain transition matrix $P$ to $Q$ having
$\pi^*$ as its stationary distribution:
\begin{align}

    P^* &= \argmin_{P: \pi^* P = \pi^*} D(P||Q)
    \label{eq:optimization_problem_in_introduction}
\end{align}
where $D(P||Q)$ is the Kullback-Leibler (KL) divergence rate between Markov
chains with transition matrices $P$ and $Q$
\cite{CsiszarCoverChoi87}~\footnote{This expression is independent of $\mu =
\mathbb{P}(X_1)$ for aperiodic, irreducible $P$.}:
\begin{equation}

    D(P||Q) \triangleq \lim_{T \rightarrow \infty} \frac{1}{T} \mathbb{E}_P
    \log\left(\frac{p_P(X_0^T)}{p_Q(X_0^T)}\right)
    \label{eq:def_KL_div_MC_discrete_time_in_introduction}
\end{equation}
where we note that $D(P||Q)$ is independent of the initial distribution $\mu$.

We think of $P$ as the ``controlled'' or ``driven'' chain and $D(P||Q)$ as the
cost of control per unit time -- the power.
We further consider the analogous question for the case of a continuous-time
Markov chain $(X_t)_{t \in [0,\infty)}$.

This setting is inspired by the following thought experiment due to Feynman
\cite{FeynmanLectures}: a person holds a heavy bag above the floor for an hour
and gets tired.  The net work done on the object is zero\footnote{Since work is
the product of force and distance displaced, the latter of which is zero.}, so
why does she get tired?  A table could hold the same bag indefinitely without
an energy source, and so could the person if she were frozen solid,
the bag hanging on her stiff, lifeless limb.
The latter observation implicates the microscopic dynamics of muscles as key to
this question.
A toy model for the motion of striated muscle fibers -- see
\cite{Huxley57,Qian00} for an extended discussion -- is of a random walk in a
periodic energy potential\footnote{The myosin protein joining the fibers
together is doing the random walking, the energy potential having
periodically-spaced minima corresponding to
discrete steps along the fiber.}.  Attachment of the heavy bag pulls on the
muscle fibers, biasing the random walk in the direction of gravity, tilting the
energy potential.  The person must use chemical energy\footnote{The hydrolysis
of ATP molecules.} to de-bias the random walk in such a way that the bag is
held at the desired height above the floor.
If the person is frozen solid,
then the underlying random walk stops\footnote{Or slows down a lot, as lowering
the temperature reduces the transition rate between the potential's energy
minima.}, and so chemical energy is no longer used to hold the bag.
We ask: what is a lower bound on the power to hold the bag?

We recast the above story as optimization problem
(\ref{eq:optimization_problem_in_introduction}).
The state space $\mathcal{X}$ corresponds to the possible configurations of the
system (the position of myosin along a fiber and its internal state).
The uncontrolled Markov chain $Q$ corresponds to the underlying fluctuations of
the myosin molecule along a filament and the controlled chain $P$ corresponds
to chemically driving the system.

The control goal is macroscopic: the net force the person exerts on the bag is
the sum of the forces due to each microscopic subsystem (myosin protein).  We
do not get to control each subsystem separately, but can control them all in
the same way, so that each subsystem corresponds to a trajectory drawn from
Markov chain $P$ independently of other subsystems.  This notion of macroscopic
control is reflected in the KL divergence cost
(\ref{eq:def_KL_div_MC_discrete_time_in_introduction}), which is stated in
terms of the probability distribution (\ref{eq:def_MC_distribution}) over
microscopic trajectories, rather than in terms of a single trajectory.

Our choice of the Kullback-Leibler divergence as the control cost function is
motivated by this quantity's appearance in thermodynamics as proportional to
the free energy difference from the equilibrium distribution over trajectories,
which in turn lower bounds the work to prepare a non-equilibrium distribution
over trajectories.
KL divergence control (the microscopic per-trajectory setting) was introduced
in the reinforcement learning literature by Todorov \cite{Todorov06,Todorov09}
and has connections to data compression; we discuss this and other motivations
for our work in section \ref{sec:KL_div_cost}. The problem of maintaining a
target nonequilibrium distribution has been studied recently by
\cite{HorowitzZhouEngland17,HorowitzEngland17}, using a different notion of
cost.  We discuss our work's relation to prior works in section \ref{subsec:
comparison}.

Minimizing the KL divergence with respect to the first argument --  computing
the I-projection --
connects this setting to large deviations theory.  Indeed, most of the
minimum-cost controlled chains we compute first appeared in the computation of
rate functions for large deviations of the empirical measure of both discrete-
and continuous-time Markov chains
\cite{DemboZeitouniBook,BaldiPiccioni98,CsiszarCoverChoi87}.  The novelty of
our work lies in relating these results to the minimum-power control setting --
showing the minimized KL divergence to be the minimum power to ``hold'' a
nonequilibrium distribution; in aggregating related problem statements -- in
continuous and discrete time and for reversible base chains; in computing some
of these minimizations more explicitly than previously available; and in
computing these minimizations for a few common examples like the
birth-and-death chain and the two-state chain.

This work is organized as follows.  Section \ref{sec:KL_div_cost} motivates the
use of KL divergence as energy cost by drawing from information theory and
optimization settings and contains definitions of this cost function in
discrete and continuous time. Section \ref{sec:minimum_power} shows how this
energy cost of holding a given target distribution may be analytically
minimized.  Section \ref{sec:examples} contains several examples that apply our
theory to calculate the minimum-power controlled chain, including a
birth-and-death chain which serves as a toy model of the muscular fiber,
addressing the motivating question of Feynman. We conclude with a summary and
outlook in section \ref{sec:discussion}.

We release code for computing the minimum-cost chains in this work at
\href{https://github.com/dmitrip/controlledMC}{https://github.com/dmitrip/controlledMC}.

\section{Kullback-Leibler divergence rate as the cost of control of Markov
chains} \label{sec:KL_div_cost}

We motivate the KL divergence rate between Markov chains in both discrete and
continuous time as the cost function lower bounding the power in the
bag-holding thought experiment.  We present a thermodynamics perspective in
subsection \ref{subsec:thermo} and an equivalent perspective due to Todorov
\cite{Todorov06,Todorov09} of a Markov decision process with log ratio cost
function in subsection \ref{subsec:log_loss_action_cost}.  We summarize known
expressions for the KL divergence between Markov chains in discrete and
continuous time in subsection \ref{sec:MC_discrete_time} and
\ref{sec:MC_continuous_time} respectively. Finally, \ref{subsec: comparison}
places this work in context with related work.

\subsection{KL divergence in thermodynamics}\label{subsec:thermo}

We summarize briefly the appearance of the KL divergence in measuring work in
statistical mechanics.
Below, let $D(p||q) = \sum_{i \in \mathcal{X}} p(i) \log(p(i)/q(i))$ denote the
Kullback-Leibler (KL) divergence between distributions $p$ and $q$ on finite
set $\mathcal{X}$ and let $H(p) = \sum_{i \in \mathcal{X}} p(i) \log(1/p(i))$ denote the
entropy of distribution $p$.
Let \begin{equation}
    q = \left(\frac{1}{Z}e^{-\beta U(i)}\right)_{i \in \mathcal{X}}

    \label{eq:Boltzmann_distribution}
\end{equation}
be the Boltzmann distribution on $\mathcal{X}$, where $U(\cdot)$
is the energy function (sometimes called the energy potential or internal
energy), $\beta$ is the inverse temperature, and $Z = \sum_{i \in \mathcal{X}} e^{-\beta U(i)}$ is the partition function.
Denote the free energy of distribution $p$ by $F(p)$
\cite{MezardMontanariBook}:
\begin{align}

    F(p) &\triangleq \mathbb{E}_p(U(X)) - \frac{1}{\beta} H(p) \\

    &= \frac{1}{\beta} \mathbb{E}_p \log
    \left(\frac{p_X}{q_X}\right) - \frac{1}{\beta} \log(Z)\\

    &= \frac{1}{\beta} D(p||q) - \frac{1}{\beta} \log(Z)
\end{align}
where the expectation $\mathbb{E}_p(\cdot)$ is over random variable $X$ with
distribution $p$.  Then $F(\cdot)$ is minimized at equilibrium $p=q$, so that
$F(q) = -\frac{1}{\beta}\log(Z)$.

In thermodynamics \cite{KittelKroemerBook,LandauLifshitz} the work $W$ to
prepare the ``controlled'' distribution $p$ starting from the ``base'',
equilibrium distribution $q$ (also known as the work on the system) is at least
the free energy difference:
\begin{equation}

    W \geq F(p) - F(q) = \frac{1}{\beta} D(p||q)
    \label{eq:work_lower_bound_free_energy}
\end{equation}
As is customary in this setting, there is in the background a notion of a
stochastic process transforming initial states into final states, and the work
$W$ in (\ref{eq:work_lower_bound_free_energy}) is an average over realizations
of this process.  In Appendix (\ref{app:physical_systems}) we provide a
physical example in the spirit of the Szilard's engine thought experiment in
thermodynamics, for which the KL divergence does emerge as the work done.
In the bag-holding thought experiment, there is a large collection of
independent myosin systems, and the total work is the sum of the works on each
system.  We imagine the number of microscopic systems to be large enough that
fluctuations about the average work $W$ per subsystem are small, so it is this
average work that's our object of study.

The KL divergence cost is familiar in data compression, where the ``energy'' of
symbol $X$ drawn from distribution $p$ is $U(X) = -\log(p_X)$.  Any compression
scheme must use at least $\mathbb{E}_p(U(X)) = H(p)$ bits  to encode a sample from
$p$ \cite{Shannon1948,CoverThomasBook} on average over draws from distribution
$p$.  If we use a compression scheme that instead uses $-\log(q_X)$ bits to
encode symbol $X$ -- a mismatched code -- then we would pay $D(p||q)$ extra
bits per symbol on average.  Section \ref{sec:MC_discrete_time} contains
analogous remarks for compressing samples drawn from Markov chain
distributions.

\subsection{Markov chains in discrete time} \label{sec:MC_discrete_time}

In the Markov chain control setting, we apply the preceding picture with
alphabet $\mathcal{X}^{T}$ (trajectories of length $T$) and with Markov chain
distributions on $\mathcal{X}^T$ with a desired marginal distribution $\pi^*$.  We
consider the continuous time setting in section \ref{sec:MC_continuous_time}.

A discrete time Markov chain distribution $p_{Q,\mu}^{(T)}$ on the set
$\mathcal{X}^{T}$ is the Boltzmann distribution with energy function
$U_{Q,\mu}^{(T)}$ parametrized by the stochastic transition matrix $Q =
(Q_{ij})_{i,j\in\mathcal{X}}$ and the initial distribution vector $\mu = (\mu)_{i
\in \mathcal{X}}$ (obtained by taking the logarithm of the rightmost quantity in
(\ref{eq:def_MC_distribution})):
\begin{equation}

    U_{Q,\mu}^{(T)}(x_1^T) = -\log(\mu_{x_1}) - \sum_{t=2}^{T}
    \log(Q_{x_{t-1},x_t})

\end{equation}
where $x_1^T \triangleq (x_1, \ldots, x_T)$.  Then $p_{Q,\mu}^{(T)}(x_1^T) = e^{-U_{Q,\mu}^{(T)}(x_1^T)}$, where $Z = \beta =
1$ (\ref{eq:Boltzmann_distribution}).  Given another transition matrix $P$ and
initial distribution $\nu$, the work to prepare Markov chain distribution
$p_{P,\nu}^{(T)}$ starting from the Boltzmann distribution $p_{Q,\mu}^{(T)}$ is
lower bounded by the free energy difference $F(p_{P,\nu}^{(T)}) -
F(p_{Q,\mu}^{(T)})$ (\ref{eq:work_lower_bound_free_energy}).  In the limit $T
\rightarrow \infty$, the work per time step -- the power -- is lower bounded by
\begin{align}
    \lim_{T \rightarrow \infty} \frac{1}{T} W(p_{P,\nu}^{(T)}) &\geq \lim_{T \rightarrow
    \infty} \frac{1}{T} \left(F(p_{P,\nu}^{(T)}) - F(p_{Q,\mu}^{(T)})\right)
    \label{eq:power_lower_bound_free_energy_mc} \\
    &= \lim_{T \rightarrow \infty} \frac{1}{T} D(p_{P,\nu}^{(T)}||p_{Q,\mu}^{(T)}) \\

    &= \lim_{T \rightarrow \infty} \frac{1}{T} \mathbb{E}_{p_{P,\nu}^{(T)}}
    \log\left(\frac{p_{P,\nu}^{(T)}(X_1,\ldots,X_T)}{p_{Q,\nu}^{(T)}(X_1,\ldots,X_T)}\right)
    \label{eq:KL_div_MC_discrete_time_expression} \\
    &= \sum_{i \in \mathcal{X}} \pi(P)_i \sum_{j \in \mathcal{X}}
    P_{ij} \log\left(\frac{P_{ij}}{Q_{ij}}\right) \\
    &\triangleq D(P||Q) \label{eq:def_KL_div_MC_discrete_time_1}
\end{align}
where $\pi(P)$ is the stationary distribution of transition matrix $P$ and the
last equality defines the KL divergence rate \cite{CsiszarCoverChoi87} between
Markov chains with transition matrices $P$ and $Q$.

In data compression, a sample $X_1^T \sim p_{Q,\mu}^{(T)}$ can be compressed on
average to at least $T H(Q)$ bits \cite{CoverThomasBook}, where $H(Q)$ is the
entropy rate of Markov chain with transition matrix $Q$:
\begin{equation}
    H(Q) = -\sum_{i \in \mathcal{X}} \pi(Q)_i \sum_{j \in \mathcal{X}} Q_{ij}
    \log(Q_{ij})
\end{equation}
Encoding samples from distribution $p_{Q,\mu}^{(T)}$ with respect to a
mismatched code based on distribution $p_{P,\nu}^{(T)}$ incurs an average cost
per unit time of at least $D(P||Q)$ extra bits.

Maximizing the lower bound on power (\ref{eq:power_lower_bound_free_energy_mc})
over transition matrices $P$ with target stationary distribution $\pi^*$ yields
the optimization problem (\ref{eq:optimization_problem_in_introduction}):
\begin{equation}
    P^* = \argmin_{P: \pi^* P = \pi^*} D(P||Q)
    \label{eq:optimization_problem_in_MC}
\end{equation}

\subsection{Log loss action cost} \label{subsec:log_loss_action_cost}

Another path to the same optimization problem
(\ref{eq:optimization_problem_in_MC}) (KL divergence as a lower bound on the
work to maintain a nonequilibrium distribution) is in terms of a Markov
decision process with log loss action cost, a setting introduced by Todorov
\cite{Todorov06,Todorov09}.  Let $Q$ be the uncontrolled chain and let $P$ be
the controlled chain.  Let $c(i, j, P, Q)$ be the \textit{microscopic cost}
paid when a transition is made from $X_t = i$ to $X_{t+1} = j$ when the
controller chooses transition probability matrix $P$.  KL divergence control
amounts to using the log likelihood ratio:
\begin{equation}

    c(i, j, P, Q) \triangleq \log\left(\frac{P_{ij}}{Q_{ij}}\right)
\end{equation}
If we view the rows of $P$ and $Q$ as Boltzmann distributions with different
energy potentials -- that is, if we choose energy functions $E_i(\cdot)$,
$E'_i(\cdot)$ such that $P_{ij} = e^{-E_i(j)}$ and $Q_{ij} = e^{-E'_i(j)}$ --
then the microscopic cost is the difference in energies: $c(i,j,P,Q) = E'_i(j)
- E_i(j)$.

If $X_t \sim \mu$, then let the \textit{cost} $D_\mu(P||Q)$
be the expected microscopic cost (the average cost paid per microscopic
system):
\begin{align}

    D_\mu(P||Q) &\triangleq \mathbb{E}_{\mu, P}(c(X_t, X_{t+1}, P, Q)) \\

    &= \sum_{i \in \mathcal{X}} \mu_i \sum_{j \in \mathcal{X}} P_{ij}\log
    \left(\frac{P_{ij}}{Q_{ij}}\right) \label{eq:cost_disc_time_KL_div}
\end{align}
Thus
$D_\mu(P||Q)$ is a $\mu$-weighted KL divergence between the rows of transition
matrices $P$ and $Q$.  We are interested in macroscopic control of $\mu_t$,
rather than microscopic control of $X_t$, so our setup differs from the setting
introduced by \cite{Todorov06,Todorov09}: we average the control cost over
$\mu_t$, so there is no randomness in our setting.  Finally for irreducible,
aperiodic transition matrix $P$, we have the identity
\begin{equation}

    D_{\pi(P)}(P||Q) = D(P||Q) \label{eq:def_KL_div_MC_discrete_time}
\end{equation}
where $\pi(P)$ is the stationary distribution of $P$.  Minimizing the cost
$D_{\pi^*}(P||Q)$ with respect to $P$ such that $\pi^* P = \pi^*$ is
optimization problem (\ref{eq:optimization_problem_in_MC}).

\subsection{Markov chains in continuous time} \label{sec:MC_continuous_time}

The setup of section \ref{sec:MC_discrete_time} has a natural counterpart for
continuous-time Markov chains.  Let $\bar{Q} = (\bar{Q}_{ij})_{i,j \in \mathcal{X}}$
denote the transition rate matrix of the uncontrolled continuous-time Markov
chain $(X_t)_{t \in \mathbb{R}_{\geq 0}}$, where henceforth the overbar notation
corresponds to rate matrices.  Let $\bar{P}$ be the controlled rate matrix and
let $X_t \sim \mu_t$.  Then
\begin{equation}
    \frac{d}{dt} \mu_t = \left(\frac{d}{dt} \mathbb{P}_{\bar{P}}(X_t = i)\right)_{i
    \in \mathcal{X}} = \mu_t \bar{P} \Rightarrow \mu_t = \mu_0 e^{\bar{P} t}
    \label{eq:controlled_state_continuous_time}
\end{equation}
where $e^{\bar{P}t}$ denotes the matrix exponential.  Note that every rate
matrix $\bar{P}$ satisfies $\bar{P}_{ij} \geq 0$ for $i \neq j$ and
$\bar{P}_{ii} = -\sum_{j\in\mathcal{X}:j\neq i} \bar{P}_{ij} \leq 0$, so the row
sums of $\bar{P}$ are 0.  Conditioned on $X_t = i$, the time until the next
jump is exponentially-distributed with a mean of $-1/\bar{P}_{ii}$, and the
probability to jump to $j$ is proportional to $\bar{P}_{ij}$ for $i \neq j$.

The natural notion of KL divergence rate $D(\bar{P}||\bar{Q})$ between
transition rate matrices $\bar{P}$ and $\bar{Q}$ is
\cite{Kesidis93,BaldiPiccioni98,deLaFortelle00, BertiniFaggionatoGabrielli14}
the limiting log likelihood ratio, analogous to
(\ref{eq:KL_div_MC_discrete_time_expression}):
\begin{align}
    D(\bar{P}\|\bar{Q}) &\triangleq \lim_{T \rightarrow \infty} \frac{1}{T}
    \mathbb{E}_{p_{\bar{P},\mu_0}^{(T)}}
    \log\left(\frac{p_{\bar{P},\mu_0}^{(T)}((X_t)_{t \in
    [0,T]})}{p_{\bar{Q},\nu_0}^{(T)}((X_t)_{t \in [0,T]})}\right) \\
    &= \sum_{i \in \mathcal{X}} \pi(\bar{P})_i \sum_{j \in \mathcal{X}: j \neq i}
    \left(\bar{Q}_{ij} - \bar{P}_{ij} + \bar{P}_{ij}
    \log\left(\frac{\bar{P}_{ij}}{\bar{Q}_{ij}}\right)\right)
    \label{eq:KL_div_MC_continuous_time_expression}
\end{align}
where $p_{\bar{P},\mu_0}^{(T)}(X_0^T)$ denotes the likelihood under rate matrix
$\bar{P}$ and initial distribution $\mu_0$, and where $\pi(\bar{P})$ is the
stationary distribution of rate matrix $\bar{P}$~\footnote{That is,
$\pi(\bar{P}) \bar{P} = (0)_{i \in \mathcal{X}}$.  Equivalently, $\pi(\bar{P})
e^{\bar{P} t} = \pi(\bar{P})$ for all $t \geq 0$.}.  The quantity in the second
summation in (\ref{eq:KL_div_MC_continuous_time_expression}) is the KL
divergence between two Poisson distributions with means $\bar{P}_{ij}$ and
$\bar{Q}_{ij}$.

The optimization problem analogous to (\ref{eq:optimization_problem_in_MC}) in
continuous time is:
\begin{equation}
    \bar{P}^* = \argmin_{\bar{P}: \pi(\bar{P}) = \pi^*} D(\bar{P}||\bar{Q})
    \label{eq:optimization_problem_in_MC_continuous_time}
\end{equation}

\subsection{Comparison to prior work} \label{subsec: comparison}

Recent work \cite{HorowitzZhouEngland17,HorowitzEngland17} considers the
question of the minimum power used to maintain a nonequilibrium state.  Their
setting uses a different notion of cost than we do and also makes some
restrictions about the base and controlled chains (they work in continuous
time, assume that the base chain $\bar{Q}$ is reversible, and only allow
controlled chains $\bar{P}$ such that $\bar{P}_{ij} \geq \bar{Q}_{ij}$ for all
$i \neq j$ -- this corresponds to the biochemical mechanism of adding transitions with non-negative rates).  \cite{HorowitzZhouEngland17,HorowitzEngland17} minimize the
entropy production rate among all controlled chains with the desired target
distribution $\pi^*$ and find that ``fast control is optimal'': there is in
general no optimally controlled chain, but given any chain $\bar{P}$ that has
the target distribution $\pi^*$, we can come arbitrarily close to the minimum
entropy production bound by speeding up $\bar{P}$ arbitrarily
much\footnote{That is, using $c \bar{P} = (c \bar{P}_{ij})_{i,j\in\mathcal{X}}$ as
the controlled chain and letting $c \rightarrow \infty$. This corresponds to the
statement in \cite{HorowitzZhouEngland17} than the ``added edges (should)
operate much faster than the equilibrium transitions''.} (while incurring an
arbitrarily large KL divergence cost $D(\bar{P}||\bar{Q})$ according to our
metric).

The difference between the two notions of cost -- KL divergence rate in this
work and entropy production rate in
\cite{HorowitzZhouEngland17,HorowitzEngland17} is the difference between total
energy used and the efficiency with which that energy is used as measured by
entropy production rate.  The very fast controlled chain of
\cite{HorowitzZhouEngland17,HorowitzEngland17} uses a lot of energy
efficiently, while our chain $\bar{P}^*$ minimizes energy use by the
controller, but not the efficiency.  A consequence of this difference is that
our optimal controlled chain $\bar{P}^*$
(\ref{eq:optimization_problem_in_MC_continuous_time}) depends on the
uncontrolled chain $\bar{Q}$ (see section \ref{sec:minimum_power}), while the
very fast close-to-optimal chain of
\cite{HorowitzZhouEngland17,HorowitzEngland17} does not, except in the
requirement that it be much faster than $\bar{Q}$.

\cite{Todorov06,Todorov09} introduced the KL divergence control setting.
\cite{Todorov06,Todorov09} uses ``microscopic'' control cost, assigning a cost
to a trajectory rather than a distribution over trajectories.  Similarly, the
control goal in \cite{Todorov06,Todorov09} is microscopic (to reach a certain
subset of the state space $\mathcal{X}$), rather than macroscopic (to maintain a
target distribution $\pi^*$ over $\mathcal{X}$).

\cite{Gopalkrishnan16} considers the problem of erasing a bit of information
encoded in the stationary distribution of a two-state continuous time Markov
chain and uses the KL divergence to measure the cost of control, as does our
work.  Whereas our control goal is to hold a target distribution $\pi^*$ and
minimize the cost per unit time, \cite{Gopalkrishnan16}'s control goal is to
have (in the notation of section \ref{sec:MC_continuous_time}) $\mu_T = \pi^*$
by a fixed time $T$ and to minimize the total cost used to achieve this.
Consequently, \cite{Gopalkrishnan16} uses a time-varying controlled chain,
while ours is constant in time.

\section{Minimum-power controlled chains} \label{sec:minimum_power}

In this section we minimize the power used to hold a desired nonequilibrium
stationary distribution $\pi^*$.

\begin{theorem} (Minimum-power chain) \label{thm:I_proj_stat_distrib} Let
    $\pi^*$ be a distribution on finite set $\mathcal{X}$ with $\pi^*_i > 0$ for all
    $i \in \mathcal{X}$.  Let $\pi(P)$ and $\pi(\bar{P})$ denote the stationary
    distributions of discrete- and continuous-time chains $P$ and $\bar{P}$,
    respectively.
    \begin{enumerate}
        \item (Discrete time)

            Let $Q$ be an irreducible, aperiodic transition probability matrix
            (the uncontrolled discrete-time chain), let $\pi = \pi(Q)$, and let
            $P^*$ denote the minimum-power controlled chain with the desired
            stationary distribution $\pi^*$:
    \begin{equation}
        P^* = \argmin_{P: \pi(P) = \pi^*} D(P||Q)
        \label{eq:def_I_proj_stat_distrib_disc_time}
    \end{equation}
    where the minimum is over all transition probability matrices with the
            desired stationary distribution and where $D(P||Q)$ is as defined
            in (\ref{eq:def_KL_div_MC_discrete_time_1}).  Then $P^*$ exists, is
            unique, and satisfies for all $i,j\in\mathcal{X}$:
    \begin{equation}
        P^*_{ij} = Q_{ij} e^{\lambda_i + \eta_j}
        \label{eq:I_proj_stat_distrib_dic_time_element}
    \end{equation}
    where $(\lambda_i)_{i \in \mathcal{X}}, (\eta_i)_{i \in \mathcal{X}}$ are real-valued
            constants satisfying the recursive relations:
    \begin{align}
        \lambda_i &= - \log\bigg(\sum_{j \in \mathcal{X}} Q_{ij} e^{\eta_j}\bigg)
        \label{eq:eta_disc_time} \\
        \eta_i &= \log(\pi^*_i) - \log\bigg(\sum_{j \in \mathcal{X}} \pi_j Q_{ji}
        e^{\lambda_j}\bigg) \label{eq:lambda_disc_time}
    \end{align}
\item (Continuous time)
    Let $\bar{Q}$ be a transition rate matrix (the uncontrolled continuous-time
            chain) with $e^{\bar{Q}}$ irreducible and aperiodic, and let
            $\bar{P}^*$ denote the minimum-power controlled chain with the
            desired stationary distribution $\pi^*$:
    \begin{equation}\label{eq:minpowerchain}
        \bar{P}^* = \argmin_{\bar{P}: \pi(\bar{P}) = \pi^*} D(\bar{P}||\bar{Q})
    \end{equation}
    where the minimum is over all transition rate matrices with the desired
            stationary distribution and where $D(\bar{P}||\bar{Q})$ is as
            defined in (\ref{eq:KL_div_MC_continuous_time_expression}).  Then
            $\bar{P}^*$ exists, is unique, and
    satisfies:
    \begin{equation}
        \bar{P}^*_{ij} = \begin{cases}
            \bar{Q}_{ij} e^{\lambda_i - \lambda_j} &: i \neq j \\
            -\sum_{j \in \mathcal{X}: j \neq i} \bar{P}^*_{ij} &: i = j
            \label{eq:def_I_proj_stat_distrib_cont_time}
        \end{cases}
    \end{equation}
    where $(\lambda_i)_{i \in \mathcal{X}}$ are real-valued constants satisfying the
            recursive relations:
    \begin{align}
        \lambda_i = \frac{1}{2} \log\left(\frac{\sum_{j \in \mathcal{X}: j \neq i}
        \pi^*_j \bar{Q}_{ji} e^{\lambda_j}}{\sum_{j \in \mathcal{X}: j \neq i}
        \pi^*_i \bar{Q}_{ij} e^{-\lambda_j}}\right) \label{eq:lambda_cont_time}
    \end{align}
    \end{enumerate}
\end{theorem}

Existence and uniqueness of $P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time})
follow as a special case of Lemma 1 of \cite{CsiszarCoverChoi87}.  Existence
and uniqueness of $\bar{P}^*$ (\ref{eq:def_I_proj_stat_distrib_cont_time}) were
shown in \cite{BaldiPiccioni98}.  We prove expressions
(\ref{eq:I_proj_stat_distrib_dic_time_element}) and
(\ref{eq:def_I_proj_stat_distrib_cont_time}) for $P^*_{ij}$ and
$\bar{P}^*_{ij}$
by setting up a Lagrange multiplier optimization problem, where $\alpha,
(\lambda_i)_{i \in \mathcal{X}}, (\eta_i)_{i \in \mathcal{X}}$ are Lagrange multipliers.
See Appendix \ref{app:proof} for proof details.

The recursive relations (\ref{eq:eta_disc_time}), (\ref{eq:lambda_disc_time}),
and (\ref{eq:lambda_cont_time}) enable an iterative computation of the chains
$P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time})
and $\bar{P}^*$ (\ref{eq:def_I_proj_stat_distrib_cont_time}).  In the
continuous-time case, for example, we initialize $(\lambda_i^{(0)})_{i \in
\mathcal{X}}$ to some value and then use relation (\ref{eq:lambda_cont_time}) to
compute $\lambda^{(t)}_i$ as a function of $(\lambda^{(t-1)}_i)_{i \in \mathcal{X}}$
at the $t$-th iteration until numerical convergence.  The example in section
\ref{sec:examples} is computed in this way.

The chains $P^*$ and $\bar{P}^*$ are the I-projections of the chains $P$ and
$\bar{P}$ on the set of discrete- and continuous-time Markov chains,
respectively, with a fixed stationary distribution.  The discrete-time case is
presented in \cite{DemboZeitouniBook,CsiszarCoverChoi87} and the
continuous-time case in \cite{BaldiPiccioni98}, where $D(P^*||Q)$ and
$D(\bar{P}^*||\bar{Q})$ arise as large deviations rate functions for the
empirical marginal distribution.
To our knowledge Theorem \ref{thm:I_proj_stat_distrib} presents the most
explicit characterization of the I-projection in terms of the Lagrange
multipliers.

We next specialize our results to the case of reversible uncontrolled Markov
chains, a case important in equilibrium thermodynamics.
Let $Q^\text{r} = (Q^\text{r}_{ij})_{i,j\in\mathcal{X}}$ denote the \textit{time-reverse} of a
transition probability matrix $Q$.  That is:
\begin{equation}
    Q^\text{r}_{ij} = \frac{\pi_j}{\pi_i} Q_{ji} \label{eq:time_rev_chain_disc_time}
\end{equation}
where $\pi = \pi(Q) = \pi(Q^\text{r})$ is the stationary distribution.  Then
$Q^\text{r}_{ij} = \lim_{t \rightarrow \infty} \mathbb{P}_Q(X_t=j|X_{t+1}=i)$; if $X_0 \sim
\pi$, then $Q^\text{r}_{ij} = \mathbb{P}_Q(X_t=j|X_{t+1}=i)$ for all $t$.   Analogously,
in continuous time, the time-reverse $\bar{Q}^\text{r}$ of a transition rate matrix
$\bar{Q}$ satisfies $\bar{Q}_{ij} = \frac{\pi_j}{\pi_i} \bar{Q}_{ji}$ for all
$i,j\in\mathcal{X}$.  A chain $Q$ is \textit{reversible} if $Q = Q^\text{r}$
(analogously, $\bar{Q} = \bar{Q}^\text{r}$ in continuous time).

\begin{theorem} \label{thm:I_proj_stat_distrib_reversible_base_chain}
    (Reversible uncontrolled chain) Let notation be as in the statement of
    Theorem \ref{thm:I_proj_stat_distrib}.

    \begin{enumerate}
        \item If the uncontrolled transition probability matrix $Q$ is
            reversible, then so is the minimum-power chain $P^*$
            (\ref{eq:def_I_proj_stat_distrib_disc_time}).  Analogously, in
            continuous time, if the uncontrolled transition rate matrix
            $\bar{Q}$ is reversible, then so is $\bar{P}^*$
            (\ref{eq:def_I_proj_stat_distrib_cont_time}).
        \item In discrete time, the Lagrange multipliers
            (\ref{eq:eta_disc_time}), (\ref{eq:lambda_disc_time}) satisfy
            $\eta_i = \lambda_i + \log(\pi^*_i/\pi_i)$ for all $i \in \mathcal{X}$.
        \item In continuous time, $\bar{P}^*$ satisfies for $i \neq j$:
            \begin{equation}
                \bar{P}^*_{ij} = \bar{Q}_{ij} \sqrt{\frac{\pi_i \pi^*_j}{\pi_j
                \pi^*_i}} \label{eq:reversible_cont_time}
            \end{equation}
            and
            \begin{equation}
                D(\bar{P}^*||\bar{Q}) = \frac{1}{2}\sum_{i,j\in\mathcal{X}:i\neq j}
                \left(\sqrt{\pi^*_i \bar{Q}_{ij}} - \sqrt{\pi^*_j
                \bar{Q}_{ji}}\right)^2 \label{eq:KL_div_rev_cont_time}
            \end{equation}
    \end{enumerate}
\end{theorem}
\begin{proof}
    1) We can check that if $Q = Q^\text{r}$, then $D(P||Q) = D(P^\text{r}||Q)$ for all
    $P$.  Suppose that $P^*$ is not reversible.  Let $P^{\leftrightarrow}
    \triangleq (P^*+P^{*\text{r}})/2$.  Then $P^{\leftrightarrow}$ is reversible
    and $\pi(P^{\leftrightarrow}) = \pi(P^*) = \pi^*$.  Since $D(P||Q)$ is
    strictly convex in $P$, we have $D(P^{\leftrightarrow}||Q) < (D(P^*||Q) +
    D(P^{*\text{r}}||Q))/2 = D(P^*||Q)$, contradicting the optimality of $P^*$.
    Another proof: suppose that $P^*$ is not reversible, then $D(P^*||Q) =
    D(P^{*\text{r}}||Q)$, contradicting the uniqueness of $P^*$ established in
    Theorem \ref{thm:I_proj_stat_distrib}.  Therefore $P^*$ is reversible.  An
    analogous argument proves $\bar{P}^*$ is reversible.  2) follows by using
    time reversal (\ref{eq:time_rev_chain_disc_time}) twice along with the
    reversibility of $Q$ and $P^*$, established in part 1):
    \begin{align}
        \pi^*_i Q_{ij}e^{\lambda_i + \eta_j} &= \pi^*_i P^*_{ij} = \pi^*_j
        P^*_{ji} = \pi^*_j Q_{ji} e^{\lambda_j + \eta_i} \\
        &= \pi^*_j \frac{\pi_i}{\pi_j} Q_{ij} e^{\lambda_j + \eta_i}
    \end{align}
    and collecting $i$- and $j$-dependent terms to separate sides of the
    equality to conclude that $\frac{\pi^*_i}{\pi_i} e^{\lambda_i - \eta_i} =
    a$ for some constant $a$ for all $i \in \mathcal{X}$.  Choosing $a=1$ yields the
    result.

    3) Expression (\ref{eq:reversible_cont_time}) for $\bar{P}^*$ is derived in
    \cite{BaldiPiccioni98} and $D(\bar{P}^*||Q)$
    (\ref{eq:KL_div_rev_cont_time}) is derived as a large deviations rate
    function in \cite{BertiniFaggionatoGabrielli14}.
\end{proof}

Part 2) of Theorem \ref{thm:I_proj_stat_distrib_reversible_base_chain} lets us
simplify computation of $P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time})
somewhat when the uncontrolled chain $Q$ is reversible.

\section{Examples} \label{sec:examples}

We conclude with numerical examples of minimum-power controlled Markov chains
with a target stationary distribution.  The first example is a two-state chain
in discrete time and the second example is a birth-and-death chain in
continuous and discrete time -- a toy model of the muscle fiber thought
experiment in the introduction (section \ref{sec:introduction}).

\subsection{Two-state chain in discrete time}

Let $Q$ be a two-state discrete-time Markov chain on set $\mathcal{X} = \{1,2\}$ and
let $\pi^* = (\pi^*_1, 1-\pi^*_1)$ be our desired nonequilibrium distribution
with $\pi^*_1 \in (0,1)$.  All two-state chains are reversible, so we apply
Theorem \ref{thm:I_proj_stat_distrib_reversible_base_chain} part 2) to compute
the minimum-power controlled chain with stationary distribution $\pi^*$.  A
computation shows the $2 \times 2$ minimum-power transition matrix $P^*$
(\ref{eq:def_I_proj_stat_distrib_disc_time}) has off-diagonal entries:
    \begin{equation}
        P^*_{ij} = \frac{1}{\pi^*_i} \cdot \frac{1-\sqrt{1-4 \pi^*_1
        (1-\pi^*_1) s}}{2 s} \ : i \neq j
    \end{equation}
    where the second factor is independent of $i,j$ and
    \begin{equation}
        s \triangleq \frac{1}{Q_{1,2}} + \frac{1}{Q_{2,1}} - \frac{1}{Q_{1,2}
        Q_{2,1}}
    \end{equation}

    The diagonal terms of $P^*$ are such that the row sums are $1$.

\subsection{Birth-and-death chain}

We next present the example of the birth-and-death chain as a toy model of
Feynman's muscle fiber thought experiment (see section \ref{sec:introduction}).
For detailed models of molecular motors see \cite{Julicher} and \cite{Qian00}.

\subsubsection{Continuous time} \label{sec:cont_time_birth_and_death_chain}

\begin{figure}[!t]
    \centering
    \includegraphics[width=3in]{figures/birthanddeathchain.pdf}
    \caption{Uncontrolled chain $Q$ (a toy model for the slack arm pulled down
    by gravity, with smaller-index states closer to the ground and rates $q_- >
    q_+$): the birth-and-death chain with 5 states in continuous time.  The
    states correspond to minima in the potential landscape experienced by a
    single protein at different positions along the muscle fiber.}
    \label{fig:birth_and_death_chain}
\end{figure}

Let $\bar{Q}$ be the continuous-time birth-and-death chain on set
$\{1,2,\ldots,K\}$ with parameters $q_+, q_- > 0$ depicted in Figure
\ref{fig:birth_and_death_chain} for $K = 5$.  The chain transitions increments
from state $i$ to $i + 1$ (resp. decrements to $i-1$) with rate $q_+$ (resp.
$q_-$).  All other transitions, as well as decrementing from state $1$ and
incrementing from state $K$, have rate $0$.  $\bar{Q}$ is reversible and its
stationary distribution is, up to normalization, \cite{LevinPeresWilmerBook}:
\begin{equation}
    \pi_i \sim \left(\frac{q_+}{q_-}\right)^i
    \label{eq:stat_distrib_birth_and_death_chain}
\end{equation}

Let our control objective be to maintain the target distribution $\pi^*$, a
geometric distribution on $\{1,\ldots,K\}$:
\begin{equation}
    \pi^*_i \sim b^i \label{eq:target_distribution_birth_and_death_chain}
\end{equation}
where $b > 0$.  Then applying Theorems \ref{thm:I_proj_stat_distrib} and
\ref{thm:I_proj_stat_distrib_reversible_base_chain} we find the minimum-power
controlled chain $\bar{P}^*$ (\ref{eq:minpowerchain}) with stationary
distribution $\pi^*$ to be another birth-and-death chain with increment and
decrement rates $p_+^*, p_-^*$:
\begin{equation}
    p_{\pm}^* = (q_+ q_-)^{\frac{1}{2}} b^{\pm \frac{1}{2}}
    \label{eq:birth_and_death_chain_cont_time_parameters}
\end{equation}
The cost per unit time of this birth-and-death chain $\bar{P}^*$ is
(\ref{eq:KL_div_rev_cont_time}):
\begin{equation}

    D(\bar{P}^*||\bar{Q}^*) = \left(\sqrt{q_+} - \sqrt{q_- b}\right)^2
    \left(\frac{1-b^{K-1}}{1-b^K}\right)
    \label{eq:birth_and_death_chain_cont_time_cost}
\end{equation}
If $b < 1$, then as $K \rightarrow \infty$, we have $D(\bar{P}^*||\bar{Q}) \rightarrow
\left(\sqrt{q_+} - \sqrt{q_- b}\right)^2$.

\begin{figure}[!t]
    \centering

    \includegraphics[width=3in]{figures/periodic_potential.pdf}
    \caption{A tilted periodic potential energy landscape.

    We consider the $i$-th local minimum of the potential as a discrete
    position state with energy $U_i$ and energy difference $B = U_{i+1} - U_i$
    between adjacent states.  The activation energy to go left (resp. right) is
    $A$ (resp. $A+B$).}
    \label{fig:energy_potential_snapshot}
\end{figure}

Recalling the motivating example of section \ref{sec:introduction}, we can
think of the birth-and-death chain as a biased random walk, where the random
walker tends to spend more time at small values of $i$ if $q_- > q_+$.  This is
a toy model of a myosin protein moving on an actin filament in muscles -- a
random walk in a tilted energy potential with periodically spaced minima
corresponding to discrete steps along the fiber.  The energy potential is
depicted in Figure \ref{fig:energy_potential_snapshot}: a transition from the
$i+1$-th energy minimum to the $i$-th energy minimum must overcome activation
energy $A$, and the reverse transition must overcome activation energy $B+A$,
where $B$ is the energy difference between adjacent states.  In the bag-holding
thought experiment, $B$ is the gravitational potential energy difference
between adjacent states, with state $1$ being closest to the ground.

In terms of these energies, the increment and decrement rates $q_+, q_-$ are:
\begin{align}
    q_- &= c e^{-\beta A} \\
    q_+ &= c e^{-\beta (A + B)}
\end{align}
for some constant $c$, where $\beta$ is the inverse temperature.  The
stationary distribution of the uncontrolled chain $\bar{Q}$ is, up to
normalization:
\begin{equation}
    \pi_i \sim \left(\frac{q_+}{q_-}\right)^i = e^{-\beta B i}
\end{equation}

Let's write the target nonequilibrium distribution as:
\begin{equation}
    \pi^*_i \sim b^i = e^{-\beta B^* i}
    \label{eq:target_distribution_cont_time_birth_and_death_chain}
\end{equation}
with energy difference $B^* = (-1/\beta) \log(b)$ between adjacent states.

Since the optimal controlled chain $\bar{P}^*$ is another birth-and-death
chain, we can write its parameters $p^*_+, p^*_-$
(\ref{eq:birth_and_death_chain_cont_time_parameters}) in terms of energies
$A^*$ and $B^*$:
\begin{align}
    p^*_- &= c e^{-\beta A^*} \\
    p^*_+ &= c e^{-\beta (A^* + B^*)}
\end{align}
where, using (\ref{eq:birth_and_death_chain_cont_time_parameters}), we find the
controlled activation energy $A^*$:
\begin{equation}
    A^* = A + \frac{1}{2}\left(B - B^*\right)
\end{equation}
The cost per unit time of this birth-and-death chain $\bar{P}^*$ is
(\ref{eq:birth_and_death_chain_cont_time_cost}):
\begin{equation}

    D(\bar{P}^*||\bar{Q}^*) = c e^{-\beta A} \left(e^{-\frac{\beta}{2} B} -
    e^{-\frac{\beta}{2} B^*}\right)^2 \left(\frac{1-e^{-\beta B^*
    (K-1)}}{1-e^{-\beta B^* K}}\right)
\end{equation}
where if $B^* > 0$, then the last factor tends to $1$ as $K \rightarrow \infty$.

In the muscle fiber thought experiment (where $B$ is the gravitational
potential energy difference between adjacent states and state $0$ is closest to
the ground) if $B^* < B$, then target distribution $\pi^*$
(\ref{eq:target_distribution_cont_time_birth_and_death_chain}) corresponds to
imposing a constant force upwards (away from state $1$) on the random-walking
myosin protein.  The control objective is macroscopic: rather than control the
microscopic trajectory of a single myosin protein, we imagine controlling a
large collection of identical, independent myosin proteins in the same way by
imposing the same controlled chain $\bar{P}^*$ on all myosins; the bag-holder's
arm position is determined by an average over the positions of this collection
of myosins.

\subsubsection{Discrete time} \label{sec:discrete_time_birth_and_death_chain}

Let $Q$ be the discrete-time birth-and-death chain on set $\mathcal{X} =
\{1,\ldots,K\}$ with transition probability $q_+$ (resp. $q_-$) to increment
(resp. decrement) the state from $i$ to $i+1$ (resp. $i-1$).  The stationary
distribution $\pi$ of $Q$ is as in
(\ref{eq:stat_distrib_birth_and_death_chain}),
the same as in the continuous time case with transition rates $q_+, q_-$.  Let
our control objective be to maintain the target distribution $\pi^*_i \sim
b^i$, a geometric distribution on $\{1,\ldots,K\}$
(\ref{eq:target_distribution_birth_and_death_chain}).

Then in contrast to the continuous time case of section
\ref{sec:cont_time_birth_and_death_chain}, the minimum-power controlled chain
$P^*$ (\ref{eq:def_I_proj_stat_distrib_disc_time}) is not in general a
birth-and-death chain.  Consider this numerical example: let the target
nonequilibrium stationary distribution be as in
(\ref{eq:target_distribution_birth_and_death_chain}) with $b = \sqrt{q_-/q_+}$:
\begin{equation}
    \pi^*_i \sim \sqrt{\pi_{K-i}} \sim \left(\frac{q_+}{q_-}\right)^{-i/2}
    \label{eq:target_distribution_birth_and_death_chain_discrete_time}
\end{equation}
$\pi^*$ (\ref{eq:target_distribution_birth_and_death_chain_discrete_time}) is
biased the other way from $\pi_i$, assigning most of its mass to large values
of $i$ if $q_+ < q_-$.  The square root in
(\ref{eq:target_distribution_birth_and_death_chain_discrete_time}) makes
$\pi^*$ look more uniform than $\pi$.

Let $q_- = 0.2$, $q_+ = 0.1$, so that the stationary distribution of $Q$ is
$\pi_i \sim 2^{-i}$, and $\pi^*_i \sim 2^{i/2}$.  Figure
\ref{fig:birth_and_death_chain_I_proj} shows the non-zero off-diagonal elements
of $P^*_{ij}$ with three-digit precision; the increment and decrement
probabilities vary with state $i$, so $P^*$ is not a birth-and-death chain.

\begin{figure}[!t]
    \centering
    \includegraphics[width=3in]{figures/birthanddeathchain_I_proj.pdf}
    \caption{The minimum-power discrete-time controlled chain $P^*$
    (\ref{eq:def_I_proj_stat_distrib_disc_time}) for the uncontrolled
    birth-and-death chain $Q$ with parameters $q_- = 0.2$, $q_+ = 0.1$ and
    target stationary distribution $\pi^*_i \sim 2^{i/2}$.  Transitions from a
    state to itself are not shown.}
    \label{fig:birth_and_death_chain_I_proj}
\end{figure}

Finally, Figure \ref{fig:birth_and_death_chain_cost} depicts the time evolution
of distribution $\mu_t = \mu_0 P^{*t}$ with $\mu_0 = \pi$, and the cost
$D_{\mu_t}(P^*||Q)$ converging to the minimum power $D(P^*||Q) \approx 0.0315$
to maintain the nonequilibrium distribution $\pi^*$.

\begin{figure}[!t]
    \centering
    \includegraphics[width=3.8in]{figures/birthanddeathchaincost.pdf}
    \caption{(top) Time evolution of the $|\mathcal{X}| = 5$ components of $\mu_t =
    \mu_0 (P^{*})^t$ with $\mu_0 = \pi$, where $P^*$ is the minimum cost chain
    (\ref{eq:def_I_proj_stat_distrib_disc_time}), showing $\mu_t \rightarrow \pi^*$ as
    $t \rightarrow \infty$.  (bottom) The cost (blue) $D_{\mu_t}(P^*||Q)$
    (\ref{eq:cost_disc_time_KL_div}) and the minimum power (red) $D(P^*||Q)$.}
    \label{fig:birth_and_death_chain_cost}
\end{figure}

Returning to the molecular motor picture, the discrete-time chain $P^*$ differs
from its continuous-time cousin of section
\ref{sec:cont_time_birth_and_death_chain} in that $P^*$ is not a
birth-and-death chain.  Thus an optimal discrete-time control policy modifies
the base birth-and-death chain $Q$ in a way that depends on state $i$, and so
can't be thought of as corresponding to a constant, state-independent force
upwards (away from state $1$) as in the continuous time case.

\section*{Acknowledgment}

The authors gratefully acknowledge Hideo Mabuchi for suggestions and insightful discussions.

\bibliographystyle{IEEEtran}