# Automated optimization of force field parameters against ensemble-averaged
measurements with Bayesian Inference of Conformational Populations

Robert M. Raddi  Vincent A. Voelz<EMAIL_ADDRESS>Department of Chemistry,
Temple University, Philadelphia, PA 19122, USA.

(August 28, 2024)

###### Abstract

Accurate force fields are essential for reliable molecular simulations. These
models are refined against quantum mechanical calculations and experimental
measurements, which are subject to random and systematic errors. Bayesian
Inference of Conformational Populations (BICePs) is a reweighting algorithm
that reconciles simulated ensembles with sparse or noisy observables by
sampling the full posterior distribution of conformational populations and
experimental uncertainty. In this method, a metric called the BICePs score is
used to perform model selection, by calculating the free energy of "turning
on" the conformational populations under experimental restraints. This
approach, when used with improved likelihood functions to deal with
experimental outliers, can be used for force field validation (Raddi et al.
2023). Here, we extend the BICePs approach to perform automated force field
refinement while simultaneously sampling the full distribution of
uncertainties, using a variational method to minimize the BICePs score. To
demonstrate the utility of this method, we refine multiple interaction
parameters for a 12-mer HP lattice model using ensemble-averaged distance
measurements as restraints. To illustrate the resilience of BICePs in the
presence of unknown random and systematic errors, we assess the performance of
our algorithm through repeated optimizations and under various extents of
experimental error. Our results suggest that variational optimization of the
BICePs score is a promising direction for robust and automatic
parameterization of molecular potentials.

††preprint: AIP/123-QED

##  I Introduction

The accuracy of physical potentials used in molecular simulations continues to
improve, especially as training data for fitting energy surfaces from quantum-
mechanical calculations becomes more available.Boothroyd  _et al._ (2023);
Smith  _et al._ (2019) Of similar importance is developing models that can
accurately predict ensemble-averaged experimental measurements. Parameterizing
models for this purpose is a nontrivial task that requires continuous
refinement of microscopic parameters through iterated sampling of the
equilibrium distribution.

Such refinement often involves global minimization of the deviation between
experimental measurements (e.g. NMR observables) and theoretical predictions
from simulated ensembles (the “forward model”), with several complications to
consider. First, the forward model is itself an approximation containing some
error. Second, ensemble-averaged experimental measurements can be sparse
and/or noisy, susceptible to random and systematic errors, which are often
unknown a priori. Hence, force field refinement requires a mechanism to
integrate these multiple sources of uncertainty to automatically discover the
optimal parameter set across multiple systems.

Other challenges arise from practical considerations. The high dimensionality
and interdependence of force field parameter space presents a significant
challenge for refinement. When parameter spaces are large, reproducibility may
also be an issue; ideally, identical experimental data should produce
parameterizations that should converge towards the same optimal values. For
these reasons, automated optimization procedures with sufficiently smooth and
differentiable objective functions are preferred.

Numerous methodsKöfinger and Hummer (2021); Wang, Martinez, and Pande (2014);
Ge and Voelz (2018); Madin  _et al._ (2022) have been developed to address
most of these challenges. However, many algorithms do not include any
treatment of uncertainty in the training data, and others lack gradients for
automatic refinement. To address these challenges, Bayesian inference methods
have been developed.Köfinger and Hummer (2021); Madin  _et al._ (2022); Ge and
Voelz (2018) These methods estimate a Bayesian posterior distribution of
conformational populations by treating molecular simulation predictions as
prior information, weighted by a likelihood function constructed from the
experimental measurements and their uncertainties. The Bayesian posterior can
then be used to optimize the prior.

In this work, we extend the Bayesian Inference of Conformational Populations
(BICePs) algorithmVoelz, Ge, and Raddi (2021); Raddi  _et al._ (2023) to
perform automated force field refinement. BICePs is a reweighting algorithm
that refines structural ensembles against sparse and/or noisy experimental
observables, and has been used in many previous applications.Voelz and Zhou
(2014); Wan  _et al._ (2020); Hurley  _et al._ (2021); Raddi, Ge, and Voelz
(2023) A key advantage of BICePs is that it does not require uncertainty
estimates for the forward model; instead, it infers the posterior distribution
of these parameters directly from the data through MCMC sampling. BICePs also
computes a free energy-like quantity called the BICePs score that can be used
for model selection.Ge and Voelz (2018); Raddi  _et al._ (2023)

More recently, we have equipped BICePs with a replica-averaging forward model,
making it a maximum-entropy (MaxEnt) reweighting method, and unique in that no
adjustable regularization parameters are required to balance experimental
information with the prior.Raddi  _et al._ (2023) With this new approach, the
BICePs score becomes a powerful objective function to select and parameterize
optimal models. Here, we show that the BICePs score, which reflects the total
evidence for a model, can be used for variational optimization of model
parameters. The BICePs score contains a form of inherent regularization, and
has specialized likelihood functions that allow for the automatic detection
and down-weighting of data points subject to systematic error.Raddi  _et al._
(2023)

To efficiently optimize complex parameter spaces, we derive the first and
second derivatives of the BICePs score. We then show how to perform automatic
force field optimization against ensemble-averaged observables for several
test systems, including a protein lattice model with adjustable bead
interaction strengths, and a von Mises-distributed polymer model.

##  II Theory

BICePs uses a Bayesian statistical framework to treat the extent of
uncertainty in experimental observables, $\sigma$, as nuisance parameters.
Previous versions of BICePs sampled conformational states $X$ and uncertainty
parameter(s) $\sigma$ from the Bayesian posterior, which takes the form

$\overbrace{p(X,\sigma|D)}^{\text{posterior}}\propto\overbrace{p(D|X,\sigma)}^{\text{likelihood}}\overbrace{p(X)p(\sigma)}^{\text{priors}}.$
(1)

Here, the prior $p(X)$ comes from a theoretical model of conformational state
populations (typically from a molecular simulation), $p(D|X,\sigma)$ is a
likelihood function quantifying how well a forward model prediction $f(X)$
agrees with the experimental data $D$, and $p(\sigma)\sim\sigma^{-1}$ is a
non-informative Jeffreys prior.

When BICePs is equipped with a replica-averaged forward model, it becomes a
MaxEnt reweighting method in the limit of large numbers of replicas Pitera and
Chodera (2012); Cavalli, Camilloni, and Vendruscolo (2013); Cesari, Reißer,
and Bussi (2018); Roux and Weare (2013); Hummer and Köfinger (2015). The
posterior takes the general form

$\begin{split}p(&\mathbf{X},\bm{\sigma}|D)\propto\\\
&\prod_{r=1}^{N_{r}}\Bigl{\\{}p(X_{r})\prod_{j=1}^{N_{j}}\frac{1}{\sqrt{2\pi\sigma_{j}^{2}}}\exp\Big{[}-\frac{(d_{j}-f_{j}(\mathbf{X}))^{2}}{2\sigma_{j}^{2}}\Big{]}p(\sigma_{j})\Bigr{\\}}\end{split}$
(2)

where $\mathbf{X}$ is a set of $N_{r}$ conformation replicas, $d_{j}$ is an
observable in the set of $N_{j}$ ensemble-averaged experimental measurements,
and $f_{j}(\mathbf{X})=\frac{1}{N_{r}}\sum_{r}^{N_{r}}f_{j}(X_{r})$ is the
replica-averaged forward model prediction of observable $j$. The $\sigma_{j}$
values are nuisance parameters that capture uncertainty in the measurements as
well as the replica-averaged forward model. In (2), a Gaussian likelihood is
used, but more sophisticated models of uncertainty can be used to capture
outliers and systematic error with fewer parameters, as discussed below.
Markov chain Monte Carlo (MCMC) is used to sample the posterior.

The replica-averaged forward model $f(\mathbf{X})$ requires a more
sophisticated treatment of uncertainty, since it is approximating the ensemble
average as a replica average.Bonomi  _et al._ (2016). Therefore, $\sigma_{j}$
becomes a combination of both the Bayesian error $\sigma^{B}_{j}$ and the
standard error of the mean sigma $\sigma^{\text{SEM}}_{j}$,
$\sigma_{j}=\sqrt{(\sigma^{\text{B}}_{j})^{2}+(\sigma^{\text{SEM}}_{j})^{2}}$.
The finite sampling error, $\sigma^{\text{SEM}}$ is estimated by taking a
windowed average over our finite sample $f(\textbf{X})$ as
$\sigma^{SEM}_{j}=\sqrt{\frac{1}{N}\sum\nolimits^{N}_{r}(f_{j}(X_{r})-\langle
f_{j}(\mathbf{X})\rangle)^{2}}$, which decreases as the square root of the
number of replicas.

###  II.1 Accounting for systematic error and outliers.

BICePs is equipped with special likelihoods that are very robust in the
presence of systematic error. The Student’s likelihood model is a data error
model that marginalizes the uncertainty parameters for individual observables,
assuming that the level of noise is mostly uniform, except for a few erratic
measurements. This limits the number of uncertainty parameters that need to be
sampled, while still capturing outliers.

Consider a model where uncertainties $\sigma_{j}$ for particular observables
$j$ are distributed about some typical uncertainty $\sigma_{0}$ according to a
conditional probability $p(\sigma_{j}|\sigma_{0})$. We derive a posterior with
a single uncertainty parameter $\sigma_{0}$ by marginalizing over all
$\sigma_{j}$. For a single replica (for simplicity), the posterior is given by

$p(X_{r},\sigma_{0}|D)\propto
p(X_{r})\prod^{N_{j}}_{j=1}\int\limits_{\sigma^{\text{SEM}}}^{\infty}p(d_{j}|\mathbf{X},\sigma_{j})p(\sigma_{j}|\sigma_{0})d\sigma_{j}$
(3)

where $\sigma_{0}=\sqrt{(\sigma^{\text{B}})^{2}+(\sigma^{\text{SEM}})^{2}}$.

Modeling the prior on uncertainties $p(\sigma_{j}|\sigma_{0})$ with a
distribution with long tails is very useful because its long tail makes it
able to tolerate outliersSivia and Skilling (2006); Bonomi  _et al._ (2016).
In most cases, however, it is unclear a priori what distribution is best for
modeling the input data. To improve the situation, we introduce a model with
an additional nuisance parameter $\beta$, that is able to tune the extent of
the distribution’s tail:

$p(\sigma_{j}|\sigma_{0},\beta)=\frac{\Gamma((\beta+1)/2)}{\Gamma(\beta/2)}\frac{2\beta^{\beta}{\sigma_{0}}^{2\beta-1}}{\sqrt{\beta}{\sigma_{j}}^{2\beta}}\exp\left(-\frac{\beta{\sigma_{0}}^{2}}{\sigma_{j}^{2}}\right).$
(4)

where $\sigma_{0}$ is defined as above, and $1\leq\beta<\infty$. When this
distribution is inserted into the posterior, and marginalized over all
$\sigma_{j}$, the result is

