††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# BAD-NEUS: Rapidly converging trajectory stratification

John Strahan  Department of Chemistry and James Franck Institute, University
of Chicago, Chicago IL 60637, USA  Chatipat Lorpaiboon  Department of
Chemistry and James Franck Institute, University of Chicago, Chicago IL 60637,
USA  Jonathan Weare  Courant Institute of Mathematical Sciences, New York
University, New York, New York 10012, USA  Aaron R. Dinner
<EMAIL_ADDRESS>Department of Chemistry and James Franck Institute,
University of Chicago, Chicago IL 60637, USA

###### Abstract

An issue for molecular dynamics simulations is that events of interest often
involve timescales that are much longer than the simulation time step, which
is set by the fastest timescales of the model. Because of this timescale
separation, direct simulation of many events is prohibitively computationally
costly. This issue can be overcome by aggregating information from many
relatively short simulations that sample segments of trajectories involving
events of interest. This is the strategy of Markov state models (MSMs) and
related approaches, but such methods suffer from approximation error because
the variables defining the states generally do not capture the dynamics fully.
By contrast, once converged, the weighted ensemble (WE) method aggregates
information from trajectory segments so as to yield unbiased estimates of both
thermodynamic and kinetic statistics. Unfortunately, errors decay no faster
than unbiased simulation in WE. Here we introduce a theoretical framework for
describing WE that shows that introduction of an element of stratification, as
in nonequilibrium umbrella sampling (NEUS), accelerates convergence. Then,
building on ideas from MSMs and related methods, we propose an improved
stratification that allows approximation error to be reduced systematically.
We show that the improved stratification can decrease simulation times
required to achieve a desired precision by orders of magnitude.

##  I Introduction

The sampling of a stochastic process can be controlled by evolving an ensemble
and splitting and merging the trajectories of its members, and various
algorithms based on this strategy have been introduced [1, 2, 3, 4, 5, 6, 7].
Because the trajectory segments between splitting events are unbiased, such
algorithms can yield dynamical statistics, such as transition probabilities
and mean first passage times (MFPTs) to selected states, in addition to
steady-state probability distributions and averages. Furthermore, splitting
algorithms generally require little communication between the members of the
ensemble, making them relatively straightforward to implement regardless of
the underlying dynamics, and community software is available [8]. As a result,
splitting algorithms are widely used.

Recent mathematical analysis of one of the oldest splitting algorithms in the
molecular simulation literature, weighted ensemble (WE) [1], shows that it is
asymptotically unbiased and can dramatically reduce the variance of statistics
when the splitting criteria, which are based on a partition of the state
space, are chosen appropriately [9]. However, the method relies on convergence
of the steady-state ensemble of trajectories, which is known to be slow. In
fact, as we argue below and was observed previously [9], it is as slow as
running direct unbiased simulations. Prior to convergence, the method is
systematically biased.

Some previous work to accelerate the convergence of WE focused on improving
the initialization of walkers using an approximation of the stationary
distribution. In one case, the approximation of the stationary distribution
was obtained from a free-energy method that biased trajectories to escape
metastable wells [10], while in another it was obtained from a machine-
learning method based on unbiased short trajectories [11]. The former approach
requires the underlying dynamics to obey detailed balance to ensure the
validity of the free energy method, while the latter approach relies on a
Markov approximation (but does not make other restrictions on the dynamics).
In both cases, however, the approximation to the stationary distribution is
used only for the first iteration, and so if the approximation is deficient,
significant bias can persist in the WE calculation for many iterations.

Others focused on improving the rate estimates from WE [12, 13]. Of particular
relevance to the present study, Copperman and Zuckerman exploited the fact
that splitting algorithms sample unbiased trajectory segments to construct
history-augmented Markov state models (haMSMs) [14] from their WE data [12].
In an haMSM, the system is divided into discrete states and the probabilities
for transitions between them are computed from data, just as in a standard
MSM, except that the path ensemble is split based on the last metastable state
visited; the rate constant is then computed as the flux from $A$ to $B$
conditioned on having last visited $A$, normalized by the steady-state
probability of having last visited $A$. The idea is that finer discretization
of the state space than the WE simulation yields improved estimates. However,
this approach does not accelerate the convergence of the sampling.

The basic approach in our work is to combine these strategies and fully
integrate them in the sampling algorithm by applying a short-trajectory based
approximation to the stationary distribution between successive iterations. As
we show below, if the approximation satisfies certain properties and is
correctly integrated in the algorithm, we can accelerate the convergence by
several orders of magnitude without introducing systematic bias. Our approach
is able to compute the full stationary distribution with minimal assumptions
about the dynamics, as well as kinetic statistics such as rates and committors
(probabilities of visiting selected states before others).

Conceptually, one can view our approach as a form of stratification, in which
the weights of not just individual trajectory segments but groups of them are
manipulated. Stratification was first introduced for controlling the sampling
of trajectory segments in nonequilibrium umbrella sampling (NEUS) [3, 15, 16,
4, 17, 18]. NEUS groups trajectories that are in the same regions of state
space, as defined by collective variables (CVs), and estimates the steady-
state probabilities of the regions by solving a global flux balance equation;
this strategy was subsequently adopted in extensions of WE[19] and exact
milestoning (EM)[6]. A unified framework for trajectory stratification that
incorporates elements of all of the above algorithms and is unbiased in the
limit that each region contains an infinite number of ensemble members is
presented in Ref. 18, which also shows that the regions can be defined in
terms of path-based quantities. Mathematical analysis of trajectory
stratification algorithms can be found in Refs. 20, 21.

Our paper is organized as follows. In Section II, we review the WE algorithm
and recast it to show why an initialization bias is slow to disappear. This
formulation also clarifies the relation of WE to the NEUS algorithm as
described in Ref. 18, and we show how software for WE can be easily extended
to enable trajectory stratification in Section III. In Section IV, we
introduce a generalization of MSMs that represents the dynamics through a
basis expansion [22, 23, 24] and in turn our strategy for accelerating
sampling, which we term Basis-Accelerated NEUS (BAD-NEUS). We show that the
NEUS algorithm is a special case of our new scheme. In Sections V and VI, we
demonstrate BAD-NEUS on a two-dimensional model with a known steady-state
distribution and a molecular example.

##  II Weighted Ensemble (WE)

As described above, in WE, the state space is partitioned into regions, and
then trajectories are copied (cloned or split) or removed (killed or pruned)
from the ensemble based on criteria. Mathematically, we denote the state of
ensemble member (henceforth, walker) $i$ at time $t$ by $X^{i}_{t}$, and we
associate with $i$ a weight $w^{i}_{t}$ such that $\sum_{i}w^{i}_{0}=1$. We
track the index of the region containing $X^{i}_{t}$ by an index process
$J^{i}_{t}$. For example, $J^{i}_{t}$ might track the element of a partition
of state space in which $X^{i}_{t}$ currently resides. Other options for the
index process are available and can be advantageous (see Refs. 18, 25 for
further discussion).

Here, we use a relatively simple procedure and represent the splitting and
pruning by resampling the ensemble within each region. Each iteration of the
sampling thus consists of two steps: evolution according to the underlying
dynamics and resampling the ensemble. In the evolution step, for each walker
$i$, we numerically integrate for a time interval $\Delta$ to obtain
$X^{i}_{t+\Delta}$ from $X^{i}_{t}$ and update $J^{i}_{t+\Delta}$ accordingly.
In this step the weights are unchanged; $w^{i}_{t+\Delta}=w^{i}_{t}$. We note
that $\Delta$ can be a random stopping time, not just a fixed time interval.
An often useful choice is to take $t+\Delta$ to be the first time after time
$t$ that the value of the index process $J^{i}$ changes. In this case, each
walker has a different value for $\Delta.$ In the resampling step, for each
value $k$ of the index process (e.g., for each region in a partition of state
space), if there is at least one walker with $J^{i}_{t+\Delta}=k$, we select a
number $N_{k}$. Then we sample from the set $\Phi=\\{i:J^{i}_{t+\Delta}=k\\}$
$N_{k}$ times with replacement according to the probabilities
$p_{i}=w^{i}_{t+\Delta}/\sum_{i\in\Phi}w^{i}_{t+\Delta}$. For each index $r$
so chosen, we append
$(X^{r}_{t+\Delta},J^{r}_{t+\Delta},\sum_{i\in\Phi}w^{i}_{t+\Delta}/N_{k})$ to
the new ensemble.

If the underlying dynamics are ergodic, in the case where $\Delta$ is a fixed
time horizon, WE can be used to compute steady-state averages of functions as

$\mathbbm{E}_{X_{0}\sim\pi}[g(X_{0},X_{1},\ldots,X_{\Delta-1})]\\\
=\lim_{T\to\infty}\frac{1}{T}\sum_{n=0}^{T}\frac{1}{N}\sum_{i=0}^{N}g(X^{i}_{n\Delta},X^{i}_{n\Delta+1},\ldots,\\\
X^{i}_{n\Delta+\Delta-1})w^{i}_{n\Delta},$ (1)

where the subscript on the expectation indicates that we draw the initial
state $X_{0}$ from the steady-state distribution $\pi,$ and $n$ indexes
successive weighted ensemble iterations. We present necessary modifications
for the case where $\Delta$ is a random variable later. Here and below, $X$
without a superscript indicates a realization of the underlying Markov process
rather than a walker in an ensemble. An important example is the steady-state
flux into a set $D$ of interest, which can be used to compute the MFPT to $D$
from another set in the right algorithmic setup [9]. This flux can be obtained
by setting
$g(X_{0},\ldots,X_{\Delta})=\mathbbm{1}_{D^{c}}(X_{0})\mathbbm{1}_{D}(X_{1})$,
where $\mathbbm{1}_{D}(x)$ is an indicator function that is 1 if $x\in D$ and
0 otherwise. This amounts to counting the number of walkers that enter $D$ in
a single time step and summing the total weight of those walkers.

Having stated the basic WE algorithm, we now present WE in a new way. While
this may initially appear to complicate the description, it reveals that WE is
slow to converge and facilitates the introduction of approaches to accelerate
convergence. In this description, which we call the distribution
representation of WE, we directly evolve distributions rather than
approximating them through individual walkers. To this end, we define the flux
distributions

$\overline{\pi}_{\ell}(dx|k)=P[X_{S(\ell)}\in dx|J_{S(\ell)}=k]$ (2)

where $dx$ is an infinitesimal volume in state space that is located at a
specific value of $x$ and $S(\ell)$ is an increasing sequence of stopping
times. For example, if we choose $S(\ell)=\ell\Delta$ for a fixed $\Delta$,
the following description corresponds to the standard version of WE already
sketched. Alternatively, we can let $S(\ell)$ be the time of the $\ell$-th
change in the value of the index process $J$:

$\displaystyle S(0)$ $\displaystyle=0$ (3) $\displaystyle S(1)$
$\displaystyle=\min\\{t>0:J_{t}\neq J_{t-1}\\}$ $\displaystyle\vdots$
$\displaystyle S(\ell)$ $\displaystyle=\min\\{t>S(\ell-1):J_{t}\neq
J_{t-1}\\}.$

In all practical implementations of the sampling algorithms that we discuss,
the conditional flux distributions take the form of empirical distributions of
$N$ walkers:

$\overline{\pi}_{\ell}(dx|k)\approx\frac{1}{\overline{z}_{\ell}^{k}}\sum_{i=1}^{N}w_{S^{i}(\ell)}^{i}\delta_{X_{S^{i}(\ell)}^{i}}(dx)\mathbbm{1}_{\\{k\\}}\left(J_{S^{i}(\ell)}^{i}\right),$
(4)

where $\delta_{x}$ is the Dirac delta function centered at position $x$, and
the normalization constant $\overline{z}_{\ell}^{k}$ is the total weight in
region $k$ at time $S(\ell)$:

$\displaystyle\overline{z}_{\ell}^{k}$
$\displaystyle=P[J_{S(\ell)}=k]\approx\sum_{i=1}^{N}w_{S^{i}(\ell)}^{i}\mathbbm{1}_{\\{k\\}}\left(J_{S^{i}(\ell)}^{i}\right).$
(5)

The outputs of the weighted ensemble iteration and the accelerated variants
that we consider here are approximations of the steady-state conditional flux
distributions, as well as estimates of the corresponding normalization
constants. The conditional flux distributions are related to joint flux
distributions through the normalization constants:

$\displaystyle\overline{\pi}_{\ell}(dx,k)$ $\displaystyle=P[X_{S(\ell)}\in
dx,J_{S(\ell)}=k]$ (6)
$\displaystyle=\overline{z}_{\ell}^{k}\overline{\pi}_{\ell}(dx|k),$ (7)

and we work with $\overline{\pi}_{\ell}(dx,k)$ below.

In particular, we now write the evolution and resampling steps in terms of
operators. For the former, we define the operator $\mathcal{U}$, which
propagates the joint flux distribution of $(X_{S(\ell)},J_{S(\ell)})$ for a
duration of $S(\ell+1)-S(\ell)$:

$\mathcal{U}\overline{\pi}_{\ell}(dx,k)=\overline{\pi}_{\ell+1}(dx,k).$ (8)

In the long-time limit, this yields the eigenequation

$\overline{\pi}(dx,k)=\mathcal{U}\overline{\pi}(dx,k),$ (9)

where $\overline{\pi}$ is the steady-state joint flux distribution. In
practice, one approximates the distributions through walkers, and we denote
the corresponding evolution by $\tilde{\mathcal{U}}$. Mathematically,

