# Computing transition path theory quantities with trajectory stratification

Bodhi P. Vani  Department of Chemistry and James Franck Institute, University
of Chicago, Chicago IL 60637  Jonathan Weare  Courant Institute of
Mathematical Sciences, New York University, New York 10012  Aaron R. Dinner∗
Department of Chemistry and James Franck Institute, University of Chicago,
Chicago IL 60637<EMAIL_ADDRESS>
###### Abstract

Transition path theory computes statistics from ensembles of reactive
trajectories. A common strategy for sampling reactive trajectories is to
control the branching and pruning of trajectories so as to enhance the
sampling of low probability segments. However, it can be challenging to apply
transition path theory to data from such methods because determining whether
configurations and trajectory segments are part of reactive trajectories
requires looking backward and forward in time. Here, we show how this issue
can be overcome efficiently by introducing simple data structures. We
illustrate the approach in the context of nonequilibrium umbrella sampling
(NEUS), but the strategy is general and can be used to obtain transition path
theory statistics from other methods that sample segments of unbiased
trajectories.

††preprint: AIP/123-QED

##  I Introduction

One of the main purposes of simulations is to learn how processes occur.
However, most processes of interest occur on timescales that are orders of
magnitude longer than the numerical integration time step. In the case of
molecular systems, on which we focus here, the time step for all-atom models
is in the femtosecond range, while conformational transitions of functional
significance typically take place in the microsecond to seconds range. As a
result, at best, only a small number of events can be observed by direct
simulation. Moreover, even if one obtains examples of events of interest, the
number of dynamical variables is generally sufficiently large that it is not
obvious which are the key steps Ma and Dinner (2005), leave alone their
likelihoods and the distributions of times that they occur.

Transition Path Theory (TPT) provides a means of computing these statistics
for transitions between two metastable states from the ensemble of reactive
trajectories (i.e., those that contain an event) Vanden-Eijnden (2006); E. and
Vanden-Eijnden (2006); Metzner, Schütte, and Vanden-Eijnden (2006). A key idea
in TPT is that these statistics can be expressed as products over the steady-
state distribution and the probabilities of coming from or going to the
metastable states (known as commitment probabilities or committors). This
factorization allows estimating these probabilities separately, which opens
the door to computing statistics of reactive trajectories from a broader range
of data.

Markov State Models (MSMs) are a popular means of computing TPT statistics Noé
_et al._ (2009); Bowman, Pande, and Noé (2013). In MSMs, dynamics are coarse-
grained to transitions between a set of discrete states, the probabilities of
which are assumed to be independent of the sequence of states visited. The
advantage relative to direct simulation comes from the fact that there is
considerable freedom in how one estimates the state-to-state transition
probabilities Wu  _et al._ (2017); Buchete and Hummer (2008). Dynamical
Galerkin Approximation (DGA) Thiede  _et al._ (2019); Strahan  _et al._ (2021)
is a generalization of MSMs for TPT statistics that accounts for the boundary
conditions of the statistics Russo  _et al._ (2021) and represents them
through expansions over sets of basis functions (which can take a form beyond
indicator functions on discrete states). Recently, TPT statistics have also
been computed with milestoning Elber  _et al._ (2017), which instead assumes
that the system reaches an equilibrium within selected parts of the
configuration space (the milestones) Faradjian and Elber (2004); Elber (2020).

Although often reasonable, the above approximation schemes can break down,
and, even if they do not, they ultimately limit the resolution with which
statistics can be computed. Methods that sample the ensemble of trajectories
that connect the metastable states without such assumptions can in principle
provide results with arbitrary resolution. One such method is Transition Path
Sampling (TPS), a Monte Carlo procedure that accepts and rejects whole
trajectories Dellago  _et al._ (1998); Dellago, Bolhuis, and Chandler (1998);
Dellago  _et al._ (2002). However, generally it is more efficient to divide
the configuration space and sample segments of trajectories that transition
between regions. This idea is at the foundation of the Weighted Ensemble (WE)
Huber and Kim (1996); Zhang, Jasnow, and Zuckerman (2007); Zuckerman and Chong
(2017), Transition Interface Sampling (TIS) Bolhuis (2003); van Erp, Moroni,
and Bolhuis (2003), Forward Flux Sampling (FFS) Allen, Warren, and Ten Wolde
(2005); Allen, Valeriani, and ten Wolde (2009), Nonequilibrium Umbrella
Sampling (NEUS) Warmflash, Bhimalapuram, and Dinner (2007); Dickson,
Warmflash, and Dinner (2009a); Vanden-Eijnden and Venturoli (2009); Dickson
and Dinner (2010); Dinner  _et al._ (2020), and Exact Milestoning (EM) Bello-
Rivas and Elber (2015) methods. With the exception of TIS, these methods allow
treatment of microscopically irreversible dynamics, which is also not
straightforward in TPS.

The differences between these methods are in the details of the trajectory
segments they select and how they track their probabilities to enable
reconstruction of the overall statistics. Generally the interfaces are
specified through collective variables that combine information from multiple
coordinates used for the underlying numerical integration. TIS and FFS require
non-intersecting interfaces because they are formulated in terms of the
probabilities of going from one to another in order. WE, NEUS, and EM do not
have this restriction and thus readily allow control of sampling in more than
one collective variable. In its traditional form, WE tracked the weight of
each trajectory segment (equivalent to a member of the ensemble) individually.
NEUS instead tracked the weight of a subpopulation within a region (termed a
stratum in statistics) and redistributed the weight to satisfy a global flux
balance condition. Later implementations of WE that incorporate this idea
Bhatt, Zhang, and Zuckerman (2010) are very similar to NEUS, as is EM, which
relaxes the distributions on the milestones. Several of these algorithms can
thus be described by a common framework, trajectory stratification Dinner  _et
al._ (2018).

Because these methods sample trajectory segments rather than whole
trajectories, care must be used to weight trajectory segments appropriately
when computing statistics. Vanden-Eijnden and Venturoli Vanden-Eijnden and
Venturoli (2009) showed that the rate can be obtained exactly from NEUS by
augmenting the dynamics with information about the last metastable state
visited, and this scheme was further refined and applied in ref. Dickson,
Warmflash, and Dinner, 2009a. However, the rate requires only tracking
transitions into the metastable states, making it easier to compute than
statistics that require information about intervening states. Vanden-Eijnden
and Venturoli Vanden-Eijnden and Venturoli (2009) provided a formula for
computing the (backward) committor from the region weights, but, to the best
of our knowledge, statistics such as committors and probability currents of
reactive trajectories (reactive currents for short), which reveal microscopic
likelihoods of reaction and reaction pathways, respectively, have not been
computed from such methods. The challenge presented by these quantities is
that each trajectory segment can contribute to an infinite number of complete
trajectories, and determining the ones that contribute to a given statistic
requires looking both forward and backward in time.

Here, we introduce a simple bookkeeping procedure and practical formulas for
computing committors and reactive currents. While we present our methods in
terms of the NEUS algorithm, they can be adapted to the other algorithms for
sampling trajectory segments described above. The paper is organized as
follows. In Section II, we define these quantities in the TPT framework, and
in an associated appendix we provide an explicit justification for an
estimator of the current. In Section III we provide a succint description of
the NEUS algorithm and show how it must be modified. We demonstrate the
algorithm for isomerization of a peptide in Section IV.

##  II Transition path theory

We are interested in learning the statistics of reactive trajectories. In TPT,
these are defined to be trajectories that start in a state $A$ and end in a
non-overlapping state $B$ without returning to $A$. To delimit their ends, we
define

$\displaystyle t_{+}(t)$ $\displaystyle=\min(t^{\prime}>t,\ X_{t^{\prime}}\in
A\cup B)$ (1) $\displaystyle t_{-}(t)$ $\displaystyle=\max(t^{\prime}<t,\
X_{t^{\prime}}\in A\cup B),$ (2)

where $X_{t^{\prime}}$ is the configuration at time $t^{\prime}$, and
$X_{t}\in(A\cup B)^{c}$. The forward committor $q_{+}(x)$ is the probability
that the system goes to $B$ before $A$ starting from state $x$:

$q_{+}(x)={\bf P}[X(t_{+}(0))\in B|X(0)=x].$ (3)

Similarly, the backward committor $q_{-}(x)$ is the probability that a system
at $x$ came from $A$ after $B$:

$q_{-}(x)={\bf P}[X(t_{-}(0))\in A|X(0)=x].$ (4)

In the present study, for visualization, we project the committors onto the
space defined by a vector of collective variables, $\theta(x)$, weighted by
$\pi(x)$:

$q_{\pm}^{\theta}(\Theta)=\int q_{\pm}(x)\pi(x)\delta(\theta(x)-\Theta)dx,$
(5)

A key element of TPT is that the probability that a trajectory is reactive can
be expressed in terms of quantities that are local in space. Specifically, the
probability that a trajectory that passes through $x$ is reactive is
proportional to $\pi(x)q_{-}(x)q_{+}(x)$, where $\pi(x)$ is the steady-state
distribution. The fact that TPT is based on local statistics aids in
interpretation, but it also can make obtaining certain statistics challenging.
We discuss how TPT of an augmented process can be used to obtain additional
statistics that require knowledge of sequences of states visited in ref.
Lorpaiboon, Weare, and Dinner, 2022. Of course, the analysis we propose here
can be complemented by direct analysis of the reactive trajectories as well.

###  II.1 Reactive current

The reactive current is a central quantity in TPT Vanden-Eijnden (2006); E.
and Vanden-Eijnden (2006); Metzner, Schütte, and Vanden-Eijnden (2009). It can
be used to define reactive pathways, and averages over reactive
trajectories—in particular the rate—can be computed from it. The reactive
current is defined by considering a surface $S$ that divides the space into
two regions: $C$ containing $A$ and its complement $C^{c}$ containing $B$.
Then, the reactive current, ${\color[rgb]{0,0,0}{I}}_{AB}(x)$, is the vector
field whose integral is the reactive flux across $S$:

$\int_{S}{\color[rgb]{0,0,0}{I}}_{AB}\cdot\hat{n}_{S}d\sigma_{S}=\lim_{\tau\to
0}\lim_{T\to\infty}\frac{1}{2T\tau}\int_{-T}^{T}[\mathbf{1}_{C}(X(t))\mathbf{1}_{C^{c}}(X(t+\tau))\\\
-\mathbf{1}_{C^{c}}(X(t))\mathbf{1}_{C}(X(t+\tau))]\mathbf{1}_{A}(X(t_{-}(t)))\mathbf{1}_{B}(X(t_{+}(t+\tau)))dt$
(6)

where $\hat{n}_{S}$ is a unit vector normal to $S$ (pointing from $C$ to
$C^{c}$), $d\sigma_{S}$ is a surface element, $\tau$ is a lag time, and

$\mathbf{1}_{D}(x)=\begin{cases}1&\text{if }x\in D\\\
0&\text{otherwise}.\end{cases}$ (7)

In general, the full state space of the system has many dimensions.
${\color[rgb]{0,0,0}{I}}_{AB}(x)$ is a vector field in the full space of
dynamical variables, of which there are many in general. However, most of
these are irrelevant for determining whether the reaction proceeds, and we
seek a small number of key dimensions. Because these are often combinations of
the original dynamical variables, we term them collective variables. We show
in the Supplemental Information that