$\begin{split}&p(X_{r},\sigma_{0},\beta|D)\propto
p(X_{r})\prod^{N_{j}}_{j=1}\int\limits_{\sigma^{\text{SEM}}}^{\infty}\frac{1}{\sqrt{2\pi}\sigma_{j}}\exp\left(-\frac{(d_{j}-f_{j}(\mathbf{X}))^{2}}{2\sigma_{j}^{2}}\right)\\\
&\quad\quad\quad\times\frac{\Gamma((\beta+1)/2)}{\Gamma(\beta/2)}\frac{2\beta^{\beta}{\sigma_{0}}^{2\beta-1}}{\sqrt{\beta}{\sigma_{j}}^{2\beta}}\exp\left(-\frac{\beta{\sigma_{0}}^{2}}{\sigma_{j}^{2}}\right)d\sigma_{j}\\\
&=p(X_{r})\prod^{N_{j}}_{j=1}\frac{\Gamma((\beta+1)/2)}{\Gamma(\beta/2)}\frac{1}{\sqrt{2\pi\beta}\sigma_{0}}\\\
&\quad\times\left[1+\frac{(d_{j}-f_{j}(X))^{2}}{2\beta\sigma_{0}^{2}}\right]^{-\beta}\gamma\left(\beta,\frac{(d_{j}-f_{j}(X))^{2}+2\beta\sigma_{0}^{2}}{2\left(\sigma^{\text{SEM}}\right)^{2}}\right),\end{split}$
(5)

Here, the marginal likelihood contains the lower incomplete gamma function,
$\gamma$. We call this the Student’s model because it is a variation of
Student’s $t$-distribution that can be interpolated between functional forms.
When $\beta=1$, the model is equivalent Metainference’s Outliers modelBonomi
_et al._ (2016). In the limit of $\beta\rightarrow\infty$, the likelihood
becomes Gaussian. When considering the full posterior, this extra nuisance
parameter is given a non-informative Jeffreys prior, $p(\beta)\sim\beta^{-1}$.
For more details regarding the Student’s model, see Figures S1-S2.

###  II.2 The BICePs score: a tool for quantitative model selection and
refinement

The BICePs score is a free energy-like quantity that rigorously characterizes
model quality. For a model of prior populations $p(X|\epsilon)$ parameterized
by $\epsilon$, with corresponding prior energy $E(X|\epsilon)=-\ln
p(X|\epsilon)$, the BICePs score $f(\epsilon)$ is calculated as the negative
logarithm of a Bayes factor,

$f(\epsilon)=-\ln\left(Z(\epsilon)\big{/}Z_{0}\right).$ (6)

$Z(\epsilon)$ is the total evidence of the specified model, marginalized over
all conformational replicas and uncertainty. By defining the BICePs energy
function as $u=-\ln p(\mathbf{X},\bm{\sigma}|D,\epsilon)$, it takes the form

$Z(\epsilon)=\int\int\exp\left(-u(\mathbf{X},\mathbf{\sigma}\mid
D,\epsilon)\right)d\mathbf{X}d\mathbf{\sigma},$ (7)

$Z_{0}$ is defined similarly, but for a uniform prior that does not depend on
$\epsilon$. $Z_{0}$ thus serves as a well-defined reference state. By
constructing a series of intermediates that scale the prior
$p_{\lambda}(X|\epsilon)\sim[p(X|\epsilon)]^{\lambda}$ as
$\lambda=0\rightarrow 1$, the BICePs score is calculated as the change in free
energy for “turning on” the prior (Figure 1).

Figure 1: Force field optimization by variationally minimizing the BICePs
score. Given ensemble-averaged experimental observables $D$, and a simulated
ensemble $p_{\lambda}(X|\epsilon)$ generated from molecular simulation using
initial force field parameters $\epsilon$ , an automated procedure in used to
find the optimal parameters $\epsilon^{*}$ that best match experiment. The
procedure is guided by the first and second derivatives of the BICePs score at
each iteration to propose new values of $\epsilon$. optimal parameters
$\epsilon$. This cycle is iterated until convergence is reached.

###  II.3 Variational minimization of the BICePs score to find optimal model
parameters.

We assume there exists some posterior distribution
$p^{*}(\mathbf{X},\mathbf{\sigma}|D,\epsilon^{*})$ with the maximal evidence,
derived from optimal parameters $\epsilon^{*}$. In this case, the free energy
$f(\epsilon^{*})=f^{*}$ is a global minimum, and any other set of parameters
$\epsilon$ has $f(\epsilon)-f^{*}\geq 0$.

Our objective is to find $\epsilon^{*}$ starting with trial parameters
$\epsilon$ and corresponding distribution
$p(\mathbf{X},\mathbf{\sigma}|D,\epsilon)$. To do this, we iteratively perform
BICePs calculations $\lambda=0\rightarrow 1$ to estimate $f(\epsilon)$, then
update the parameters to further minimize this quantity, until the stopping
criteria is reached. Given a differentiable prior, $f(\epsilon)$ is also
differentiable, enabling application of any number of gradient-based convex
optimization schemes, with a stopping criteria of $\frac{\partial
f}{\partial\epsilon}\big{|}_{\epsilon=\epsilon^{*}}=0$, within some tolerance.

The derivative of the BICePs score with respect to parameters $\epsilon$
reduces to the ensemble-averaged value of $\partial u/\partial\epsilon_{i}$,

$\begin{split}\frac{\partial
f}{\partial\epsilon_{i}}&=\iint\frac{1}{Z(\epsilon)}\left[\frac{\partial
u}{\partial\epsilon_{i}}\right]\exp\left(-u\right)d\mathbf{X}d\mathbf{\sigma}=\bigg{\langle}\frac{\partial
u}{\partial\epsilon_{i}}\bigg{\rangle},\end{split}$ (8)

where

$\begin{split}\frac{\partial
u}{\partial\epsilon_{i}}=\sum^{N}_{r=1}\biggl{\\{}\frac{\partial
E(X_{r}|\epsilon)}{\partial\epsilon_{i}}-\bigg{\langle}\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}}\bigg{\rangle}\biggr{\\}}\end{split}$
(9)

The reference ensemble $Z_{0}$ does not affect the gradient of $f$, since a
uniform prior is independent of $\epsilon$, and $\langle\partial
u_{0}/\partial\epsilon\rangle$ is zero (see section Automated optimization of
force field parameters against ensemble-averaged measurements with Bayesian
Inference of Conformational Populations for more details).

The second partial derivatives of the BICePs score with respect to parameters
$\epsilon_{i}$ and $\epsilon_{j}$ are:

$\displaystyle\frac{\partial^{2}f(\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}=\bigg{\langle}\frac{\partial^{2}u}{\partial\epsilon_{i}\partial\epsilon_{j}}\bigg{\rangle}-\biggl{(}$
$\displaystyle\bigg{\langle}\left(\frac{\partial
u}{\partial\epsilon_{i}}\right)\left(\frac{\partial
u}{\partial\epsilon_{j}}\right)\bigg{\rangle}$ (10)
$\displaystyle-\bigg{\langle}\left(\frac{\partial
u}{\partial\epsilon_{i}}\right)\left(\frac{\partial
u}{\partial\epsilon_{j}}\right)\bigg{\rangle}^{\prime}\biggr{)},$

where the notation $\langle\cdot\rangle^{\prime}$ denotes the ensemble average
with respect to $\left(\frac{1}{Z}\exp(-u)\right)^{2}$. The first term on the
right is the ensemble-averaged second derivative of the energy function $u$
with respect parameters $\epsilon_{i}$ and $\epsilon_{j}$:

$\begin{split}\frac{\partial^{2}u}{\partial\epsilon_{i}\partial\epsilon_{j}}=\sum^{N}_{r=1}\biggl{\\{}&\frac{\partial^{2}E(X_{r}|\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}-\bigg{\langle}\frac{\partial^{2}E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}\bigg{\rangle}\\\
&+\bigg{\langle}\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}}\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{j}}\right)\bigg{\rangle}\\\
&+\frac{\sum_{X}\exp\left(-2E(X|\epsilon)\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}}\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{j}}\right)}{\left(\sum_{X}\exp\left(-E(X|\epsilon)\right)\right)^{2}}\biggr{\\}}.\end{split}$
(11)

When $i$ = $j$, the second partial derivative of the BICePs score reduces to
the difference between the ensemble-averaged second derivative of the energy
$u$ and the variance of its first partial derivative:

$\begin{split}\frac{\partial^{2}f}{\partial\epsilon^{2}}&=\bigg{\langle}\frac{\partial^{2}u}{\partial\epsilon^{2}}\bigg{\rangle}-\left(\bigg{\langle}\left(\frac{\partial
u}{\partial\epsilon}\right)^{2}\bigg{\rangle}-\bigg{\langle}\frac{\partial
u}{\partial\epsilon}\bigg{\rangle}^{2}\right)\end{split}$ (12)

In practice, all of these quantities are estimated as ensemble-averaged
observables at the same time $f(\epsilon)$ is estimated, using the MBAR free
energy estimator. MCMC sampling at several intermediates $\lambda=0\rightarrow
1$ enables accurate estimates of all quantities (Figure 1).

##  III Methods

###  III.1 A toy protein lattice model

The HP lattice model of Dill and ChanDill  _et al._ (1995) is a simplified
protein folding model that represents a protein sequence as a chain of beads
on a 2-D square lattice. In this model, each amino acid residue is classified
as either a hydrophobic (H) or polar (P) bead. Each microscopic chain
conformation $X$ has an energy $V(X)$ proportional to the number of H-H
contacts $n(X)$,

$V(X)=-\epsilon\cdot n(X),$ (13)

A “contact” is defined as a non-sequential pair of hydrophobic residues at
adjacent lattice sites. As a test system, we consider a 12-mer protein
sequence HPHPHPHPPHPH (Figure 2a) . When $\epsilon$ is sufficiently large,
there is a single lowest-energy microstate that corresponds to the global free
energy minimum; such sequences are deemed foldable. For the 12-mer
HPHPHPHPPHPH sequence, there are 15037 unique microstate conformations that
can be binned into 72 macrostates, each with a unique set of contacts.

Figure 2: (a) The folded state of the 12-mer HP lattice model protein with
sequence HPHPHPHPPHPH. Highlighted in yellow are the favorable non-bonded
contacts between the hydrophobic residues (black beads). (b) Extended diagram
of the eight distance measurements used in this work.

To explore a system with multiple parameters, we modify the HP model to
include multiple interactions. Each bead with sequence position $i$ is
assigned a unique interaction strength, $\epsilon_{i}$, and a geometric-
average combination rule is used to compute the energy $V$ of each microstate:

$V=\sum_{{i,j}\text{ in contact}}-(\epsilon_{i}\epsilon_{j})^{1/2},$ (14)

where the sum is taken over all pairs of residues $(i,j)$ that are in contact.

The first derivative of the contact energy with respect to $\epsilon_{i}$ is
given by