$\tilde{\mathcal{U}}\left[\frac{1}{N}\sum_{i=1}^{N}w_{S^{i}(\ell)}^{i}\delta_{X_{S^{i}(\ell)}^{i}}(dx)\mathbbm{1}_{\\{k\\}}\left(J_{S^{i}(\ell)}^{i}\right)\right]\\\
=\frac{1}{N}\sum_{i=1}^{N}w_{S^{i}(\ell)}^{i}\delta_{X_{S^{i}(\ell+1)}^{i}}(dx)\mathbbm{1}_{\\{k\\}}\left(J_{S^{i}(\ell+1)}^{i}\right).$
(10)

The subscript on the weight does not change because the evolution step only
affects the state and index of a walker, not its weight. $\tilde{\mathcal{U}}$
as defined in (10) is an unbiased stochastic approximation of $\mathcal{U}$ in
the sense that, for any function $g$ and any distribution $\rho$,

$\sum_{k}\int g(x,k)[\mathcal{U}\rho](dx,k)=\\\ \mathbbm{E}\left[\sum_{k}\int
g(x,k)[\tilde{\mathcal{U}}\rho](dx,k)\right].$ (11)

If distributions are represented with some ansatz such as a neural network,
one can also apply an approximate propagator based on a variational method as
outlined in Refs. 26, 27. Such a scheme would not require new trajectories to
be run at each iteration, nor would it require a resampling step.

As noted above, we also represent resampling through an operator,
$\mathcal{R}_{k}$. We define it such that, for any distribution $\rho$ and any
function $g$

$\sum_{j}\int
g(x,j)\left(\frac{1}{\zeta_{k}}\mathbbm{1}_{\\{k\\}}(j)\rho(dx,j)\right)=\\\
\mathbbm{E}\left[\sum_{j}\int g(x,j)\mathcal{R}_{k}\rho\right]$ (12)

and

$\zeta_{k}=\sum_{j}\int\mathbbm{1}_{\\{k\\}}(j)\rho(dx,j).$ (13)

Above, the LHS selects the portion of the distribution $\rho$ in region $k$
and then renormalizes it, while the RHS corresponds to resampling from the
distribution. A simple choice which is commonly employed for the operator
$\mathcal{R}_{k}$ is to sample the set of walkers in region $k$ at the end of
the evolution (denoted $K$) from a multinomial distribution with probability
proportional to
$\\{w_{S^{i}(\ell)}^{i}\mathbbm{1}_{\\{k\\}}(J_{S^{i}(\ell)}^{i})\\}_{i\in K}$
with $N_{k}$ trials and then return the distribution

$\mathcal{R}_{k}\left[\frac{1}{N}\sum_{i=1}^{N}w_{S^{i}(\ell)}^{i}\delta_{X_{S^{i}(\ell)}^{i}}(dx)\mathbbm{1}_{\\{k\\}}(J_{S^{i}(\ell)}^{i})\right]\\\
=\frac{1}{N_{k}}\sum_{i\in K}\delta_{X_{S^{i}(\ell)}^{i}}(dx)$ (14)

and normalization

$\overline{z}_{\ell}^{k}\approx\sum_{i}\mathbbm{1}_{\\{k\\}}(J_{S^{i}(\ell)}^{i})w_{S^{i}(\ell)}^{i}.$
(15)

In words, (14) shows that ${\cal R}_{k}$ constructs the state distribution
from a sum over the resampled walkers in each region. There are other possible
ways of resampling in the walker representation, such as stratified or pivotal
resampling [28], and these may in practice be better. However, any choice that
satisfies (12) suffices for our discussion.

Finally, we connect the flux distributions $\overline{\pi}(dx|k)$ to the
steady-state distribution $\pi(x,k)$ for the process $(X_{t},J_{t})$, where
the absence of superscripts again indicates the underlying Markov process
rather than a walker in an ensemble. We do this with a key identity, which we
now state. For a function $g$ of a length $\tau$ trajectory,

$\mathbbm{E}_{(x,k)\sim\pi}[g(X_{0},J_{0},\dots,X_{t},J_{t})]=\frac{\sum_{k}\overline{z}^{k}\int\overline{\pi}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}\left[\sum_{t=0}^{S(1)-1}g(X_{t},J_{t},\dots,X_{t+\tau},J_{t+\tau})\right]}{\sum_{k}\overline{z}^{k}\int\overline{\pi}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}\left[S(1)\right]},$
(16)

where now $\pi$ is the steady-state joint distribution of $(X_{t},J_{t})$.
This identity says that, for each region $k$, we draw initial walker states
$X_{0}$ from $\overline{\pi}(dx|k)$ and set $J_{0}=k$; we evolve those walkers
until they leave region $k$; then we compute averages over them and weight the
contribution from the walkers that started in region $k$ by
$\overline{z}^{k}$, which is the steady-state limit of
$\overline{z}_{\ell}^{k}$. We derive (16) in Appendix A. The distribution
representation of WE is summarized in Algorithm 1.

Algorithm 1 WE (distribution representation)

1: Approximate Markov propagator $\tilde{\mathcal{U}}$, resampling operator
$\mathcal{R}_{k}$, starting flux distributions
$\\{\overline{\pi}_{0}(dx|k)\\}_{k=1}^{n}$, initial region weights
$\\{\overline{z}_{0}^{k}\\}_{k=1}^{n}$ with $\sum_{k}\overline{z}_{0}^{k}=1$.

2:for $\ell=0,\dots,L$ do

3:
$\overline{\pi}_{\ell}(dx,k)\leftarrow\overline{z}_{\ell}^{k}\overline{\pi}_{\ell}(dx|k)$

4:
$\tilde{\overline{\pi}}\leftarrow[\tilde{\mathcal{U}}\overline{\pi}_{\ell}](dx,k)$
$\triangleright$ One Power Iteration Step

5: for $k=1,\dots,n$ do

6:
$\overline{\pi}_{\ell+1}(dx|k),\overline{z}_{\ell+1}^{k}\leftarrow\mathcal{R}_{k}[\tilde{\overline{\pi}}]$

7: end for

8:end for

9:return $\overline{\pi}_{L}(dx|k)$ and $\overline{z}_{L}^{k}$

We see that each step of WE applies the operator $\mathcal{U}$ to the previous
distribution $\overline{\pi}_{\ell}(dx,k)$. Therefore, WE can be seen as
performing a power iteration in ${\cal U}$, starting from some initial
distribution. The resampling step in WE plays a role analogous to the
normalization step in standard power iteration in the following sense. In
standard power iteration, the normalization step serves to prevent the iterate
from becoming too small or too large and introducing numerical instability. In
WE, the resampling step serves to prevent any of the region distributions from
becoming poorly sampled, which would increase variance. However, the
normalization step in power iteration and, in turn, the resampling step in WE
does not accelerate decay of initialization bias, nor is the WE
autocorrelation time reduced relative to direct sampling.

Algorithm 1 suggests that, if $\tilde{\overline{\pi}}=\overline{\pi}$ in step
3, convergence would be achieved. Therefore, our strategy for accelerating WE
is to replace the power iteration step with a more accurate approximation.
Specifically, we compute changes of measure $\nu_{\ell}$ such that

$\overline{\pi}(dx,k)\approx\nu_{\ell}(x,k)\overline{\pi}_{\ell}(dx,k).$ (17)

If we set up our approximation such that whenever $\pi_{\ell}=\pi$, we get
$\nu_{\ell}=1$, then the resulting modification of WE will have a fixed point
at the true steady-state distribution, and there will be no approximation
error resulting from deficiencies in the approximation algorithm’s specific
ansatz. Our scheme has the advantage that it can potentially approach steady
state much faster than WE. This is the idea that we develop in Sections III
and IV.

##  III Nonequilibrium umbrellla sampling (NEUS)

We now present NEUS as a simple extension to WE that accelerates convergence
in the way suggested in (17); EM [6] can be formulated similarly. Our
development is based on the algorithm in Ref. 18, which corrects a small
systematic bias in earlier NEUS papers [3, 15, 16, 4, 17]. We begin by making
the observation that, at steady state,

$\mathbbm{E}_{X_{0}\sim\overline{\pi}}[\mathbbm{1}_{\\{k\\}}(J_{0})]=\mathbbm{E}_{X_{0}\sim\overline{\pi}}[\mathbbm{1}_{\\{k\\}}(J_{S(1)})],$
(18)

or, in terms of the density we wish to compute,

$\int\overline{\pi}(dx,k)=\sum_{j}\int\overline{\pi}(dx,j)\mathbbm{E}_{X_{0}=x,J_{0}=j}[\mathbbm{1}_{\\{k\\}}(J_{S(1)})].$
(19)

This observation introduces one equation per possible value of the index $j$.
Therefore, we can parameterize the change of measure with as many free
parameters:

$\overline{\pi}(dx,k)\approx
c_{\ell}^{k}\overline{z}_{\ell}^{k}\overline{\pi}_{\ell}(dx|k)$ (20)

Substituting this ansatz into (19) and simplifying yields the matrix equation

$c_{\ell}^{k}\overline{z}_{\ell}^{k}=\sum_{j}c_{\ell}^{j}\overline{z}_{\ell}^{j}\overline{G}_{jk},$
(21)

where

$\overline{G}_{jk}=\int\overline{\pi}(dx|j)\mathbbm{E}_{X_{0}=x,J_{0}=j}[\mathbbm{1}_{\\{k\\}}(J_{S(1)})].$
(22)

We summarize NEUS in Algorithm 2.

Algorithm 2 NEUS

1: Approximate Markov propagator $\tilde{\mathcal{U}}$, Resampling operator
$\mathcal{R}_{k}$, starting flux distributions
$\\{\overline{\pi}_{0}(dx|k)\\}_{k=1}^{n}$, initial region weights
$\\{\overline{z}_{0}^{k}\\}_{k=1}^{n}$ with $\sum_{k}\overline{z}_{0}^{k}=1$.

2:for $\ell=0,\dots,L$ do

3:
$\overline{G}_{jk}\leftarrow\int\overline{\pi}_{\ell}(dx|j)\mathbbm{E}_{X_{0}=x,J_{0}=j}[\mathbbm{1}_{\\{k\\}}(J_{S(1)})]$

4: Solve $c_{\ell}\overline{z}_{\ell}\overline{G}=c_{\ell}\overline{z}_{\ell}$
for $c_{\ell}\overline{z}_{\ell}$

5: $\overline{\pi}_{\ell}(dx,k)\leftarrow
c_{\ell}^{k}\mathbbm{1}_{\\{k\\}}\overline{z}_{\ell}^{k}\overline{\pi}_{\ell}(dx|k)$

6:
$\tilde{\overline{\pi}}\leftarrow\tilde{\mathcal{U}}\overline{\pi}_{\ell}(dx,k)$

7: for $k=1,\dots,n$ do

8:
$\overline{\pi}_{\ell+1}(dx|k),\overline{z}_{\ell+1}^{k}\leftarrow\mathcal{R}_{k}[\tilde{\overline{\pi}}]$

9: end for

10:end for

11:return $\overline{\pi}_{L}(dx|k)$ and $\overline{z}_{L}^{k}$

We thus see that NEUS is distinguished from WE by steps 2 and 3. In practice,
walkers are drawn and their dynamics are simulated to determine $S(\ell+1)$ in
the loop over $\ell$ prior to step 2 (the resulting state distribution is used
later in step 5); then, $\overline{G}_{jk}$ is computed in step 2 from walkers
at iteration $\ell$ as

$\overline{G}_{jk}\approx\frac{\sum_{i}\mathbbm{1}_{\\{j\\}}\left(J^{i}_{S^{i}(\ell)}\right)\mathbbm{1}_{\\{k\\}}\left(J^{i}_{S^{i}(\ell+1)}\right)}{\sum_{i}\mathbbm{1}_{\\{j\\}}\left(J^{i}_{S^{i}(\ell)}\right)}.$
(23)

We solve for the product $c_{\ell}\overline{z}_{\ell}$ directly and update the
approximation for $\overline{\pi}$ accordingly. We then proceed as in WE.

We note that, if the flux distributions $\overline{\pi}(dx|k)$ and
normalization constants $\overline{z}^{k}$ assume their steady-state values,
then (21) is satisfied by $c^{k}=1$, and NEUS has the correct fixed point.

##  IV BAD-NEUS

To improve on NEUS, it is necessary to improve the approximation of the change
of measure in (17). Here we do so by introducing a basis expansion. This
allows us to vary the expressivity of the approximation through the number of
basis functions, which can exceed the number of regions, in contrast to the
steady-state condition in (18).

To this end, we note that, for any lag time $\tau$ and any function $f(x,k)$,

$\mathbbm{E}_{(X_{0},J_{0})\sim\pi}[f(X_{0},J_{0})-f(X_{\tau},J_{\tau})]=0.$
(24)

We stress that this relation holds for any fixed time $\tau$. Expanding this
expectation using (16), with
$g(X_{0},J_{0},\dots,X_{\tau},J_{\tau})=f(X_{0},J_{0})-f(X_{\tau},J_{\tau})$
and multiplying through by the normalization in the denominator, we obtain

$0=\sum_{k}\overline{z}^{k}\int\overline{\pi}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}\left[\sum_{t=0}^{S(1)-1}f(X_{t},J_{t})-f(X_{t+\tau},J_{t+\tau})\right].$
(25)

We then introduce the basis set $\\{\phi_{p}(x,k)\\}_{p}$, and write

$\overline{\pi}(dx,k)=\overline{z}^{k}\overline{\pi}(dx|k)\approx\sum_{p}c_{\ell}^{p}\phi_{p}(x,k)\overline{z}_{\ell}^{k}\overline{\pi}_{\ell}(dx|k).$
(26)