${\color[rgb]{0,0,0}{I}}_{AB}\cdot\nabla\theta(x)=\pi(x)q_{-}(x)q_{+}(x)\\\
\times\lim_{\tau\to 0}\frac{1}{2\tau}{\bf
E}[(\theta(X(\tau))-\theta(X(-\tau))|\\\ X(t_{-}(0))\in A,X(t_{+}(0))\in
B,X(0)=x].$ (8)

We use $\bf E$ to denote expectations over the steady-state trajectory
ensemble. This formula is simply the product of the probability of being on a
reactive trajectory and the average increment per unit time of $\theta$ along
a reactive trajectory. We project this current onto the space of collective
variables:

$\displaystyle{\color[rgb]{0,0,0}{I}}_{AB}^{\theta}(\Theta)$
$\displaystyle=\int{\color[rgb]{0,0,0}{I}}_{AB}(x)\cdot\nabla\theta(x)\delta(\theta(x)-\Theta)dx.$
(9) $\displaystyle=\lim_{\lvert d\Theta\rvert\rightarrow 0}\frac{1}{\lvert
d\Theta\rvert}\int_{\\{\theta(x)\in
d\Theta\\}}{\color[rgb]{0,0,0}{I}}_{AB}(x)\cdot\nabla\theta(x)dx.$ (10)

We recently showed that this projected reactive current can also be used to
compute the reactive flux through an integral similar in form to (6) but in
the $\theta$ space; the projected current
${\color[rgb]{0,0,0}{I}}_{AB}^{\theta}(\Theta)$ can in turn be used to compute
the rate Strahan  _et al._ (2021). The expressions in (8) and (9) allow us to
obtain the reactive current in the space of collective variables from reactive
trajectories.

In the next section, we discuss how reactive trajectories can be sampled
efficiently. Although we treat time as continuous above, with a view toward
developing a practical algorithm, we henceforth assume a dynamics with
discrete time steps; the lag time $\tau$ is thus an integer multiple of the
time step in practice.

##  III Trajectory stratification

The idea of stratification in general is to obtain better overall statistics
by controlling the sampling of subpopulations (strata). The most familiar form
of stratification in molecular simulations is umbrella sampling Torrie and
Valleau (1977); Pangali, Rao, and Berne (1979); Frenkel and Smit (2001);
Thiede  _et al._ (2016); Dinner  _et al._ (2020), which is frequently used to
compute free energies. In this case, the strata are defined by regions in a
space of collective variables. Each of the strata is sampled independently by
a copy of the system, and the information is then combined to obtain unbiased
averages. Trajectory stratification extends this strategy from states to
trajectories and thus enables treating microscopically irreversible systems
and estimating dynamical statistics (for microscopically reversible or
irreversible systems) Warmflash, Bhimalapuram, and Dinner (2007); Dickson,
Warmflash, and Dinner (2009a); Vanden-Eijnden and Venturoli (2009); Dickson
and Dinner (2010); Dinner  _et al._ (2020).

Here, we focus on the case of steady-state, time-independent quantities. To
formulate the algorithm mathematically, we associate with $X(t)$ a variable
$J(t)$ that reports the index of the stratum at time $t$. For example, a
trajectory that spent one time step in stratum 1, two time steps in stratum 2,
went back to stratum 1 for one time step, and then spent three time steps in
stratum 5 would have $J(t=1)=1$, $J(t=2)=2$, $J(t=3)=2$, $J(t=4)=1$,
$J(t=5)=5$, $J(t=6)=5$, and $J(t=7)=5$. Then, denoting the time that a
trajectory first exits a stratum as $s=\min\\{t:J(t)\neq J(0)\\}$, we can
write expectations as

${\bf E}[f]=\frac{1}{\Omega}\sum_{i}z_{i}f_{i},$ (11)

where $\Omega$ is a normalization factor and

$f_{i}=\int_{x}\sum_{t}f(X(t)){\bf P}_{x,i}[X(t),t<s]\pi_{i}(dx).$ (12)

${\bf P}_{x,i}$ indicates that we are conditioning the probability on starting
at configuration $x$ in stratum $i$; $\pi_{i}(dx)\equiv\pi_{i}(x)dx$ is the
probability of being in the differential volume element $dx$ conditioned on
entering stratum $i$ at that state, and $z_{i}$ is the normalization factor
for $\pi_{i}(dx)$. Operationally, (12) corresponds to initiating trajectories
from steady-state entry points to stratum $i$, terminating them when they
leave the stratum, and averaging over the sampled points. We show in ref.
Dinner  _et al._ , 2018 that, for ergodic averages, the vector $z$ can be
obtained by solving an eigenequation:

${z}^{\text{\tiny T}}G={z}^{\text{\tiny T}}.$ (13)

The matrix $G$ tracks the number of transitions between strata; we define it
precisely below in (17). Physically, (13) can be viewed as a statement of
global flux balance.

The key idea in the present work is that we define the index process so as to
track the last metastable state visited, in the spirit of refs. Dickson,
Warmflash, and Dinner, 2009a and Vanden-Eijnden and Venturoli, 2009. In other
words, if there are $K$ strata with domains $D_{j}$ defined in terms of a set
of order parameters, then $J(t)$ runs from 1 to $2K$:

$J(t)=K\mathbf{1}[X(t_{-}(t))\in B]+\sum_{j=1}^{K}j\mathbf{1}[X(t)\in D_{j}].$
(14)

Thus strata with indices from $1$ to $K$ capture data for trajectories that
last came from state $A$, and strata with indices $K+1$ to $2K$ represent
trajectories that last came from state $B$.

###  III.1 A self-consistent iteration

Given the index process, we can now define $z_{i}$, $\pi_{i}$, and $G_{ij}$
more precisely (we note that these quantities correspond to the steady-state
averages of $\bar{z}_{i}$, $\bar{\pi}_{i}$, and $\bar{G}_{ij}$ defined in
(24), (31), and (34), respectively in ref. Dinner  _et al._ , 2018; we
suppress the overbars here for simplicity of notation):

$\pi_{i}(dx)=\lim_{T\to\infty}\frac{1}{T}\frac{\sum_{{\color[rgb]{0,0,0}{t=0}}}^{T}{\bf
P}[J(t)=i,J(t-1)\neq i,X(t)\in dx]}{z_{i}}$ (15)

and

$z_{i}=\lim_{T\to\infty}\frac{1}{T}\sum_{{\color[rgb]{0,0,0}{t=0}}}^{T}{\bf
P}[J(t)=i,J(t-1)\neq i].$ (16)

To obtain $G_{ij}$, we draw states from $\pi_{i}(dx)$ and evolve them until
they leave stratum $i$:

$\displaystyle G_{ij}$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\frac{\sum_{{\color[rgb]{0,0,0}{t=0}}}^{T}{\bf
P}[J(t)=j,J(t-1)=i]}{z_{i}}$ (17) $\displaystyle=\int{\bf
P}_{x,i}[J(s)=j]\pi_{i}(dx).$ (18)

Furthermore, we add to the distribution of entry points to stratum $j$ that
come from stratum $i$:

$\displaystyle\gamma_{ij}(dx)$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\frac{\sum_{{\color[rgb]{0,0,0}{t=0}}}^{T}{\bf
P}[J(t)=j,J(t-1)=i,X(t)\in dx]}{G_{ij}z_{i}}$ (19) $\displaystyle=\int{\bf
P}_{y,i}[J(s)=j,X(s)\in dx]\pi_{i}(dy).$ (20)

Thus, given $\pi_{i}(dx)$, we can compute $G_{ij}$ and $\gamma_{ij}(dx)$, and
in turn $z_{i}$ from (13). Because the entries of $G$ depend on the dynamics
within the strata, and those are initialized using $G$, the algorithm is
iterative in nature. To complete the iteration, we start from (15) and write
$\pi_{j}(dx)$ in terms of the other quantities:

$\displaystyle\pi_{j}(dx)$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\frac{\sum_{{\color[rgb]{0,0,0}{t=0}}}^{T}\mathbf{P}[J(t)=j,J(t-1)\neq
j,X(t)\in dx]}{z_{j}}$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\frac{\sum_{i\neq
j}\sum_{{\color[rgb]{0,0,0}{t=0}}}^{T}\mathbf{P}\left[J(t)=j,J(t-1)=i,X(t)\in
dx\right]}{z_{j}}$ $\displaystyle=\frac{1}{z_{j}}\sum_{i\neq
j}z_{i}G_{ij}\gamma_{ij}(dx).$ (21)

###  III.2 Algorithm

If one had a long trajectory that passed between the metastable states many
times, one could compute (8) by identifying the reactive trajectories from $A$
to $B$ and accumulating the increment per unit time in $\theta$ as a function
of $\theta$. However, in the stratification scheme above, we generally do not
know whether a copy of the system (random walker, or walker for short) is part
of a reactive trajectory as we are evolving it.

The key to addressing this issue is based on the fact that, as described
previously Warmflash, Bhimalapuram, and Dinner (2007); Dickson, Warmflash, and
Dinner (2009b, a); Dinner  _et al._ (2018), in practice, $\gamma_{ij}(dx)$ is
a list of configurations that are saved when walkers enter stratum $j$ from
stratum $i$. Consequently, for each walker, we can save the $\theta$ points
visited, the increments in $\theta$ at those points, and a pointer to the
configuration in $\gamma_{ij}(dx)$ from which the walker originated; then,
when a walker enters $A$ or $B$, we reconstruct the full sequence of walkers
connecting the metastable states and in turn the sequence of $\theta$ values
visited and the associated increments, and we add these data to the
appropriate sums. The algorithm is as follows.

  1. 1.

Initialize simulation quantities as follows.

    1. (a)

Define the $2K$ strata as in (14).

    2. (b)

Use preliminary trajectories (obtained from either unbiased simulation or an
enhanced sampling method such as STePS Guttenberg, Dinner, and Weare (2012))
to populate $\gamma_{ij}(dx)$. In practice, $\gamma_{ij}(dx)$ is a list of
configurations used to represent the distribution of entry points to stratum
$j$ from stratum $i$, and each element contains a full configuration and a
pointer (set to null for the configurations from the preliminary
trajectories). In the present study, we allow the list to grow throughout the
simulation (i.e., we do not overwrite configurations, as doing so may bias the
statistics as discussed in ref. Dickson, Warmflash, and Dinner, 2009a).

    3. (c)

Estimate $G_{ij}$ as the number of transitions from stratum $i$ to stratum $j$
in the preliminary trajectories, normalized by the number of transitions out
of stratum $i$.

    4. (d)

Initialize $z$ by solving (13).

    5. (e)

Create a grid over the collective variables and denote the volume of each grid
element by $\Delta\Theta$. Denoting the grid points by $\Theta$, initialize
all elements of the following arrays to zero:

      1. i.

$p^{+}_{A}(\Theta)$ (resp. $p^{+}_{B}(\Theta)$), which accumulates the weight
of trajectories entering state $A$ (resp. $B$);

      2. ii.

$p^{-}_{A}(\Theta)$ (resp. $p^{-}_{B}(\Theta)$), which accumulates the weight
of trajectories exiting state $A$ (resp. $B$);

      3. iii.

$v_{AB}(\Theta)$ (resp. $v_{BA}(\Theta)$), which accumulates the weighted CV
increments of trajectories exiting state $A$ and entering state $B$ (resp.
exiting state $B$ and entering state $A$).

  2. 2.

Run NEUS essentially as described in Section 3.3 of ref. Dinner  _et al._ ,
2018. Namely, at each iteration $l$:

    1. (a)

Initialize all elements of a matrix $g^{(l)}$ with the same dimensions as $G$
(i.e., $2K\times 2K$) to zero.

    2. (b)

For each stratum $j$, draw $N$ walkers from $\pi_{j}$. To be precise, for each
walker, select an element of $\gamma_{ij}(dx)$ with probability $z_{i}G_{ij}$
and use the saved configuration to initialize the walker; this procedure
represents (21). Associate with the walker a pointer to the element of
$\gamma_{ij}(dx)$.

    3. (c)

Evolve each walker by unbiased dynamics until it exits its initial stratum
($j$ in step 2b).

    4. (d)

When a walker exits from stratum $j$ to stratum $k$, save its final
configuration and associated pointer to $\gamma_{jk}(dx)$ and add 1 to
$g^{(l)}_{jk}$. Also save the sequence of $\theta$ values visited by the
walker trajectory.

    5. (e)

Compute $G_{ij}$ as

$G_{ij}=\frac{\sum_{l^{\prime}=L}^{l}g^{(l^{\prime})}_{ij}}{(l-L+1)N},$

where $L=\max[1,\min(l-L^{\prime},L^{\prime})]$, and $L^{\prime}$ is chosen to
minimize bias from the preliminary trajectories (i.e., to provide a burn-in
period for the simulation).

    6. (f)

Update $z$ by solving (13).

    7. (g)

Check convergence criteria and go to step 2a if not satisfied. We examine both
the change in the $z$ vector and the number of complete reactive trajectories.

  3. 3.

When a walker enters $A\cup B$, reconstruct the sequence of walkers leading to
that event from the last exit of $A\cup B$ and in turn the sequence of
$\theta$ values and strata visited. For each stratum $i$ and each point
$\theta(t)$ in the sequence, determine the nearest grid point $\Theta$ and

    1. (a)

increment $p_{A}^{+}(\Theta)$ (resp. $p_{B}^{+}(\Theta)$) by $z_{i}$ if the
sequence terminated in $A$ (resp. $B$);

    2. (b)

increment $p_{A}^{-}(\Theta)$ (resp. $p_{B}^{-}(\Theta)$) by $z_{i}$ if the
sequence originated in $A$ (resp. $B$);

    3. (c)

increment $v_{AB}(\Theta)$ (resp. $v_{BA}$) by
$z_{i}[\theta(t+\tau)-\theta(t-\tau)]$ if the sequence originated in $A$ and
terminated in $B$ (resp. originated in $B$ and terminated in $A$).

Above, $z_{i}$ accounts for the weight of each trajectory segment, which is
necessary to correct for the enhanced sampling.

  4. 4.

Construct the TPT quantities for reactive trajectories from $A$ to $B$ as
follows (exchange $A$ and $B$ below for reactive trajectories from $B$ to
$A$).

    1. (a)

Compute the forward committor to $B$ as

$q_{+}^{\theta}(\Theta)=\frac{p^{+}_{B}(\Theta)}{p^{+}_{A}(\Theta)+p^{+}_{B}(\Theta)}.$

    2. (b)

Compute the backward committor from $A$ as

$q_{-}^{\theta}(\Theta)=\frac{p^{-}_{A}(\Theta)}{p^{+}_{A}(\Theta)+p^{+}_{B}(\Theta)}.$

Note that $p_{A}^{-}+p_{B}^{-}=p_{A}^{+}+p_{B}^{+}$ is the total weight
associated with all trajectories that start and end in $A\cup B$.

    3. (c)

Compute the current from $A$ to $B$ as

${\color[rgb]{0,0,0}{I}}_{AB}^{\theta}(\Theta)=\frac{v_{AB}(\Theta)}{2\tau\Delta\Theta\sum_{\Theta^{\prime}}[p^{+}_{A}(\Theta^{\prime})+p^{+}_{B}(\Theta^{\prime})]}.$

Figure 1 shows how the algorithm samples the ensemble of reactive trajectories
in a piecewise fashion and how they are reconstructed.

Figure 1: Schematic of algorithm. (upper left) Define the metastable states
and the strata (black grid) and initialize the data structures (step 1). In
the example shown, trajectories are initialized from the upper left metastable
state boundary, so the last metastable state visited is known. (upper right)
In each stratum, draw configurations from $\gamma_{ij}(dx)$, save pointers to
the associated $\gamma_{ij}(dx)$ elements, and then run unbiased dynamics from
these configurations until the trajectories exit the stratum (steps 2b and
2c). A single example stratum is shown. (lower right) As the simulation
progresses, each point in $\gamma_{ij}(dx)$ may give rise to segments of
reactive or unreactive trajectories. Save the collective variable values
associated with these trajectory segments (step 2d). (lower left) When a
reactive or unreactive trajectory is realized (i.e., a trajectory segment
enters a metastable state), trace back from pointer to pointer to determine
the sequences of strata and collective variable grid sites (gray grid) visited
(step 3). Use these to increment unnormalized statistics with appropriate
weights (steps 3a-c). Note that the collective variable grid sites and the
strata need not coincide.

##  IV Numerical Results

We demonstrate our method by analyzing the C${}_{\text{7ax}}$ $\to$
C${}_{\text{7eq}}$ transition of the alanine dipeptide (N-acetyl-
alanyl-N’-methylamide) in vacuum. The peptide is represented by the CHARMM36m
force field MacKerell  _et al._ (1998); Best  _et al._ (2012); Huang  _et al._
(2017). The simulations are performed at 300 K with the SD (stochastic
dynamics) integratorGoga  _et al._ (2012) in GROMACS 5.1.4 Abraham  _et al._
(2015); LINCS was used to constrain the lengths of bonds to hydrogen atoms;
the step size was 1 fs, and the friction coefficient was $10\text{ ps}^{-1}$.
We used PLUMED 2.3 Bonomi  _et al._ (2009); Tribello  _et al._ (2014); Bonomi
_et al._ (2019) to extract collective variables; simulations were terminated
when walkers left their strata using the PLUMED COMMITTOR function.

The metastable states and strata are defined in terms of the $\phi$ and $\psi$
dihedral angles. The metastable states are taken to be circles of radius
$20^{\circ}$ around $(\phi,\psi)=(-83^{\circ},75^{\circ})$ and
$(70^{\circ},-70^{\circ})$, which correspond to the minima of the potential of
mean force determined below. The system was initialized from preliminary
trajectories that were generated by using Steered Transition Path Sampling
(STePS) Guttenberg, Dinner, and Weare (2012) to drive the system from each
metastable state to the other. STePS builds up trajectories by repeatedly
shooting bursts of short segments and then preferentially selecting those that
make forward progress for continuation. Trajectories were initialized from
points in the equilibrium ensemble within each metastable state, and progress
was measured by the distance to the center of the other metastable state. The
segment length was 1 ps, and the bias threshold was 0.75 (i.e., forward
trajectories are selected with a probability of 0.75 or their probability
estimated from the burst, whichever is higher). Initially, 10 trajectories
were released in each burst; subsequently, we varied the number with the
progress distance such as to expect two forward segments in each burst based
on statistics of the previous run. We do not use the statistics to correct for
the STePS bias—i.e., each reactive trajectory is weighted equally.

The strata form a partition of unity; their centers are at $(-7+2n)30^{\circ}$
for integer $n\in[1,6]$, and their widths are $100^{\circ}$, so they overlap
$40^{\circ}$ both ways, and the strata wrap around to enforce periodicity. We
release $10$ walkers in every stratum in each iteration and run each iteration
until all 10 walkers exit their stratum. The collective variables $\phi$ and
$\psi$ and their increments are saved every 100 fs. The lag time used to
compute the increments ($\tau$ in (8)) was 2 ps. The simulation is run until
we see at least $1000$ crossings in each direction between the metastable
states, i.e., walkers that cross from $j\leq K$ to $j>K$ and the other way. We
used a burn-in period of $L^{\prime}=100$ iterations. The full simulation was
approximately 140 million steps. Though we stratify in $\phi$ and $\psi$, we
can project the results onto other variables as well, and we also consider
$\omega$ (Figure 2). Statistics are accumulated for 50 values of each dihedral
angle.

Figure 2: Alanine dipeptide with the three dihedral angles labeled: $\phi$
(C5-N7-C9($\alpha$)-C15), $\psi$ (N7-C9($\alpha$)-C15-N17), and
(C1-C5-N7-C9($\alpha$)). Colors are cyan for carbon, blue for nitrogen, white
for hydrogen, and red for oxygen.

Figure 3: Potentials of mean force computed from the NEUS simulations:
projections on (top) $\phi$ and $\psi$ and (bottom) $\phi$ and $\omega$. The
stable states are marked with white circles on the top plot.

Potentials of mean force (PMF) computed by NEUS are shown in Figure 3. In
Figure 4 we show the forward committor and the $A\to B$ and $B\to A$ currents
in the $\phi\psi$-plane. We can see the clear existence of two main pathways,
one of which cuts diagonally across the middle of the figure and one of which
cuts across the upper left and lower right corners. Three representative
trajectories consistent with the pathways are shown in Figure 5. In Figure 6,
we show the commitor and currents in the $\phi\omega$-plane. We observe a
diagonal structure in these plots that is consistent with coupling between
distortion of the peptide plane and $\phi$ noted previously Bolhuis, Dellago,
and Chandler (2000).

To assess the accuracy of our committor estimates, we compute the committor as
a function of all three dihedral angles (not shown), select configurations
predicted to have $0.49<q_{+}^{\theta}<0.51$, and evaluate their committors by
shooting 20 independent simulations from each configuration. The resulting
histogram of values is peaked at $q_{+}\approx 0.5$, as desired. To
characterize the convergence of the committor estimates, we compute the
backwards committor and plot $q_{+}^{\theta}+q_{-}^{\theta}-1$, which should
be close to zero for this system (Figures 9 and 10); the deviation provides an
indication of the numerical error. We see that after 250 crossings between the
metastable states the error overall is small, with the largest deviations in
regions with high free energy close to where $q_{+}^{\theta}=0.5$. These
deviations decrease with additional crossings.

Figure 4: Committor (top) and forward (middle) and backward (bottom) reactive
currents for the C7ax $\to$ C7eq transition of the alanine dipeptide, plotted
on the $\phi\psi$-plane.  Figure 5:  Example trajectories plotted on the
$\phi\psi$-plane. Dots mark the centers of the metastable states.

Figure 6: Committor (top) and forward (middle) and backward (bottom) reactive
currents for the C7ax $\to$ C7eq transition of the alanine dipeptide, plotted
on the $\phi\omega$-plane. For the reactive currents, only grid points with
data from at least 5 trajectories are shown.  Figure 7: Histogram of
committors computed by shooting for 137 structures predicted to have
$0.49<q_{+}^{\theta}<0.51$ for reaction from $A$ to $B$. 20 independent
simulations were used for each configuration.

To assess the reactive current estimates, we compute the rate as a function of
the number of crossings two ways. First we directly sum the flux into
$B$Vanden-Eijnden and Venturoli (2009); Dickson, Warmflash, and Dinner
(2009a):

$R_{\rm NEUS}=\frac{1}{\sum_{j=1}^{2K}z_{j}}\sum_{j\leq
K}z_{j}\frac{n^{B}_{j}}{T_{j}},$ (22)

where $n^{B}_{j}$ is the number of reactive trajectories that enter $B$ from
stratum $j$, $T_{j}$ is the total time simulated in that stratum, and the sum
over $j\leq K$ selects for contributions from trajectories that were last in
$A$. Second, we sum the reactive currents that cross a dividing surface $S$:

$R_{\rm
TPT}=\sum_{S}{\color[rgb]{0,0,0}{I}}_{AB}^{\theta}(\Theta)\cdot\hat{n}_{S}\Delta_{S},$
(23)

where $\hat{n}_{S}$ is the unit vector normal to $S$, and $\Delta_{S}$ is a
differential element with the same dimension as $S$. Since the collective
variable space is periodic, there are two pathways, and we average over a set
of surfaces for each. Specifically, we use vertical lines in the $\phi\psi$
plots: (i) $\phi=(-57.6+7.2n)^{\circ}$ for integer $n$ in $[0,15]$, which cut
through the middle of Figure 4, and (ii) $\phi=(-172.8+7.2n)^{\circ}$ for
integer $n$ in $[0,9]$ and $[37,49]$, which line the left and right sides of
Figure 4.

The two estimates are about a factor of two different from each other (Figure
8); we also obtain a direct estimate of $1.4\times 10^{-6}$ ps-1 from an
unbiased trajectory of length 2.5 $\mu$s. We consider the three estimates to
be good agreement given that the rate can be a very challenging statistic to
converge to even its order of magnitude. We expect (22) to be most precise (as
evidenced by the scale of fluctuations in Figure 8) and recommend its use; we
consider (23) only as a means of validating the reactive currents. That said,
we can use the currents to obtain the fluxes associated with each of the
pathways. We find them to be $R_{\rm TPT}=4.6\times 10^{-7}$ ps-1 for the
pathway that crosses $\phi=0^{\circ}$ and $R_{\rm TPT}=8.1\times 10^{-7}$ ps-1
for the pathway that crosses $\phi=180^{\circ}$.

To further assess the reactive currents, we also compute
${\color[rgb]{0,0,0}{I}}_{AB}^{\theta}(\Theta)+J_{BA}^{\theta}(\Theta)$
(Figure 11). We expect this sum to be close to but not exactly equal to zero
given that the integrator is underdamped. We see that this is the case, with
the only significant deviations close to the metastable states. Put together,
these results validate the method.

Figure 8:  Rates computed directly from the flux into $B$ (top) and by
integrating the reactive currents (bottom). An equal number of crossings
between metastable states each way is used for each point on the horizontal
axis.

##  V Conclusions

In this paper, we investigated the computation of committors and reactive
currents with trajectory stratification. The fundamental challenge is that
these quantities require knowledge of the metastable states in which
trajectories begin and end, but this information cannot be obtained directly
from most walkers within the simulation. To address this issue, we store the
collective variable values visited by each walker and the associated
collective variable increments sampled, as well as a pointer to the previous
walker. This enables us to reconstruct the sequence of walkers that gave rise
to a trajectory once it reaches a metastable state, despite the fact that the
number of trajectories to which a single walker can contribute grows
exponentially with the number of walkers that follow it. The current work
builds on previous trajectory stratification studies that showed that
computation of the rate requires separating the ensemble of trajectories based
on the metastable state in which they originate Dickson, Warmflash, and Dinner
(2009a); Vanden-Eijnden and Venturoli (2009). As in previous trajectory
stratification studies, the results are exact in the limit of infinite
sampling, in contrast to MSMs and DGA.

Although we demonstrated the approach with NEUS, it is general and can be
applied to other algorithms that sample trajectory segments Huber and Kim
(1996); Zhang, Jasnow, and Zuckerman (2007); Zuckerman and Chong (2017);
Bolhuis (2003); van Erp, Moroni, and Bolhuis (2003); Allen, Warren, and Ten
Wolde (2005); Allen, Valeriani, and ten Wolde (2009); Bello-Rivas and Elber
(2015). By the same token, we focused on a simple numerical example that
permits validation of the results. Extending the approach to more complex
systems should be straightforward, so long as collective variables that enable
an efficient stratification and the accumulation of statistics can be
identified. This challenge is system specific, though methods for facilitating
the selection of collective variables based on simulation data exist (e.g.,
refs. Ma and Dinner, 2005; Strahan  _et al._ , 2021). In the case of NEUS, the
most significant challenge may be that the lists representing
$\gamma_{ij}(dx)$ can become large (up to 5500 elements in the present
example). Although this issue has not been limiting for complex systems
treated to date Dickson  _et al._ (2011), it is an important practical
consideration that we leave for future work.

##  VI Supplementary Material

Justification of (8), committor convergence plots, sums of forward and
backward reactive currents.

##  VII Acknowledgments

The authors wish to thank Chatipat Lorpaiboon and John Strahan for useful
discussions and critical readings of the manuscript. This work was supported
by National Institutes of Health award R35 GM136381 and National Science
Foundation award DMS-2054306.

## References

  * Ma and Dinner (2005) A. Ma and A. R. Dinner, “Automatic method for identifying reaction coordinates in complex systems,” Journal of Physical Chemistry B 109, 6769–6779 (2005).
  * Vanden-Eijnden (2006) E. Vanden-Eijnden, “Transition path theory,” in  _Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1_, edited by M. Ferrario, G. Ciccotti, and K. Binder (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 453–493.
  * E. and Vanden-Eijnden (2006) W. E. and E. Vanden-Eijnden, “Towards a theory of transition paths,” Journal of Statistical Physics 123, 503 (2006).
  * Metzner, Schütte, and Vanden-Eijnden (2006) P. Metzner, C. Schütte, and E. Vanden-Eijnden, “Illustration of transition path theory on a collection of simple examples,” Journal of Chemical Physics 125, 084110 (2006).
  * Noé  _et al._ (2009) F. Noé, C. Schütte, E. Vanden-Eijnden, L. Reich, and T. R. Weikl, “Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations,” Proceedings of the National Academy of Sciences 106, 19011–19016 (2009).
  * Bowman, Pande, and Noé (2013) G. R. Bowman, V. S. Pande, and F. Noé,  _An Introduction to Markov State Models and their Application to Long Timescale Molecular Simulation_ , Vol. 797 (Springer Science & Business Media, 2013).
  * Wu  _et al._ (2017) H. Wu, F. Nüske, F. Paul, S. Klus, P. Koltai, and F. Noé, “Variational Koopman models: Slow collective variables and molecular kinetics from short off-equilibrium simulations,” Journal of Chemical Physics 146, 154104 (2017).
  * Buchete and Hummer (2008) N.-V. Buchete and G. Hummer, “Coarse master equations for peptide folding dynamics,” Journal of Physical Chemistry B 112, 6057–6069 (2008).
  * Thiede  _et al._ (2019) E. H. Thiede, D. Giannakis, A. R. Dinner, and J. Weare, “Galerkin approximation of dynamical quantities using trajectory data,” Journal of Chemical Physics 150, 244111 (2019).
  * Strahan  _et al._ (2021) J. Strahan, A. Antoszewski, C. Lorpaiboon, B. P. Vani, J. Weare, and A. R. Dinner, “Long-time-scale predictions from short-trajectory data: A benchmark analysis of the trp-cage miniprotein,” Journal of chemical theory and computation 17, 2948–2963 (2021).
  * Russo  _et al._ (2021) J. D. Russo, J. Copperman, D. Aristoff, G. Simpson, and D. M. Zuckerman, “Unbiased estimation of equilibrium, rates, and committors from Markov state model analysis,” arXiv:2105.13402 (2021).
  * Elber  _et al._ (2017) R. Elber, J. M. Bello-Rivas, P. Ma, A. E. Cardenas, and A. Fathizadeh, “Calculating iso-committor surfaces as optimal reaction coordinates with milestoning,” Entropy 19, 219 (2017).
  * Faradjian and Elber (2004) A. K. Faradjian and R. Elber, “Computing time scales from reaction coordinates by milestoning,” Journal of Chemical Physics 120, 10880–10889 (2004).
  * Elber (2020) R. Elber, “Milestoning: An efficient approach for atomically detailed simulations of kinetics in biophysics,” Annual Review of Biophysics 49, 69–85 (2020).
  * Dellago  _et al._ (1998) C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, “Transition path sampling and the calculation of rate constants,” Journal of Chemical Physics 108, 1964–1977 (1998).
  * Dellago, Bolhuis, and Chandler (1998) C. Dellago, P. G. Bolhuis, and D. Chandler, “Efficient transition path sampling: Application to Lennard-Jones cluster rearrangements,” Journal of Chemical Physics 108, 9236–9245 (1998).
  * Dellago  _et al._ (2002) C. Dellago, P. Bolhuis, P. L. Geissler, _et al._ , “Transition path sampling,” Advances in Chemical Physics 123 (2002).
  * Huber and Kim (1996) G. A. Huber and S. Kim, “Weighted-ensemble Brownian dynamics simulations for protein association reactions.” Biophysical Journal 70, 97–110 (1996).
  * Zhang, Jasnow, and Zuckerman (2007) B. W. Zhang, D. Jasnow, and D. M. Zuckerman, “Efficient and verified simulation of a path ensemble for conformational change in a united-residue model of calmodulin,” Proceedings of the National Academy of Sciences 104, 18043–18048 (2007).
  * Zuckerman and Chong (2017) D. M. Zuckerman and L. T. Chong, “Weighted ensemble simulation: review of methodology, applications, and software,” Annual Review of Biophysics 46, 43–57 (2017).
  * Bolhuis (2003) P. G. Bolhuis, “Transition-path sampling of $\beta$-hairpin folding,” Proceedings of the National Academy of Sciences 100, 12129–12134 (2003).
  * van Erp, Moroni, and Bolhuis (2003) T. S. van Erp, D. Moroni, and P. G. Bolhuis, “A novel path sampling method for the calculation of rate constants,” Journal of Chemical Physics 118, 7762–7774 (2003).
  * Allen, Warren, and Ten Wolde (2005) R. J. Allen, P. B. Warren, and P. R. Ten Wolde, “Sampling rare switching events in biochemical networks,” Physical Review Letters 94, 018104 (2005).
  * Allen, Valeriani, and ten Wolde (2009) R. J. Allen, C. Valeriani, and P. R. ten Wolde, “Forward flux sampling for rare event simulations,” Journal of Physics: Condensed Matter 21, 463102 (2009).
  * Warmflash, Bhimalapuram, and Dinner (2007) A. Warmflash, P. Bhimalapuram, and A. R. Dinner, “Umbrella sampling for nonequilibrium processes,” Journal of Chemical Physics 127, 114109 (2007).
  * Dickson, Warmflash, and Dinner (2009a) A. Dickson, A. Warmflash, and A. R. Dinner, “Separating forward and backward pathways in nonequilibrium umbrella sampling,” Journal of Chemical Physics 131, 154104 (2009a).
  * Vanden-Eijnden and Venturoli (2009) E. Vanden-Eijnden and M. Venturoli, “Exact rate calculations by trajectory parallelization and tilting,” Journal of Chemical Physics 131, 044120 (2009).
  * Dickson and Dinner (2010) A. Dickson and A. R. Dinner, “Enhanced sampling of nonequilibrium steady states,” Annual Review of Physical Chemistry 61, 441–459 (2010).
  * Dinner  _et al._ (2020) A. R. Dinner, E. H. Thiede, B. Van Koten, and J. Weare, “Stratification as a general variance reduction method for Markov chain Monte Carlo,” SIAM/ASA Journal on Uncertainty Quantification 8, 1139–1188 (2020).
  * Bello-Rivas and Elber (2015) J. M. Bello-Rivas and R. Elber, “Exact milestoning,” Journal of Chemical Physics 142, 03B602_1 (2015).
  * Bhatt, Zhang, and Zuckerman (2010) D. Bhatt, B. W. Zhang, and D. M. Zuckerman, “Steady-state simulations using weighted ensemble path sampling,” Journal of Chemical Physics 133, 014110 (2010).
  * Dinner  _et al._ (2018) A. R. Dinner, J. C. Mattingly, J. O. Tempkin, B. Van Koten, and J. Weare, “Trajectory stratification of stochastic dynamics,” SIAM Review 60, 909–938 (2018).
  * Lorpaiboon, Weare, and Dinner (2022) C. Lorpaiboon, J. Weare, and A. R. Dinner, “Augmented transition path theory for sequences of events,” in preparation (2022).
  * Metzner, Schütte, and Vanden-Eijnden (2009) P. Metzner, C. Schütte, and E. Vanden-Eijnden, “Transition path theory for Markov jump processes,” Multiscale Modeling & Simulation 7, 1192–1219 (2009).
  * Torrie and Valleau (1977) G. M. Torrie and J. P. Valleau, “Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling,” Journal of Computational Physics 23, 187–199 (1977).
  * Pangali, Rao, and Berne (1979) C. Pangali, M. Rao, and B. Berne, “A Monte Carlo simulation of the hydrophobic interaction,” Journal of Chemical Physics 71, 2975–2981 (1979).
  * Frenkel and Smit (2001) D. Frenkel and B. Smit,  _Understanding Molecular Simulation: from Algorithms to Applications_ (Elsevier, 2001).
  * Thiede  _et al._ (2016) E. H. Thiede, B. Van Koten, J. Weare, and A. R. Dinner, “Eigenvector method for umbrella sampling enables error analysis,” Journal of Chemical Physics 145, 084115 (2016).
  * Dickson, Warmflash, and Dinner (2009b) A. Dickson, A. Warmflash, and A. R. Dinner, “Nonequilibrium umbrella sampling in spaces of many order parameters,” Journal of Chemical Physics 130, 02B605 (2009b).
  * Guttenberg, Dinner, and Weare (2012) N. Guttenberg, A. R. Dinner, and J. Weare, “Steered transition path sampling,” Journal of Chemical Physics 136, 234103 (2012).
  * MacKerell  _et al._ (1998) A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, and M. Karplus, “All-atom empirical potential for molecular modeling and dynamics studies of proteins,” Journal of Physical Chemistry B 102, 3586–3616 (1998).
  * Best  _et al._ (2012) R. B. Best, X. Zhu, J. Shim, P. E. M. Lopes, J. Mittal, M. Feig, and A. D. MacKerell, “Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone $\phi$, $\psi$ and side-chain $\chi_{1}$ and $\chi_{2}$ dihedral angles,” Journal of Chemical Theory and Computation 8, 3257–3273 (2012).
  * Huang  _et al._ (2017) J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B. L. de Groot, H. Grubmüller, and A. D. M. Jr, “CHARMM36m: An improved force field for folded and intrinsically disordered proteins,” Nature Methods 14, 71–73 (2017).
  * Goga  _et al._ (2012) N. Goga, A. Rzepiela, A. De Vries, S. Marrink, and H. Berendsen, “Efficient algorithms for Langevin and DPD dynamics,” Journal of Chemical Theory and Computation 8, 3637–3649 (2012).
  * Abraham  _et al._ (2015) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindah, “GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers,” SoftwareX 1-2, 19–25 (2015).
  * Bonomi  _et al._ (2009) M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R. A. Broglia, and M. Parrinello, “PLUMED: A portable plugin for free-energy calculations with molecular dynamics,” Computer Physics Communications 180, 1961–1972 (2009).
  * Tribello  _et al._ (2014) G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and G. Bussi, “PLUMED 2: New feathers for an old bird,” Computer Physics Communications 185, 604–613 (2014).
  * Bonomi  _et al._ (2019) M. Bonomi, G. Bussi, C. Camilloni, G. A. Tribello, P. Banáš, A. Barducci, M. Bernetti, P. G. Bolhuis, S. Bottaro, D. Branduardi, R. Capelli, and P. Carloni, “Promoting transparency and reproducibility in enhanced molecular simulations,” Nature Methods 16, 670–673 (2019).
  * Bolhuis, Dellago, and Chandler (2000) P. G. Bolhuis, C. Dellago, and D. Chandler, “Reaction coordinates of biomolecular isomerization,” Proceedings of the National Academy of Sciences 97, 5877–5882 (2000).
  * Dickson  _et al._ (2011) A. Dickson, M. Maienschein-Cline, A. Tovo-Dwyer, J. R. Hammond, and A. R. Dinner, “Flow-dependent unfolding and refolding of an RNA by nonequilibrium umbrella sampling,” Journal of Chemical Theory and Computation 7, 2710–2720 (2011).

##  Appendix A Formulas for the reactive current

In this section, we derive an expression for the reactive current of a
diffusion and use it to motivate (7). To this end, we introduce the generator,
which describes the evolution of the expectation of functions over an
infinitesimal time:

$Lf(x)=\lim_{\tau\to 0}\frac{{\bf E}[f(X(\tau))|X(0)=x]-f(x)}{\tau}.$ (24)

Backward-in-time statistics such as $q_{-}(x)$ can be accessed through the
adjoint of $L$ in a $\pi$-weighted inner product (see Strahan  _et al._ (2021)
for further discussion), and in this sense we can write

$L_{\pi}^{\dagger}f(x)=\lim_{\tau\to 0}\frac{{\bf
E}[f(X(-\tau))|X(0)=x]-f(x)}{\tau}.$ (25)

For the diffusive dynamics

$\dot{x}(t)=b(x(t))+\sqrt{2}\sigma(x(t))\eta(t)$ (26)

with $a=\sigma\sigma^{\top}$ and $\eta$ is a white noise (i.e., the
$\pi$-weighted inner product is
$\langle\eta_{i}(t)\eta_{j}(s)\rangle=\delta_{ij}\delta(t-s)$), these
operators take the forms

$\displaystyle{L}f(x)$ $\displaystyle=\sum_{j}b_{j}(x)\frac{\partial
f}{\partial x_{j}}(x)+\sum_{ij}a_{ij}(x)\frac{\partial^{2}f}{\partial
x_{i}\partial x_{j}}(x)$ (27) $\displaystyle L^{\dagger}_{\pi}f(x)$
$\displaystyle=-\sum_{i}b_{i}(x)\frac{\partial f(x)}{\partial
x_{i}}+\frac{2}{\pi(x)}\sum_{j}\frac{\partial}{\partial
x_{i}}(a_{ij}(x)\pi(x))\frac{\partial f(x)}{\partial
x_{i}}+\sum_{ij}a_{ij}(x)\frac{\partial^{2}f(x)}{\partial x_{i}\partial
x_{j}}.$ (28)

###  A.1 Reactive current for diffusive dynamics

Here, we use the forms in (27) and (28) to show

${\color[rgb]{0,0,0}{I}}_{AB}\cdot\nabla\theta(x)=\frac{\pi(x)}{2}\left(q_{-}(x)\mathcal{L}[q_{+}\theta](x)-q_{+}(x)\mathcal{L}_{\pi}^{\dagger}[q_{-}\theta](x)\right)$
(29)

for $x\in(A\cup B)^{c}$. For clarity, we suppress the $x$ arguments of
functions in the rest of this section.

We start from the $i$-th component of the reactive current Vanden-Eijnden
(2006):

${\color[rgb]{0,0,0}{I}}_{AB,i}=q_{+}q_{-}{\color[rgb]{0,0,0}{I}}_{i}+\pi\bigg{(}q_{-}\sum_{j}a_{ij}\frac{\partial
q_{+}}{\partial x_{j}}-q_{+}\sum_{j}a_{ij}\frac{\partial q_{-}}{\partial
x_{j}}\bigg{)}$ (30)

where ${\color[rgb]{0,0,0}{I}}_{i}$ is the $i$-th component of the equilibrium
probability current:

${\color[rgb]{0,0,0}{I}}_{i}=\pi b_{i}-\sum_{j}\frac{\partial}{\partial
x_{j}}\big{(}\pi a_{ij}\big{)}.$ (31)

Substituting (31) into (30), summing over components, rearranging, and dotting
with $\nabla\theta$,

$\displaystyle{\color[rgb]{0,0,0}{I}}_{AB}\cdot\nabla\theta=\frac{\pi}{2}\sum_{i}\Bigg{[}q_{-}\bigg{(}$
$\displaystyle q_{+}b_{i}-\frac{1}{\pi}\sum_{j}q_{+}\frac{\partial}{\partial
x_{j}}\big{(}\pi a_{ij}\big{)}+2\sum_{j}a_{ij}\frac{\partial q_{+}}{\partial
x_{j}}\bigg{)}$ $\displaystyle+q_{+}\bigg{(}$ $\displaystyle
q_{-}b_{i}-\frac{1}{\pi}\sum_{j}q_{-}\frac{\partial}{\partial x_{j}}\big{(}\pi
a_{ij}\big{)}-2\sum_{j}a_{ij}\frac{\partial q_{-}}{\partial
x_{j}}\bigg{)}\Bigg{]}\frac{\partial\theta}{\partial x_{i}}.$ (32)

Adding the identity $\theta Lq^{+}=0$ to the first set of parentheses and
rearranging,

$\sum_{i}b_{i}\frac{\partial
q_{+}\theta}{dx_{i}}+\sum_{i,j}a_{ij}\frac{\partial^{2}q_{+}\theta}{\partial
x_{i}\partial x_{j}}-\frac{1}{\pi}q_{+}\sum_{i,j}\frac{\partial}{\partial
x_{j}}\big{(}\pi a_{ij}\big{)}\frac{\partial\theta}{\partial
x_{i}}+q_{+}\sum_{i,j}a_{ij}\frac{\partial^{2}\theta}{\partial x_{i}\partial
x_{j}}\\\ =L[q_{+}\theta]-\frac{1}{\pi}q_{+}\sum_{j}\frac{\partial}{\partial
x_{j}}\big{(}\pi a_{ij}\big{)}\frac{\partial\theta}{\partial
x_{i}}+q_{+}\sum_{i,j}a_{ij}\frac{\partial^{2}\theta}{\partial x_{i}\partial
x_{j}}.$ (33)

Similarly, adding the identity $\theta L^{\dagger}_{\pi}q_{-}=0$ to the second
set of parentheses in (32) and rearranging,

$\sum_{i}\bigg{(}q_{-}b_{i}-\frac{1}{\pi}\sum_{j}q_{-}\frac{\partial}{\partial
x_{j}}\big{(}\pi a_{ij}\big{)}-2\sum_{j}a_{ij}\frac{\partial q_{-}}{\partial
x_{j}}\bigg{)}\frac{\partial\theta}{\partial x_{i}}\\\
+\sum_{i}\theta\bigg{(}b_{i}\frac{\partial q_{-}}{\partial
x_{i}}-\frac{2}{\pi}\frac{\partial}{\partial x_{i}}(a_{ij}\pi)\frac{\partial
q_{-}}{\partial x_{i}}-\sum_{j}a_{ij}\frac{\partial^{2}q_{-}}{\partial
x_{i}\partial x_{j}}\bigg{)}\\\
=-L_{\pi,i}^{\dagger}[q_{-}\theta]+\frac{1}{\pi}q_{-}\sum_{j}\frac{\partial}{\partial
x_{j}}\big{(}\pi a_{ij}\big{)}\frac{\partial\theta}{\partial
x_{i}}-q_{-}\sum_{i,j}a_{ij}\frac{\partial^{2}\theta}{\partial x_{i}\partial
x_{j}}.$ (34)

Substituting these expressions into (32) gives (29).

###  A.2 General form for the reactive current

Given (29), we can use (24) and (25) to write

$\displaystyle{\color[rgb]{0,0,0}{I}}_{AB}\cdot\nabla\theta(x)=\lim_{\tau\to
0}\frac{\pi(x)}{2\tau}\\{$ $\displaystyle q_{-}(x){\bf
E}[\theta(X(\tau))q_{+}(X(\tau))|X(0)=x]$ $\displaystyle-$ $\displaystyle
q_{+}(x){\bf E}[\theta(X(-\tau))q_{-}(X(-\tau))|X(0)=x]\\}.$ (35)

We can re-write the first term on the right hand side as

$\displaystyle{\bf E}[\theta(X(\tau))q_{+}(X(\tau))|X(0)=x]$
$\displaystyle={\bf E}[{\bf
E}[\theta(X(\tau))\mathbf{1}_{B}(X(t_{+}(0)))\,|\,X(\tau),\,X(0)]|X(0)=x]$
(36) $\displaystyle={\bf
E}[\theta(X(\tau))\mathbf{1}_{B}(X(t_{+}(0)))\,|\,X(0)=x]$ (37)
$\displaystyle=q_{+}(x){\bf E}[\theta(X(\tau))\,|\,X(t_{+}(0))\in B,\
X(0)=x].$ (38)

The second equality follows from the tower property of conditional
expectations. The third equality follows from the definition of the
conditional probability and the definition of the committor, $q_{+}(x)={\bf
P}[X(t_{+}(0))\in B|X(0)=x]$. Applying the same logic to the backward-in-time
process,

${\bf E}[\theta(X(-\tau))q_{-}(X(-\tau))|X(0)=x]=q_{-}(x){\bf
E}[\theta(X(-\tau))|X(t_{-}(0))\in A,X(0)=x].$ (39)

Substituting into (35) and combining terms,

${\color[rgb]{0,0,0}{I}}_{AB}\cdot\nabla\theta(x)=\pi(x)q_{-}(x)q_{+}(x)\\\
\times\lim_{\tau\to 0}\frac{{\bf
E}[(\theta(X(\tau))-\theta(X(-\tau))\,|\,X(t_{-}(0))\in A,\,X(t_{+}(0))\in
B,\,X(0)=x]}{2\tau}$ (40)

for $x\in(A\cup B)^{c}$. In combining terms into a single expectation, we have
exploited the fact that, for a fixed $X(0)$ of a Markov process, $X(-\tau)$ is
independent of $X(t_{+}(0))$, and $X(\tau)$ is independent of $X(t_{-}(0))$.

##  Appendix B Forward and backward reaction quantities

To quantify the convergence of our results, we compute the backward committor
and plot $q^{+}-q^{-}-1$, which should be zero for a reversible system. We
show these data as a function of the number of completed crossings between the
metastable states in Figs. 9 and 10. We also plot the sum of the reactive
currents from A to B and from B to A. Because the dynamics are underdamped,
this quantity should not be exactly zero, though we would expect it to be
close in general. We see that this is the case in Fig. 11, with the main
deviations close to the boundaries of the metastable states.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 9: Sum of forward committor and backward committors offset by 1 (i.e., $q^{+}+q^{-}-1$) projected onto $\phi$ and $\psi$ after harvesting (a) 100, (b) 250, (c) 500, (d) 1000 metastable state crossings. |
---|---
(a) | (b)
|
(c) | (d)
Figure 10: Sum of forward committor and backward committors offset by 1 (i.e.,
$q^{+}+q^{-}-1$) projected onto $\phi$ and $\omega$ after harvesting (a) 100,
(b) 250, (c) 500, (d) 1000 metastable state crossings.

Figure 11:  Sum of forward and backward reactive currents for indicated axes.


Computer Science and Artificial Intelligence Lab, Massachusetts Institute of
Technology, USA and
https://erikdemaine.org/edemaine@mit.eduhttps://orcid.org/0000-0003-3803-5703
Department of Mathematics, Massachusetts Institute of Technology, USA and
https://yaelkirk.github.io/yaelkirk@mit.eduhttps://orcid.org/0009-0007-6718-7390NSF
Graduate Research Fellowship under Grant No. 2141064. Computer Science and
Artificial Intelligence Lab, Massachusetts Institute of Technology, USA and
https://rebeccayelin.github.io/ryelin@mit.eduhttps://orcid.org/0000-0003-4747-4978MIT
Stata Family Presidential Fellowship. Erik D. Demaine, Yael Kirkpatrick, and
Rebecca Lin <ccs2012> <concept>
<concept_id>10003752.10003777.10003779</concept_id> <concept_desc>Theory of
computation Problems, reductions and completeness</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10002950.10003624.10003633.10010917</concept_id>
<concept_desc>Mathematics of computing Graph algorithms</concept_desc>
<concept_significance>300</concept_significance> </concept> <concept>
</ccs2012> [500]Theory of computation Problems, reductions and completeness
[300]Mathematics of computing Graph algorithms

###### Acknowledgements.

We thank Lily Chung and Timothy Gomez for their helpful comments. 1

# Graph Threading with Turn Costs

Erik D. Demaine  Yael Kirkpatrick  Rebecca Lin

###### Abstract

How should we thread a single string through a set of tubes so that pulling
the string taut self-assembles the tubes into a desired graph? While prior
work [ITCS 2024] solves this problem with the goal of minimizing the length of
string, we study here the objective of minimizing the total turn cost. The
frictional force required to pull the string through the tubes grows
exponentially with the total absolute turn angles (by the Capstan equation),
so this metric often dominates the friction in real-world applications such as
deployable structures. We show that minimum-turn threading is NP-hard, even
for graphs of maximum degree $4$, and even when restricted to some special
cases of threading. On the other hand, we show that these special cases can in
fact be solved efficiently for graphs of maximum degree $3$, thereby fully
characterizing their dependence on maximum degree. We further provide
polynomial-time exact and approximation algorithms for variants of turn-cost
threading: restricting to threading each edge exactly twice, and on
rectangular grid graphs.

###### keywords:

Deployable structures, reconfiguration, turn costs, Eulerian tour, Hamiltonian
cycle, traveling salesperson, computational complexity, approximation
algorithms

###### category:

##  1 Introduction

Deployable structures are dynamic systems capable of transitioning between
multiple geometric configurations, with applications to compact storage, ease
of transportation, and versatility. One intriguing approach to building
deployable structures is to thread a collection of tubes with string so that
different configurations result from adjustments to string tension. For
example, this idea is at the heart of the “push puppet” family of toys.
Interest in this approach is growing within various applied fields, including
fabrication [Lin2024push] and robotics [bern2022contact], with the goal of
combining the advantages in adaptability and safety of soft materials with the
precision and integrity of rigid structures [bern2023fabrication,
patterson2023safe].

Threading was formalized into theoretical computer science in recent work by
Demaine, Kirkpatrick, and Lin [demaine2024graph]. Specifically, a threading of
a graph is a closed walk that, at every vertex, induces a connected graph on
the incident edges. (See Section 2 for precise definitions.) Thus routing a
single string along this walk, and pulling the string taut, self-assembles the
tubes (edges) to mimic the connectivity of the given graph. Past work
[demaine2024graph] gave a polynomial-time algorithm for finding the _minimum-
length_ threading of a given graph.

In this paper, we analyze threading according to a _turn_ metric, motivated by
challenges that arise when constructing threadings in practice [Lin2024push].
While minimizing length makes sense to save on material cost and manufacturing
time, the feasibility of deployable structures is limited by the force
required to pull the string, which is usually dominated by the amount of
friction in the system. According to the Capstan equation, friction increases
_exponentially_ with the sum of the absolute values of turn angles in the
threading route [Kyosev2015braiding]. This motivates us to find threadings of
minimum total “turn cost”, where we allow each vertex to specify the
(symmetric) cost of turning between every pair of incident edges. (Again see
Section 2 for precise definitions.)

###  1.1 Our Results

In this paper, we analyze the Threading problem under a turn metric, where the
goal is to find the threading of a graph with minimum total turn cost. We
present several new results:

  * •

In Section 3, we prove NP-completeness of Threading under a turn metric in
various settings. These results provide a stark contrast to the length metric
where Threading can be solved in polynomial time [demaine2024graph].

    * –

We give a simple proof that Threading is NP-complete, even if the objective is
simply to minimize the number of turns, that is, the turn costs are either $0$
or $1$. This reduction shows how just one linear-degree vertex can effectively
simulate an arbitrary Hamiltonian Cycle instance, motivating the restriction
to bounded degree.

    * –

We prove that Threading is NP-complete even for graphs of maximum degree $4$,
via a reduction from 1-in-3 SAT.

    * –

Perfect Threading is a special case where we require using the smallest
possible turn cost when restricted to every vertex. Such a threading is not
always possible, and we prove that it is NP-complete to decide even for graphs
of maximum degree $4$. The analogous problem for the length metric (visiting
each vertex the fewest possible times) is of course polynomial, and was a
helpful stepping stone toward the full polynomial-time algorithm
[demaine2024graph].

    * –

It follows from our reductions that both Threading and Perfect Threading are
NP-complete even when we further constrain the threading to at most two visits
per edge, a problem we call Double Threading.

  * •

In Section 4, we provide polynomial-time algorithms for Perfect Threading and
Double Threading in graphs of maximum degree $\leq 3$. As we show these
problems NP-hard for graphs of maximum degree $4$, our results give a complete
characterization of the complexity of these two problems. (It remains open
whether the more general Threading problem is hard in graphs of maximum degree
$\leq 3$.)

  * •

In Section 5, we further analyze interesting special cases of Threading under
our turn metric, in particular developing several more polynomial-time
algorithms:

    * –

When each edge must be traversed _exactly_ twice, Exactly-Double Threading, we
show that the problem is equivalent to solving a sequence of Traveling
Salesman problems, one at each vertex. We use this fact to prove the
intractability of Exactly-Double Threading for unbounded-degree graphs, to
develop polynomial-time approximation algorithms, and to develop exact
algorithms for graphs of logarithmic degree or for minimizing the number of
turns (turn costs $\in\\{0,1\\}$).

    * –

We present a strategy for threading on rectangular grid graphs.

    * –

If the ratio between the maximum and minimum turn costs is $r$, we provide a
$2r$-approximation algorithm.

Minimizing turn cost is a generalization of minimizing length (e.g., define
the $e\to f$ turn cost to be $\frac{|e|+|f|}{2}$), but minimum-turn threadings
can have arbitrarily longer-than-optimal lengths (see Figure 2).

###  1.2 Related Work

Turn costs have been studied for several other problems. In the angular-metric
traveling salesman problem [aggarwal2000angular], the goal is to minimize the
total absolute turn angle of the tour. In minimum-turn milling
[arkin2005optimal], the goal is to cover a region with a moving tool while
minimizing the total turn cost or number of turns, motivated by machinery such
as mills which operate more efficiently moving in straight lines rather than
turning. In minimum-turn Eulerian tour [ellis2015dna], the goal is to minimize
the total turn cost; in surprising contrast to regular Eulerian tours, this
problem is NP-hard even for minimizing the number of turns (turn costs
$\in\\{0,1\\}$) in graphs of bounded degree.

##  2 Preliminaries

In this section, we discuss the necessary background for Threading and
introduce several special variants of the problem. We assume all graphs in
question to be connected with a minimum degree of at least $2$.

######  Definition 2.1 (Threading [demaine2024graph]).

Let $G=(V,E)$ be a graph with $n=|V|$ vertices and $m=|E|$ edges, where each
edge $e\in E(G)$ represents a tube and each vertex $v\in V(G)$ represents the
junction where the tubes incident to $v$ meet. A

Figure 1:  Two distinct threadings of a tetrahedron, where each junction graph
is either a tree (a) or a cycle (b) ([demaine2024graph, Figure 2]) Their
threading graphs are shown in the corners.

While prior work [demaine2024graph] focused on minimum-_length_ threading, we
shift our attention to minimum-_turn_ threading. Let us first define turn
costs formally:

######  Definition 2.2 (Turn Cost).

A

Now we can state the problem of finding a minimum-turn threading as:

######  Problem 2.3 (Threading).

Given $(G,\alpha)$, find a threading $T$ of $G$ that minimizes $\alpha(T)$.

The decision version of Threading is as follows:

######  Problem 2.4 (Decision Threading).

Given $(G,\alpha)$ and a nonnegative real number $c$, does there exist a
threading $T$ of $G$ such that $\alpha(T)\leq c$?

#### Naïve Double-Threading Solution.

We can naïvely compute a threading. First, assign each junction graph to be a
cycle; see the example in Figure 1b. Then, compute an Euler tour in $O(m)$
time through the union of these junction graphs (Figure 1) that prohibits
U-turns [demaine2024graph, Section 2.2]. This tour is a threading that
“double-threads” each edge.

######  Definition 2.5 (Minimum-Turn Perfect Threading).

A

By definition, each junction graph induced by a threading is connected, and so
contains a spanning tree. Hence, the cost incurred by the threading at each
junction must be at least that of a minimum spanning tree at that junction.
Therefore, if a minimum-turn perfect threading exists, it is a minimum-turn
threading. It follows that finding a minimum-turn perfect threading is a
subproblem of Threading.

######  Problem 2.6 (Perfect Threading).

Given $(G,\alpha)$, find a minimum-turn perfect threading of $G$, if it
exists.

At times, it may be pragmatic to limit the number of traversals through a
tube, which motivates our final problem definition:

######  Problem 2.7 (Double Threading).

Given $(G,\alpha)$, find a minimum-turn threading of $G$ that traverses each
edge at most twice.

We prove our first result: Decision Threading is in NP. It is straightforward
to verify that the same conclusion holds for the decision versions of Perfect
Threading and Double Threading.

######  Theorem 2.8.

Decision Threading is in NP.

######  Proof 2.9.

A purported solution $T$ to Decision Threading can be verified by checking
that $T$ is indeed a threading of $G$ and that $\alpha(T)\leq c$. Each of
these checks can be performed in $O(|T|)$ time.

We now prove that, for any threading $T$, there exists a threading
$T^{\prime}$ with size polynomial in $m$ such that
$\alpha(T^{\prime})\leq\alpha(T)$, thereby showing that a polynomial-size
witness for the solution exists. We can construct $T^{\prime}$ by repeatedly
deleting cycles in $T$ that do not disconnect any junction graphs until no
such cycles remain. We next show that $T^{\prime}$ traverses each edge at most
$O(m^{2})$ times.

Assume for a contradiction that $T^{\prime}$ traverses an edge $uv\in E(G)$
more than $|T(G)|+1$ times, where $T(G)$ denotes the set of all turns in $G$.
Then $T^{\prime}$ contains at least $|T(G)|+1$ cycles, each formed via a
revisit to $uv$. By the pigeonhole principle, there exists at least one cycle
that fails to traverse a turn not already covered by the other cycles. This
redundant cycle may be removed to produce a shorter walk for which junction
graphs are still connected. This contradicts our construction of $T^{\prime}$.
Hence, it follows that $uv$ is traversed $\leq|T(G)|+1$ times. Since
$|T(G)|=\sum_{v\in V(G)}d(v)^{2}\in O(m^{2})$ [de1998upper], our conclusion
follows.

While the proof of Theorem 2.8 demonstrates that an edge is traversed
$O(m^{2})$ times, Figure 2 shows an example where an edge is traversed
$\Omega(m)$ times.

Figure 2:  A tree-like graph with degree-3 junctions as marked by the circles.
A green line connecting an edge pair indicates a turn cost of $0$, whereas a
red line indicates a turn cost of $1$. All other turns have zero cost. The
minimum-turn threading of this graph deviates far from a perfect threading —
it traverses the “root” triangle as many times as there are “leaf” triangles.

##  3 Intractability of Threading

In Section 3.1, we show that Threading is NP-complete, even with turn costs
$\in\left\\{0,1\right\\}$, via a simple reduction from the Hamiltonian Cycle
problem. The reduction uses the turn costs at a single linear-degree vertex to
encode the connectivity of an arbitrary graph. We later utilize this technique
to give positive results in Section 5.1. In Section 3.2, we show that
Threading remains NP-complete for graphs of maximum degree $4$ by simulating
1-in-3 SAT with vertices of degree $\leq 4$.

###  3.1 General Graphs

######  Theorem 3.1.

Threading is NP-complete, even if only minimizing the number of turns.

We prove the hardness of Threading by reducing from Hamiltonian Cycle, which
is NP-hard [garey1979computers]. The Hamiltonian Cycle problem asks whether a
graph $G$ contains a Hamiltonian cycle, a cycle that visits every vertex
exactly once.

Given an instance $G$ of Hamiltonian Cycle with $n$ vertices and $m$ edges, we
construct an instance $(H,\alpha)$ of Threading as follows (Figure 3): First,
create $n$ copies of the cycle graph $C_{3}$, one per vertex in $G$, and
connect each cycle via a bridge to a central vertex $c$. We refer to each
cycle and its incident bridge as an “arm” and denote the arm corresponding to
vertex $v_{i}$ as $a_{i}$. Next, at vertex $c$, assign each pair of incident
edges from distinct arms $a_{i}$ and $a_{j}$ a turn cost of $0$ if
$v_{i}v_{j}\in E(G)$. Assign all remaining turns a cost of $1$. By design,
traversing each arm incurs a cost of $4$. The time to construct $(H,\alpha)$
is $O(m+n)$.

Figure 3:  Constructing an instance $H$ of Threading from an instance $G$ of
Hamiltonian Cycle. The turn costs at the central vertex in $H$ are depicted
inside the right-most circle, where a green line connecting a pair of edges
indicates a turn cost of $0$ and a red line indicates a turn cost of $1$. Turn
costs for unlabeled turns are assumed to be $1$. The edges of the junction
graph induced by the threading in orange are in bold.

######  Lemma 3.2.

Graph $G$ contains a Hamiltonian cycle if and only if its corresponding
instance $(H,\alpha)$ to Threading has a threading of turn cost $4n$.

######  Proof 3.3.

Suppose $G$ has a Hamiltonian cycle $C$. Label the vertices of $G$ such that
$C$ is given by the sequence of vertices $v_{1},v_{2},...,v_{n},v_{1}$.
Construct a walk $T$ through $H$ that traverses the arms the same order:
$a_{1},a_{2},\dots,a_{n},a_{1}$. Here, $a_{i}a_{j}a_{k}$ indicates the walk
entering arm $a_{j}$ directly from arm $a_{i}$, traversing all edges of
$a_{j}$ in either a clockwise or counter-clockwise fashion, and then exiting
$a_{j}$ directly to arm $a_{k}$. This walk creates a connected junction graph
at each vertex. Most notably, at the central vertex $c$, the junction graph is
given by the Hamiltonian cycle $C$ where each vertex $v_{i}$ maps to the
bridge from arm $a_{i}$. It follows that $T$ is a threading of $H$. This
threading incurs a turn cost of $4n$: Each arm of cost $4$ is traversed
exactly once because $C$ is a Hamiltonian cycle, and moving between
consecutive arms $a_{i}$ and $a_{j}$ in $H$ incurs zero cost because
$v_{i}v_{j}$ is in the Hamiltonian cycle and is hence in $E(G)$.

Now assume there exists a threading $T$ of $H$ of turn cost $4n$. By
definition, $T$ visits each arm at least once, and upon entering an arm, it
fully traverses the arm because U-turns are prohibited. Given that each arm
incurs a cost of $4$, turn cost $\alpha(T)=4n$ implies two facts: (1) $T$
traverses each arm exactly once and (2) incurs zero turn cost at vertex $c$.
Labeling the sequence of arms traversed by $T$ as $a_{1},a_{2},\dots,a_{n}$,
we map $T$ to vertices $v_{1},v_{2},\dots,v_{n},v_{1}$ in $G$. Consecutive
vertices in this sequence are adjacent by (2), hence this sequence forms a
cycle. Furthermore, each vertex is visited once per (1), so this cycle is
Hamiltonian.

It follows from Theorem 2.8 and Lemma 3.2 that Threading is NP-complete, even
with turn costs $\in\left\\{0,1\right\\}$.

###  3.2 Maximum-Degree-4 Graphs

Next we consider bounded degree graphs and show that even the special case of
Perfect Threading is NP-complete in this setting.

######  Theorem 3.4.

Perfect Threading in graphs of maximum degree $4$ is NP-complete.

Since Perfect Threading is a special case of Threading (Section 2), we
conclude:

######  Corollary 3.5.

Threading in graphs of maximum degree $4$ is NP-complete.

We prove Theorem 3.4 via a reduction from 1-in-3 SAT, which is NP-hard
[garey1979computers].

######  Definition 3.6 (1-in-3 SAT).

Given $m$ clauses consisting of three literals each, over a set of $n$ boolean
variables, the 1-in-3 SAT problem asks to determine whether there exists a
truth assignment to the variables such that each clause contains exactly one
true literal.

For example, the formula

$\varphi=\mathbbm{1}(w,x,\lnot y)\wedge\mathbbm{1}(\lnot w,x,\lnot
z)\wedge\mathbbm{1}(x,y,z),$ (1)

where $\mathbbm{1}(x,y,z)$ represents that exactly one of $x,y,z$ is true, has
a satisfying assignment of $x=\textsc{false},y=\textsc{true},z=\textsc{false}$
and $w=\textsc{true}$.

Given an arbitrary instance $\varphi$ of 1-in-3 SAT, we construct an instance
of Perfect Threading. This graph comprises several gadgets to simulate
$\varphi$: a “variable gadget” for each variable, a “clause gadget” for each
clause, and “split gadgets” to handle the occurrence of the same literal in
multiple clauses.

#### Variable Gadget.

For each variable $x$, construct a degree-$3$ junction as shown in Figure 4a.
A turn connected by a green line indicates a turn cost of 1 and a red line
indicates a turn cost of 2. The left and right edges incident to the junction
represent the literals $x$ and $\lnot x$, respectively. A key observation is
that any threading that forms a minimum spanning tree (MST) at this junction
will single-thread the top edge and select one of the bottom edges to single-
thread and the other to double-thread. We use this double-threading to
represent the assignment of $x$: If the bottom left edge is double-threaded,
then $x$ is assigned true; conversely, if the bottom right edge is double-
threaded, then $x$ is assigned false, that is, $\lnot x$ is assigned true
(Figure 4).

Figure 4:  (a) The junction representing variable $x$, with green indicating a
cost of $1$ and red indicating a cost of $2$. (b) Two possible minimum-cost
threadings of the junction and their corresponding junction graphs; double-
threading the left and right edges of the junction corresponds to setting $x$
and $\lnot x$ to true, respectively.

#### Clause Gadget.

To simulate a clause $\mathbbm{1}(x,y,z)$, we construct a junction that is
adjacent to variable gadgets $x$, $y$, and $z$ or their copies (as can arise
from the use of the splitting gadget), with all turn costs being $1$ (Figure
5). Any threading that induces an MST on this junction will double-thread
exactly one of the incident edges and single-thread the remaining two. By
using the double-threading to represent a true assignment, this corresponds to
having exactly one literal out of $x,y,z$ assigned true.

Figure 5:  (a) The gadget for clause $\mathbbm{1}(x,y,z)$ and (b) an example
threading that corresponds to assigning $x$ to true.

#### Split Gadget.

This final and most involved gadget allows us to “split” an edge representing
a literal $x$. It transforms a single edge simulating $x$ into three edges,
all simulating $x$ and thus receiving the same assignment. The gadget is
depicted in Figure 6, where the green, red, and blue lines indicate turns with
costs of $1$, $2$, and $3$, respectively.

Figure 6:  A split gadget. All unmarked turns have unit cost.

######  Claim 8.

Any Perfect Threading of the split gadget (Figure 6) threads edges $a,b,c,d$
the same number of times.

######  Proof 3.7.

First, note that the junction graph at vertex $x$ must be one of the two MSTs
shown in Figure 6b. Therefore, under any threading that induces an MST for
$J(x)$, either edges $a$ and $b$ are single-threaded while $e$ and $f$ are
double-threaded, or vice versa. Next, consider vertex $w$. To induce an MST on
$w$, edges $g$ and $h$ must be single-threaded. Hence, by simply counting the
number of traversals through vertex $y$, if edge $f$ is single-threaded, then
edge $c$ must be threaded twice, and vice versa. The same conclusion applies
symmetrically to edges $e$ and $d$. We conclude that a perfect threading of
the split gadget follows one of two options depicted in Figure 6b. In both
scenarios, the edges $a,b,c$, and $d$ are all threaded the same number of
times.

#### Full Construction.

Given a 1-in-3 SAT instance
$\varphi=\mathbbm{1}(\ell_{1,1},\ell_{1,2},\ell_{1,3})\land\mathbbm{1}(\ell_{2,1},\ell_{2,2},\ell_{2,3})\land\dots\land\mathbbm{1}(\ell_{m,1},\ell_{m,2},\ell_{m,3})$
using variables $x_{1},\ldots,x_{n}$, construct the graph $G(\varphi)$ as
follows:

Algorithm 1 Constructing a (multi)graph $G(\varphi)$ with turn costs from a
1-in-3 SAT formula $\varphi$.

  1. 1.

For every variable $x_{i}$, construct a variable gadget as shown in Figure 4.
We refer to the center-top edge as the “edge simulating variable $x_{i}$”, and
the bottom-left edge as an “edge simulating _literal_ $x_{i}$” and the bottom-
right edge as an “edge simulating _literal_ $\lnot x_{i}$”.

  2. 2.

If literal $x_{i}$ appears more than once in $\varphi$, construct a split
gadget and connect one of its outgoing edges to the edge simulating $x_{i}$.
Treat the remaining three outgoing edges of the gadget as “edges simulating
literal $x_{i}$”. Repeat this process until the number of edges simulating
$x_{i}$ is at least the number of appearances of $x_{i}$ in $\varphi$.

  3. 3.

For every clause $\mathbbm{1}(\ell_{i,1},\ell_{i,2},\ell_{i,3})$, construct a
clause gadget and connect it to edges simulating $\ell_{i,1},\ell_{i,2}$, and
$\ell_{i,3}$.

  4. 4.

Duplicate our entire construction up to this point. Connect each dangling edge
simulating a variable or a literal directly to its copy.

This construction runs in $O(n+m)$ time. An example construction for Equation
1 is provided in Figure 7a.

Figure 7:  (a) Construction of multigraph $G(\varphi)$ for the formula
$\varphi=\mathbbm{1}(x,y,z)\land\mathbbm{1}(x,\lnot
y,w)\land\mathbbm{1}(x,\lnot z,\lnot w)$. By directly connecting the dangling
edges simulating the (unused) literal $\lnot x$ (dashed) in Step 4 of
Algorithm 2, we create multiedges. (b) We can instead connect the edges via a
de-multigraph gadget (Figure 8) to obtain a simple graph
$G^{\prime}(\varphi)$.

######  Lemma 3.8.

$G(\varphi)$ has a minimum-turn perfect threading if and only if $\varphi$ has
a satisfying assignment.

######  Proof 3.9.

Let’s first construct a minimum-turn perfect threading of $G(\varphi)$ from a
satisfying assignment of $\varphi$ as follows: For each variable $x$, thread
its gadget according to Figure 4b, double-threading the edge simulating $x$ if
$x=\textsc{true}$ (left) and the edge simulating $\lnot x$ if
$x=\textsc{false}$ (right). For each literal $\ell$, thread its split
gadget(s) according to Figure 6b so that the edges $a,b,c,d$ are either all
double-threaded if $\ell=\textsc{true}$ (right) or all single-threaded if
$\ell=\textsc{false}$ (left). For each clause, thread its gadget as shown in
Figure 5b so that the edge simulating its single true literal is double-
threaded, while the other two edges are single-threaded. By design, this
threading induces an MST at each junction of $G(\varphi)$ so it is a minimum-
turn perfect threading.

Now let’s demonstrate the converse. Suppose we have a minimum-turn threading
of $G(\varphi)$. For each variable $x$, if the edge simulating $x$ in its
gadget is double-threaded, then assign $x$ to true, and false, otherwise. We
claim that this assignment satisfies $\varphi$. Consider a clause
$\mathbbm{1}(\ell_{i,1},\ell_{i,2},\ell_{i,3})$ in $\varphi$. Under a perfect
threading, the gadget representing this clause must have an edge double-
threaded and two single-threaded. Since these edges represent the literals
$\ell_{i,1},\ell_{i,2},\ell_{i,3}$, this means that exactly one of these
literals is assigned true, thereby satisfying
$\mathbbm{1}(\ell_{i,1},\ell_{i,2},\ell_{i,3})$.

Figure 8:  A de-multigraph gadget. All turns are of unit cost.

#### De-Multigraph.

As shown in Figure 7a, if a literal does not appear in $\varphi$, it creates a
multiedge with the edge simulating its corresponding variable. To address this
issue and transform the multigraph $G(\varphi)$ into a simple graph
$G^{\prime}(\varphi)$, we modify Step 4 of the construction by connecting the
dangling edge simulating the literal to its copy using the “de-multigraph
gadget” shown in Figure 8, rather than directly (Figure 7b). This gadget can
be threaded in a way that induces an MST at each of its junctions, regardless
of whether the incoming literal edges are threaded once or twice; see Figure
8b. Hence, we can adapt the proof of Lemma 3.8 and conclude that
$G^{\prime}(\varphi)$, with the specified turn costs, qualifies as an instance
of Perfect Threading and has a minimum-turn perfect threading if and only if
the 1-in-3 SAT instance $\varphi$ has a satisfying assignment.

This result concludes our proof of Theorem 3.4. As each edge in graph
$G^{\prime}(\varphi)$ is traversed no more than twice, our reduction further
implies the hardness of Double Threading. Hence:

######  Corollary 3.10.

Double Threading is NP-complete for graphs of maximum degree $4$.

##  4 Maximum-Degree-3 Graphs

Here we present polynomial-time algorithms for solving Perfect Threading and
Double Threading on graphs of maximum degree $3$, thereby fully characterizing
the complexity of these two problems.

###  4.1 Perfect Threading

######  Theorem 4.1.

Perfect Threading in graphs of maximum degree $3$ is in P.

Let us consider an instance $(G,\alpha)$ of Perfect Threading with graph $G$
of maximum degree $3$. We begin our proof of the theorem by observing that a
spanning tree at a degree-$3$ junction $v\in V(G)$ is a path of length $2$.
Hence, one of the edges incident to $v$ must be double-threaded while the
remaining two edges are single-threaded. From this observation, we deduce that
finding a perfect threading in $G$ is a matching problem: Selecting one
incident edge per degree-$3$ vertex to double-thread is equivalent to finding
a set of edges such that each degree-$3$ vertex is incident to one edge from
this set.

However, not every spanning tree qualifies as a _minimum_ spanning tree (MST).
Hence, to achieve a minimum-turn perfect threading, we must double-thread only
those edges that could induce MSTs. That is, we seek matchings that select for
edges $uv$ corresponding to degree-$2$ vertices in MSTs for both junction
graphs $J(u)$ and $J(v)$. See our proposed algorithm below:

Algorithm 2 Perfect Threading on maximum-degree-3 graphs.

  1. 1.

Delete all degree-$2$ vertices in $G$ and their incident edges.

  2. 2.

For every remaining vertex $v\in V$, consider all MSTs for $J(v)$ and delete
any edge incident to $v$ of degree $1$ in all MSTs.

  3. 3.

Compute a perfect matching of the remaining graph $G^{\prime}$. Return the
matching if it exists.

#### Correctness.

The correctness the algorithm follows from the claim below:

######  Claim 9.

Graph $G$ has a perfect threading if and only if the graph $G^{\prime}$
produced by Algorithm 2 has a perfect matching.

######  Proof 4.2.

Suppose $G$ has a minimum-turn perfect threading $T$. Consider the set $M$ of
edges that are double-threaded by $T$. Since $T$ induces an MST at every
junction, none of the edges in $M$ could be removed in Step 2 of Algorithm 2.
Hence, every edge of $M$ exists in $G^{\prime}$. Furthermore, since every
vertex in $G^{\prime}$ corresponds to a degree-$3$ vertex in $G$, it must have
exactly one incident edge that is double-threaded by $T$, implying one
incident edge in $M$. Therefore, $M$ is a perfect matching in $G^{\prime}$.

Conversely, given a perfect matching $M$ of $G^{\prime}$, we can construct
junction graphs that result in a minimum-turn perfect threading: For every
degree-2 vertex, we take the trivial junction graph comprising a single edge.
For every degree-$3$ vertex $v$ with its matching $u$ in $M$, we set $J(v)$ to
be the path of length $2$, with the central vertex corresponding to the edge
$uv$. This path forms an MST by the construction of $G^{\prime}$ (Step 2).

#### Running-Time.

Step 1 runs in $O(n)$ time. Step 2 also runs in $O(n)$ time given that each
degree-$3$ vertex has at most three MSTs. Finally, since $|G^{\prime}|\in
O(n)$, Step 3 runs in $O(n^{2}\log{n})$ time via the perfect matching
algorithm by Galil, Micali, and Gabow [galil1986ev].

###  4.2 Double Threading

While perfect threadings may not exist for some graphs, double-threadings are
always possible — refer to our naïve algorithm in Section 2 for details. Here
we present an efficient solution to Double Threading in maximum-degree-$3$
graphs.

######  Theorem 4.3.

Double Threading in graphs of maximum degree 3 can be solved in time
$O(m^{2}\log m)$.

Suppose $G$ is a graph of maximum degree $3$. Let $\mathcal{C}$ be a
collection of vertex-disjoint simple cycles in $G$. Prior work
[demaine2024graph, Lemma 9] achieves a threading $T$ that single-threads every
edge in $\mathcal{C}$ and double-threads every other edge. Because $G$ has
maximum degree $3$, the turn cost of this threading is:

$\alpha(T)=\underbrace{\sum_{(uv,vw)\in T(G)}\alpha(uv,vw)}_{\text{total cost
of all turns in $G$}}-\sum_{C\in\mathcal{C}}\alpha(C).$

Therefore, solving Double Threading on $G$ is equivalent to finding a
collection of vertex-disjoint simple cycles in $G$ whose total turn cost is
maximized. We give an algorithm to solve this problem; see Algorithm 3 below.

Algorithm 3 Computing maximum-turn vertex-disjoint simple cycles.

  1. 1.

Construct a weighted graph $G^{\prime}$ from $(G,\alpha)$ (Figure 9) as
follows:

    1. (a)

For each edge $uv$, create a zero-weight path of length $3$ comprising
vertices $uv^{+}$, $uv^{-}$, $vu^{-}$, and $vu^{+}$.

    2. (b)

For each turn $(uv,vw)$ in $T(G)$, add an edge of weight $\alpha(uv,vw)$
connecting the vertices $vu^{+}$ and $vw^{+}$.

  2. 2.

Compute a maximum weight perfect matching $M$ in $G^{\prime}$.

  3. 3.

Return the set of edges $S=\left\\{uv\in E(G):(uv^{-},vu^{-})\in M\right\\}$.

Figure 9:  Construction of $G^{\prime}$ from $G$ in Algorithm 3.

#### Correctness.

We briefly explain why the matching $M$ in $G^{\prime}$ corresponds to a
collection of vertex-disjoint simple cycles in $G$. Consider two possible
cases for each vertex $v\in V(G)$:

  1. 1.

Suppose $M$ contains a weighted edge $(vu^{+},vw^{+})$ created in Step 1b.
Only one such edge is possible for $v$ given that the set of vertices
$\\{v{x}^{+}\\}_{x\in N(v)}$ are pairwise connected and have size at most $3$.
It follows that the remaining vertices $\\{v{x}^{+}\\}_{x\in
N(v)\setminus\\{u,w\\}}$ are matched with their counterparts in
$\\{vx^{-}\\}_{x\in N(v)\setminus\\{u,w\\}}$, and so vertices $vu^{-}$ and
$vw^{-}$ must be saturated by edges $(vu^{-},uv^{-})$ and $(vw^{-},wv^{-})$,
respectively. These edges correspond to some cycle traversing edges $uv$ and
$vw$ through $u$ in $G$.

  2. 2.

Otherwise, all vertices in $\\{vx^{+}\\}_{x\in N(v)}$ are matched to vertices
in $\\{v{x}^{-}\\}_{x\in N(v)}$. This indicates that $v$ does not participate
in any cycle.

#### Running-Time.

Graph $G^{\prime}$ has $O(m)$ vertices and edges. We can adapt the minimum-
weight perfect matching algorithm by Galil, Micali, and Gabow [galil1986ev]
for Step 2. This procedure runs in
$O(|V(G^{\prime})||E(G^{\prime})|\log|V(G^{\prime})|)=O(m^{2}\log m)$ time,
which dominates the $O(m)$ construction time of $G^{\prime}$ in Step 1. Hence,
the overall running time of Algorithm 3 is $O(m^{2}\log m)$.

Finding a double-threading of $G$ from the collection of maximum-turn vertex
disjoint simple cycles given by Algorithm 3 can be done in $O(m)$ time
[demaine2024graph, Section 2.2].

It remains an open question whether the general Threading problem can be
solved efficiently in graphs with maximum degree $3$. Recall that Figure 2
depicts a graph whose minimum-turn threadings are far from perfect.

##  5 Other Special Cases

In this section, we give polynomial-time algorithms for several additional
variants of Threading. We show that Exactly-Double Threading is equivalent to
the Traveling Salesperson problem. This result paves the way for several
insights: a fixed-parameter tractable algorithm contingent on the maximum
degree of the graph, approximation algorithms, and a polynomial-time algorithm
if the goal is simply to minimize the number of turns. Next, we give a
polynomial-time algorithm to thread rectangular grids. Finally, we conclude
our paper with a polynomial-time approximation for Threading given bounded
positive turn costs.

###  5.1 Exactly-Double Threading

Consider $n$ cities, with distance $d(i,j)$ between cities $i$ and $j$. The
Traveling Salesperson problem (TSP) asks if there exists a tour of length at
most some constant $C\geq 0$ that visits each city exactly once.

######  Theorem 5.1.

Exactly-Double Threading is equivalent to Traveling Salesperson.

First, we sketch a reduction from TSP to Exactly-Double Threading. Given an
instance of TSP, we construct an instance $(G,\alpha)$ of Exactly-Double
Threading: Create $n$ distinct arms incident to a central vertex $c$ as in the
proof of Lemma 3.2. Set the turn costs between any two arms $a_{i}$ and
$a_{j}$ about vertex $c$ to $d(i,j)$. Set all other turn costs to $1$. It
follows that there exists a tour of length at most $C$ that visits each city
exactly once if and only if there exists a threading of $G$ of cost at most
$C+4n$.

Let’s now consider the reverse direction. We reduce Exactly-Double Threading
to TSP: Leveraging the key observation that each junction graph must be
precisely a cycle, we find the minimum-cost cycle at each junction $v$ by
solving TSP on $v$ — we assign a city to each edge incident to $v$, with
distance $\alpha(uv,vw)$ between cities $uv$ and $vw$. Subsequently, we apply
the polynomial-time procedure by Demaine, Kirkpatrick, and Lin
[demaine2024graph] to construct a threading, using these minimum-cost cycles
as junction graphs.

A fixed-parameter tractable algorithm for Exactly-Double Threading follows
directly from the proof of Theorem 5.1:

######  Corollary 5.2.

Exactly-Double Threading can be solved in polynomial time if the maximum
vertex degree $\Delta$ of the graph is in $O(\log n)$.

Solving TSP on each junction can be achieved in $O(d(v)^{2}+2^{d(v)})$ time
via the Bellman-Held algorithm [bellman1962dynamic, held1962dynamic], so the
runtime is polynomial as long as $d(v)\leq\Delta$ is logarithmic in $n$.

We can further inherit results from TSP: If the turn costs at each junction
form a metric space, that is, if they obey the triangle inequality, then we
can approximate Exactly-Double Threading in polynomial time by solving metric
TSP approximately at each junction. For instance, if all turn costs are either
$1$ or $2$, then we can utilize the approximation algorithm for (1, 2)-TSP by
Adamaszek, Mnich, and Paluch [adamaszek2018new] to obtain an
$8/7$-approximation scheme that runs in $O(\sum_{v\in V(G)}d(v)^{3}+m)$ time.
More generally, we can apply the Christofides-Serdyukov algorithm
[christofides2022worst] to obtain a $3/2$-approximation scheme that achieves
the same time complexity.

######  Corollary 5.3.

Exactly-Double Threading is in P if turn costs are in $\\{0,1\\}$.

This result follows from the fact that we can solve TSP in a _complete_ graph
$G$ with costs in $\\{0,1\\}$ in polynomial time.

The hardness of Double Threading given turn costs in $\left\\{0,1\right\\}$
remains open.

###  5.2 Rectangular Grid Graphs

A

######  Definition 10.

rectangular grid graph $G$ is a $w\times h$ lattice graph consisting of $w$
rows and $h$ columns of vertices, referred to as