$\frac{\partial V}{\partial\epsilon_{i}}=\sum_{{i,j}\text{ in
contact}}-\frac{1}{2}(\epsilon_{i}\epsilon_{j})^{-1/2}\epsilon_{j},$ (15)

and the second derivative of the contact energy with respect to $\epsilon_{i}$
and $\epsilon_{j}$ is given by:

$\frac{\partial^{2}V}{\partial\epsilon_{i}\partial\epsilon_{j}}=\sum_{{i,j}\text{
in
contact}}\frac{1}{4}(\epsilon_{i}\epsilon_{j})^{-3/2}\epsilon_{j}\epsilon_{i}.$
(16)

We use the 72 macrostates with unique contacts as a set of prior populations
for use with BICePs. For each macrostate $k$, the populations $p_{k}$ are
computed as

$p_{k}=\frac{g_{k}\exp(-V_{k})}{\sum_{k}g_{k}\exp(-V_{k})},$ (17)

where $V_{k}$ is the reduced energy (in units $k_{B}T$) of each microstate in
macrostate $k$, and $g_{k}$ is the multiplicity (the number of microstates
belonging to macrostate $k$).

The ensemble-averaged distance observables for each macrostate (and their
variance over all microstates belonging to it) are computed as follows. Let
$(d_{ij})_{k}$ represent the average distance between residues $i$ and $j$ in
the $k^{th}$ macrostate. Then, the ensemble-averaged distance observables are
given by:

$\langle d_{ij}\rangle=\sum_{k}p_{k}(d_{ij})_{k},$ (18)

and the variance of the distance observables is given by:

$\sigma^{2}(\langle d_{ij}\rangle)=\sum_{k}p_{k}^{2}\sigma^{2}((d_{ij})_{k}),$
(19)

where $\sigma^{2}((d_{ij})_{k})$ is the variance of the average distance
between residues $i$ and $j$ in the $k^{th}$ macrostate.

For ease of interpretation, we assume that the distance units of the square
lattice are Ångstroms (1Å=$10^{-10}$ m).

##  IV Results and Discussion

###  IV.1 The BICePs score landscape correctly identifies optimal parameters
for an HP lattice model.

To understand BICePs’ ability to learn optimal parameters, we start by
scanning the BICePs score landscape over a single $\epsilon$ parameter. For
this purpose, we use an HP model where all six hydrophobic beads share the
same value of $\epsilon_{i}=\epsilon$ .

To generate a set of experimental distance observables for testing, we set the
"true" value of $\epsilon$ to 1.0., and compute average distances and
variances of eight distances for each macrostate (Figure 2b). The BICePs score
and its first and second derivatives are then evaluated across a range of
$\epsilon$ values, spanning from 0.0 to 5.5, in increments of 0.125, where
prior energies are unique for each $\epsilon$ value. BICePs scores and
derivatives are averaged over five independent calculations using 100k steps
of MCMC with 8 replicas. The associated error bars represent the standard
error of the mean across these five instances.

The resulting $f(\epsilon)$ landscape pinpoints the optimal parameter as a
minimum, and correctly computes the first and second derivatives (Figure 3).
Despite the stochastic sampling inherent to BICePs, the BICePs score landscape
is remarkably smooth and precise, revealing the optimal epsilon to be
$\epsilon^{*}=1.0$, the value with the lowest BICePs score,
$f(\epsilon^{*})$(Figure 3a.), with its derivative $df/d\epsilon$ at this
location equal to 0.0 (Figure 3b.). This optimal model is achieved when the
overlap between the likelihood and prior distributions is maximized (Figure
1). When $\epsilon$ is significantly larger than the true value, the
uncertainty in $f$ and $\partial f/\partial\epsilon$ increases. In this case,
the physical interpretation is that for large values of $\epsilon$ the 12-mer
primarily populates the folded macrostate, which results in poor overlap
between the likelihood and prior, and more finite sampling error.

Figure 3: (a) 1-D scans over $\epsilon$ to reveal the landscape of the BICePs
score (blue dots). The green tangent lines are the computed derivatives at
each $\epsilon$ value. Uncertainties in the BICePs scores and their
derivatives come from five independent scans along $\epsilon$. (b) The
derivative of the BICePs score (blue dots) and the second derivative of the
BICePs score (yellow lines) at each epsilon value. The dotted black line at
$\epsilon=1.0$ shows the true value, which is where the derivative of the
BICePs score equals zero. BICePs calculations are run using the Student’s
model for 100k steps with 8 replicas.

###  IV.2 Automated optimization schemes.

With the ability to compute first and second derivatives of the BICePs score,
we next explored the performance of automated optimization methods L-BFGS-
BByrd  _et al._ (1995) and Trust-NCGWright (2006) on the single-$\epsilon$
system. L-BFGS-B is a first-order approach that uses information from
gradients, while Trust-NCG is a second-order approach that additionally makes
use of information from the Hessian. All default optimization parameters from
SciPy were used except a value of ftol=1e-8 (default: ftol=1e-9), and the
bounds for each parameter to be optimized was set to be in the range [0, 10].
With only a single parameter, automatic refinement using both optimization
methods quickly finds the minimum in a single iteration (Figure S3). However,
use of the Hessian (second-order) enables faster convergence, and tends to be
slightly more robust for this particular example.

###  IV.3 The Student’s model confers resiliency when dealing with
experimental error.

We next examined the robustness of the BICePs score to errors and
uncertainties present in experimental data, by introducing random and
systematic errors to the 2–11 and 0–11 distance observables (Figure S4). When
we used BICePs with a standard Gaussian likelihood model in which a single
uncertainty parameter $\sigma$ is used for all experimental observables,
introducing random and systematic errors significantly reduces the accuracy of
the model. When 2–11 and 0–11 distances are perturbed by +3 to +5 Å, the
BICePs score predicts a value of $\epsilon$ smaller than the “true” value,
predicting structures that are primarily "unfolded". This in turn causes a
leftward shift in the gradient plot (Figure 3b).

In contrast, the Student’s likelihood model deals with the outliers much more
handily. The minimum BICePs score correctly identifies the optimal value of
$\epsilon$, and the leftward shift in the gradient plot is absent. The
Student’s model is able to automatically identify outliers directly from the
input data, and mitigate their influence, even when 25% of the data (2 of the
8 distances) are corrupted.

To more critically assess the performance of the Gaussian and Student’s models
in dealing with systematic errors, we examine the extent to which the gradient
of the BICePs score deviates from zero at the “true” value of
$\epsilon^{*}=1$. We introduced random perturbations to the 2–11 and 0–11
distances of various standard deviations, $\sigma_{data}$ (Figure 4). When the
Gaussian likelihood model was used, the computed gradient at $\epsilon^{*}=1$
increases proportionally with the magnitude of $\sigma_{data}$., and becomes
particularly unreliable when the error surpasses 0.75 Å. When the Student’s
model was used, however, we find that such deviations are reduced. As
$\sigma_{data}$ increases from zero, the gradient at $\epsilon^{*}=1$ begins
to increase, but then begins to decrease beyond $\sigma_{data}$ = 0.6 Å. This
happens because the perturbations in the distance observables become large
enough to be detected as outliers by the Student’s model.

Figure 4: Comparative analysis in performance of the Student’s likelihood
model (blue) and a Gaussian likelihood model (red), when random and systematic
error of varying magnitude ($\sigma_{data}$) is introduced to the 2-11 and
0-11 distances. The vertical axis shows the derivative of the BICePs score
evaluated at the “true” target value, $\epsilon=1$. This value corresponds to
the "True" minima, illustrated in Figure 3. The Gaussian likelihood’s
derivative becomes notably less dependable when data incorporates errors,
especially when $\sigma_{data}$ exceeds 0.75. surpassing one standard
deviation. The values shown were calculated using 5000 random perturbations to
the 2-11 and 0-11 distances, and represent the average of 300 BICePs
calculations. Error bars represent the standard deviation.

###  IV.4 Multiparameter optimization

While our toy protein lattice model is very simple, the inference problem it
presents is challenging, especially as we increase the number of epsilon
parameters to refine. Our 12-mer system has only eight experimental distances
as experimental restraints, and these are significantly influenced by
ensemble-averaging. Moreover, the protein lattice model has a highly
cooperative folding landscape, which can lead to correlated distance
restraints that don’t independently provide informative insights.

With this in mind, we further tested multi-parameter minimization of the
BICePs score, using the modified HP lattice model with tunable interaction
strengths $\epsilon_{i}$. Ideally, our algorithm should be able to efficiently
find the global minimum of $f(\epsilon_{1},\epsilon_{2},...)$ starting from
many different initial parameter sets, and be resilient against random and
systematic errors. For two or more $\epsilon_{i}$ parameters, it becomes
computationally expensive to perform multidimensional scans of the BICePs
score landscape, necessitating automated parameter optimization.

We first generated test data for a model where we set the “true” values to
$\epsilon_{2}=1.25,\epsilon_{4}=1.5$ and the remaining $\epsilon_{i}$ set to
3.0. We then tested two-parameter $(\epsilon_{2},\epsilon_{4})$ optimizations,
starting from three different initial parameter sets
$(\epsilon_{2}^{0},\epsilon_{4}^{0})$ = $\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$,
using both L-BFGS-B and ‘trust-ncg’ optimization for BICePs scores calculated
using 8 replicas. Optimizations using L-BFGS-B are unable to fully reach the
“true” values of the interaction strength parameters
($\epsilon_{2}^{\text{True}}=1.25$, $\epsilon_{4}^{\text{True}}=1.5$) within
10 iterations (Figure S5). The ‘trust-ncg’ optimization, however, successfully
found the correct minima within 10 iterations, demonstrating the notable
advantage of second-order methods (Figure S6).

To visualize the two-parameter BICePs score landscape
$f(\epsilon_{2},\epsilon_{4})$, we constructed a smooth 2-D landscape by
averaging five repeated evaluations over a (6$\times$6) set of grid points,
and trained a Gaussian process regression on the grid values using an RBF
kernel. The landscape matches the computed BICePs scores, and shows that
trajectories of successive ($\epsilon_{2},\epsilon_{4}$) values correctly
follow gradients, normal to the landscape contour lines.

The reported uncertainties in the optimized $(\epsilon_{2},\epsilon_{4})$
values comes from an estimation of the covariance through inversion of the
Hessian (the matrix of second partial derivatives of the BICePs score):

$\sigma_{\epsilon}=\sqrt{\mathrm{diag}\left[\mathbf{H}^{-1}\right]}.$ (20)