Making the choice $f=\phi_{r}$ and inserting the basis expansion into (25)
gives

$0=\sum_{kp}c_{\ell}^{p}\overline{z}_{\ell}^{k}\int\phi_{p}(x,k)\overline{\pi}_{\ell}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}\left[\sum_{t=0}^{S(1)-1}\phi_{r}(X_{t},J_{t})-\phi_{r}(X_{t+\tau},J_{t+\tau})\right].$
(27)

This is a linear system which can be solved for the expansion coefficients:

$0=\sum_{p}c_{\ell}^{p}M_{pr}.$ (28)

with the matrix entries given by

$M_{pr}=\sum_{k}\overline{z}_{\ell}^{k}\int\phi_{p}(x,k)\overline{\pi}_{\ell}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}\left[\sum_{t=0}^{S(1)-1}\phi_{r}(X_{t},J_{t})-\phi_{r}(X_{t+\tau},J_{t+\tau})\right].$
(29)

The basis set must include the constant function in its span so that the
system has a unique nontrivial solution. The NEUS algorithm corresponds to the
choice

$\phi_{p}(x,k)=\mathbbm{1}_{\\{p\\}}(k)$ (30)

and $\tau=1$. We summarize the BAD-NEUS algorithm in Algorithm 3.

Consistent with WE and NEUS, in practice, we estimate the integral $M_{pr}$
through a sum over walkers:

$M_{pr}\approx\sum_{ki}\overline{z}_{\ell}^{k}\phi_{p}(X^{i}_{S^{i}(\ell)},k)\mathbbm{1}_{\\{k\\}}(J^{i}_{0})\left[\sum_{t=S^{i}(\ell)}^{S^{i}(\ell+1)-1}(\phi_{r}(X^{i}_{t},J_{t}^{i})-\phi_{r}(X^{i}_{t+\tau},J_{t+\tau}^{i}))\right].$
(31)

Given $M$ and, in turn, the estimated expansion coefficients, we update the
weights as

$w_{S^{i}(\ell)}^{i}=\frac{\overline{z}_{\ell}^{J_{S^{i}(\ell)}^{i}}}{Z}\sum_{p}c_{\ell}^{p}\phi_{p}\left(X_{S^{i}(\ell)}^{i},J_{S^{i}(\ell)}^{i}\right),$
(32)

where $Z$ normalizes the total of the walker weights to one. If the basis set
contains the constant function in its span, and the flux distributions
$\overline{\pi}_{\ell}(dx|k)$ and region weights $\overline{z}_{\ell}^{k}$
take their steady-state values, then the linear system (27) is solved by
$\sum_{p}c_{p}\phi(x,k)=1$, and BAD-NEUS has the correct fixed point. To
reduce variance in our estimate of the matrix $M$ to ensure a stable solve, it
is often useful to work with the mixture distribution for the last $h$
iterations:

$\overline{\pi}_{\ell,h}(dx|k)=\frac{1}{h}\sum_{t=\ell-h+1}^{\ell}\overline{\pi}_{t}(dx|k).$
(33)

When working with walkers, this mixture distribution corresponds to
concatenating the walkers from $h$ prior iterations, and therefore contains
more data than using a single iteration. We simply substitute
$\overline{\pi}_{\ell}$ for $\overline{\pi}_{\ell,h}$ in all of our
algorithms. A unified walker-based algorithm that we use in practice is given
in Algorithm 4, with the additional steps needed for NEUS and BAD-NEUS in
Algorithm 5.

Finally, we note that the general strategy of improving the approximation of
the change of measure is not limited to using a basis expansion. Just as one
can use various means to solve the equations of the operator that encodes the
statistics of the dynamics [22, 23, 29, 27], one can represent the change of
measure here by either a basis expansion or a nonlinear representation. In
particular, our tests based on the neural-network approach in Ref. 27 show
significant promise (unpublished results).

Algorithm 3 BAD-NEUS

1: Approximate Markov propagator $\tilde{\mathcal{U}}$, Resampling operator
$\mathcal{R}_{k}$, starting flux distributions
$\\{\overline{\pi}_{0}(dx|k)\\}_{k=1}^{n}$, initial region weights
$\\{\overline{z}_{0}^{k}\\}_{k=1}^{n}$ with $\sum_{k}\overline{z}_{0}^{k}=1$,
basis set $\\{\phi_{p}(x,k)\\}_{p}$, lag time $\tau$

2:for $\ell=0,\dots,L$ do

3: $\overline{\pi}_{\ell}\leftarrow\overline{\pi}_{\ell,h}$

4:
$\Delta\phi\leftarrow\sum_{t=0}^{S(1)-1}\phi_{r}(X_{t},J_{t})-\phi_{r}(X_{t+\tau},J_{t+\tau})$

5:
$M_{pr}\leftarrow\sum_{k}\int\overline{z}_{\ell}^{k}\phi_{p}(x,k)\overline{\pi}_{\ell}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}\left[\Delta\phi\right]$

6: Solve $cM=0$

7: $\overline{\pi}_{\ell}(dx,k)\leftarrow\sum
c_{\ell}^{p}\phi_{p}(x,k)\overline{z}_{\ell}^{k}\overline{\pi}_{\ell}(dx|k)$

8:
$\tilde{\overline{\pi}}\leftarrow\tilde{\mathcal{U}}\overline{\pi}_{\ell}(dx,j)$

9: for $k=1,\dots,n$ do

10:
$\overline{\pi}_{\ell+1}(dx|k),\overline{z}_{\ell+1}^{k}\leftarrow\mathcal{R}_{k}[\tilde{\overline{\pi}}]$

11: end for

12:end for

13:return $\overline{\pi}_{L}(dx|k)$ and $\overline{z}_{L}^{k}$

Algorithm 4 Trajectory Stratification (walker representation)

1: Starting ensemble of walkers
$E_{0}=\\{(X_{0}^{i},J_{0}^{i},w_{0}^{i})\\}_{i=1}^{N},$ with
$\sum_{i=1}^{N}w_{0}^{i}=1$, number of past iterations to retain $h$, total
iterations $L,$ number of strata $n.$ lag time $\tau$

2:for $\ell=0,\dots,L$ do

3: $N\leftarrow\text{length}(E_{\ell})$

4: Generate a list of trajectories
$T_{\ell}\leftarrow\\{[(X_{S^{i}(\ell)}^{i},J_{S^{i}(\ell)}^{i},w_{S^{i}(\ell)}^{i}),(X_{S^{i}(\ell)+1}^{i},J_{S^{i}(\ell)+1}^{i},w_{S^{i}(\ell)}^{i}),\dots,$$(X_{S^{i}(\ell+1)+\tau}^{i},J_{S^{i}(\ell+1)+\tau}^{i},w_{S^{i}(\ell)}^{i})]\\}_{i=1}^{N}$,
where
$(X_{S^{i}(\ell)}^{i},J_{S^{i}(\ell)}^{i},w_{S^{i}(\ell)}^{i})=E_{\ell}[i]$

5: $T_{\ell}\leftarrow\text{concatenate}(\\{T_{\ell},T_{\ell-1},\dots
T_{\ell-h+1}\\})$

6: $N\leftarrow\text{length}(T_{\ell})$

7: Renormalize the weights in $T_{\ell}$ by dividing each by $h$.

8: Update weights $\\{w^{i}_{S^{i}(\ell)}\\}_{i=1}^{N}$ using Algorithm 5.
(Omit for Weighted Ensemble)

9: Initialize an empty list $E_{\ell+1}=\\{\\}$

10: for $k=1,\dots,n$ do

11:
$\overline{z}^{k}_{\ell+1}\leftarrow\sum_{i=1}^{N}w_{S^{i}(\ell)}^{i}\mathbbm{1}_{\\{k\\}}(J^{i}_{S^{i}(\ell+1)})$

12: for $r=1,\dots,N_{k}$ do

13: Sample an index $b$ with probability proportional to
$w_{S^{b}(\ell)}^{b}\mathbbm{1}_{\\{k\\}}(J^{b}_{S^{b}(\ell+1)})$

14: Append to $E_{\ell}$ the configuration
$(X_{S^{b}(\ell+1)}^{b},J_{S^{b}(\ell+1)}^{b},\overline{z}^{k}_{\ell+1}/N_{k})$

15: end for

16: end for

17:end for

18:return $T_{L}$ and $\overline{z}_{L}^{k}$

Algorithm 5 Approximating the steady-state flux distribution (walker
representation)

1: List of trajectories $T_{\ell}$, lag time $\tau$, basis functions
$\phi_{p}$, bin weights $\overline{z}_{\ell}$

2:$M_{pr}\leftarrow\sum_{ki}\overline{z}_{\ell}^{k}\phi_{p}(X^{i}_{S^{i}(\ell)},k)\mathbbm{1}_{\\{k\\}}(J^{i}_{0})\left[\sum_{t=S^{i}(\ell)}^{S^{i}(\ell+1)-1}(\phi_{r}(X^{i}_{t},J_{t}^{i})-\phi_{r}(X^{i}_{t+\tau},J_{t+\tau}^{i}))\right]$

3:Solve $c_{\ell}M=0$

4:$w_{S^{i}(\ell)}^{i}\leftarrow(\overline{z}_{\ell}^{J_{S^{i}(\ell)}^{i}}/Z)\sum_{p}c_{\ell}^{p}\phi_{p}\left(X_{S^{i}(\ell)}^{i},J_{S^{i}(\ell)}^{i}\right)$

5:return $\\{w_{S^{i}(\ell)}^{i}\\}_{i=1}^{N}$

##  V Two-dimensional test system

We first illustrate our approach by sampling a two-dimensional system, which
enables us to compare estimates of the steady-state distribution to the
Boltzmann probability to ensure that simulations are run until convergence.
Specifically, the system $x=(u,v)$ is defined by the Müller-Brown potential
[30], which is a sum of four Gaussian functions:

$V(u,v)=\frac{1}{20}\sum_{i=1}^{4}C_{i}\exp[a_{i}(u-v_{i})^{2}\\\
+b_{i}(u-u_{i})(v-v_{i})+c_{i}(v-v_{i})^{2}].$ (34)

For all results shown, we use $C_{i}=\\{-200,-100,-170,15\\}$,
$a_{i}=\\{-1,-1,-6.5,0.7\\}$, $b_{i}=\\{0,0,11,0.6\\}$,
$c_{i}=\\{-10,-10,-6.5,0.7\\}$, $u_{i}=\\{1,-0.27,-0.5,-1\\}$,
$v_{i}=\\{0,0.5,1.5,1\\}$. The potential with these parameter choices is shown
in Figure 1. The presence of multiple metastable states and a minimum energy
pathway that does not parallel the axes of the variables used for the
numerical integration make this system representative of difficulties commonly
encountered in molecular simulations. Because the model is two-dimensional, we
can readily visualize results and compare them with statistics independently
computed using the grid-based scheme in Ref. 27.

We evolve the system with overdamped Langevin dynamics, discretized with the
BAOAB algorithm [31]:

$X_{t+dt}=X_{t}-\nabla V(X_{t})dt+\sqrt{\frac{dt}{2\beta}}(Z_{t}+Z_{t-dt}),$
(35)

where $dt$ is the time step, $\beta$ is the inverse temperature, and
$Z_{t}\sim N(0,I)$ is a random vector with components drawn from the normal
distribution with zero mean and unit standard deviation ($I$ is the two-
dimensional identity matrix). For all results shown, we use $dt=0.001$ and
$\beta=2$.

Figure 1:  Müller-Brown potential. Orange and red ovals indicate states $A$
and $B$, respectively.

###  V.1 Comparison of algorithms

A key point of our theoretical development is that NEUS and BAD-NEUS are
simple elaborations of WE. This enables us to define a unified walker-based
algorithm that we use in practice (Algorithm 4, with the added operations
needed for NEUS and BAD-NEUS in Algorithm 5). Splitting and stratification are
defined through the rule for switching the index process. The sequence of
stopping times $S(\ell)$ is then determined by the index process through (3).
The update rule for $J_{t}$ that we employ for all three algorithms is as
follows. Let $\\{\psi_{k}(x)\\}_{k=1,\ldots,n}$ be a set of nonnegative
functions. If $\psi_{J_{t}}(X_{t+dt})>0$, then $J_{t+dt}=J_{t}$; otherwise,
$\mathbbm{P}[J_{t+dt}=k]=\psi_{k}(X_{t+dt})/\sum_{k}\psi_{k}(X_{t+dt})$. That
is, the value of the index process remains the same until the walker leaves
the support of $\psi_{J_{t}}$, and then a new index value $k$ is drawn with
likelihood proportional to $\psi_{k}(X_{t+dt})$. For the Müller-Brown model,
we use

$\psi_{k}(u,v)=\begin{cases}1&\text{if}\ |v-v_{k}|<\varepsilon_{k}\
\text{and}\ 1<k<n\\\ 1&\text{if}\ v-v_{k}<\varepsilon_{k}\ \text{and}\ k=1\\\
1&\text{if}\ v_{k}-v>\varepsilon_{k}\ \text{and}\ k=n\\\
0&\text{otherwise}.\end{cases}$ (36)

Unless otherwise indicated, we set $n=10$, space the $v_{k}$ uniformly in the
interval $[-0.2,1.8]$, and set $\varepsilon_{k}=0.6(v_{k+1}-v_{k})$, so that
the regions of support (strata) overlap. This choice prevents walkers in
barrier regions from rapidly switching back and forth between index values,
limiting the overhead of the algorithms. Unless otherwise indicated, we use
2000 walkers per region and run them until their index processes switch
values. Initial configurations for $J_{0}=k$ are generated by uniformly
sampling the support of $\psi_{k}$.