Conceptually, this means that the estimated uncertainties in the best-fit
values represent the widths of the basins on the BICePs score landscape. This
can sometimes lead to an overestimation of the uncertainty with trust-region
optimization schemes, especially with flatter basins like those we observe in
the BICePs score landscape. Conversely, the sharper the curvature (larger
eigenvalues of the Hessian), the more certain we are about the location of the
minimum in that direction, resulting in smaller computed uncertainties.
Inspection of Equations (9), (12) and (39) show that $\sigma_{\epsilon}\propto
1/\sqrt{N}$, where $N$ is the number of BICePs replicas. This suggests that
one can obtain smaller uncertainties in the optimized parameters by using more
replicas (albeit at greater computational expense). Indeed, repeating the
‘trust-ncg’ optimization with 32 BICePs replicas results in average optimized
parameters of $\epsilon_{2}$ = 1.12 $\pm$ 0.39, and $\epsilon_{4}$ = 1.43
$\pm$ 0.34, which have smaller uncertainty estimates compared to the
calculations made with 16 replicas ($\sigma_{\epsilon_{2}}$ = 0.59, and
$\sigma_{\epsilon_{4}}$ = 0.47). (Figure S7).

In all our multiparameter optimizations, we use a standardized convergence
criterion, determined by assessing the average relative change between the old
parameters $\epsilon^{\text{old}}$ and the new parameters
$\epsilon_{i}^{\text{new}}$ at each iteration, calculated as

$\frac{1}{N_{\epsilon}}\sum_{i=1}^{N_{\epsilon}}\left|\frac{\epsilon_{i}^{\text{old}}-\epsilon_{i}^{\text{new}}}{\epsilon_{i}^{\text{old}}}\right|<0.05,$
(21)

where convergence is considered to have been reached when the average change
in parameters is less than 0.05. In addition, we strictly require that the
convergence condition is satisfied three times in succession.

###  IV.5 Multiparameter optimization in the presence of systematic error

To assess the resiliency of multiparameter optimization to systematic error,
we added extra complexity to the optimization problem by introducing
systematic error of +3 Å and +3.5 Å to the 2–11 and 4–9 distances observables,
respectively. This makes total error in the data $\sigma_{data}=1.63$ Å, a
significant obstacle to the refinement of parameters $\epsilon_{2}$ and
$\epsilon_{4}$ , considering these beads each participate in only two of
measured distances (Figure 2b).

We performed ‘trust-ncg’ optimization using the Student’s model with 200k MCMC
steps and 32 replicas, starting from three different initial parameters:
$(\epsilon_{2}^{0},\epsilon_{4}^{0})$ = $\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$
(Figure 5). Indeed, repeating the ‘trust-ncg’ optimization with 32 BICePs
replicas results in average optimized parameters of $\epsilon_{2}$ = 1.24
$\pm$ 0.41, and $\epsilon_{4}$ = 1.31 $\pm$ 0.33, which have smaller
uncertainty estimates compared to the calculations made with 16 replicas
($\epsilon_{2}$ = 1.28 $\pm$ 0.59, and $\epsilon_{4}$ = 1.31 $\pm$ 0.47).
Notably, this was achieved in only a few iterations, showing the algorithm’s
proficiency in finding the optimal parameters, irrespective of the initial
reference parameters.

Figure 5: Average traces over a total of 25 independent rounds of parameter
$(\epsilon_{2},\epsilon_{4})$ optimizations using second-order (trust-ncg)
method with BICePs, for a maximum of ten iterations. Optimizations converge to
the same parameters ($\epsilon_{2}^{\text{True}}=1.25$,
$\epsilon_{4}^{\text{True}}=1.5$) when starting from different initial
parameters $(\epsilon_{2}^{0},\epsilon_{4}^{0})$ =
$\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$. Average BICePs optimized parameters
were determined to be $\epsilon_{2}$ = 1.24 $\pm$ 0.41, and $\epsilon_{4}$ =
1.31 $\pm$ 0.33., where the uncertainties are esimated from the inverse
Hessian. The BICePs score landscape was generated from the average values of
five scans over $\epsilon_{2}$ and $\epsilon_{4}$. All calculations used the
Student’s model with 200k MCMC steps and 32 replicas. The experimental data is
corrupted with systematic error in the 2-11 and 4-9 distances for +3 and +3.5Å
shift, respectively. The total error in the data is $\sigma_{data}=1.63$.

In adherence to best practicesBeiranvand, Hare, and Lucet (2017), we compute
accuracy profiles for two-parameter ‘trust-ncg’ optimizations using 32 BICePs
replicas with no systematic error (Figure S8) and systematic error (Figure
S9). The accuracy profiles show the relative error
$||\vec{\epsilon}-\vec{\epsilon^{*}}||_{2}=\sqrt{\sum_{i}(\epsilon_{i}-\epsilon^{*}_{i})^{2}}$
as a function of number of iterations. Optimizations with no systematic error
were repeated using only 16 BICePs replicas (Figure S10). A comparison of the
16-replica and 32-replica results shows that the variance between independent
optimization traces decreases as the number of replicas increases. The error
bars plotted in the accuracy profiles report uncertainties estimated from the
inverse Hessian. In each case, it can be seen that the actual variation in the
relative error across many optimizations is well within the Hessian-estimated
uncertainty.

As we did with the one-parameter optimizations, we also examined the
performance of two-parameter optimizations when using a standard Gaussian
likelihood model, in the absence of systematic error (Figure S11a) and in the
presence of systematic error (Figure S11b). As before, we find that the
accuracy of the refined parameters drops significantly when using a standard
Gaussian likelihood model. With only a single parameter to estimate the
uncertainty in the data set of distances, the Gaussian likelihood model is
unable to properly detect the two outliers, and as a result predicts optimal
$\epsilon_{2}$ and $\epsilon_{4}$ values that are much smaller than the “true”
values. Conversely, the accuracy of the Student’s model is very robust in the
presence of systematic error (Figure S12). These results yet again highlight
the ability of BICePs to handle uncertainties and systematic error
comprehensively when used with a likelihood model that can deal with outliers.

###  IV.6 Performance in higher-dimensional parameter landscapes

We next attempted to optimize a three-parameter BICePs score landscape, by
introducing an additional parameter, $\epsilon_{6}$. This is a challenging
task that may take many iterations to reach the convergence criteria. Our
objective was to examine the typical percentage of optimization traces that
converge within a designated number of iterations. For this test, we refrained
from adding systematic error to the data.

Five initial starting points were selected:
$(\epsilon_{2}^{0},\epsilon_{4}^{0},\epsilon_{6}^{0})$ = $(4.0,5.0,3.0),$
$(0.5,5.0,0.5)$, $(2.0,4.0,4.0)$, $(1.0,0.5,3.0)$, $(3.0,0.5,0.5)$.
Optimization targeted the “true” interaction strength parameters:
$\epsilon_{2}^{\text{True}}=1.25$, $\epsilon_{4}^{\text{True}}=1.5$, and
$\epsilon_{6}^{\text{True}}=1.5$, with an acceptable average relative
deviation of 0.2 (refer to Figure 6). Twenty-five independent optimizations
were performed per starting point, using 32 BICePs replicas and 200k MCMC
steps per iteration.

For these tests, we required that the convergence condition must be satisfied
two times in succession, which is a typical stopping point. Traces of the
relative error throughout the optimization show a heterogeneous collection of
convergence trajectories, depending on the starting point (Figure 6). At each
iteration, we track the percentage of traces that remain unconverged. For some
starting points, all optimization trajectories converge within ten iterations,
while for others, up to 28% remain unconverged going into the tenth iteration.
The average optimized parameters from the converged trajectories were
determined to be $(\epsilon_{2}=1.17\pm 0.38,\epsilon_{4}=1.37\pm
0.32,\epsilon_{6}=1.51\pm 0.37)$, where the uncertainties are calculated from
the inverse Hessian. An overlay of average traces on the BICePs score
landscape is shown in Figure S13.

A robust correlation ($R^{2}=0.96$) between the converged BICePs score and the
relative error in the parameters further confirms the validity of our approach
in higher dimensions (see Figure S14). This relationship illustrates that an
increased distance from the “true” parameters (i.e. a higher relative error)
corresponds to an enlargement of the BICePs score $f$. A comparative analysis
of average BICePs scores across optimization traces (Figure S15) and the
landscape contours (Figure S13) indicates that our Gaussian process regression
does a fine job of approximating the actual BICePs score landscape in higher
dimensions.

Figure 6: Accuracy profile of optimization traces, demonstrating convergence
with a relative error of 0.2. A total of 25 independent rounds of parameter
$(\epsilon_{2},\epsilon_{4},\epsilon_{6})$ optimizations were conducted for
each of the five sets of initial parameters
$(\epsilon_{2}^{0},\epsilon_{4}^{0},\epsilon_{6}^{0})$: $(4.0,5.0,3.0)$
(gray), $(0.5,5.0,0.5)$ (yellow), $(2.0,4.0,4.0)$ (green), $(1.0,0.5,3.0)$
(magenta), $(3.0,0.5,0.5)$ (cyan). Color-coded percentages report the fraction
of trajectories that remain unconverged at each iteration. Error bars (black)
represent the uncertainties in predicted parameters for each iteration,
calculated using the diagonalized inverse Hessian (Equation (20)) and averaged
over the 25 independent optimizations. Average optimized parameter values from
converged trajectories were $(\epsilon_{2}=1.17\pm 0.38,\epsilon_{4}=1.37\pm
0.32,\epsilon_{6}=1.51\pm 0.37)$.

###  IV.7 Köfinger & Hummer’s toy polymer model

Köfinger and Hummer (KH) recently developed a method to optimize molecular
simulation force fields using a Bayesian inference maximum-entropy reweighting
method called BioFF.Köfinger and Hummer (2021) To validate their method, they
introduced a simple parameterizable 2-D polymer model, and demonstrated that
BioFF can efficiently discover the optimal parameters that reproduce ensemble-
averaged observables. Here, we demonstrate that BICePs optimization can also
reproduce this result.

The KH polymer model is a chain of $N$ beads, each separated by unit length,
with orientations defined by $N-1$ angles $\Delta\phi$ between contiguous bond
vectors. Each bond angle is distributed according a von Mises distribution
parameterized by $\kappa$, which controls the polymer’s stiffness. Full
details about this model are given in Appendix B. Unlike the HP lattice model,
which possesses an exact solution for the prior, the KH polymer model requires
estimation of the prior through statistical sampling, which is subject to
finite sampling error.

To implement KH model optimization with BICePs, we derived the first and
second partial derivatives of the energy function required for computation of
the BICePs score (see SI). Following Köfinger and Hummer, we consider a
100-bead chain, and use the end-to-middle distance (between bead 1 and bead
50) as the ensemble-averaged observable, which we estimate as the mean of
10000 samples from the von Mises distribution with the “true” value of
$\kappa$ set to 10.0.

We then performed ‘trust-ncg’ optimization using 200 BICePs replicas and 100k
MCMC steps, starting from an initial parameter value of $\kappa$ = 20.0. At
the first optimization step, draw 500 samples from the von Mises distribution,
and use their energies to compute their Boltzmann populations; each sample is
a separate conformational state. In each subsequent optimization step, we draw
500 more samples to add to the collection, and compute energies and Boltzmann
populations for the combined collection. The result of this procedure is
highly efficient discovery of the optimal value of $\kappa^{*}$ = 10, reaching
convergence within two iterations (Figure S16).

While BICePs optimization achieves similar results to BioFF, it is worth
mentioning some of the key differences of the two methods. BioFF is a maximum-
entropy method that reweights a set of reference populations by minimizing a
negative log-likelihood that uses a regularization parameter, $\theta$, to
balance two terms: the relative entropy with respect to the reference, and a
$\chi^{2}$ term representing the deviation between model predictions of
ensemble-averaged quantities and experiment. The BioFF-based optimization
algorithm iteratively (1) reweights the populations given a current set of
force field parameters, then (2) optimizes the force field parameters given
the updated populations, repeating this cycle until convergence.

One issue with this approach is that it requires a coarse-grained description
of the relevant conformational states and their populations, which need to be
revised as the force field parameters are optimized. Köfinger and Hummer deal
with this by keeping track of each new thermodynamic ensemble being sampled,
and using the running collection of all the samples obtained to estimate each
new set of coarse-grained state populations with the binless WHAMShirts and
Chodera (2008); Rosta  _et al._ (2011) free energy estimator. In our BICePs
approach, we avoid this expense by simply adding new samples as additional
conformational states, and recalculating populations as their Boltzmann
weights.

Another issue with maximum-entropy methods like BioFFKöfinger, Różycki, and
Hummer (2019); Köfinger and Hummer (2021) and BioEnBottaro, Bengtsen, and
Lindorff-Larsen (2020) is the need to choose an appropriate value for the
regularization parameter $\theta$. BICePs does not require a regulation
parameter, because the BICePs score is computed by sampling the entire
posterior distribution of uncertainty parameters $\sigma$. (In BICePs, the
$\sigma$ are the analogous quantities that balance prior information against
the experimental restraints.)

Very recently, Köfinger & Hummer put the regularization parameter $\theta$ on
a firmer theoretical footing by relating it to the variance of the reduced
energy difference, thus enabling estimation of $\theta$ a priori.Köfinger and
Hummer (2023). Another recent work, Gilardoni et al. Gilardoni, Fröhlking, and
Bussi (2023) introduced an approach similar to BICePs, in that ensemble
refinement and force field parameter refinement are integrated into a singular
process. This method bears resemblance to both BioEn and BioFF, as it employs
a regularization parameter to scale the Kullback-Leibler divergence. Gilardoni
et al.’s technique is distinguished by the incorporation of two distinct
confidence hyperparameters. The first parameter is designed to adjust the
reliability of the prior, while the second modulates the trustworthiness of
the force field corrections.

###  IV.8 Further Discussion

Our work here greatly improves upon previous work by Ge et al., which used an
HP lattice model as an example system to perform force field parameterization
and model selection using BICePs.Ge and Voelz (2018). That work approached the
problem through brute force, scanning parameter space in search of the minimum
BICePs score, which would be costly and impractical for multiple parameters.
Another issue with that work was the use of a previous version of BICePs
(called BICePs v1.0) that does not use replica-averaging, and therefore was
not a maximum-entropy method. Because of this, applying BICePs in practical
cases, such as refining conformational populations against NMR observables for
beta-hairpins trploop 2A and 2B, resulted in ensembles dominated by the
state(s) that best agree with experiment. Since that work, we have developed
the current version of BICePs that performs replica averaging and is equipped
with specialized likelihood functions that can account for systematic
error.Raddi, Ge, and Voelz (2023) The HP lattice model continues to serve as
an excellent test system for force field refinement as it grants the option to
use multiple parameters and provides exact solutions for generated ensembles.

While our algorithm has the advantage of being able to incorporate priors on
force field parameters, this is a possibility we do not pursue in this work.
Obtaining priors for force field parameters is a more involved approach, but
we leave room for this in future development.

###  IV.9 Optimization across multiple systems

An eventual application of our BICePs optimization method is transferable
force field parameterization across many different molecular systems (e.g. a
set of proteins) against ensemble-averaged experimental measurements.

Our approach can easily be adapted for this purpose, and performed efficiently
so that each molecular system can be evaluated in a separate BICePs
calculation. To do this, suppose we have a set of molecular systems indexed by
$s$, of total number $N_{s}$. We make use of the partial derivatives of the
BICePs score with respect to each $\epsilon_{i}$ to determine a weight for
each system $s$. The normalized weight for the $i^{th}$ parameter of system
$s$ is

$w_{i,s}=\frac{\left(\frac{\partial
f}{\partial\epsilon_{i}}\right)_{s}\left(\epsilon_{i,s}-\epsilon_{i}^{\text{old}}\right)}{\sum_{s=1}^{N_{s}}\left(\frac{\partial
f}{\partial\epsilon_{i}}\right)_{s}\left(\epsilon_{i,s}-\epsilon_{i}^{\text{old}}\right)}.$
(22)

At each optimization iteration, we update the $i^{th}$ parameter,
$\epsilon_{i}^{\text{new}}$ according to a weighted average across all of the
systems:

$\epsilon_{i}^{\text{new}}\leftarrow\sum^{N_{s}}_{s=1}w_{i,s}\epsilon_{i,s}.$
(23)

In this scheme, the derivatives play a pivotal role, where the weights are
determined by how significantly each parameter change influences the BICePs
score. This provides a balance, ensuring parameters from a system with larger
effects on the score have more influence in determining the next set of
parameters $\epsilon_{\text{new}}$.

###  IV.10 Outlook

Looking ahead, there is tremendous potential to extend the scope of the BICePs
optimization against ensemble averaged observables to many applications. One
area of application would be to optimize physics-based molecular mechanics
force fields for proteins against the many published solution-NMR data that
exists.Cavender  _et al._ (2023) Another application would be to optimize
parameters of general-purpose force fields for bespoke applications. Examples
of this may include refining existing all-atom microscopic models to better
predict macroscopic properties (e.g. ensemble-average properties of disordered
peptides like radius of gyration), or developing classes of models that can
better predict the experimental observables for peptidomimetics without having
to perform expensive fitting to quantum mechanical energy surfaces. In
principle, any parameterizable model for the prior can be used with BICePs,
including neural-network based potentials, or even generative models obtained
from deep learning methods.

##  V Conclusion

In this work, we showed how the BICePs approach can be used for automatic and
robust optimization of force field parameters against ensemble-averaged
experimental measurements. This is done by variational minimization of the
BICePs score, a free energy-like quantity characterizing how well a model
agrees with the experimental data, calculated by sampling over the posterior
distribution of possible uncertainties. Using information from first and
second partial derivatives of the BICePs score with respect to model
parameters, we show how robust multi-parameter optimization of the BICePs
score can be performed. Demonstrations in simple toy protein lattice models
and polymer models show the utility of the method, and how specialized
likelihood functions to deal with outliers are especially good at dealing
systematic error. These results suggest automated variational optimization of
the BICePs score is a powerful approach for the parameterization of molecular
potentials.

## Conflicts of interest

Authors declare no conflicts of interest.

## Data and software

The BICePs algorithm is freely available at github.com/vvoelz/biceps and can
be installed using `pip install biceps`. The HP lattice model is open source
software, located at: github.com/vvoelz/HPSandbox. All data, scripts, Jupyter
notebook examples and analysis can be found here:
github.com/robraddi/biceps_FFOpt. Detailed documentation and tutorials can be
found here: https://biceps. readthedocs.io. For any issues or questions,
please submit the request on GitHub.

## Acknowledgements

RMR and VAV are supported by National Institutes of Health grant R01GM123296.