For BAD-NEUS, we use a basis set consisting of 10 indicator functions per
stratum. The indicator functions are determined by clustering all the samples
in a stratum with $k$-means clustering and partitioning it into Voronoi
polyhedra based on the cluster means. Unless otherwise indicated, we compute
statistics using data from the last $h=3$ iterations and use a lag time of
$\tau=10$ time steps.

We measure convergence by computing the root mean square (RMS) difference in
$\ln(\tilde{\pi}(x))$ from $-\beta V$, where $\tilde{\pi}(x)$ is an estimate
of the steady-state distribution from the normalized histogram of samples. For
the histogram, we use a uniform $50\times 50$ grid on the rectangle defined by
the minimum and maximum values in the NEUS dataset. Since some grid regions
are empty because they correspond to energies too high for NEUS to sample, we
only compute errors over bins with $V<7$. We stop each simulation when the RMS
difference drops below one. Results are shown in Figure 2. We see that WE
requires about 73 times more iterations than NEUS, which in turn requires
about 13 times more iterations than BAD-NEUS to reach the convergence
criterion.

Figure 2:  Convergence of the steady-state distribution of the Müller-Brown
system for different algorithms as indicated.

###  V.2 Choice of hyperparameters

Having demonstrated that trajectory stratification outperforms WE for this
system, we now examine the effect of key hyperparameters in NEUS and BAD-NEUS.

Figure 3:  Effect of the number of strata on convergence of BAD-NEUS applied
to the Müller-Brown system. The number of iterations required (top) and the
total number of time steps (bottom)—a proxy for computational effort—until
convergence. The total number of walkers is held fixed. Error bars are the
standard deviation over 10 replicate simulations. For these simulations, there
are 40,000 walkers total, and the lag time is $\tau=1$ time step.

First, we investigate the effect of the number of strata. To compare
simulations that require essentially the same resources, we hold the total
number of walkers fixed at 40,000, and we divide the walkers uniformly across
the strata. For these simulations, we use a lag time of $\tau=1$ and retain
$h=3$ past iterations. In Figure 3 we plot both the number of iterations to
reach convergence and the total number of time steps needed for convergence as
we vary the number of strata. We see that the number of BAD-NEUS iterations
required for convergence increases linearly with the number of strata.
However, the lengths of trajectories decrease because the strata are spaced
more closely. Due to this interplay, the total effort to achieve convergence
decreases until about 12 strata and then levels off. While actual performance
will depend on the overhead incurred by stopping and starting the dynamics
engine and computing the BAD-NEUS weights, these results indicate that finer
stratification is better.

Figure 4:  Effect of retaining past iterations on convergence of NEUS (top)
and BAD-NEUS (bottom) applied to the Müller-Brown system.  Figure 5:
Dependence of convergence on the lag time for BAD-NEUS applied to the Müller-
Brown system. Error bars are the standard deviation over 10 replicate
simulations.

Next, we examine the dependence of convergence on the number of iterations
used to compute statistics, $h$ (other hyperparameters are set to the default
values given above). Using data from a larger number of iterations allows for
more averaging, but it can also contaminate the statistics with samples
obtained before the steady-state distribution has converged. While the trends
are clearer for NEUS than BAD-NEUS because the former converges more slowly,
Figure 4 indicates that retaining fewer iterations is better for both NEUS and
BAD-NEUS. It is therefore important to use enough walkers per iteration that
data from past iterations need not be retained. If one has access to a large
number of processing elements, and many walkers can be simulated in parallel,
this is not a severe restriction. Once convergence is achieved one can begin
accumulating data to reduce variance.

Finally, we investigate the impact of the lag time in Figure 5 (other
hyperparameters are set to the default values given above). We find that using
longer lag times weakly decreases the number of iterations required for
convergence.

###  V.3 Kinetic statistics

A common problem in molecular dynamics is that while equilibrium averages are
relatively straightforward to calculate from biased simulations (via, for
example, the various methods reviewed in Ref. 32), dynamical averages are much
harder. In this section, we use BAD-NEUS to compute dynamical averages
efficiently. To do so, we split the ensemble of transition paths based on the
last metastable state visited as previously for NEUS [33, 16, 25]. We
specifically compute the backward committor, $q_{-}$, which is the probability
that the system last visited a reactant state $A$ rather than a product state
$B$, and the transition path theory (TPT) rate, defined as the mean number of
$A$ to $B$ transitions per unit time divided by the fraction of time spent
last in $A$. The former can be obtained without saving additional information
once we split the path ensemble. The forward committor, $q_{+}$, which is the
probability that the system next visits state $B$ rather than state $A$, can
also be obtained from NEUS but requires saving additional information [25].
Here, both of the systems that we consider are microscopically reversible, so
we can obtain $q_{+}$ from $q_{+}=1-q_{-}$.

For the Müller-Brown system, we consider the following history dependent
stratification. Let $\\{\psi_{k}^{A}(x)\\}_{k=1,\ldots,n}$ and
$\\{\psi_{k}^{B}(x)\\}_{k=1,\ldots,m}$ be non-negative functions, let
$\psi_{k}=\psi_{k}^{A}$ if $k\leq n$; otherwise, let
$\psi_{k}=\psi_{k-n}^{B},$ and let $\mathbbm{1}_{A}$ and $\mathbbm{1}_{B}$ be
indicator functions on sets $A$ and $B$, respectively. Let $T^{-}_{A\cup
B}(t)$ be the time the system sampled at time $t$ was last in $A$ or $B$, so
that $\mathbbm{1}_{A}(X_{T^{-}_{A\cup
B}(t)})=1-\mathbbm{1}_{B}(X_{T^{-}_{A\cup B}(t)})$ reports whether $A$ rather
than $B$ was last visited. Define the index process by the following update
rule:

$\mathbbm{P}[J_{t+dt}=k]=\begin{cases}\psi_{k}(X_{t+dt})\mathbbm{1}_{A}(X_{T^{-}_{A\cup
B}(t+dt)})/\sum_{i}\psi_{i}^{A}(X_{t+dt})&k\leq
n,\psi_{J_{t}}(X_{t+dt})=0\text{ or }\\\
&\quad\quad\quad\mathbbm{1}_{A}(X_{T^{-}_{A\cup
B}(t+dt)})\neq\mathbbm{1}_{A}(X_{T^{-}_{A\cup B}(t)})\\\
\psi_{k}(X_{t+dt})\mathbbm{1}_{B}(X_{T^{-}_{A\cup
B}(t+dt)})/\sum_{i}\psi_{i}^{B}(X_{t+dt})&k>n,\psi_{J_{t}}(X_{t+dt})=0\text{
or }\\\ &\quad\quad\quad\mathbbm{1}_{A}(X_{T^{-}_{A\cup
B}(t+dt)})\neq\mathbbm{1}_{A}(X_{T^{-}_{A\cup B}(t)})\\\
\mathbbm{1}_{\\{k\\}}(J_{t})&\text{otherwise}\end{cases}$ (37)

Thus there are $n+m$ total strata, the first $n$ of which contain walkers that
last visited $A$, and the remainder of which contain walkers that last visited
$B$. The steady-state distribution of walkers which last visited $A$ is
proportional to $\pi(x)q_{-}(x)$, and that which last visited $B$ is
proportional to $\pi(1-q_{-}(x))$. This ensemble lets us compute dynamical
averages. The backward committor projected onto a space of CVs can be computed
using

$q_{-}^{\theta}(s)=\frac{\mathbbm{E}[\mathbbm{1}_{\\{0,\ldots,n\\}}(J_{t})\mathbbm{1}_{\\{ds\\}}(\theta(X_{t}))]}{\mathbbm{E}[\mathbbm{1}_{\\{ds\\}}(\theta(X_{t}))]},$
(38)

where $\theta$ is a vector with components that are the CVs and $ds$ is a bin
in the space. The TPT rate constant is defined as

$k_{AB}=\frac{\operatorname{\mathbbm{E}}_{\pi}[q_{-}(x)\mathbbm{1}_{B^{c}}(x)\mathcal{T}_{1}\mathbbm{1}_{B}(x)]}{\operatorname{\mathbbm{E}}_{\pi}[q_{-}(x)]},$
(39)

where $\mathcal{T}_{t}$ is the operator that describes the evolution of
expectations of functions:

$\mathcal{T}_{t}g(x,j)=E_{X_{0}=x,J_{0}=j}[f(X_{t},J_{t})].$ (40)

In practice, we compute (39) from (16) by choosing

$g(X_{0},J_{0}...,X_{t},J_{t})=\\\
\mathbbm{1}_{B^{c}}(X_{0})\mathbbm{1}_{B}(X_{1})\mathbbm{1}_{\\{0,\ldots,n\\}}(J_{0})$
(41)

for the numerator and

$g(X_{0},J_{0}...,X_{t},J_{t})=\mathbbm{1}_{\\{0,\ldots,n\\}}(J_{0})$ (42)

for the denominator. This corresponds to the steady-state flux into $B$ from
trajectories that originated in $A$.

Figure 6:  Convergence of the TPT inverse rate estimate for the Müller-Brown
system. The green line is the grid-based reference obtained from 44. For these
simulations, there are 2000 walkers per stratum, $h=4$ iterations are
retained, and the lag time is $\tau=10$ time steps. Figure 7:  Illustration of
the backward committor computed from BAD-NEUS. (top) Reference. (bottom) BAD-
NEUS results. (left) Backward committor. (right) $\ln(q_{-}/(1-q_{-}))$, which
emphasizes states near $A$ and $B$. These results are obtained from the same
simulations as Fig. 6.

For the Müller-Brown system, we define states $A$ and $B$ as

$\displaystyle A=\\{y,v:6.5(u+0.5)^{2}-11(u+0.5)(v-1.5)+6.5(u-1.5)^{2}<0.3\\}$
(43) $\displaystyle B=\\{u,v:(u-0.6)^{2}+0.5(v-0.02)^{2}<0.2\\}.$

We use a similar index process construction as for the equilibrium
calculations, except there are 10 $v_{k}$ that are evenly spaced in the
interval $[0,1.6]$ for $\psi_{k}^{A}$ and 10 $v_{k}$ that are evenly spaced in
the interval $[-0.3,0.8]$ for $\psi_{k}^{B}$, for a total of 20 strata. We use
different definitions for $\psi_{k}^{A}$ and $\psi_{k}^{B}$ because, for
walkers originating from state $A$ ($B$), a stratum located below (above)
state $B$ ($A$) is unlikely be populated. Both NEUS and BAD-NEUS use 2000
walkers per stratum, and we compute statistics using data from the last $h=3$
iterations. For BAD-NEUS, we use a basis set consisting of 10 indicator
functions per stratum, again based on $k$-means clustering; we use a lag time
of $\tau=1$.

To compute reference values for these kinetic statistics, we use the grid-
based approximation to the generator in Ref. 27, with the same grid
parameters. Since the dynamics are microscopically reversible, we obtain the
backward committor by solving for the forward committor using the approach in
Ref. 27, then setting $q_{-}=1-q_{+}$. We solve for the reaction rate using

$k_{AB}=\frac{2(\vec{\pi
q}_{-})^{\top}(P\vec{\mathbbm{1}}_{B})}{(\vec{\pi}^{\top}\vec{q}_{-})(\beta\epsilon_{x})},$
(44)

where $P$ is the discretized transition matrix defined on the grid,
$\epsilon_{x}$ is the grid spacing, and arrows indicate vectors of function
values on the grid points ($\vec{\pi q}_{-}$ is a single vector of product
values).

Figure 6 shows estimates for the TPT rate. We see that BAD-NEUS converges
several fold faster than NEUS. Each BAD-NEUS iteration requires an average of
407 time units of sampling (arising from 2000 walkers in each of 20 strata,
with an average time before leaving the stratum of 10.2 time steps). Each
iteration is therefore significantly less total computational cost than
generating a single $A$ to $B$ transition on average (1200 time units, as
evidenced by $k_{AB}^{-1}$) and is amenable to parallelization.

Figure 7 shows the backward committor computed at BAD-NEUS iteration 20 using
(38) with the variables of numerical integration as the CVs. We represent the
results in two ways: the committor itself and its logit function. The former
emphasizes the transition region ($q_{-}\approx 0.5$) while the latter
emphasizes values close to states $A$ and $B$ ($q_{-}\approx 0$ and
$q_{-}\approx 1$). Both show excellent agreement with the reference.

##  VI Molecular Test System

As a more challenging test of the method, we estimate the folding rate of the
first 39 amino acids of the N-terminal domain of ribosomal protein L9 (NTL9).
The native secondary structure of NTL9 is $\beta\beta\alpha\beta\alpha$
sequentially. In the native state, the $\beta$-strands form an anti-parallel
$\beta$-sheet, and the $\alpha$-helices pack on either side of it. NTL9(1-39)
lacks the C-terminal helix (Fig. 8). We choose this system because it enables
direct comparison with an earlier computational study [34] (described further
below), which estimated the folding to be on the millisecond timescale,
consistent with experimental studies of (full) NTL9 [35, 36] and earlier
simulations of NTL9(1-39) [37].

Figure 8:  Representative NTL9(1-39) structures drawn from the (left) native
and (right) denatured states in the BAD-NEUS simulation. There is residual
helical structure in the denatured state in both our simulations and those to
which we compare [34].

All molecular dynamics simulations were performed with OpenMM [38]. We compute
CVs including backbone RMSDs and fractions of native contacts using MDtraj
[39]. To enable direct comparison with Ref. 34, we model NTL9(1-39) using the
AMBER FF14SB force field [40] and the Hawkins, Cramer and Truhlar generalized
Born implicit solvent model [41, 42]. We use a Langevin thermostat with a time
step of 2 fs, a temperature of 300. K, and a friction constant of $\gamma=80$
ps-1, which corresponds to the high-viscosity case considered in Ref. 34.

We run the molecular dynamics for each walker in intervals of 20 ps and
compute the CVs to update the index process at the end of each interval. For
BAD-NEUS, after an index process changes values (($S(1)$ in (3)), we run the
walker $\tau$ additional molecular dynamics steps to ensure that we have a
total of $\tau$ steps beyond $S(1)$. While this procedure allows walkers to
run beyond the time that they exit their stratum, it does not bias the
results. Since our algorithm only requires us to stop and start the molecular
dynamics engine while running unbiased dynamics, we do not need to modify
OpenMM in any way.

The stratification is similar to that for the Müller-Brown system. Namely, we
set

$\psi_{k}^{A}(X)=\psi_{k}^{B}(X)=\begin{cases}1&\text{if}\
|Q(X)-Q_{k}|<\varepsilon_{k}\\\ 0&\text{otherwise},\end{cases}$ (45)

where $Q(X)$ is the fraction of native contacts using the definition of Ref.
43. We determine the native contacts from the crystal structure, PDB ID
2HBB[44]. We space 20 $Q_{k}$ values uniformly in the interval $[0.35,0.85]$
and set $\varepsilon_{k}=0.6(Q_{k+1}-Q_{k})$, so that the regions overlap.
Thus there are 40 total strata. We define the folded state as full backbone
RMSD $<0.28$ nm and $Q>0.85$; the denatured state has full backbone RMSD
$>0.80$ nm and $Q<0.35$.

To initialize the simulation, we minimize the energy starting from the PDB
structure, draw velocities from a Maxwell-Boltzmann distribution for 380 K,
and simulate for 80 ns at that temperature; this denatures the protein (Fig.
8). We sample 4000 frames in each of the 40 strata from this trajectory and
weight them uniformly; this set forms our starting walkers,
$\\{X^{i}_{0},J^{i}_{0},1/N\\}_{i=1}^{N}$, where $N=40$. We use a lag time of
10 ps, and we retain a history of $h=4$ iterations. The initialization
pipeline, integrator settings, and stratification choice are the same between
NEUS and BAD-NEUS. The hyperparameter choices are summarized in Table 1.

Table 1: Hyperparameter choices used for NTL9(1-39) simulations Hyperparameter | Kinetic statistics
---|---
Stratification CV |  $Q$ and $\mathbbm{1}_{A}(T^{-}_{A\cup B}(t))$
Stopping rule | Stratum exit
Walkers per stratum | 4000
Number of strata | 40
Number of iterations retained, $h$ | 4
Type of basis set |  $k$-means indicator
Number of basis functions per stratum | 10
Lag time, $\tau$ | 10 ps

To construct the basis set for BAD-NEUS, we define the following seven CVs:
(1) the fraction of native contacts, (2) the full backbone RMSD, (3-5) the
backbone RMSD for each of the three $\beta$-strands taken individually
(residues 1-4, 17-20, and 36-38), (6) the backbone RMSD for the $\alpha$-helix
(residues 22-29), and (7) the backbone RMSD for the three $\beta$-strands
together. Our basis set consists of 10 indicator functions on each stratum
constructed by $k$-means clustering on the seven-dimensional CV space. At each
iteration (i.e., one pass through the outer loop in Algorithm 4), we determine
the basis functions for stratum $j$ by taking the associated cluster centers
from the prior iteration and refining them with 10 iterations of Lloyd’s
algorithm [45] using all the samples with index process $J_{t}^{i}=j$ from the
last $h$ iterations.

Here we focus on the convergence of the folding rate; a detailed analysis of
the NTL9 folding mechanism (with and without the C-terminal helix) based on
potentials of mean force and committors will be presented elsewhere (see also
Ref. 46). As mentioned above, we compare our results to those of Ref. 34,
which estimates rates using haMSMs applied to data from WE (haMSM-WE). As
discussed in the Introduction, BAD-NEUS goes beyond haMSM-WE by using the
basis set (Markov model) to accelerate the convergence of the sampling rather
than just the rate estimates.

Figure 10 shows the inverse rate constant obtained from BAD-NEUS and haMSM-WE
[34]. The average inverse rate over the final 15 iterations of BAD NEUS is
4.74 ms, with a standard deviation of 0.99 ms, while the authors of Ref. 34
used a Bayesian bootstrapping approach to estimate a confidence interval of
0.17–1.9 ms, which they state is likely an underestimate of the true 95
percent confidence interval owing to the limited number of independent samples
used in the analysis. Given that we have only a single BAD NEUS run, the
standard deviation of the final 15 (highly correlated) inverse rate estimates
is also likely a significant underestimate of our statistical error. We thus
view the results in Figure 10 as in agreement.

We examine the weights in more detail in Figure 9. The increase in the inverse
rate around iteration 35 in Figure 10 (and concommitant decrease in below 5 in
Figure 9(bottom)) corresponds to a shift in the weights of strata with indices
close to 25 (compare iterations 30 and 45 in Figure 9(top)). These strata
contain walkers that are nearly denatured from trajectories that last visited
the native state, and the shift in weights in that region stabilizes the
denatured state relative to the native state.

Both BAD-NEUS and haMSM-WE required significantly less aggregate simulation to
compute the rate constant (an average) than the millisecond timescale of a
single folding event. The haMSM-WE calculation used an aggregate simulation
time of 252 $\mu$s, while we use an average of 4 $\mu$s per iteration, for a
total of about 204 $\mu$s. The speedup is actually more than it appears
because the haMSM-WE calculation is for the forward (folding) direction only;
computing potentials of mean force and committors would also require the
corresponding calculation for the backward (unfolding) direction. By contrast,
the BAD-NEUS simulation provides these statistics in its present form.

Figure 9:  Comparison of inverse rates for NTL9(1-39) folding from NEUS, BAD-
NEUS and haMSM-WE. The haMSM-WE results are from Ref. 34 The shaded area shows
the Bayesian $95\%$ credible interval reported in that publication.  Figure
10:  (top) BAD NEUS bin weights for the NTL9(1-39) system. Strata with indices
less than 20 correspond to the folding direction (trajectories last in the
denatured state), and strat with indices greater than or equal to 20
correspond to the unfolding direction (trajectories last in the native state).
(bottom) Deviation of weights ($\ln(\overline{z})$) from those in the last
iteration.

##  VII Conclusions

Trajectory stratification methods enable enhanced sampling for estimating both
thermodynamic (equilibrium) and kinetic (nonequilibrium) statistics. Here, we
introduced a new trajectory stratification method, BAD-NEUS, which converges
faster than existing ones. We show that our method is a natural generalization
of WE as originally formulated [1] but that, without modification, WE
converges no faster than unbiased dynamics. The key modification that we
introduce is the insertion of an approximation to the steady-state
distribution before the resampling step. This approximation algorithm is
designed so that it preserves the steady-state distribution of the dynamics,
and therefore does not introduce systematic errors into the algorithm in the
large data limit.

In the present study, we use a simple basis expansion to model the steady-
state distribution, but our strategy is general, and alternatives are also
possible. One can use neural networks to learn the steady-state distribution
and/or the basis functions [26, 27], and one can incorporate memory into a
basis expansion [47]. These alternatives, which can be used separately or
together, alleviate the need to identify basis sets that describe the dynamics
well and, in the case of memory, may reduce the sampling required. We thus
expect our method to be a significant be step toward accurate estimates of
rates for protein folding and similarly complex molecular conformational
changes.

###### Acknowledgements.

This work was supported by National Institutes of Health award R35 GM136381
and National Science Foundation award DMS-2054306. This work was completed in
part with resources provided by the University of Chicago Research Computing
Center. “Beagle-3: A Shared GPU Cluster for Biomolecular Sciences” is
supported by the National Institutes of Health (NIH) under the High-End
Instrumentation (HEI) grant program award 1S10OD028655-0.

##  Appendix A Derivation of (16)

Here we derive (16) by modifying the proof from Ref. 48. Let $Z_{t}$ be a
Markov process, $\mathcal{T}_{t}$ be operator defined in (40), and
$S_{n}=\sum_{t=0}^{n-1}g(Z_{t})$ be the partial sum of a function $g$. Let
$\omega(z)$ satisfy the linear problem

$(\mathcal{T}_{1}-\mathcal{I})\omega(z)=-(g(z)-\mathbbm{E}_{\pi}[g(z)])$ (46)

subject to the constraint $\mathbbm{E}_{\pi}[\omega]=0$.

Consider the process:

$U_{n}=S_{n}-\mathbbm{E}_{\pi}[g]n+\omega(Z_{n})$ (47)

Then $U_{t}$ is a martingale with respect to $Z_{0},...,Z_{t}$. To see this,
we compute:

$\displaystyle\mathbbm{E}[U_{t+1}-U_{t}|Z_{0},...Z_{t}]$
$\displaystyle=\mathbbm{E}[S_{t+1}-\operatorname{\mathbbm{E}}_{\pi}[g](t+1)+\omega(Z_{t+1})-\left(S_{t}-\mathbbm{E}_{\pi}[g]t+\omega(Z_{t})\right)|Z_{0},...Z_{t}]$
$\displaystyle=\mathbbm{E}[g(Z_{t})-\operatorname{\mathbbm{E}}_{\pi}[g]+\omega(Z_{t+1})-\omega(Z_{t})|Z_{0},...Z_{t}]$
$\displaystyle=g(Z_{t})-\operatorname{\mathbbm{E}}_{\pi}[g]+\mathbbm{E}[\omega(Z_{t+1})-\omega(Z_{t})|Z_{0},...Z_{t}]$
$\displaystyle=g(Z_{t})-\operatorname{\mathbbm{E}}_{\pi}[g]+(\mathcal{T}_{1}-\mathcal{I})\omega(Z_{t})$
$\displaystyle=g(Z_{t})-\operatorname{\mathbbm{E}}_{\pi}[g]-g(Z_{t})+\mathbbm{E}_{\pi}[g]$
$\displaystyle=0.$ (48)

We refer the reader to Ref. 48 for a more technical discussion and conditions
under which the relevant expectations are bounded such that we may apply the
optional stopping theorem. Then, for any stopping time $T$ with
$\operatorname{\mathbbm{E}}[T]<\infty$,

$0=\operatorname{\mathbbm{E}}[U_{T}-U_{0}]=\operatorname{\mathbbm{E}}[S_{T}]-\operatorname{\mathbbm{E}}_{\pi}[g]\operatorname{\mathbbm{E}}[T]+\operatorname{\mathbbm{E}}[\omega(Z_{T})]-\operatorname{\mathbbm{E}}[\omega(Z_{0})],$
(49)

which yields the result:

$\frac{\operatorname{\mathbbm{E}}[\sum_{t=0}^{T-1}g(Z_{t})]}{\operatorname{\mathbbm{E}}[T]}=\operatorname{\mathbbm{E}}_{\pi}[g]+\frac{\operatorname{\mathbbm{E}}[\omega(Z_{0})]-\operatorname{\mathbbm{E}}[\omega(Z_{T})]}{\operatorname{\mathbbm{E}}[T]}.$
(50)

Now take $Z_{t}=(X_{t},J_{t},...,X_{t+\tau},J_{t+\tau})$ and take the stopping
time to be $S(1)$. In the case where $X_{0},J_{0}$ is distributed according to
the steady state flux distribution $\overline{\pi}(dx,k),$ (9) implies that
$(X_{T},J_{T})\sim\overline{\pi}(dx,k)$, and hence the distribution of $Z_{0}$
and $Z_{T}$ are the same, and so the residual term in (50) is zero. In this
case, we note that the distribution of $(X_{0},J_{0})$ is given by:

$\displaystyle\pi(dx)$ $\displaystyle=\mathbbm{P}[X_{0}\in dx]$ (51)
$\displaystyle=\sum_{k}\mathbbm{P}[X_{0}\in dx\mid
J_{0}=k]\mathbbm{P}[J_{0}=k]$ (52)
$\displaystyle=\sum_{k}\overline{z}^{k}\overline{\pi}(dx|k).$ (53)

The joint distribution for $(X,J)$ is
$\overline{\pi}(dx,k)=\overline{z}^{k}\overline{\pi}(dx|k)$. Therefore,

$\mathbbm{E}\left[\sum_{t=0}^{S(1)-1}g(X_{t},J_{t},...,X_{t+\tau},J_{t+\tau})\right]=\sum_{k}\overline{z}^{k}\int\overline{\pi}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}\left[\sum_{t=0}^{S(1)-1}g(X_{t},J_{t},...,X_{t+\tau},J_{t+\tau})\right]$
(54)

and

$\mathbbm{E}[S(1)]=\sum_{k}\overline{z}^{k}\int\overline{\pi}(dx|k)\mathbbm{E}_{X_{0}=x,J_{0}=k}[S(1)].$
(55)

Substituting these into (50) and remembering that the residual term is zero
yields (16).