## References

  * Boothroyd  _et al._ (2023) S. Boothroyd, P. K. Behara, O. C. Madin, D. F. Hahn, H. Jang, V. Gapsys, J. R. Wagner, J. T. Horton, D. L. Dotson, M. W. Thompson, _et al._ , “Development and benchmarking of open force field 2.0. 0: The sage small molecule force field,” Journal of Chemical Theory and Computation (2023).
  * Smith  _et al._ (2019) J. S. Smith, B. T. Nebgen, R. Zubatyuk, N. Lubbers, C. Devereux, K. Barros, S. Tretiak, O. Isayev, and A. E. Roitberg, “Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning,” Nature communications 10, 2903 (2019).
  * Köfinger and Hummer (2021) J. Köfinger and G. Hummer, “Empirical Optimization of Molecular Simulation Force Fields by Bayesian Inference,” The European Physical Journal B 94, 245 (2021).
  * Wang, Martinez, and Pande (2014) L.-P. Wang, T. J. Martinez, and V. S. Pande, “Building Force Fields: An Automatic, Systematic, and Reproducible Approach,” The journal of physical chemistry letters 5, 1885–1891 (2014).
  * Ge and Voelz (2018) Y. Ge and V. A. Voelz, “Model Selection Using BICePs: a Bayesian Approach for Force Field Validation and Parameterization,” J. Phys. Chem. B 122, 5610–5622 (2018).
  * Madin  _et al._ (2022) O. C. Madin, S. Boothroyd, R. A. Messerly, J. Fass, J. D. Chodera, and M. R. Shirts, “Bayesian-Inference-Driven Model Parametrization and Model Selection for 2CLJQ Fluid Models,” J. Chem. Inf. Model. 62, 874–889 (2022).
  * Voelz, Ge, and Raddi (2021) V. A. Voelz, Y. Ge, and R. M. Raddi, “Reconciling Simulations and Experiments with BICePs: a Review,” Front. Mol. Biosci. 8, 661520 (2021).
  * Raddi  _et al._ (2023) R. M. Raddi, T. Marshall, Y. Ge, and V. Voelz, “Model Selection Using Replica Averaging with Bayesian Inference of Conformational Populations,” (2023).
  * Voelz and Zhou (2014) V. a. Voelz and G. Zhou, “Bayesian Inference of Conformational State Populations from Computational Models and Sparse Experimental Observables,” J. Comput. Chem. 35, 2215–2224 (2014).
  * Wan  _et al._ (2020) H. Wan, Y. Ge, A. Razavi, and V. A. Voelz, “Reconciling Simulated Ensembles of Apomyoglobin with Experimental Hydrogen/Deuterium Exchange Data Using Bayesian Inference and Multiensemble Markov State Models.” J. Chem. Theory Comput. 16, 1333–1348 (2020).
  * Hurley  _et al._ (2021) M. F. D. Hurley, J. D. Northrup, Y. Ge, C. E. Schafmeister, and V. A. Voelz, “Metal Cation-Binding Mechanisms of Q-Proline Peptoid Macrocycles in Solution,” J. Chem. Inf. Model. 61, 2818–2828 (2021).
  * Raddi, Ge, and Voelz (2023) R. M. Raddi, Y. Ge, and V. A. Voelz, “BICePs V2. 0: Software for Ensemble Reweighting Using Bayesian Inference of Conformational Populations,” Journal of chemical information and modeling 63, 2370–2381 (2023).
  * Pitera and Chodera (2012) J. W. Pitera and J. D. Chodera, “On the Use of Experimental Observations to Bias Simulated Ensembles,” Journal of chemical theory and computation 8, 3445–3451 (2012).
  * Cavalli, Camilloni, and Vendruscolo (2013) A. Cavalli, C. Camilloni, and M. Vendruscolo, “Molecular Dynamics Simulations with Replica-Averaged Structural Restraints Generate Structural Ensembles According to the Maximum Entropy Principle,” The Journal of chemical physics 138, 03B603 (2013).
  * Cesari, Reißer, and Bussi (2018) A. Cesari, S. Reißer, and G. Bussi, “Using the Maximum Entropy Principle to Combine Simulations and Solution Experiments,” Computation 6, 15 (2018).
  * Roux and Weare (2013) B. Roux and J. Weare, “On the Statistical Equivalence of Restrained-Ensemble Simulations with the Maximum Entropy Method,” The Journal of chemical physics 138, 02B616 (2013).
  * Hummer and Köfinger (2015) G. Hummer and J. Köfinger, “Bayesian Ensemble Refinement by Replica Simulations and Reweighting,” The Journal of chemical physics 143, 12B634_1 (2015).
  * Bonomi  _et al._ (2016) M. Bonomi, C. Camilloni, A. Cavalli, and M. Vendruscolo, “Metainference: A Bayesian Inference Method for Heterogeneous Systems,” Sci. Adv. 2, e1501177 (2016).
  * Sivia and Skilling (2006) D. Sivia and J. Skilling,  _Data Analysis: a Bayesian Tutorial_ (OUP Oxford, 2006).
  * Dill  _et al._ (1995) K. A. Dill, S. Bromberg, K. Yue, H. S. Chan, K. M. Ftebig, D. P. Yee, and P. D. Thomas, “Principles of protein folding—a perspective from simple exact models,” Protein science 4, 561–602 (1995).
  * Byrd  _et al._ (1995) R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM Journal on scientific computing 16, 1190–1208 (1995).
  * Wright (2006) S. J. Wright,  _Numerical optimization_ (2006).
  * Beiranvand, Hare, and Lucet (2017) V. Beiranvand, W. Hare, and Y. Lucet, “Best Practices for Comparing Optimization Algorithms,” Optim. Eng. 18, 815–848 (2017).
  * Shirts and Chodera (2008) M. R. Shirts and J. D. Chodera, “Statistically Optimal Analysis of Samples from Multiple Equilibrium States,” J. Chem. Phys. 129, 124105–11 (2008).
  * Rosta  _et al._ (2011) E. Rosta, M. Nowotny, W. Yang, and G. Hummer, “Catalytic mechanism of rna backbone cleavage by ribonuclease h from quantum mechanics/molecular mechanics simulations,” Journal of the American Chemical Society 133, 8934–8941 (2011).
  * Köfinger, Różycki, and Hummer (2019) J. Köfinger, B. Różycki, and G. Hummer, “Inferring Structural Ensembles of Flexible and Dynamic Macromolecules Using Bayesian, Maximum Entropy, and Minimal-Ensemble Refinement Methods,” Biomolecular Simulations: Methods and Protocols , 341–352 (2019).
  * Bottaro, Bengtsen, and Lindorff-Larsen (2020) S. Bottaro, T. Bengtsen, and K. Lindorff-Larsen, “Integrating Molecular Simulation and Experimental Data: a Bayesian/Maximum Entropy Reweighting Approach,” in  _Structural Bioinformatics_ (Springer, 2020) pp. 219–240.
  * Köfinger and Hummer (2023) J. Köfinger and G. Hummer, “Encoding prior knowledge in ensemble refinement,” (2023).
  * Gilardoni, Fröhlking, and Bussi (2023) I. Gilardoni, T. Fröhlking, and G. Bussi, “Boosting ensemble refinement with transferable force field corrections: synergistic optimization for molecular simulations,” arXiv preprint arXiv:2311.13532 (2023).
  * Cavender  _et al._ (2023) C. E. Cavender, D. A. Case, J. C.-H. Chen, L. T. Chong, D. A. Keedy, K. Lindorff-Larsen, D. L. Mobley, O. Ollila, C. Oostenbrink, P. Robustelli, V. A. Voelz, M. E. Wall, D. C. Wych, and M. K. Gilson, “Structure-based experimental datasets for benchmarking of protein simulation force fields,” arXiv preprint arXiv:2303.11056 (2023).

##  Appendix A Deriving the first and second derivatives of the BICePs score.

The BICePs score is defined as the negative logarithm of the Bayes factor,
which is the quotient of two model evidences:

$f(\epsilon)=-\ln\left(Z(\epsilon)\big{/}Z_{0}\right),$ (24)

where the evidence for model with force field parameters $\epsilon$ is

$Z(\epsilon)=\int\int\exp\left(-u(\mathbf{X},\mathbf{\sigma}\mid
D,\epsilon)\right)d\mathbf{X}d\mathbf{\sigma}$ (25)

and $Z_{0}$ is the model evidence for the reference distribution, which is set
to be a uniform ensemble (see Figure 1) and so, does not depend on $\epsilon$.
The partial derivatives of the BICePs score, $f(\epsilon)$ , with respect to
the model parameters $\epsilon_{i}$ , is:

$\frac{\partial
f(\epsilon)}{\partial\epsilon_{i}}=\frac{\partial}{\partial\epsilon_{i}}\left[-\ln\frac{Z(\epsilon)}{Z_{0}}\right]=-\frac{1}{Z}\frac{\partial
Z(\epsilon)}{\partial\epsilon_{i}},$ (26)

where

$\frac{\partial Z}{\partial\epsilon_{i}}=\int\int\left[-\frac{\partial
u}{\partial\epsilon_{i}}\right]\exp(-u)dXd\sigma$ (27)

This allows us to express the derivative of the BICePs score with respect to
the force field parameter $\epsilon_{i}$ as:

$\begin{split}\frac{\partial
f(\epsilon)}{\partial\epsilon_{i}}&=\int\int\frac{1}{Z(\epsilon)}\left[\frac{\partial
u}{\partial\epsilon_{i}}\right]\exp\left(-u\right)d\mathbf{X}d\mathbf{\sigma}\\\
\end{split}$ (28)

which is the difference of Boltzmann averaged values of $\partial
u/\partial\epsilon_{i}$. This can also be explained in terms of prior
energies. BICePs uses $\lambda$ to scale the prior energies, where the
reference ensemble is given $\lambda=0$. Therefore, the average energy of the
reference ensemble is zero and

$\begin{split}\frac{\partial
f(\epsilon)}{\partial\epsilon_{i}}=\bigg{\langle}\frac{\partial
u}{\partial\epsilon_{i}}\bigg{\rangle}\end{split}$ (29)

The Hessian matrix and corresponding second partial derivatives of the BICePs
score is required for computing uncertainties of the inferred parameters by
the covariance matrix (see equation 20), and when using a trust-region methods
for optimization (see methods like ‘trust-ncg‘, ‘trust-exact‘ in
‘scipy.optimize.minimize‘). Trust-region methods can be more robust for noisy
gradients; they may perform better with noisy objective functions because it
does not rely solely on the gradient information to determine the search
direction. Instead, it uses a combination of the gradient and second-order
information with the local quadratic model to estimate the direction of
steepest descent. The Hessian matrix contains elements of second-order partial
derivatives with respect to each parameter and gives the curvature of the
objective function.