## References

  * [1] G. A. Huber and S. Kim.  Weighted-ensemble Brownian dynamics simulations for protein association reactions.  Biophysical Journal, 70(1):97–110, January 1996.
  * [2] Rosalind J. Allen, Chantal Valeriani, and Pieter Rein ten Wolde.  Forward flux sampling for rare event simulations.  Journal of Physics: Condensed Matter, 21(46):463102, October 2009\.
  * [3] Aryeh Warmflash, Prabhakar Bhimalapuram, and Aaron R. Dinner.  Umbrella sampling for nonequilibrium processes.  Journal of Chemical Physics, 127(15):154112, October 2007.
  * [4] Alex Dickson and Aaron R Dinner.  Enhanced sampling of nonequilibrium steady states.  Annual review of physical chemistry, 61:441–459, 2010.
  * [5] Nicholas Guttenberg, Aaron R Dinner, and Jonathan Weare.  Steered transition path sampling.  The Journal of Chemical Physics, 136(23), 2012.
  * [6] Juan M. Bello-Rivas and Ron Elber.  Exact milestoning.  Journal of Chemical Physics, 142(9):094102, March 2015.
  * [7] Daniel M Zuckerman and Lillian T Chong.  Weighted ensemble simulation: review of methodology, applications, and software.  Annual review of biophysics, 46:43–57, 2017.
  * [8] John D. Russo, She Zhang, Jeremy M. G. Leung, Anthony T. Bogetti, Jeff P. Thompson, Alex J. DeGrave, Paul A. Torrillo, A. J. Pratt, Kim F. Wong, Junchao Xia, Jeremy Copperman, Joshua L. Adelman, Matthew C. Zwier, David N. LeBard, Daniel M. Zuckerman, and Lillian T. Chong.  WESTPA 2.0: High-Performance Upgrades for Weighted Ensemble Simulations and Analysis of Longer-Timescale Applications.  Journal of Chemical Theory and Computation, 18(2):638–649, February 2022.  Publisher: American Chemical Society.
  * [9] D. Aristoff, J. Copperman, G. Simpson, R. J. Webber, and D. M. Zuckerman.  Weighted ensemble: Recent mathematical developments.  Journal of Chemical Physics, 158(1):014108, January 2023.
  * [10] Surl-Hee Ahn, Anupam A. Ojha, Rommie E. Amaro, and J. Andrew McCammon.  Gaussian-Accelerated Molecular Dynamics with the Weighted Ensemble Method: A Hybrid Method Improves Thermodynamic and Kinetic Sampling.  Journal of Chemical Theory and Computation, 17(12):7938–7951, December 2021.  Publisher: American Chemical Society.
  * [11] Anupam Anand Ojha, Saumya Thakur, Surl-Hee Ahn, and Rommie E. Amaro.  DeepWEST: Deep Learning of Kinetic Models with the Weighted Ensemble Simulation Toolkit for Enhanced Sampling.  Journal of Chemical Theory and Computation, 19(4):1342–1359, February 2023.  Publisher: American Chemical Society.
  * [12] Jeremy Copperman and Daniel M. Zuckerman.  Accelerated Estimation of Long-Timescale Kinetics from Weighted Ensemble Simulation via Non-Markovian “Microbin” Analysis.  Journal of Chemical Theory and Computation, 16(11):6763–6775, November 2020.  Publisher: American Chemical Society.
  * [13] Alex J. DeGrave, Anthony T. Bogetti, and Lillian T. Chong.  The RED scheme: Rate-constant estimation from pre-steady state weighted ensemble simulations.  The Journal of Chemical Physics, 154(11):114111, March 2021.
  * [14] Ernesto Suárez, Steven Lettieri, Metthew C Zwier, Sundar Raman Subramanian, Lillian T Chong, and Daniel M Zuckerman.  Simultaneous computation of dynamical and equilibrium information using a weighted ensemble of trajectories.  Biophysical Journal, 106(2):406a, 2014.
  * [15] Alex Dickson, Aryeh Warmflash, and Aaron R Dinner.  Nonequilibrium umbrella sampling in spaces of many order parameters.  Journal of chemical physics, 130(7), 2009.
  * [16] Alex Dickson, Aryeh Warmflash, and Aaron R Dinner.  Separating forward and backward pathways in nonequilibrium umbrella sampling.  Journal of chemical physics, 131(15), 2009.
  * [17] Alex Dickson, Mark Maienschein-Cline, Allison Tovo-Dwyer, Jeff R Hammond, and Aaron R Dinner.  Flow-dependent unfolding and refolding of an rna by nonequilibrium umbrella sampling.  Journal of Chemical Theory and Computation, 7(9):2710–2720, 2011\.
  * [18] Aaron R. Dinner, Jonathan C. Mattingly, Jeremy O. B. Tempkin, Brian Van Koten, and Jonathan Weare.  Trajectory Stratification of Stochastic Dynamics.  SIAM Review, 60(4):909–938, January 2018.  Publisher: Society for Industrial and Applied Mathematics.
  * [19] Divesh Bhatt, Bin W. Zhang, and Daniel M. Zuckerman.  Steady-state simulations using weighted ensemble path sampling.  Journal of Chemical Physics, 133(1):014110, July 2010.
  * [20] Gabriel Earle and Jonathan Mattingly.  Convergence of Stratified MCMC Sampling of Non-Reversible Dynamics, February 2022.  arXiv:2111.05838 [math].
  * [21] Gabriel Earle and Brian Van Koten.  Aggregation methods for computing steady states in statistical physics.  Multiscale Modeling & Simulation, 21(3):1170–1209, 2023.
  * [22] Erik H. Thiede, Dimitrios Giannakis, Aaron R. Dinner, and Jonathan Weare.  Galerkin Approximation of Dynamical Quantities using Trajectory Data.  Journal of Chemical Physics, 150(24):244111, 2019.
  * [23] John Strahan, Adam Antoszewski, Chatipat Lorpaiboon, Bodhi P. Vani, Jonathan Weare, and Aaron 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(5):2948–2963, 2021\.
  * [24] Justin Finkel, Robert J. Webber, Dorian S. Abbot, Edwin P. Gerber, and Jonathan Weare.  Learning forecasts of rare stratospheric transitions from short simulations.  Monthly Weather Review, 149(11):3647–3669, 2021.
  * [25] Bodhi P. Vani, Jonathan Weare, and Aaron R. Dinner.  Computing transition path theory quantities with trajectory stratification.  Journal of Chemical Physics, 157(3):034106, July 2022.
  * [26] Junfeng Wen, Bo Dai, Lihong Li, and Dale Schuurmans.  Batch stationary distribution estimation.  In Proceedings of the 37th International Conference on Machine Learning, ICML’20, pages 10203–10213. JMLR.org, 2020.
  * [27] John Strahan, Spencer C. Guo, Chatipat Lorpaiboon, Aaron R. Dinner, and Jonathan Weare.  Inexact iterative numerical linear algebra for neural network-based spectral estimation and rare-event prediction.  Journal of Chemical Physics, 2023.
  * [28] Randal Douc, Olivier Cappé, and Eric Moulines.  Comparison of Resampling Schemes for Particle Filtering.  ISSN: 1845-5921, pages 64–69, 2005.
  * [29] John Strahan, Justin Finkel, Aaron R. Dinner, and Jonathan Weare.  Predicting rare events using neural networks and short-trajectory data, March 2023.  arXiv:2208.01717 [physics].
  * [30] Klaus Müller and Leo D. Brown.  Location of saddle points and minimum energy paths by a constrained simplex optimization procedure.  Theoretica chimica acta, 53(1):75–93, 1979.
  * [31] Benedict Leimkuhler and Charles Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Applied Mathematics Research eXpress, 2013(1):34–56, 2013.
  * [32] Jérôme Hénin, Tony Lelièvre, Michael Shirts, Omar Valsson, and Lucie Delemotte.  Enhanced sampling methods for molecular dynamics simulations [article v1. 0].  Living Journal of Computational Molecular Science, 4(1):1583–1583, 2022.
  * [33] Eric Vanden-Eijnden and Maddalena Venturoli.  Exact rate calculations by trajectory parallelization and tilting.  The Journal of chemical physics, 131(4), 2009.
  * [34] Upendra Adhikari, Barmak Mostofian, Jeremy Copperman, Sundar Raman Subramanian, Andrew A. Petersen, and Daniel M. Zuckerman.  Computational Estimation of Microsecond to Second Atomistic Folding Times.  Journal of the American Chemical Society, 141(16):6519–6526, April 2019.  Publisher: American Chemical Society.
  * [35] Donna L Luisi, Brian Kuhlman, Kostandinos Sideras, Philip A Evans, and Daniel P Raleigh.  Effects of varying the local propensity to form secondary structure on the stability and folding kinetics of a rapid folding mixed $\alpha/\beta$ protein: characterization of a truncation mutant of the N-terminal domain of the ribosomal protein L911Edited by P. E. Wright.  Journal of Molecular Biology, 289(1):167–174, May 1999.
  * [36] Satoshi Sato, Jae-Hyun Cho, Ivan Peran, Rengin G. Soydaner-Azeloglu, and Daniel P. Raleigh.  The N-Terminal Domain of Ribosomal Protein L9 Folds via a Diffuse and Delocalized Transition State.  Biophysical Journal, 112(9):1797–1806, May 2017.
  * [37] Vincent A Voelz, Gregory R Bowman, Kyle Beauchamp, and Vijay S Pande.  Molecular simulation of ab initio protein folding for a millisecond folder ntl9 (1- 39).  Journal of the American Chemical Society, 132(5):1526–1528, 2010\.
  * [38] Peter Eastman, Jason Swails, John D. Chodera, Robert T. McGibbon, Yutong Zhao, Kyle A. Beauchamp, Lee-Ping Wang, Andrew C. Simmonett, Matthew P. Harrigan, Chaya D. Stern, Rafal P. Wiewiora, Bernard R. Brooks, and Vijay S. Pande.  OpenMM 7: Rapid development of high performance algorithms for molecular dynamics.  PLOS Computational Biology, 13(7):e1005659, July 2017.  Publisher: Public Library of Science.
  * [39] Robert T. McGibbon, Kyle A. Beauchamp, Matthew P. Harrigan, Christoph Klein, Jason M. Swails, Carlos X. Hernández, Christian R. Schwantes, Lee-Ping Wang, Thomas J. Lane, and Vijay S. Pande.  Mdtraj: A modern open library for the analysis of molecular dynamics trajectories.  Biophysical Journal, 109(8):1528 – 1532, 2015.
  * [40] James A. Maier, Carmenza Martinez, Koushik Kasavajhala, Lauren Wickstrom, Kevin E. Hauser, and Carlos Simmerling.  ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB.  Journal of Chemical Theory and Computation, 11(8):3696–3713, August 2015.  Publisher: American Chemical Society.
  * [41] Gregory D. Hawkins, Christopher J. Cramer, and Donald G. Truhlar.  Parametrized Models of Aqueous Free Energies of Solvation Based on Pairwise Descreening of Solute Atomic Charges from a Dielectric Medium.  Journal of Physical Chemistry, 100(51):19824–19839, January 1996\.  Publisher: American Chemical Society.
  * [42] Gregory D. Hawkins, Christopher J. Cramer, and Donald G. Truhlar.  Pairwise solute descreening of solute charges from a dielectric medium.  Chemical Physics Letters, 246(1):122–129, November 1995.
  * [43] Robert B. Best, Gerhard Hummer, and William A. Eaton.  Native contacts determine protein folding mechanisms in atomistic simulations.  Proceedings of the National Academy of Sciences, 110(44):17874–17879, October 2013.  Publisher: Proceedings of the National Academy of Sciences.
  * [44] Jae-Hyun Cho, Wenli Meng, Satoshi Sato, Eun Young Kim, Hermann Schindelin, and Daniel P. Raleigh.  Energetically significant networks of coupled interactions within an unfolded protein.  Proceedings of the National Academy of Sciences, 111(33):12079–12084, August 2014.  Publisher: Proceedings of the National Academy of Sciences.
  * [45] S. Lloyd.  Least squares quantization in PCM.  IEEE Transactions on Information Theory, 28(2):129–137, March 1982\.  Conference Name: IEEE Transactions on Information Theory.
  * [46] John Champlin Strahan.  Short trajectory methods for rare event analysis and sampling.  2024\.
  * [47] Chatipat Lorpaiboon, Spencer C Guo, John Strahan, Jonathan Weare, and Aaron R Dinner.  Accurate estimates of dynamical statistics using memory.  The Journal of Chemical Physics, 160(8), 2024.
  * [48] George V. Moustakides.  Extension of Wald’s First Lemma to Markov Processes.  Journal of Applied Probability, 36(1):48–59, 1999.  Publisher: Applied Probability Trust.


remarkRemark hypothesisHypothesis claimClaim Stability estimatesF. Triki and
M. Karamehmedović and K. Linder-Steinlein

# Stability estimates for the inverse source problem with passive
measurements††thanks: Submitted to the editors . This work was funded by the
by The Villum Foundation (grant no. 25893).

Kristoffer Linder-Steinlein Department of Applied Mathematics, Technical
University of Denmark, Kgs. Lyngby, Denmark (,
https://orbit.dtu.dk/en/persons/kristoffer-linder-steinlein-2<EMAIL_ADDRESS>Mirza Karamehmedović Department of Applied Mathematics, Technical University
of Denmark, Kgs. Lyngby, Denmark (, https://orbit.dtu.dk/en/persons/mirza-
karamehmedovic<EMAIL_ADDRESS>Faouzi Triki  Laboratoire Jean Kuntzmann,
Université Grenoble-Alpes, Grenoble, France (, https://membres-
ljk.imag.fr/Faouzi.Triki/contact.html<EMAIL_ADDRESS>
###### Abstract

We consider the multi-frequency inverse source problem in the presence of a
non-homogeneous medium using passive measurements. Precisely, we derive
stability estimates for determining the source from the knowledge of only the
imaginary part of the radiated field on the boundary for multiple frequencies.
The proof combines a spectral decomposition with a quantification of the
unique continuation of the resolvent as a holomorphic function of the
frequency. The obtained results show that the inverse problem is well posed
when the frequency band is larger than the spatial frequency of the source.

###### keywords:

Passive-imaging, Stability estimate, Helmholtz equation, Inverse source
problem.

34A55

##  1 Introduction

This paper is concerned with the stability estimate for determining a one-
dimensional source in the presence of a medium using passive measurements.
That is, reconstruction of the source having access to only the imaginary part
of the resulting wave-field at the domain’s boundary. This is a type of
analysis of the time-reversal experiment [11]. Time reversal has been well
studied and is one of the most commonly used reconstruction methods in direct
imaging for inverse source problems. It was first proposed for energy focusing
in physics. Here we apply the Helmholtz-Kirchhoff identity to study time
reversal [2, 11].

Multiple results already exist in the cases where the entire wave-field at the
boundary is obtainable [5, 8, 4, 13, 9, 6, 1]. We here prove for the one-
dimensional Helmholtz equation that it is possible to complete the missing
real part of the resulting wave-field using the measurable quantities, and to
recover the source. These results can be seen in Theorem 1.1 for a large
frequency band and in Theorem 1.2 for a short frequency band.

The remaining part of this section introduces the mathematical model and
spaces from which the sources and media originate. This is followed by a
statement regarding the available data and the considered inverse problem with
passive measurements. In section 2 the proof of the first main result is
provided. It is based on a spectral decomposition of the source in a specific
orthonormal basis of eigenfunctions.

In the last section, the second main stability estimate for the inverse
problem is established using a quantification of the unique continuation for
the resolvent of the Helmholtz operator as a holomorphic function of the
frequency and making use of the results of the previous section 2.

###  1.1 Mathematical model

We here focus on the one dimensional Helmholtz equation, which can be
expressed as:

(1) $\phi^{\prime\prime}(x,k)+k^{2}(1+q(x))\phi(x,k)=f(x),\qquad
x\in\mathbb{R},$

where $1+q$ and $f$ are real-valued functions called respectively the
refractive index and the source; the resulting solution $\phi$ to (1) is the
field generated by the source $f$ in the presence of the medium $q$ both
supported in the interval $[0,1]$. The coefficient $k$ is any positive number
and referred to as the frequency. We furthermore impose the Sommerfeld
radiation condition on the field:

(2) $\begin{split}&\phi^{\prime}(0,k)+ik\phi(0,k)=0,\\\
&\phi^{\prime}(1,k)-ik\phi(1,k)=0.\end{split}$

The medium and source functions are assumed to belong to the spaces of real-
valued functions

$\mathcal{M}(q_{0},M,m,n_{0}):=\\{q\in
C^{m+1}_{0}\left([0,1]\right):\>\|q-q_{0}\|_{C^{m+1}([0,1])}\leq M,n_{0}\leq
1+q\\},$

where $m\geq 1$, $M>0$, $q_{0}\in C^{m+1}_{0}\left([0,1]\right)$ and
$1+q_{0}\geq n_{0}$ for some $n_{0}\in(0,1)$ and

$\mathcal{F}(L):=\\{f\in H^{1}(0,1):\>\|f\|_{H^{1}(0,1)}\leq
L,\>\operatorname{supp}f\subset(0,1)\\},$

for $L>0$.

###  1.2 Passive measurements model and main results

The data assumed accessible in this work is passive measurements. We consider
conventional full apparatus-based passive imaging, that is, the sensors
completely surround the domain imaged. The described setup is an application
of the analysis of the time-reversal experiment [11]. Time reversal has been
well studied and is one of the most commonly used reconstruction methods in
direct imaging for inverse source problems. It was first proposed for energy
focusing in physics. Here we apply the Helmholtz-Kirchhoff identity to study
time reversal [2, 11].

The analysis is based on an integral representation of the solution of the
Helmholtz equation (1), given by

(3) $\phi(x,k)=\int_{0}^{1}{G(x,z,k)}f(z)dz,$

where $G$ denotes the Green function of the Helmholtz equation (1), satisfying

(4)
$\begin{cases}G^{\prime\prime}(x,z,k)+k^{2}(1+q(x))G(x,z,k)=\delta_{z}(x),\qquad
x\in\mathbb{R},\\\ G^{\prime}(0,z,k)+ikG(0,z,k)=0,\\\
G^{\prime}(1,z,k)-ikG(1,z,k)=0.\end{cases}$

The time-reversal imaging functional in the multidimensional case given by [2]

$\mathcal{I}(x,k)=\int_{\Gamma}\overline{G(x,z,k)}\phi(z,k)ds(z).$

where $\Gamma$ is a closed curve where the measurements will be taken. We here
define the integral along the end-points of an interval as the Lebesgue
integral w.r.t. the signed measure $\gamma(E)=\mathbbm{1}_{b\in
E}-\mathbbm{1}_{a\in E}$ leading to

$\mathcal{I}(x,k)=\int_{\\{0,1\\}}\overline{G(x,z,k)}\phi(z,k)d\gamma(z).$

Based on [2, Corollary 2.1] or [11, Theorem 2.2] the above imaging functional
is equivalent to

$\mathcal{I}(x,k)=-\frac{1}{k}\int_{0}^{1}\Im G(x,y,k)f(y)dy,$

where the one-dimensional version of the Helmholtz-Kirchhoff identity is used
reading

$k\int_{\\{0,1\\}}\overline{G(x,z,k)}G(y,z,k)d\gamma(z)=-\Im G(x,y,k).$

Therefore

$\mathcal{I}(x,k)=-\frac{1}{k}\Im\phi(x,k),\;x\in\\{0,1\\}.$

In [11] it is proven that the Helmholtz-Kirchhoff identity is closely related
to the cross-correlation. The Helmholtz-Kirchhoff identity not only holds for
the Green’s functions but for the total field as well, and it is shown in
[10], in higher dimensions, the cross-correlation matrix is given by the
imaginary near field, however, in one dimension this amounts to the cross-
correlation being given by the imaginary field at the boundary.

Let $K$ be a fixed positive constant, and let $I=(0,K)$ be the index set for
the wave-numbers used for multifrequency measurements. We will only work with
measurements at the boundary $x\in\\{0,1\\}$, and the inverse problem can be
stated as:

_From measurements $\Im\phi(0,k)$ and $\Im\phi(1,k),\;k\in I$, reconstruct the
source $f$._

This is solved in two steps, first reconstructing the source $f$ using the
system satisfied by $\Re\phi$ on a large band of frequencies, and second
employing unique continuation techniques for holomorphic functions to recover
the required boundary data from a small band of frequencies. We give here the
main results of the paper.

######  Theorem 1.1.

Let $f\in\mathcal{F}$. Then there exist constants $c_{0}>0$ and $C_{0}>0$ that
only depend on $\mathcal{M}$ such that

(5) $\|f\|_{H^{-1}}\leq
C_{0}(\mu_{f}+1)\left(\|\Im\phi(0,\cdot)\|_{L^{\infty}(0,c_{0}\mu_{f})}+\|\Im\phi(1,\cdot)\|_{L^{\infty}(0,c_{0}\mu_{f})}\right),$

where $\mu_{f}=\|f\|_{L^{2}}/\|f\|_{H^{-1}}$.

######  Theorem 1.2.

Assume that

$\varepsilon:=C_{0}(\mu_{f}+1)\left(\|\Im\phi(0,\cdot)\|_{L^{\infty}(0,K)}+|\Im\phi(1,\cdot)\|_{L^{\infty}(0,K)}\right)<1.$

Then there exist constants $C=C(\mathcal{M},L)>0$ and
$n_{h}=n_{h}(\mathcal{M})\in\mathbb{N}^{*}$ such that the inequality

(6) $\|f\|_{H^{-1}}\leq
C^{1-\eta(c_{0}\mu_{f},K)}\varepsilon^{\eta(c_{0}\mu_{f},K)},$

holds, where $\mu_{f}=\|f\|_{L^{2}}/\|f\|_{H^{-1}}$, and

$\eta(s,K):=\frac{2}{\pi}\arctan\left(\frac{(e^{K}-1)^{n_{h}}}{\left((e^{s}-1)^{{2n_{h}}}-(e^{K}-1)^{{2n_{h}}}\right)^{\frac{1}{2}}}\right),\quad\textrm{for
}s>K.$

######  Remark 1.3.

i) The number $\mu_{f}$ represents the frequency of the source and
characterizes its spatial oscillations [3, 16].

ii) The stability estimate (5) shows that the inverse problem with passive
measurements is well posed if the frequency band is large enough and covers
the spatial frequency of the source. The Lipschitz constant grows linearly
with respect to the frequency of the source which is in agreement with the
known physical resolution in source imaging.

iii) The stability estimate (6) indicates that the inverse problem with
passive measurements becomes ill-posed if the frequency band shrinks to zero.
Notice that when $K=c_{0}\mu_{f}$ we recover the stability estimate (5) from
(6).

iv) The integer $n_{h}$ is related to $h=\frac{2\pi}{n_{h}}$ which is the
width of the complex strip $S$ around the real axis (see Proposition 3.1)
where the system (1) is free from scattering resonances. It can be shown that
$\eta(c_{0}\mu_{f},K)$ decreases exponentially to zero when $h$ tends to zero.
This shows that recovering the source from a small band of frequencies becomes
ill-posed when the imaginary part of the resonances are closer to the real
axis (trapped modes).

##  2 Proof of Theorem 1.1

Let $w=\Re(\phi)$. We deduce from (1) and (2) that $w$ satisfies the following
system

(7)
$\displaystyle\begin{cases}w^{\prime\prime}(x,k)+k^{2}(1+q(x))w(x,k)=f(x),\qquad
x\in\mathbb{R},\\\ w^{\prime}(0,k)=k\Im\phi(0,k),\\\
w^{\prime}(1,k)=-k\Im\phi(1,k).\end{cases}$

Based on the passive measurements the data at hand is equivalent to the
Neumann boundary data in the system (7). Recall that this later may not have a
unique solution for all $k\in I$.

The strategy next is to recover $f(x),\,x\in(0,1)$ from $\Im(\phi(x,k)),\,k\in
I,\,x\in\\{0,1\\}$ using the system (7). The approach is based on a spectral
decomposition of the source in an orthonormal basis formed by the
eigenfunctions of the system. Taking the frequency within the set of the
associated real eigenvalues we are able to determine the coefficient of the
expansion of the source in terms of passive boundary measurements. This method
has been applied previously in multifrequency inverse source problems with
full data [1, 13].

###  2.1 Spectral decomposition

Since $q\in\mathcal{M}$, the weighted space $L^{2}_{q}(0,1)$ endowed with the
inner product

$\langle\varphi,\psi\rangle_{L^{2}_{q}}=\int_{0}^{1}(1+q(x))\varphi(x)\psi(x)dx,$

is well defined, and its norm $\|\cdot\|_{L^{2}_{q}}$ is equivalent to the
classical $L^{2}$ norm.

Let $\\{\mu_{j},\phi_{j}\\}_{j=1}^{\infty}$ be all the pairs solutions to the
Neumann spectral problem:

(8)
$\begin{cases}-\phi^{\prime\prime}_{j}(x)=\mu_{j}^{2}(1+q(x))\phi_{j}(x)&x\in(0,1),\\\
\phi_{j}^{\prime}=0&x\in\\{0,1\\},\\\ \|\phi_{j}\|_{L^{2}_{q}}=1.\end{cases}$

Since $q\in\mathcal{M}$, we have $\mu_{1}=0$, $(\mu_{j})_{j\in\mathbb{N}^{*}}$
is an increasing real sequence, $\phi_{1}=\|1+q\|_{L^{1}}^{-1/2}$, and
$(\phi_{j})_{j\in\mathbb{N}^{*}}$ is an orthonormal basis of $L^{2}_{q}(0,1)$
[14].

Therefore

$f(x)(1+q(x))^{-1}=\sum_{j=1}^{\infty}f_{j}\phi_{j}(x),\quad x\in(0,1),$

where $f_{j}=\langle
f(1+q),\phi_{j}\rangle_{L^{2}_{q}}=\int_{0}^{1}f(x){\phi}_{j}dx$.

Multiplying the Helmholtz equation (7) taken at $k=\mu_{j}$ by $\phi_{j}$, and
integrating by parts, lead to

(9)
$f_{j}=\mu_{j}\left(\phi_{j}(0)\Im(\phi)(0,\mu_{j})-\phi_{j}(1)\Im(\phi)(1,\mu_{j})\right).$

######  Proposition 2.1.

Let $(\mu_{j},\phi_{j})$ be a pair of eigenelements of the Neumann spectral
problem (8). Then there exists a constant $C=C(\mathcal{M})>0$ such that

(10) $|\phi_{j}(0)|+|\phi_{j}(1)|\leq C(\mu_{j}+1),\quad\forall
j\in\mathbb{N}^{*}.$

######  Proof 2.2.

We have

$\phi_{j}(0)=\phi_{j}(x)-\int_{0}^{x}\phi_{j}^{\prime}(t)dt,\quad\forall
x\in(0,1).$

Hence

$|\phi_{j}(0)|\leq|\phi_{j}(x)|+\|\phi_{j}^{\prime}\|_{L^{2}}.$

Integrating both sides over $(0,1)$, yields

$|\phi_{j}(0)|\leq\|(1+q)^{-\frac{1}{2}}\|_{L^{2}}+\|\phi_{j}^{\prime}\|_{L^{2}}.$

On the other hand, we deduce from (8) that
$\|\phi_{j}^{\prime}\|_{L^{2}}=\mu_{j}$. By taking
$C=\max(1,\|(1+q)^{-\frac{1}{2}}\|_{L^{2}})$, we then obtain the final
estimate for $\phi_{j}(0)$. The estimate for $\phi_{j}(1)$ can be derived by
following the same steps of the previous proof for $\phi_{j}(0)$.