$\displaystyle\mathbf{H}=\begin{bmatrix}\frac{\partial^{2}f}{\partial\epsilon_{1}^{2}}&\frac{\partial^{2}f}{\partial\epsilon_{1}\partial\epsilon_{2}}&\cdots&\frac{\partial^{2}f}{\partial\epsilon_{1}\partial\epsilon_{n}}\\\
\frac{\partial^{2}f}{\partial\epsilon_{2}\partial\epsilon_{1}}&\frac{\partial^{2}f}{\partial\epsilon_{2}^{2}}&\cdots&\frac{\partial^{2}f}{\partial\epsilon_{2}\partial\epsilon_{n}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{\partial^{2}f}{\partial\epsilon_{n}\partial\epsilon_{1}}&\frac{\partial^{2}f}{\partial\epsilon_{n}\partial\epsilon_{2}}&\cdots&\frac{\partial^{2}f}{\partial\epsilon_{n}^{2}}\end{bmatrix}$
(30)

The second-order partial derivative of the BICePs score for parameters $i$ and
$j$, where $i\neq j$ is given as

$\begin{split}\frac{\partial^{2}f}{\partial\epsilon_{i}\partial\epsilon_{j}}&=\left[-\frac{1}{Z(\epsilon)}\frac{\partial^{2}Z(\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}+\frac{1}{Z(\epsilon)^{2}}\left(\frac{\partial
Z(\epsilon)}{\partial\epsilon_{i}}\right)\left(\frac{\partial
Z(\epsilon)}{\partial\epsilon_{j}}\right)\right],\end{split}$ (31)

where

$\begin{split}\frac{\partial^{2}Z}{\partial\epsilon_{i}\partial\epsilon_{j}}&=\frac{\partial}{\partial\epsilon_{j}}\left[\iint\left(-\frac{\partial
u}{\partial\epsilon_{i}}\right)\exp(-u)dXd\sigma\right]\\\
&=\iint\left[\left(\frac{\partial
u}{\partial\epsilon_{i}}\right)\left(\frac{\partial
u}{\partial\epsilon_{j}}\right)-\frac{\partial^{2}u}{\partial\epsilon_{i}\partial\epsilon_{j}}\right]\exp\left(-u\right)dXd\sigma\end{split}$
(32)

Putting this altogether and solving for equation 31 gives

$\begin{split}\frac{\partial^{2}f(\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}=&-\iint\frac{1}{Z(\epsilon)}\left[\left(\frac{\partial
u}{\partial\epsilon_{i}}\right)\left(\frac{\partial
u}{\partial\epsilon_{j}}\right)\right]\exp\left(-u\right)dXd\sigma\\\
&+\iint\frac{1}{Z(\epsilon)}\left[\frac{\partial^{2}u}{\partial\epsilon_{i}\partial\epsilon_{j}}\right]\exp\left(-u\right)dXd\sigma\\\
&+\iint\frac{1}{Z(\epsilon)^{2}}\left(\frac{\partial
u}{\partial\epsilon_{i}}\right)\left(\frac{\partial
u}{\partial\epsilon_{j}}\right)\exp\left(-2u\right)dXd\sigma\\\ \end{split}$
(33)

which reduces to

$\displaystyle\frac{\partial^{2}f(\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}=\bigg{\langle}\frac{\partial^{2}u}{\partial\epsilon_{i}\partial\epsilon_{j}}\bigg{\rangle}-\biggl{(}$
$\displaystyle\bigg{\langle}\left(\frac{\partial
u}{\partial\epsilon_{i}}\right)\left(\frac{\partial
u}{\partial\epsilon_{j}}\right)\bigg{\rangle}$ (34)
$\displaystyle-\bigg{\langle}\left(\frac{\partial
u}{\partial\epsilon_{i}}\right)\left(\frac{\partial
u}{\partial\epsilon_{j}}\right)\bigg{\rangle}^{\prime}\biggr{)},$

where the notation $\langle\cdot\rangle^{\prime}$ denotes the average with
respect to $\left(\frac{1}{Z}\exp(-u)\right)^{2}$. When $i=j$ (diagonal
elements of the Hessian), the second-order partial derivatives of the BICePs
score can be seen as the difference of the Boltzmann-averaged second-order
partial derivative of the energy and the variance in its first-order partial
derivative:

$\begin{split}\frac{\partial^{2}f(\epsilon)}{\partial\epsilon^{2}}=\bigg{\langle}\frac{\partial^{2}u}{\partial\epsilon^{2}}\bigg{\rangle}-\left(\bigg{\langle}\left(\frac{\partial
u}{\partial\epsilon}\right)^{2}\bigg{\rangle}-\bigg{\langle}\frac{\partial
u}{\partial\epsilon}\bigg{\rangle}^{2}\right)\end{split}$ (35)

To compute the derivatives of the BICePs score, we need the first and second
derivatives of the BICePs energy function $u$, which consists of taking
derivatives of the negative logarithm of the posterior, since $u=-\log
p(X,\bm{\sigma}|D,\epsilon)$. The derivative of this function with respect to
some force field parameter $\epsilon$ is

$\begin{split}\frac{\partial u}{\partial\epsilon_{i}}&=\frac{\partial
E}{\partial\epsilon_{i}}\left(\sum^{N}_{r=1}-\log\left[\frac{\exp(-E(X_{r}|\epsilon))}{\sum_{X}\exp\left(-E(X|\epsilon)\right)}\right]\right)\\\
&=-\sum^{N}_{r=1}\frac{\sum_{X}\exp\left(-E(X|\epsilon)\right)}{\exp\left(-E(X|\epsilon)\right)}\frac{\partial}{\partial\epsilon_{i}}\left[-\frac{\frac{\partial
E(X|\epsilon)}{\partial\epsilon_{i}}\exp(-E(X|\epsilon))}{\sum_{X}\exp\left(-E(X|\epsilon)\right)}\right.\\\
&\left.\quad\quad-\frac{\sum_{X}-\exp\left(-E(X|\epsilon)\right)\frac{\partial
E(X|\epsilon)}{\partial\epsilon_{i}}\exp(-E(X|\epsilon))}{\left(\sum_{X}\exp\left(-E(X|\epsilon)\right)\right)^{2}}\right]\\\
&=\sum^{N}_{r=1}\biggl{\\{}\frac{\partial
E(X_{r}|\epsilon)}{\partial\epsilon_{i}}-\frac{\sum_{X}\exp\left(-E(X|\epsilon)\right)\frac{\partial
E(X|\epsilon)}{\partial\epsilon_{i}}}{\sum_{X}\exp\left(-E(X|\epsilon)\right)}\biggr{\\}}\\\
&=\sum^{N}_{r=1}\biggl{\\{}\frac{\partial
E(X_{r}|\epsilon)}{\partial\epsilon_{i}}-\big{\langle}\frac{\partial
E(X|\epsilon)}{\partial\epsilon_{i}}\big{\rangle}\biggr{\\}}\end{split}$ (36)

The second derivative is

$\begin{split}\frac{\partial^{2}u}{\partial\epsilon_{i}\partial\epsilon_{j}}&=\sum^{N}_{r=1}\left[\frac{\partial^{2}E(X_{r}|\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}-\frac{\sum_{X}\exp\left(-E(X|\epsilon)\right)\frac{\partial^{2}E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}}{\sum_{X}\exp\left(-E(X|\epsilon)\right)}\right.\\\
&\left.+\frac{\sum_{X}\exp\left(-E(X|\epsilon)\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}}\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{j}}\right)}{\sum_{X}\exp\left(-E(X|\epsilon)\right)}\right.\\\
&\left.-\frac{\sum_{X}\exp\left(-2E(X|\epsilon)\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}}\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{j}}\right)}{\left(\sum_{X}\exp\left(-E(X|\epsilon)\right)\right)^{2}}\right]\\\
\end{split}$ (37)

which can be reduced to

$\begin{split}\frac{\partial^{2}u}{\partial\epsilon_{i}\partial\epsilon_{j}}=\sum^{N}_{r=1}\biggl{\\{}&\frac{\partial^{2}E(X_{r}|\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}-\bigg{\langle}\frac{\partial^{2}E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}\partial\epsilon_{j}}\bigg{\rangle}\\\
&+\bigg{\langle}\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}}\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{j}}\right)\bigg{\rangle}\\\
&-\frac{\sum_{X}\exp\left(-2E(X|\epsilon)\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{i}}\right)\left(\frac{\partial
E(\mathbf{X}|\epsilon)}{\partial\epsilon_{j}}\right)}{\left(\sum_{X}\exp\left(-E(X|\epsilon)\right)\right)^{2}}\biggr{\\}}\end{split}$
(38)

and when $i=j$, the second-order partial derivative of the energy function
becomes

$\begin{split}\frac{\partial^{2}u}{\partial\epsilon^{2}}=\sum^{N}_{r=1}\biggl{\\{}&\frac{\partial^{2}E(X_{r}|\epsilon)}{\partial\epsilon^{2}}-\bigg{\langle}\frac{\partial^{2}E(X|\epsilon)}{\partial\epsilon^{2}}\bigg{\rangle}\\\
&+\bigg{\langle}\left(\frac{\partial
E(X|\epsilon)}{\partial\epsilon}\right)^{2}\bigg{\rangle}-\bigg{\langle}\frac{\partial
E(X|\epsilon)}{\partial\epsilon}\bigg{\rangle}^{2}\biggr{\\}}\end{split}$ (39)

##  Appendix B Kofinger & Hummer’s toy polymer model

As mentioned in the main text, the two dimensional polymer modelKöfinger and
Hummer (2021) consists of a string of $n=100$ beads separated by a unit of
length, where $\Delta\phi_{i}$ is the angle between the two bond vectors
$v_{i}$ and $v_{i+1}$. The associated probability of a polymer with $n$ beads
is given by

$\begin{split}p(\bm{\Delta}\bm{\phi}\mid\bm{\mu},\bm{\kappa})&=\prod_{i=1}^{n-1}p\left(\Delta\phi_{i}\mid\mu_{i},\kappa_{i}\right)\\\
&=\prod_{i=1}^{n-1}\frac{\exp(\kappa_{i}\cos\left(\Delta\phi_{i}-\mu_{i}\right))}{2\pi
I_{0}\left(\kappa_{i}\right)}\end{split}$ (40)

where $I_{0}(\kappa)$ is the modified Bessel function of order zero and
$\mu_{i}$ introduces directional biases. The parameters $\kappa_{i}>0$
determine the stiffness of the polymer at positions $i$. Prior energy, $E$ of
the polymer configuration follow the Boltzmann distribution, which corresponds
to

$E=-\log\prod_{i=1}^{n-1}p\left(\Delta\phi_{i}\mid\mu_{i},\kappa_{i}\right)$
(41)

The authors did not provide derivatives of the toy model energy function, so
we will briefly derive derivatives of the prior energies that will be used in
our algorithm. The first derivative of the prior energy $E$ with respect to
$\kappa$ is

$\begin{split}\frac{\partial
E}{\partial\kappa}&=-\sum_{i=1}^{n-1}\frac{\partial}{\partial\kappa}\log
p\left(\Delta\phi_{i}\mid\mu_{i},\kappa_{i}\right)\\\
&=-\sum_{i=1}^{n-1}\frac{1}{p\left(\Delta\phi_{i}\mid\mu_{i},\kappa_{i}\right)}\frac{\partial
p\left(\Delta\phi_{i}\mid\mu_{i},\kappa_{i}\right)}{\partial\kappa},\end{split}$
(42)

where the derivative of the von Mises distribution with respect to $\kappa$ is

$\begin{split}\frac{\partial
p}{\partial\kappa}=&\frac{\exp\left(\kappa_{i}\cos(\mu_{i}-\Delta\phi_{i})\right)\cos(\mu_{i}-\Delta\phi_{i})}{2\pi
I_{0}\left(\kappa_{i}\right)}\\\
&\quad-\frac{\exp\left(\kappa_{i}\cos(\mu_{i}-\Delta\phi_{i})\right)I_{1}\left(\kappa_{i}\right)}{2\pi{I_{0}}^{2}\left(\kappa_{i}\right)}\end{split}$
(43)

which is relatively straightforward since the derivative of the modified
Bessel function merely increases it’s order. Next, the second derivative of
$E$ with respect to $\kappa_{i}$:

$\frac{\partial^{2}E}{\partial\kappa_{i}^{2}}=\sum\left\\{-\frac{1}{p}\frac{\partial^{2}p}{\partial\kappa_{i}^{2}}+\left(\frac{1}{p}\frac{\partial
p}{\partial\kappa_{i}}\right)^{2}\right\\},$ (44)

where the second derivative of the von Mises distribution w.r.t $\kappa_{i}$
is

$\displaystyle\frac{\partial^{2}p}{\partial\kappa_{i}^{2}}=\biggl{[}$
$\displaystyle\frac{\exp\left(\kappa_{i}\cos(\mu_{i}-\Delta\phi_{i})\right)\left(\frac{I_{0}\left(\kappa_{i}\right)}{2}+\frac{I_{2}\left(\kappa_{i}\right)}{2}\right)}{2\pi{I_{0}}^{2}\left(\kappa_{i}\right)}$
(45)
$\displaystyle+\frac{\exp\left(\kappa_{i}\cos(\mu_{i}-\Delta\phi_{i})\right)\cos^{2}(\mu_{i}-\Delta\phi_{i})}{2\pi
I_{0}\left(\kappa_{i}\right)}$
$\displaystyle-\frac{\exp\left(\kappa_{i}\cos(\mu_{i}-\Delta\phi_{i})\right)\cos(\mu_{i}-\Delta\phi_{i})I_{1}\left(\kappa_{i}\right)}{\pi{I_{0}}^{2}\left(\kappa_{i}\right)}$
$\displaystyle+\frac{\exp\left(\kappa_{i}\cos(\mu_{i}-\Delta\phi_{i})\right){I_{1}}^{2}\left(\kappa_{i}\right)}{\pi{I_{0}}^{3}\left(\kappa_{i}\right)}\biggr{]}$

##  Appendix C Supplemental Information

Figure S1: The probability density function (a) and energy landscape (b) of
the marginal likelihood for the Student’s model with respect to the replica-
averaged forward model data $f(\mathbf{X})$. The colored curves are different
values of nuisance parameter $\beta$. The Student’s model is equivalent to the
Outliers model (Cauchy) when $\beta=1$, and as $\beta$ goes to infinity, the
pdf becomes Gaussian and the energy becomes harmonic. Figure S2: The
probability density function (a) computed as
$p(D|\mathbf{X},\sigma_{0},\beta)=\prod^{N_{j}}_{j}p(d_{j}|\mathbf{X},\sigma_{0},\beta)$
and energy landscape (b) of the marginal likelihood for the Student’s model
with respect to the replica-averaged forward model data $f(\mathbf{X})$ using
multiple data points. Shown here, are three data points, two good data points
{0.25, 0.3} and one outlier {1.0}. The Cauchy likelihood ($\beta=1$) is
centered about the mean of the two good data points, demonstrating that this
model can distinguish the good and bad data. The standard Gaussian likelihood
($\beta=300$) is centered about the mean of all three data points. The colored
curves are different values of nuisance parameter $\beta$.  Figure S3: Trust-
region optimization with BICePs is found to converge faster and be more robust
than L-BFGS-B. Here, we compare first-order L-BFGS-B method (blue) to the
second-order Trust-NCG optimization method (red). A total of 25 independent
rounds of parameter ($\epsilon$) optimizations using the BICePs score for a
maximum of ten iterations. Data points are placed at the mean value for each
iteration and error bars are the standard deviation. It is typical that the
trust-ncg method reaches convergence slightly faster (the 4th iteration), and
the L-BFGS-B method tends to slightly overshoot the "True" $\epsilon$ value on
the first iteration. Figure S4: The BICePs score landscape using a standard
Gaussian likelihood versus a Student’s likelihood in the presence of
systematic error. (a) 1-D scans over $\epsilon$ to reveal the landscape of the
BICePs score (scatter dots). The green tangent lines are the derivatives at
each $\epsilon$ value. Uncertainties in the BICePs scores and derivatives come
from five independent scans along $\epsilon$. (b) The derivative of the BICePs
score (scatter dots) and the second derivative of the BICePs score (yellow
lines) at each epsilon value. The dotted black line at $\epsilon^{*}=1.0$
shows the true value, which is where the derivative of the BICePs score equals
zero. BICePs calculations are run using the Student’s model for 100k steps
with 8 replicas. Systematic error (+4Å shift) was induced to the 2-11
distance, with a total experimental error of $\sigma_{data}=1.41$. Figure S5:
Average traces over a total of 25 independent rounds of parameter
$(\epsilon_{2},\epsilon_{4})$ optimizations using first-order optimization
method (L-BFGS-B) with BICePs, for a maximum of ten iterations. Traces are
unable to fully reach the "True" iteration strength parameters
($\epsilon_{2}^{*}=1.25$, $\epsilon_{4}^{*}=1.5$) within 10 iterations when
starting from different initial parameters
$(\epsilon_{2}^{0},\epsilon_{4}^{0})$ = $\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$.
The BICePs score landscape was generated from the average values of five scans
over $\epsilon_{2}$ and $\epsilon_{4}$. All calculations used the Student’s
model with 200k MCMC steps and 8 replicas.  Figure S6: Average traces over a
total of 25 independent rounds of parameter $(\epsilon_{2},\epsilon_{4})$
optimizations using second-order (trust-ncg) method with BICePs, for a maximum
of ten iterations. Optimizations converge to the "True" optimal iteration
strength parameters ($\epsilon_{2}^{*}=1.25$, $\epsilon_{4}^{*}=1.5$) when
starting from different initial parameters
$(\epsilon_{2}^{0},\epsilon_{4}^{0})$ = $\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$.
The BICePs score landscape was generated from the average values of five scans
over $\epsilon_{2}$ and $\epsilon_{4}$. All calculations used the Student’s
model with 200k MCMC steps and 8 replicas.  Figure S7: Average traces over a
total of 25 independent rounds of parameter $(\epsilon_{2},\epsilon_{4})$
optimizations using second-order (trust-ncg) method with BICePs, for a maximum
of ten iterations. In these tests, no error was added to the experimental
data. Optimizations converge to the "True" parameters
($\epsilon_{2}^{*}=1.25$, $\epsilon_{4}^{*}=1.5$) when starting from different
initial parameters $(\epsilon_{2}^{0},\epsilon_{4}^{0})$ =
$\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$. The BICePs score landscape was
generated from the average values of five scans over $\epsilon_{2}$ and
$\epsilon_{4}$. All calculations used the Student’s model with 200k MCMC steps
and 32 replicas. Average optimized parameter values were ($\epsilon_{2}$ =
1.12 $\pm$ 0.39, and $\epsilon_{4}$ = 1.43 $\pm$ 0.34), where the
uncertainties are estimated from inverse Hessian.  Figure S8: Accuracy
profiles for $(\epsilon_{2},\epsilon_{4})$ optimizations using second-order
(trust-ncg) method with BICePs, for a total of 25 independent rounds and a
maximum of ten iterations. In these tests, no error was added to the
experimental data. Optimizations converge to the "True" parameters
($\epsilon_{2}^{*}=1.25$, $\epsilon_{4}^{*}=1.5$) with a relative error less
than 1.0 when starting from different initial parameters
$(\epsilon_{2}^{0},\epsilon_{4}^{0})$ = $\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$.
All calculations used the Student’s model with 200k MCMC steps and 32
replicas. Average optimized parameter values were ($\epsilon_{2}$ = 1.12 $\pm$
0.39, and $\epsilon_{4}$ = 1.43 $\pm$ 0.34), where the uncertainties are
estimated from inverse Hessian. Error bars represent uncertainties estimated
from the inverse Hessian, averaged over the 25 independent optimizations.
Figure S9: Accuracy profiles for $(\epsilon_{2},\epsilon_{4})$ optimizations
using second-order (trust-ncg) method with BICePs, for a total of 25
independent rounds and a maximum of ten iterations. In these tests, the
experimental data was corrupted with systematic error in the 2–11 and 4–9
distances for +3 and +3.5Å shift, respectively. The total error in the data is
$\sigma_{data}=1.63$. Optimizations converge to the "True" optimal iteration
strength parameters ($\epsilon_{2}^{*}=1.25$, $\epsilon_{4}^{*}=1.5$) with a
relative error less than 1.0 when starting from different initial parameters
$(\epsilon_{2}^{0},\epsilon_{4}^{0})$ = $\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$.
Error bars represent uncertainties estimated from the inverse Hessian,
averaged over the 25 independent optimizations. All calculations used the
Student’s model with 200k MCMC steps and 32 replicas.  Figure S10: Accuracy
profiles for $(\epsilon_{2},\epsilon_{4})$ optimizations using a second-order
(trust-ncg) method with BICePs, for a total of 25 independent rounds and a
maximum of ten iterations. In these tests, the experimental data was corrupted
with systematic error in the 2–11 and 4–9 distances for +3 and +3.5Å shift,
respectively. The total error in the data is $\sigma_{data}=1.63$.
Optimizations converge to the "True" parameters ($\epsilon_{2}^{*}=1.25$,
$\epsilon_{4}^{*}=1.5$) with a relative error less than 1.0 when starting from
different initial parameters $(\epsilon_{2}^{0},\epsilon_{4}^{0})$ =
$\\{(0.5,5.0),(4.0,5.0),(4.0,0.5)\\}$. Error bars represent uncertainties
estimated from the inverse Hessian, averaged over the 25 independent
optimizations. For the optimization traces that remain unconverged for 10
iterations, the relative error was determined to be $0.24\pm 0.75$, where the
uncertainty is estimated from the inverse Hessian. All calculations used the
Student’s model with 200k MCMC steps and 16 replicas.  Figure S11: The
accuracy of a standard Gaussian likelihood drops in the presence of systematic
error. Demonstrating optimization accuracy when the data is corrupted with and
without systematic error using $\epsilon_{2}=4.0$ and $\epsilon_{4}=5.0$ as
initial parameters, same conditions as in Figure 5, where the "True"
parameters are $(\epsilon_{2}^{*}=1.25,\epsilon_{4}^{*}=1.5)$. (a) The
Gaussian likelihood performs great when using experimental data not affected
by systematic error. (b) The accuracy suffers when the data is corrupted with
systematic error in the 2-11 and 4-9 distances for +3 and +3.5Å shift,
respectively. The total error in the data is $\sigma_{data}=1.63$. All results
shown here are averages over 25 independent optimizations, and error bars are
computed as standard error of the mean for each iteration. Figure S12: The
accuracy of the Student’s model is very robust in the presence of systematic
error. (a) Average traces over 25 independent optimizations when the data is
corrupted with systematic error in the 2–11 and 4–9 distances for +3 and +3.5Å
shift, respectively. The total error in the data is $\sigma_{data}=1.63$. (b)
A single representative trace demonstrating convergence to the "True"
parameters $(\epsilon_{2}^{*}=1.25,\epsilon_{4}^{*}=1.5)$ Figure S13: Average
traces over a total of 25 independent rounds of parameter
$(\epsilon_{2},\epsilon_{4},\epsilon_{6})$ optimizations using second-order
(trust-ncg) method with BICePs, for a maximum of ten iterations. Optimizations
converge to the "True" optimal iteration strength parameters
($\epsilon_{2}^{*}=1.25$, $\epsilon_{4}^{*}=1.5$, $\epsilon_{6}^{*}=1.5$) when
starting from different initial parameters. All calculations used the
Student’s model with 200k MCMC steps and 32 replicas. Figure S14: Sensitivity
analysis of parameter $(\epsilon_{2},\epsilon_{4},\epsilon_{6})$ optimizations
from various initial conditions. Regardless of different starting parameters,
optimizations universally converge to "True" optimal iteration strength
parameters $(\epsilon_{2}^{*}=1.25$, $\epsilon_{4}^{*}=1.5$,
$\epsilon_{6}^{*}=1.5$) and succeed in diminishing the relative error. All
BICePs score optimizations used 32 BICePs replicas with the Student’s
likelilhood model, and 200k MCMC steps per iteration. The dashed line is a
linear fitting with equation $y=0.36x-0.32$, resulting in a strong correlation
between the BICePs score and relative model error ($R^{2}=0.96$).  Figure S15:
Minimization of the BICePs score, $f$ is demonstrated when averaging over 25
independent rounds of parameter $(\epsilon_{2},\epsilon_{4},\epsilon_{6})$
optimizations using second-order (trust-ncg) method, for a maximum of ten
iterations. Error bars represent the standard error over the 25 optimizations.
Figure S16: Refinement of the $\kappa$, the stiffness parameter in the KH
polymer model. An initial parameter value, $\kappa^{0}=20$, was established.
Employing second-order optimization via the ‘trust-ncg‘ method, convergence to
the “true” parameter, $\kappa^{*}=10$, was achieved within a few iterations.
The computations used 200 replicas and 100k MCMC steps, encompassing 500
configurations. Parameter uncertainties throughout each iteration were
estimated from the inverse Hessian. These results successfully replicate the
outcome of BioFF by Kofinger and Hummer.Köfinger and Hummer (2021).