Combining the identity (9) with the estimate (10), we get

(11) $|f_{j}|\leq
C(1+\mu_{j})\left(|\Im(\phi)(0,\mu_{j})|+|\Im(\phi)(1,\mu_{j})|\right),\quad\forall
j\in\mathbb{N}^{*}.$

For $\psi\in H^{s}(0,1)$ with $s\in\mathbb{R}$, we further define the norm

$\|\psi\|_{s}:=\left(\sum_{j=1}^{\infty}(1+\mu_{j}^{2})^{s}\psi_{j}^{2}\right)^{\frac{1}{2}},$

where $\psi_{j}=\langle\psi,\phi_{j}\rangle_{L^{2}_{q}}.$ For
$\mu\in(0,+\infty)$, we then have

$\|f(1+q)^{-1}\|_{-1}^{2}=\sum_{j=1}^{\infty}(1+\mu_{j}^{2})^{-1}f_{j}^{2}\leq\sum_{\mu_{j}\leq\mu}(1+\mu_{j}^{2})^{-1}f_{j}^{2}+\frac{1}{\mu^{2}}\sum_{\mu_{j}>\mu}f_{j}^{2}.$

Consequently

$(1-\frac{\tilde{\mu}_{f}^{2}}{\mu^{2}})\|f(1+q)\|_{-1}^{2}\leq\sum_{\mu_{j}\leq\mu}(1+\mu_{j}^{2})^{-1}f_{j}^{2},$

where

$\tilde{\mu}_{f}=\frac{\|f(1+q)^{-1}\|_{0}}{\|f(1+q)^{-1}\|_{-1}}.$

Taking $\mu=\sqrt{2}\tilde{\mu}_{f}$, we obtain

(12) $\|f(1+q)^{-1}\|_{-1}^{2}\leq
2\sum_{\mu_{j}\leq\sqrt{2}\tilde{\mu}_{f}}(1+\mu_{j}^{2})^{-1}f_{j}^{2}.$

On the other hand we deduce from estimate (9)

(13) $|f_{j}|^{2}\leq
4C(1+\mu_{j}^{2})\left(|\Im(\phi)(0,\mu_{j})|^{2}+|\Im(\phi)(1,\mu_{j})|^{2}\right),\quad\forall
j\in\mathbb{N}^{*}.$

Combining inequalities (12) and (13) yields

(14) $\|f(1+q)^{-1}\|_{-1}^{2}\leq
8C\left(\sum_{\mu_{j}\leq\sqrt{2}\tilde{\mu}_{f}}1\right)\sup_{k\in(0,\sqrt{2}\tilde{\mu}_{f})}\left(|\Im(\phi)(0,k)|^{2}+|\Im(\phi)(1,k)|^{2}\right).$

######  Proposition 2.3.

Let $\mu_{j}$ be an eigenvalue of the Neumann spectral problem (8). Then there
exist constants $c_{i}=c_{i}(\mathcal{M})>0,\,i=1,2,$ such that

(15) $c_{1}(j-1)\leq\mu_{j}\leq c_{2}(j-1),\quad\forall j\in\mathbb{N}^{*}.$

######  Proof 2.4.

Using the Min-max principle we have [12]

$\mu_{j}=\min_{E_{j}\subset H^{1}(0,1),\,\textrm{dim}(E_{j})=j}\max_{\phi\in
E_{j}\setminus\\{0\\}}\frac{\|\phi^{\prime}\|_{L^{2}}^{2}}{\|\phi\|^{2}_{L^{2}_{q}}}.$

Since

$\min_{x\in[0,1]}(1+q(x))\|\phi\|_{L^{2}}\leq\|\phi\|_{L^{2}_{q}}\leq\max_{x\in[0,1]}(1+q(x))\|\phi\|_{L^{2}},$

we obtain

$\left(\max_{x\in[0,1]}(1+q(x))\right)^{-1}\pi(j-1)\leq\mu_{j}\leq\left(\min_{x\in[0,1]}(1+q(x))\right)^{-1}\pi(j-1).$

Consequently there exist constants $c_{i}>0,\,i=1,2,$ that depend only on
$\mathcal{M}$ such that the inequalities (15) hold.

We deduce from estimates (15) that

$\sum_{\mu_{j}\leq\sqrt{2}\tilde{\mu}_{f}}1\leq c_{3}\tilde{\mu}_{f}+1,$

for some constant $c_{3}>0$ that depend only on $\mathcal{M}$. Then it follows
from inequality (14) that

(16) $\|f(1+q)^{-1}\|_{-1}^{2}\leq
8C(c_{3}\tilde{\mu}_{f}+1)\sup_{k\in(0,\sqrt{2}\tilde{\mu}_{f})}\left(|\Im(\phi)(0,k)|^{2}+|\Im(\phi)(1,k)|^{2}\right).$

Since $q\in\mathcal{M}$, one can easily show that for $s\in\mathbb{R}$, the
norm $\|\cdot\|_{H^{s}}$ is equivalent to $\|\cdot\|_{s}$. Hence there exists
a constant $c_{4}>0$ that only depend on $\mathcal{M}$ such that
$\tilde{\mu}_{f}\leq c_{4}\mu_{f}$. Finally the main result of the Theorem
follows directly from the fact that $q\in\mathcal{M}$ and the last inequality
(16).

##  3 Proof of Theorem 1.2

In this section, we adopt the method used in [13] to the problem considered
here. The significant differences hinge on the fact that no Black-Box operator
theory is used, and the region where the resolvent is holomorphic is given by
[7][Proposition 2.3], which follows from Gel’fand-Levitan techniques that
convert the Helmholtz equation into a Schrödinger equation. Indeed in the
obtained Schrödinger equation, the refractive index and the frequency are
separated, which allows a better understanding of the behavior of the
solutions as functions of the frequency. For the sake of completeness,
[7][Proposition 2.3] is restated:

######  Proposition 3.1.

There exists a constant $h=h(\mathcal{M})>0$ such that the strip

$S=\\{k\in\mathbb{C};\>-h\leq\Im(k)\leq h\\},$

is free from resonances of the system (1).

We deduce from Proposition 3.1 that the Green function defined in (4) has no
poles in the strip $S$. Therefore $\phi(x,\cdot)$ which satisfies the integral
representation (3) is holomorphic and bounded in $S$. Let

(17) $M_{f}=\max_{k\in S}\left(|\phi(0,\cdot)|+|\phi(1,\cdot)|\right).$

Notice that $M_{f}>0$ depends only on $\mathcal{M}$ and $L$. We also remark
that $\overline{\phi(\cdot,k)}$ is a solution to the system (1) with radiation
conditions (2) when substituting $k$ by $-k$. We deduce from the uniqueness of
the system that $\overline{\phi(\cdot,k)}=\phi(\cdot,-k)$. Therefore for fixed
$x\in\mathbb{R}$,
$k\rightarrow\Im(\phi(x,k))=\frac{1}{2i}\left(\phi(\cdot,k)-\phi(x,-k)\right)$,
is also a holomorphic function in the strip $S$.

Consequently

(18) $F(z):=\Im(\phi(0,z))+i\Im(\phi(1,z)),$

is holomorphic in $S$. In addition it satisfies $F(-z)=F(z)$ for
$z\in\mathbb{R}$, and $|F(z)|\leq 2M_{f},\;\forall z\in S$.

Next, we aim to estimate $F(z)$ within the strip $S$ in terms of its values on
a segment $F(z),z\in(-K,K)$.

Without loss of generality we can assume that $h=\frac{\pi}{2n_{h}},$ where
$n_{h}\in\mathbb{N}^{*}$. Let
$S_{+}=\\{k\in\mathbb{C};\Re(k)>0,\,|\Im(k)|<h\\}$, be half a strip, and
$w(k,K)$ be the harmonic measure of the complex open domain
$S_{+}\setminus[0,K]\times\\{0\\}$. It is the unique solution to the system:

(22) $\displaystyle\left\\{\begin{array}[]{lllccc}\Delta w(k,K)&=&0\quad k\in
S_{+}\setminus[0,K]\times\\{0\\},\\\ w(k,K)&=&0\quad k\in\partial S_{+},\\\
w(k,K)&=&1\quad k\in(0,K]\times\\{0\\}.\end{array}\right.$

We infer from the maximum principle that $w_{0}(k,K)\in[0,1]$ for all $k\in
S_{+}$. Moreover the holomorphic unique continuation of the functions $F$
using the Two constants Theorem ([15, Chap. III, Section 2.1], [17]), gives:

######  Lemma 3.2.

Let $F$ be defined by (18), let $M_{f}$ be given by (17) and let $w$ be the
same as in (22). Then, we have

$\displaystyle|F(k)|\leq(2M_{f})^{1-w_{0}(k,K)}\|F\|_{L^{\infty}(0,K)}^{w_{0}(k,K)},\qquad\forall
k\in(K,+\infty).$

The following estimate is needed [7, Proposition 5.1].

######  Proposition 3.3.

Let $w_{0}(k,K)$ be the harmonic measure of
$S_{+}\setminus(0,K]\times\\{0\\}$. Then

(23) $\displaystyle\
w_{0}(k,K)\geq\frac{2}{\pi}\arctan(\frac{(e^{K}-1)^{n_{h}}}{\left((e^{k}-1)^{{2n_{h}}}-(e^{K}-1)^{{2n_{h}}}\right)^{\frac{1}{2}}}),$

for all $k\geq K$.

Without loss of generality we can assume that $2M_{f}\geq 1$. Since
$\|F\|_{L^{\infty}(0,c_{0}\mu_{f})}<1$, we deduce from Lemma 3.2 that

$\displaystyle\|F\|_{L^{\infty}(0,c_{0}\mu_{f})}\leq(2M_{f})^{1-w_{0}(c_{0}\mu_{f},K)}\|F\|_{L^{\infty}(0,K)}^{w_{0}(c_{0}\mu_{f},K)}.$

Applying now the derived lower bound in (23) on the last inequality, we get

$\displaystyle\|F\|_{L^{\infty}(0,c_{0}\mu_{f})}\leq(2M_{f})^{1-\eta(c_{0}\mu_{f},K)}\|F\|_{L^{\infty}(0,K)}^{\eta(c_{0}\mu_{f},K)}.$

Combining this inequality with the estimate in Lemma 3.2 achieves the proof of
the theorem.

## References

  * [1] S. Acosta, S. Chow, J. Taylor, and V. Villamizar, On the multi-frequency inverse source problem in heterogeneous media, Inverse problems, 28 (2012), p. 075013.
  * [2] H. Ammari and H. Zhang, Super-resolution in high-contrast media, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471 (2015), p. 20140946, https://doi.org/10.1098/rspa.2014.0946.
  * [3] K. Ammari and F. Triki, On weak observability for evolution systems with skew-adjoint generators, SIAM Journal on Mathematical Analysis, 52 (2020), pp. 1884–1902.
  * [4] G. Bao, P. Li, J. Lin, and F. Triki, Inverse scattering problems with multi-frequencies, Inverse Problems, 31 (2015), p. 093001.
  * [5] G. Bao, J. Lin, and F. Triki, A multi-frequency inverse source problem, Journal of Differential Equations, 249 (2010), pp. 3443–3465.
  * [6] G. Bao, J. Lin, F. Triki, et al., Numerical solution of the inverse source problem for the helmholtz equation with multiple frequency data, Contemp. Math, 548 (2011), pp. 45–60.
  * [7] G. Bao and F. Triki, Stability for the multifrequency inverse medium problem, Journal of Differential Equations, 269 (2020), pp. 7106–7128, https://doi.org/10.1016/J.JDE.2020.05.021.
  * [8] J. Cheng, V. Isakov, and S. Lu, Increasing stability in the inverse source problem with many frequencies, Journal of Differential Equations, 260 (2016), pp. 4786–4804.
  * [9] M. N. Entekhabi and V. Isakov, On increasing stability in the two dimensional inverse source scattering problem with many frequencies, Inverse Problems, 34 (2018), p. 055005.
  * [10] J. Garnier, H. Haddar, and H. Montanelli, The linear sampling method for random sources, (2022), https://doi.org/10.48550/ARXIV.2210.15560, https://arxiv.org/abs/2210.15560.
  * [11] J. Garnier and G. Papanicolaou, Passive Imaging with Ambient Noise, Cambridge University Press, 2016, https://doi.org/10.1017/CBO9781316471807.
  * [12] A. Henrot, Eigenvalues of elliptic operators, Extremum Problems for Eigenvalues of Elliptic Operators, (2006), pp. 1–16.
  * [13] P. Li, J. Zhai, and Y. Zhao, Stability for the acoustic inverse source problem in inhomogeneous media, SIAM Journal on Applied Mathematics, 80 (2020), pp. 2547–2559.
  * [14] W. C. H. McLean, Strongly elliptic systems and boundary integral equations, Cambridge university press, 2000.
  * [15] R. Nevanlinna, H. Behnke, H. Grauert, L. V. Ahlfors, D. C. Spencer, L. Bers, K. Kodaira, M. Heins, and J. A. Jenkins, Analytic functions, vol. 11, Springer, 1970.
  * [16] A. Osses and F. Triki, An improved spectral inequality for sums of eigenfunctions, arXiv e-prints, (2023), pp. arXiv–2312.
  * [17] F. Triki and C.-H. Tsou, Inverse inclusion problem: A stable method to determine disks, Journal of Differential Equations, 269 (2020), pp. 3259–3281.