# Adding alchemical variables to metadynamics to enhance sampling in free
energy calculations

Wei-Tse Hsu  Pascal T. Merz  Department of Chemical and Biological
Engineering, University of Colorado Boulder, Boulder, CO 80305  Giovanni Bussi
Scuola Internazionale Superiore di Studi Avanzati, via Bonomea 265, 34136
Trieste, Italy  Michael R. Shirts  Department of Chemical and Biological
Engineering, University of Colorado Boulder, Boulder, CO 80305
<EMAIL_ADDRESS>
###### Abstract

Performing alchemical transformations, in which one molecular system is
nonphysically changed to another, is a popular approach adopted in free energy
calculations associated with various biophysical processes, such as protein-
ligand binding, or the transfer of a molecule between environments. While the
sampling of alchemical intermediate states in either parallel (e.g.
Hamiltonian replica exchange) or serial manner (e.g. expanded ensemble) can
bridge the high-probability regions in the configurational space between two
end states of interest, alchemical methods can fail in the scenarios where the
most important slow degrees of freedom in the configurational space are almost
orthogonal to the alchemical variable, or if the system gets trapped in a deep
basin expanding in both the configurational and alchemical space.

To alleviate these issues, we propose to use alchemical variables as an
additional dimension in metadynamics, augmenting the ability both to sample
collective variables and to enhance sampling in free energy calculations. In
this study, we validate our implementation of “alchemical metadynamics” in
PLUMED with test systems and alchemical processes with varying complexities
and dimensions of collective variable space, including the binding and
unbinding of CB8-G3 host-guest binding complex. We show that multi-dimensional
alchemical metadynamics can address the challenges mentioned above challenges
of traditional alchemical free energy methods and further accelerate sampling
by introducing configurational collective variables. The method can trivially
be combined with other metadynamics-based algorithms implemented in PLUMED.
The necessary PLUMED code changes have already been released for general use
in PLUMED 2.8.

###### keywords:

American Chemical Society, LaTeX

##  1 Introduction

With the fast advent of high-performance computing (HPC) and parallel
computing over the years, molecular simulations such as molecular dynamics
(MD) and Monte Carlo (MC) simulations have become increasingly useful in
elucidating the dynamics of transformation processes of condensed matter
systems. They are most useful when the system can visit all the energetically
relevant conformations, in which case we can extract valuable thermodynamic
and structural information about the system, such as the solvation free energy
of a molecule, or the binding ensemble of a binding complex. However,
comprehensive sampling in the phase space is generally challenging in
traditional MD simulations because the system must rely on very rare
fluctuations to cross the free energy barriers that separate metastable states
of interest. This low transition probability between metastable states makes
the time scale required for achieving system ergodicity in unbiased sampling
impractically long in most systems of interest.

To address this challenge, a wide variety of advanced sampling methods have
been developed in the past decades 1. One particular flavor of advanced
sampling methods involves sampling along a set of predefined coarse-grained
descriptors of the system, or collective variables (CVs). Traditionally, CVs
could be any function of the atomic coordinates of the system, but the optimal
ones should correspond to the slowest degrees of freedoms that distinguish
different metastable states. Methods relying on the use of CVs include
umbrella sampling 2, adaptive biasing force 3, metadynamics 4 and their
variations 5, 6, 7.

Another category of advanced sampling methods is known as generalized ensemble
algorithms, which includes simulated tempering 8, replica exchange 9, 10,
expanded ensemble 11, $\lambda$-dynamics 12, and their variations 13, 14.
These methods do not require any predefined CVs, but a series of intermediate
or auxiliary states with modified probability of the coordinates of the
system. These intermediates are chosen to remove or lessen the kinetic
barriers between end states of interest. These states are defined by varying
the temperature or the Hamiltonian of the system to allow an additional
sampling dimension, which created nonphysical pathways for the system to
circumvent free energy barriers in the configurational space. Increasing the
probability overlap between adjacent intermediate states can often enhance the
diffusion in not only the temperature/Hamiltonian direction, but also the
configurational space. In replica exchange, these states are sampled with the
ensembles progressing forwards in time in parallel, while in serial
approaches, a single simulations can move between states in either discrete
(expanded ensemble) or continuous space ($\lambda$-dynamics). In order to have
even sampling between the states, one must add biases or weights to the higher
free energy states so they can be sampled. These weights, which are absent in
replica exchange, similarly modify the underlying free energy surface as the
biasing potentials in metadynamics do. Notably, these intermediate or
auxiliary states can be useful for reasons other than simply accelerating
sampling. For example, the sampling in the temperature space allows us to
determine thermodynamic observables of interest as a function of temperature,
while simulations with alchemical intermediate states are useful for
calculating the free energy difference between the end states of alchemical
processes. In free energy calculations where the temperature is constant,
sampling in alchemical intermediate states obviates the need of defining CVs,
which could be non-trivial in systems where the slowest-relaxing coordinates
are not intuitive, such as the escape of a ligand from a buried binding pocket
15, 16.

Although sampling in the alchemical variable can create new ensembles where
the slowest physical collective variables are not longer so slow, it cannot
necessarily enhance sampling where the configurational barriers are almost
orthogonal to the alchemical direction. For example, a configurational free
energy barrier can be present for all the alchemical states (scenario A in
Figure 1) so that the system could remain stuck even if it is able to drift to
other alchemical states or cross free energy barriers along the alchemical
direction. Another scenario that could easily trap the system is the presence
of large free energy basins in both the configurational and alchemical
directions (scenario B in Figure 1). In both scenarios, state-dependent
weights can only shift the whole free energy profile as a function of relevant
configurational CVs with the same amount. They do not change the partition
function ratio between points at the same alchemical state but with different
values of the configurational CVs, so the sampling in the configurational
space can still be limited in these cases.

Figure 1: Scenarios where pure alchemical bias fails to accelerate
configurational sampling. In scenario A, the free energy barrier extending
across all $\lambda$ states prevents the system from sampling both metastable
states at $\lambda$ being 0. In scenario B, the free energy basins in both the
alchemical and configurational direction can trap the system during the
adaptive build-up of weights.

Currently, the challenges mentioned above can be addressed to varying extents.
For example, the alchemical flying Gaussian method 17 or the conveyor belt
method combined with specific biasing 18 does not bias the alchemical but only
configurational space in alchemical sampling. Although these methods might
have difficulties in getting the system out of free energy basins in scenario
B in Figure 1, they are sufficient to overcome the extended configurational
free energy barrier in scenario A in Figure 1. Methods such as $\lambda$ local
elevation umbrella sampling ($\lambda$-LEUS) 19, 20, orthogonal space random
walk (OSRW) 21, double-integration orthogonal space tempering (DI-OST) 22, and
adaptive landscape flattening (ALF) 23, which all work with continuous
alchemical space in their proposals, are able to apply biasing potentials in
both the alchemical or configurational directions. Theoretically, they are
able to facilitate the escape of the system from the deep free energy basins
shown in scenario B in Figure 1, but each of them was implemented in different
contexts and algorithmic designs and the allowed forms of the configurational
bias tend to be specific rather than general.

In light of the need for a more generalized approach that can address the
issues in both scenarios in Figure 1, especially in cases where
multidimensional biases in the configurational space are needed, we propose to
use alchemical variables as an additional dimension in metadynamics. Although
it is not a collective variable of coordinates whose values divide the
coordinate space into disjoint sets, such as a center of mass distance or
radius of gyration, it has an associated free energy and dynamics, and thus
can fit into the same formalism. We term this approach alchemical metadynamics
and have implemented it in PLUMED 2.8 24 (initially, only for GROMACS). Given
the well-developed library and flexible syntax that PLUMED has for defining
configurational CVs and various types of restraints, the implementation of
alchemical metadynamics within PLUMED is particularly useful for smoothly
flattening highly-dimensional free energy landscapes in a more general way
compared to methods such as $\lambda$-LEUS, OSRW, or DI-OST. To demonstrate
the usage of such an algorithm and validate its implementation in PLUMED, we
employed this method to estimate the free energy difference of different
alchemical processes as elaborated in the Methods section, from decoupling an
argon atom or a molecule composed 4 interaction sites, to the binding and
unbinding between a guest molecule and a host molecule in solution.

##  2 Theory

###  2.1 Metadynamics

As one of the most popular CV-based advanced sampling methods, metadynamics 4
accelerates the sampling by depositing Gaussian biasing potentials to the
underlying free energy surface of the system. The biasing potential is a
function of the vector of collective variables of interest $\xi$, which can be
regarded as a reduced dimensional space calculated from a configuration $x$ by
the mapping $\Phi({\mbox{\boldmath{$x$}}})={\mbox{\boldmath{$\xi$}}}$. During
the simulation, biasing potentials are deposited to seek roughly equal
sampling across the reduced dimensional space of interest. Let the CV vector
$\xi$ be $d$-dimensional, i.e.
${\mbox{\boldmath{$\xi$}}}=(\xi_{1}({\mbox{\boldmath{$x$}}}),\xi_{2}({\mbox{\boldmath{$x$}}}),...,\xi_{d}({\mbox{\boldmath{$x$}}}))$.
The total biasing potential added after a period of time $t$ can be expressed
as

$V({\mbox{\boldmath{$\xi$}}},t)=W\sum_{t^{\prime}=k\tau,k\in
N}^{t^{\prime}<t}\exp\left(-\sum_{i=1}^{d}\frac{\left(\xi_{i}-\xi_{i}({\mbox{\boldmath{$x$}}}(k\tau))\right)^{2}}{2\sigma_{i}^{2}}\right)$
(1)

where $W$ is the height of the Gaussian, $k$ is the number of Gaussian
depositions, $\tau$ is the deposition stride and $\sigma_{i}$ is the width of
the Gaussian along the $i$-th dimension. Notably, the Gaussian height $W$ can
be either constant (in standard metadynamics) or time-dependent (in well-
tempered metadynamics 25) during the course of the simulation, with the latter
more commonly adopted for a smoother convergence and better concentration on
the physically relevant regions of the configurational space. Specifically,
the time-dependent Gaussian height $W(k\tau)$ can be written as:

$W(k\tau)=W_{0}\exp\left(-\frac{V(\vec{\xi}({\mbox{\boldmath{$x$}}}(k\tau)),k\tau)}{k_{B}\Delta
T}\right)$ (2)

where $W_{0}$ is the initial Gaussian height and $\Delta T$ is a temperature
parameter that incorporates a user-defined bias factor $\gamma=(T+\Delta T)/T$
for adjusting the decay rate of the bias. In well-tempered metadynamics, the
free energy surface as a function of the multi-dimensional CV $\xi$ can be
estimated by the following relationship:

$V({\mbox{\boldmath{$\xi$}}},t\rightarrow\infty)=-\frac{\Delta T}{T+\Delta
T}F({\mbox{\boldmath{$\xi$}}})=-(1-\frac{1}{\gamma})F({\mbox{\boldmath{$\xi$}}})$
(3)

To efficiently obtain a reasonable estimate of the free energy difference of
interest, the dimensions of the chosen set of CVs must be as low as possible
while still capturing the slowest degrees of freedom of the system, as the
space to be explored increases exponentially with the number of CVs, leading
to prohibitive time to converge weights.

###  2.2 Alchemical metadynamics

In alchemical metadynamics, the alchemical variable $\lambda$ is introduced in
the generalized CV vector
${\mbox{\boldmath{$\xi^{\prime}$}}}=(\lambda,\xi_{1}({\mbox{\boldmath{$x$}}}),\xi_{2}({\mbox{\boldmath{$x$}}}),...,\xi_{d}({\mbox{\boldmath{$x$}}}))$
such that the joint space of $\lambda$ and $\xi$ is sampled with the aid of
the biasing potential $V({\mbox{\boldmath{$\xi^{\prime}$}}})$. Unlike the
configurationally defined CVs, the alchemical variable is not a function of
atomic coordinates. In the current implementation, we assume that the
alchemical variable takes discrete values, i.e. the state index that can be
mapped to a vector of coupling parameters for decoupling/switching different
interactions, such as van der Waals interactions, electrostatic interactions,
or any kinds of restraints. Just like expanded ensemble, alchemical
metadynamics alternates the sampling along the alchemical direction and the
coordinate direction. The sampling in the discrete alchemical space can be
done by Monte Carlo sampling just as the alchemical sampling in expanded
ensemble, while the coordinate direction is sampled by molecular dynamics just
as in any other type of metadynamics. Currently, our implementation of
alchemical metadynamics is available in the combination of PLUMED 2.8
interfaced with GROMACS 2020, or a more recent versions of each. When using
alchemical metadynamics, the state index of the alchemical or coupling
parameter $\lambda$ is passed from GROMACS into PLUMED along with the system
configuration required to calculate configurational CVs. PLUMED uses these
alchemical indices and any other CV present to track the visited states of the
system and calculate the metadynamics bias. GROMACS updates the alchemical
state via MC. When calculating the energy of the current and the candidate
$\lambda$ state, GROMACS adds the metadynamics bias provided by PLUMED. This
approach is compatible with all MC schemes in alchemical space offered by
GROMACS (Metropolis-Hastings algorithm 26, Barker transition method 27, Gibbs
sampling 28, 29, and Metropolized-Gibbs sampling 30, 31).

Theoretically, one-dimensional alchemical metadynamics, which does not apply
configurational but only alchemical bias, is effectively equivalent to
expanded ensemble with a different weight updating procedure for allowing
roughly equal sampling across alchemical states. For example, in an expanded
ensemble where the Wang-Landau algorithm 32, 33, 34 is used for weight
updating, the reduced potential of the system is incremented by a Wang-Landau
incrementor whenever a move across alchemical states is attempted. This is
analogous to the periodic deposition of Gaussian potentials in 1D alchemical
metadynamics, especially when the Gaussian deposition stride is the same as
the number of integration steps between attempted moves in the alchemical
space. For converging the free energy surface in the alchemical space, both
algorithms have mechanisms for decreasing the bias along the course of the
simulation. In expanded ensemble, the Wang-Landau incrementor is modified by a
scaling factor whenever the state visitation reaches a certain flatness
criterion. In 1D well-tempered alchemical metadynamics, a bias factor is
applied to enforce a continuous exponential decay of the Gaussian height,
which would lead to marginally smoother convergence compared to expanded
ensemble or other similar free energy methods. There are also alternative
strategies proposed to this aim. 35.

In multidimensional alchemical metadynamics, introducing additional
configurational CVs could further enhance the sampling of the metastable
states that might have been missed by methods not applying configurational
biases. As the multidimensional biasing potentials can flatten out the free
energy landscape in both configurational and alchemical space, the system
would not get stuck in the phase space like the scenarios shown in Figure 1.
This approach can be easily generalized to continuous alchemical space, but
such a generalization is not explored in this study because methods such as
$\lambda$-dynamics are not currently implemented in GROMACS.

###  2.3 Free energy calculations

Theoretically, the free energy estimator for alchemical metadynamics is the
same as the one used in any other metadynamics except that the CV vector is
generalized with the introduction of the alchemical variable. Upon the
deposition of the biasing potential $V({\mbox{\boldmath{$\xi^{\prime}$}}})$ in
alchemical metadynamics, the probability distribution sampled during the
simulation is
$\tilde{P}({\mbox{\boldmath{$\xi^{\prime}$}}})\propto\exp(-\beta(F({\mbox{\boldmath{$\xi^{\prime}$}}})+V({\mbox{\boldmath{$\xi^{\prime}$}}})))$,
where $\beta=1/k_{\text{B}}T$ is the inverse temperature. To recover the
underlying free energy landscape
$F({\mbox{\boldmath{$\xi^{\prime}$}}})=-k_{\text{B}}T\ln
P({\mbox{\boldmath{$\xi^{\prime}$}}})$, we need to reweight the histogram by
assigning an unbiasing weight $w({\mbox{\boldmath{$\xi^{\prime}$}}})$ to each
sample with the CV $\xi^{\prime}$. 36 Such an unbiasing weight can be
expressed as

$w({\mbox{\boldmath{$\xi^{\prime}$}}})\propto\exp\left(\frac{V({\mbox{\boldmath{$\xi^{\prime}$}}},t_{f})}{k_{\text{B}}T}\right)$
(4)

where $t_{f}$ is the simulation length and
$V({\mbox{\boldmath{$\xi^{\prime}$}}},t_{f})$ is the total bias. The maximum
of $V({\mbox{\boldmath{$\xi^{\prime}$}}},t_{f})$ over $\xi^{\prime}$ is
usually subtracted before taking the exponential to avoid overflow, which does
not affect the normalized weights. More frequently,
$V({\mbox{\boldmath{$\xi^{\prime}$}}},t_{f})$ is replaced with
$\bar{V}({\mbox{\boldmath{$\xi^{\prime}$}}},t_{0})$, the total bias averaged
over the time period from $t_{0}=(1-f_{a})t_{f}$ to $t_{f}$ 37, where $f_{a}$
is the average fraction. Given that $t_{f}=t_{0}+N\tau$,
$\bar{V}({\mbox{\boldmath{$\xi^{\prime}$}}},t_{0})$ can be written as

$\bar{V}({\mbox{\boldmath{$\xi^{\prime}$}}},t_{0})=\frac{1}{N+1}\sum_{i=0}^{N}V({\mbox{\boldmath{$\xi^{\prime}$}}},t_{0}+i\tau)$
(5)

where $N$ is the number of Gaussians deposited from $t_{0}$ to $t_{f}$. To
calculate $P({\mbox{\boldmath{$\xi^{\prime}$}}})$ and its uncertainty (hence
$F({\mbox{\boldmath{$\xi^{\prime}$}}})$ and its uncertainty), we first discard
the equilibrium phase during which the major free energy basins were being
filled 38, with a truncation fraction of $f_{tr}$, which we here set to
$1-f_{a}$. Then, we divide the remaining part of the simulation into blocks,
for each of which we construct a weighted histogram of the CVs. Lastly, we
calculate the free energy from the probability averaged over all the blocks,
with the error of the free energy determined as the standard deviation of all
bootstrap iterations in block bootstrapping. In practice, the uncertainty of
the free energy is dependent on the number of blocks. We hence analyzed the
data with different numbers of blocks ranging 20 up to 2000 and we report the
maximum uncertainty, with its corresponding number of blocks/block size
adopted.

##  3 Methods

We validated our implementation of alchemical metadynamics with free energy
calculations for different test systems/alchemical processes with varying
complexities and dimensions of the CV space. These range from decoupling an
argon atom (System 1) from water, or a model molecule composed of 4
interaction sites (System 2, as shown in Figure 2A) from water, to the binding
and unbinding event of CB8-G3 host-guest binding complex (System 3, as shown
in Figure 2B). In the following subsections, we describe the simulation
methods of different test systems, along with the details of the corresponding
free energy calculations. All simulations were performed at 298K using GROMACS
2021.4. For metadynamics simulations, either a testing branch based on PLUMED
2.7 or PLUMED 2.8 was used, with no difference between the two versions except
for the way of specifying relevant parameters. All files relevant to our
simulations and test systems can be found in our project repository.

Figure 2: (A) The two torsional metastable states of System 2: A molecule
composed of 4 interaction sites. The first and last atom have charges of +0.2e
and -0.2e, respectively, with the other two atoms uncharged. (B) The bound
state of System 3: CB8-G3 host-guest binding complex, with the host molecule
(CB8, blue) and guest molecule (G3, red) colored in red and blue,
respectively. The fictitious sphere shown in the host molecule was used for
estimating the number of water molecules that replace the space of the guest
molecule in the binding cavity.

###  3.1 System 1: Hydration of an argon atom

As a sanity check for our implementation of alchemical metadynamics, we used
1D well-tempered alchemical metadynamics to calculate the solvation free
energy of an argon atom, which was then compared with the result obtained from
an expanded ensemble given the similarity of the two methods. With System 1,
the goal is to check if 1D alchemical metadynamics can accurately reproduce
the free energy of a simple system calculated by expanded ensemble.

####  3.1.1 Preparation of simulation inputs

The argon atom was solvated in a cubic box of length 2.4 nm and was energy-
minimized by the steepest descent algorithm until the maximum force was lower
than 100 kJ/mol/nm. Argon is modeled as a Lennard-Jones sphere with
$\epsilon=0.996$ kJ/mol and $\sigma=0.341$ nm. The system was then
equilibrated in the NVT and then NPT ensembles, both for 200 ps. The reference
temperature and pressure were maintained at 298 K and 1 bar by the velocity
rescaling method 39 and a Berendsen barostat 40, respectively. Lastly, 5 ns of
NPT MD simulation with Parrinello-Rahman barostat 41, 42 keeping the pressure
at 1 bar was performed, in which the cutoff distance for van der Waals
interactions was specified as 0.9 nm. The PME (particle mesh Ewald) method 43
was used with a switching function between 0.89 nm and 0.9 nm for efficient
calculations of long-range electrostatic interactions. A spacing of 0.10 nm
was used for the PME grids. The LINCS 44 algorithm was used to constrain bonds
involving hydrogens, with the highest order in the expansion of the constraint
coupling matrix set as 12 and the number of iterative corrections set as 2.
The structure whose box volume was the closest to the average volume of the MD
trajectory was extracted to serve as the input configuration for the expanded
ensemble and alchemical metadynamics simulations, which were both performed in
an NVT ensemble to avoid any potential issues with $\lambda$ dependence of
pressure. Although this could lead to a slightly different estimate of the
solvation free energy as compared to the one solved in the NPT ensemble, the
objective is simply to compare two methods with the same alchemical processes
in the same ensemble.

####  3.1.2 Expanded ensemble simulation

We divided our expanded ensemble calculations into two separate stages: an
equilibration and a production stage. Both stages of simulations were
performed in the NVT ensemble with 6 states for decoupling the van der Waals
interactions between the argon atom and the water molecules. In the
equilibration stage, we employed the Wang-Landau algorithm to adaptively
estimate the weight for each alchemical state, with the initial Wang-Landau
incrementor set as 0.5 $\text{k}_{\text{B}}\text{T}$. The histogram that kept
track of the state visitation was updated using the Metropolized-Gibbs Monte
Carlo moves between all alchemical states every 10 integration steps. We
adopted the default value of 0.8 for the cutoff for the flatness ratio $R$,
which means that the histogram was considered flat only if all intermediate
states had an $R$ value and its reciprocal larger than 0.8. ($R$ is defined as
the ratio between the count of a state and the average count of all states.)
Whenever the histogram was considered flat, the state counts were all reset to
0 and the Wang-Landau incrementor was scaled by a factor of 0.8. This process
for updating weights was stopped when the incrementor fell below 0.001
$\text{k}_{\text{B}}\text{T}$. The equilibrated weights were then used as the
frozen weights in the production simulation for 100 ns. For data analysis, we
only considered the time series of the Hamiltonian derivative obtained from
the production stage. After truncating the non-equilibrium regime 31 and
decorrelating the time series, we ran MBAR 45 to compute the free energy
difference between the coupled and uncoupled states, which is the solvation
free energy of the system.

####  3.1.3 1D alchemical metadynamics

To compare with expanded ensemble, we adopted the same simulation length (100
ns), starting configuration, state transition scheme (Metropolized-Gibbs
sampling), and coupling parameters of the same 6 states in 1D alchemical
metadynamics. We set the initial height of the Gaussian biasing potential as
0.5 $\text{k}_{\text{B}}\text{T}$, which is the same as the initial Wang-
Landau incrementor despite different weight updating schemes. While in
alchemical metadynamics, the strides for Gaussian depositions and MC moves are
decoupled and do not need to be the same, we set both strides as 10
integration steps to better compare with expanded ensemble, where a weight is
always applied when a move is proposed. To accommodate such a fast pace for
applying Gaussian biases, we set the bias factor as 50 to avoid an excessively
fast decay in the Gaussian height, which could potentially slow down the
compensation of the underlying free energy surface. The width of the Gaussian
was set as 0.01 to avoid any overlap between the Gaussians deposited at
different $\lambda$ values. Note that having such an overlap would not
invalidate the method, but might have made the comparison to expanded ensemble
less straightforward. For the solvation free energy calculation, we set the
truncation fraction and the average fraction as 0.25 and 0.75, respectively.
1000 blocks were used in the histogram construction, which led to a block size
of 75 ps and the largest uncertainty among the considered numbers of blocks in
the estimation of the solvation free energy. 200 bootstrap iterations were
used in block bootstrapping.

###  3.2 System 2: A 4-site system

To demonstrate the advantages of introducing configurational CVs in alchemical
metadynamics, we designed a fictitious molecule composed of 4 linearly bonded
interaction sites as the second system (see Figure 2A). We placed opposite
charges on the first and last atoms (+0.2e and -0.2e) and set the force
constant for the only torsional angle as 60 kJ/mol so that the two torsional
metastable states of the system were separated by a large free energy barrier
(around 48.5 $\text{k}_{\text{B}}\text{T}$, see Figure S3A). Such a large free
energy barrier poses a challenge of sampling both torsional metastable states
to alchemical free energy methods that do not apply any configurational bias,
such as expanded ensemble or 1D alchemical metadynamics. Such a challenge is
useful for us to highlight the difference between methods with or without the
application of configurational biases in free energy calculations. With this
system, we performed a 1D and 2D well-tempered alchemical metadynamics
starting from either metastable state to calculate the solvation free energy
of the system. The goal is to show that 2D alchemical metadynamics starting
from different torsional metastable states leads to statistically consistent
estimates of the solvation free energy, which cannot be accomplished by 1D
alchemical metadynamics due to restricted configurational sampling.

####  3.2.1 Preparation of simulation inputs

After solvating System 2 in a dodecahedral box with 1 nm between the solute
and box edges, we processed the system with the same procedure of solvation,
energy minimization, NVT and NPT equilibration, and NPT MD simulation like the
one used in System 1. A structure that had a volume closest to the average
over the NPT volume was taken as the input of a 5 ns NVT metadynamics that
only biased the torsional angle of the system. This torsional metadynamics
applied a Gaussian biasing potential every 500 integration steps, with a bias
factor of 10. The width and initial height of the Gaussians were set as 0.5
rad and 1 $\text{k}_{\text{B}}\text{T}$, respectively. With this setup, the
system was able to sample both torsional metastable states frequently in the
torsional metadynamics (see Figure S3B), from which we generated the
structures corresponding to the trans isomer (State A) and cis isomer (State
B) for starting subsequent alchemical metadynamics elaborated in later
sections.

####  3.2.2 Alchemical metadynamics

For each torsional metastable state of the system, we started a 1D and 2D
well-tempered alchemical metadynamics in an NVT ensemble, where the only
torsional angle of the system was introduced in the latter as the second CV.
The MC moves between alchemical states were proposed every 10 integration
steps using the Metropolized-Gibbs MC scheme. All four simulations were
performed for 200 ns and adopted a deposition stride of 500 steps. 20
alchemical intermediate states were used to decouple the van der Waals
interactions and Coulombic interactions. 1D alchemical metadynamics started
with a Gaussian height of 1 $\text{k}_{\text{B}}\text{T}$, with the bias
factor set as 60, while 2D alchemical metadynamics utilized a Gaussian height
of 2 $\text{k}_{\text{B}}\text{T}$, and a bias factor of 120 for flattening a
much larger phase space. The widths of the Gaussian along the alchemical
direction (for 1D and 2D simulations) were set as 0.001, while the Gaussian
width in the torsional dimension (for the 2D simulations) was set as 0.5 rad.
For all four simulations, we specified a truncation fraction of 0.25 and an
average fraction of 0.75. 1000 blocks with a block size of 150 ps were used
for all simulations, except for the 1D alchemical metadynamics starting from
state A, which used 2000 blocks for histogramming in free energy calculations.
200 bootstrap iterations were used in block bootstrapping.

###  3.3 System 3: CB8-G3 host-guest binding complex

To further exemplify the usage of alchemical metadynamics, we considered the
CB8-G3 host-guest binding complex (see Figure 2B) as the third test system,
which is composed of a cucurbit[8]uril (CB8, host) and a quinine (G3, guest)
and is one of the binding complexes in SAMPL6 SAMPLing challenge 46. The
binding free energy calculation of this system is considered particularly
challenging due to the long time scale of the water motion, especially the
exclusion of water molecules from the binding cavity during the binding event
(the so-called water-exclusion problem) 47. The past literature 47 shows that
expanded ensemble was not able to equilibrate the Wang-Landau weights within
the same time used for other binding complexes in the same challenge. Methods
such as OpenMM/SOMD and NAMD/BAR replicate calculations could not converge the
average free energy to uncertainties under 1 kcal/mol. OpenMM/HREX and
AMBER/APR could reach sub-kcal/mol uncertainties, but both of them displayed a
significant and slowly decaying bias, requiring a very long time to converge
the free energy. Here, our goal is to compare 2D alchemical metadynamics with
its 1D analog to show its superiority in accelerating the water motion and
accurate estimation of the binding free energy of the binding complex.

####  3.3.1 Preparation of simulation inputs

We started from the initial files provided by SAMPL6 SAMPLing challenge, which
included the coordinates and topology files of the solvated binding complex
and the solvated ligand alone. We re-equilibrated both systems in an NPT MD
simulation with the same settings as used in Systems 1 and 2. We then
extracted the configurations corresponding to the average NPT volumes as the
input structures of later simulations.

####  3.3.2 The double decoupling scheme for the binding free energy
calculation

We used the double decoupling method 48 with the thermodynamic cycle shown in
Figure 3 to decompose the binding free energy calculation of CB8-G3 binding
complex. Starting from state A, we split the binding event into two alchemical
processes (from A to B and from C to D) bridged by the transfer of the ligand
molecule into the binding cavity of the host molecule (from B to C). The first
alchemical process (from A to B) decoupled the ligand from its surroundings
and was termed as the solvent simulation. The second alchemical process (from
C to D) was characterized in the complex simulation, where we started from
state D and gradually turn off the non-bonded interactions between the host
and guest molecules to calculate $-\Delta A^{\text{solv}}_{\text{(elec +
vdW)}_{\text{off}}}$ rather than $\Delta A^{\text{solv}}_{\text{(elec +
vdW)}_{\text{off}}}$. Frequently, in the complex simulation, a harmonic
restraint is applied between the host and guest molecules to prevent the guest
molecule from drifting away and slowing down the convergence. However, the
periodic boundary conditions may be easily broken by a known bug in the pull
code working with expanded ensemble, so we did not apply any restraint in the
system but assumed that the ligand was able to sample to the entire simulation
box in the decoupled state. The assumption turned out to be valid (see Figure
S1) and allowed us to calculate the box volume correction term $\Delta
A_{\text{corr}}$ analytically to bring the affinity in units of the standard
concentration. Namely,

$\Delta A^{\circ}_{\text{bind}}=\Delta A_{\text{bind}}+\Delta
A_{\text{corr}}=\Delta
A_{\text{bind}}-k_{\text{B}}T\ln\left(\frac{V_{\text{box}}}{V_{0}}\right)$ (6)

where $\Delta A^{\circ}_{\text{bind}}$ is the binding free energy under
standard temperature and pressure, $V_{\text{box}}$ is the box volume and
$V_{0}$ is the molecular volume (1.6605 nm3) corresponding to the 1 mol/L
reference concentration.

Figure 3: The double decoupling thermodynamic cycle for the binding free
energy calculation of CB8-G3 host-guest binding complex. The host molecule was
colored in blue, while the guest ligand was colored in either white (fully
decoupled from the host) or red (fully coupled with the host).

####  3.3.3 Alchemical metadynamics

We used 1D well-tempered alchemical metadynamics in the NVT ensemble as the
solvation simulation. 40 alchemical intermediate states were used to gradually
decouple the non-bonded interactions between the ligand and the solvent, with
Metropolized-Gibbs moves between them attempted every 10 integration steps.
Gaussians with a width of 0.01 and the initial height of 2
$\text{k}_{\text{B}}\text{T}$ were deposited with the same frequency as the MC
moves, with its decaying rate characterized by a bias factor of 100. We
computed $-\Delta A^{\text{solv}}_{\text{(elec + vdW)}_{\text{off}}}$ from 50
block histograms (block size of 1.5 ns) constructed upon the CV time series
processed with a truncation fraction of 0.25 and an average fraction of 0.75
for reweighting. 200 bootstrap iterations were carried out in block
bootstrapping to estimate the error of the free energy difference.

For the complex simulation, we performed and compared both 1D and 2D well-
tempered alchemical metadynamics in the NVT ensemble. The two simulations
utilized the same parameters except the ones in 2D alchemical metadynamics for
sampling the additional CV. Specifically, both simulations applied 40
alchemical intermediate states to decouple the non-bonded interactions between
the host and guest molecules. The strides for MC moves and Gaussian
depositions were both 10 integration steps. The Gaussian biasing potentials
had a width of 0.01 and the initial height of 5 $\text{k}_{\text{B}}\text{T}$.
The bias factor was set as 150 to overcome the significant probability
imbalance between the coupled and uncoupled states. Note that this is common
in both the solvent and complex simulations and it necessitated the relatively
frequent deposition of the biasing potentials.

In 2D alchemical metadynamics, the number of water molecules $N$ in the
binding cavity that replaced the space of the guest molecule was introduced as
the additional configurational CV. The value of $N$ was estimated as the
coordination number $S$ between any host atom (set A) and any oxygen atom of
the water molecules (set B) in the fictitious sphere shown in Figure 2B. The
sphere was centered at the COM of the host molecule and had a radius of
$R_{0}=0.35$ nm. Mathematically, this is done by a switching function (see
Figure S4) to ensure continuous derivatives:

$S=\sum_{i\in A}\sum_{j\in
B}\frac{1}{1+\left(\frac{r_{ij}}{R_{0}}\right)^{6}}$ (7)

where $r_{ij}$ is the inter-particle distance between atoms in sets $A$ and
$B$. Note that we did not aim to exclude all the water molecules because the
bound state of the binding complex could still have water molecules in the
binding cavity. To avoid wasting time in configurations that were likely to be
unphysical, we set an upper and lower harmonic wall potential (with the force
constant equal to 200 kJ/mol) at $N=10.5$ and $N=0.7$, respectively. These
cutoffs were decided from the time series obtained from 10 standard MD
simulations of the binding complex (see the section of Standard MD simulations
of System 3 in the SI). All the standard MD simulations were performed in an
NVT ensemble for 200 ns and utilized the same protocol and parameters as the
ones used in the NVT equilibration during the preparation of Systems 1 or 2.

For the calculation of the free energy difference $\Delta
A^{\text{complex}}_{\text{(elec+vdW)}_{\text{on}}}$, we adopted a truncation
fraction of 0.25 and an average fraction of 0.75 for both 1D and 2D alchemical
metadynamics. 100 blocks (block size: 1.5 ns) and 50 blocks (block size: 3 ns)
were respectively used in the two cases for data histogramming, which was
followed by 200 bootstrap iterations in block bootstrapping. While the
protocol for data analysis of the 1D simulation was not different from the one
previously described, for the 2D alchemical metadynamics, we defined the
unbiasing weight for each sample differently as

$w({\mbox{\boldmath{$\xi^{\prime}$}}})=\exp\left(\frac{\bar{V}_{\text{tot}}({\mbox{\boldmath{$\xi^{\prime}$}}},t_{0})}{k_{\text{B}}T}\right)$
(8)

where
$\bar{V}_{\text{tot}}({\mbox{\boldmath{$\xi^{\prime}$}}},t_{0})=\bar{V}({\mbox{\boldmath{$\xi^{\prime}$}}},t_{0})+V_{\text{walls}}(N)$
is the total bias that includes the biases from metadynamics and the sum of
both wall potentials ($V_{\text{walls}}$). Such unbiasing weights could remove
the influence coming from the wall potentials during the reweighting. With the
free energy difference calculated for each alchemical process, we calculate
$\Delta A^{\circ}_{\text{bind}}$ by simply summing them up, along with the
correction term $\Delta A_{\text{corr}}$.

##  4 Results and discussion

###  4.1 System 1: An argon atom

With the weights fixed at the values equilibrated by the Wang-Landau
algorithm, the solvation free energy of System 1 estimated by expanded
ensemble with MBAR was -3.302 $\text{k}_{\text{B}}\text{T}$, with an
uncertainty as small as 0.027 $\text{k}_{\text{B}}\text{T}$ owing to
sufficiently even state visitation (see Figure S5A). This estimation is
statistically consistent with the one obtained from the 1D alchemical
metadynamics, which was -3.296 $\pm$ 0.017 $\text{k}_{\text{B}}\text{T}$.
Notably, the essential difference in the weight updating approaches between
the two methods makes it infeasible to directly compare the performance of the
methods. In expanded ensemble, the incrementor decreases in a step-wise manner
and the same value is applied to all intermediate states. On the other hand,
the Gaussian height in well-tempered metadynamics continuously decays with the
amount of biases that have been deposited in the state being sampled, which
means that the potential energies of different states are elevated with
different amounts depending on how frequently the states have been visited. In
our case, the state visitation of expanded ensemble could have been more even
if we had adopted a stricter, rather than a typical criterion for histogram
flatness, which could also lead to an even smaller uncertainty. With the
chosen parameters in our case, though, the expanded ensemble still reached a
low uncertainty. More importantly, the data collected from the two simulations
is sufficient to indicate that 1D alchemical metadynamics free energy
calculations yield results that are equivalent within statistical significance
to the expanded ensemble results.

###  4.2 System 2: A molecule composed of 4 interaction sites

The purpose of System 2 is to demonstrate the difference in the
configurational sampling between alchemical metadynamics with or without the
introduction of the configurational bias. As a consequence, in both 1D
alchemical metadynamics simulations, the accumulation of one-dimensional
alchemical biases allowed the system to freely sample all the intermediate
states (see Figure S6). However, such biases did not facilitate the
compensation of the free energy wells along the torsional direction. With the
lack of direct biases in the torsional direction, the interconversion of the
two isomers became the slowest degree of freedom that trapped the system.
Accordingly, it can be seen in Figure 4A that 1D alchemical metadynamics
failed to sample both metastable states regardless of which torsional state
the simulation was initialized at. This insufficient sampling of the torsional
space caused the dependence of the estimated solvation free energy on the
starting torsional metastable state. Specifically, the solvation free energies
estimated by 1D alchemical metadynamics starting from State A and State B were
0.707 $\pm$ 0.032 $\text{k}_{\text{B}}\text{T}$ and -0.232 $\pm$ 0.033
$\text{k}_{\text{B}}\text{T}$ (see Figure 5), respectively. As either
simulation failed to account for the potential energy contribution of the
other metastable state, the free energy estimates were statistically
inconsistent with each other.

Figure 4: The sampled torsional angle in (A) 1D alchemical metadynamics and
(B) 2D alchemical metadynamics as a function of time. As can be seen in the
figure, the sampling in the torsional space is restricted in 1D alchemical
metadynamics, but essentially complete in its 2D analog.

By contrast, the 2D Gaussian biasing potentials in both 2D alchemical
metadynamics simulations flattened out the free energy surface both along the
alchemical and torsional directions simultaneously, so both simulations were
able to sample the alchemical and torsional space exhaustively (see Figure S7
and Figure 4B). This comprehensive and even sampling along all the slow
degrees of freedom led to statistically consistent binding free energy
estimations from State A (0.643 $\pm$ 0.031 $\text{k}_{\text{B}}\text{T}$) and
State B (0.697 $\pm$ 0.032 $\text{k}_{\text{B}}\text{T}$), as shown in Figure
5. In addition, the 2D free energy surface of the system can be accurately
recovered from either of the two simulations. Figure 6 shows the average of
the 2D free energy surfaces obtained from the two 2D alchemical metadynamics
simulations and its corresponding contour plot. The difference plot of the two
surfaces can be found in Figure S8. From the 2D free energy surface, it is
clear that the alchemical variable is orthogonal to the torsional angle of the
system, which explains the ineffectiveness of the one-dimensional alchemical
bias in the torsional sampling. Notably, this is exactly Scenario A in Figure
1 that can fail alchemical free energy methods without the application of
configurational biases, so the success in System 2 just verified the usage of
2D alchemical metadynamics in overcoming the extensive free energy barrier
present in a certain configurational CV direction.

Figure 5: The solvation free energies estimated by 1D and 2D alchemical
metadynamics starting from either states. The two 1D alchemical metadynamics
simulations starting from different torsional metastable states led to
statistically different estimations of the binding free energy, while the
values estimated by the 2D simulations are statistically consistent with each
other. Figure 6: (A) The average of the 2D free energy surfaces obtained from
the two 2D alchemical metadynamics simulations starting from State A and State
B. (B) The average of the contour plots obtained from the 2D alchemical
metadynamics simulations starting from State A and State B.

###  4.3 System 3: CB8-G3 host-guest binding complex

In System 3, the impediment to frequent sampling of both the bound and unbound
states is the long time scale of the water molecules entering and exiting the
binding cavity. Our hypothesis is that the 2D alchemical metadynamics can
accelerate the slow water motion and thus reduce the time needed to converge
the sampling.

To examine how the dynamics of the water molecules was changed upon the
introduction of the configurational CV in the complex simulation, we plotted
the time series of the number of water molecules in the fictitious sphere (see
Figure 2) representing the volume occupied by the guest ligand (see Figure 7).
As a result, in the course of the 200 ns 1D alchemical metadynamics, there
were only 3 rounds of binding/unbinding events. More specifically, Figure 7A
shows that it only took 2 to 9 ns for the water to enter the binding cavity
after the ligand unbound, while the water exclusion always took more than 45
ns. Notably, the alchemical variable is not strictly orthogonal to the number
of water molecules, which implies that the alchemical biases could also
indirectly bias the motion of the water. This explains the shortened time
scale of water inclusion in the 1D alchemical metadynamics as compared to the
standard MD simulations of System 3 (see Figure 7A and the section of Standard
MD simulations of the binding complex in the SI), as the alchemical biases
allowed the system to drift to weakly coupled states where the ligand could
unbind easily. This indirect influence from alchemical biases is also
evidenced by the several transitions between the two metastable states
occurred in the 1D alchemical metadynamics, as contrasted with the lack of
transitions in the standard MD simulations of the System 3 (see Figure 7A and
the section of Standard MD simulations of System 3 in the SI). It also
explains the gradually increasing time scale of the water motion in both
directions, which can be ascribed to the fact that decreasing alchemical
biases were less efficient in pushing the water molecules into or out of the
binding cavity. Unlike the 1D alchemical metadynamics in System 2, the 1D
simulation here did sample both configurational metastable states, although
the event was too rare to provide reliable statistics in free energy
reconstructions, yielding a free energy difference between the coupled and
uncoupled states ($-\Delta A^{\text{complex}}_{\text{(elec +
vdW)}_{\text{on}}}$) of 133.87 $\pm$ 0.15 $\text{k}_{\text{B}}\text{T}$.

Figure 7: The number of water molecules as a function of time in (A) 1D
alchemical metadynamics and (B) 2D alchemical metadynamics. The two figures
are using the same scale for the y-axis. As can be seen from the figure, there
were much more transitions in 2D alchemical metadynamics compared with the 1D
simulation.

In the 2D alchemical metadynamics, the acceleration of the water motion no
longer merely relied on the marginal contribution from the alchemical bias,
but was directly aided by the bias in the configurational direction. As such,
the time scale of water exclusion was largely shortened, allowing over 20
transitions between the two metastable states during the simulation. Notably,
the water molecules were squeezed in the binding cavity by the configurational
bias, so the maximum number of water molecules (13.4) was larger than that in
the standard MD simulation (around 10.5) even under the influence of the upper
wall potential set at 10.5. Similar to 1D alchemical metadynamics, the time
scale of the water motion gradually increased with the decreasing Gaussian
height. Given the well-sampled configurational space, the free energy
difference between the coupled and uncoupled states ($-\Delta
A^{\text{complex}}_{\text{(elec + vdW)}_{\text{on}}}$) was estimated as 132.92
$\pm$ 0.15 $\text{k}_{\text{B}}\text{T}$, which was not far but statistically
different from the one obtained in the 1D simulation.

With the free energy difference $-\Delta A^{\text{complex}}_{\text{(elec +
vdW)}_{\text{on}}}$ estimated, we could finally estimate the binding free
energy of the binding complex $\Delta A^{\circ}_{\text{bind}}=\Delta
A^{\text{solv}}_{\text{(elec + vdW)}_{\text{off}}}+\Delta
A^{\text{complex}}_{\text{(elec +
vdW)}_{\text{on}}}-k_{\text{B}}T\ln\left(\frac{V_{\text{box}}}{V_{0}}\right)$.
Specifically, $\Delta A^{\text{solv}}_{\text{(elec + vdW)}_{\text{off}}}$ was
estimated to be 131.89 $\pm$ 0.16 $\text{k}_{\text{B}}\text{T}$ by 1D
alchemical metadynamics as the solvent simulation, which in theory should be
accurate given the comprehensive alchemical sampling (see Figure S9) and the
lack of configurational metastable states that the simulation could have
missed. The correction term $\Delta A_{\text{corr}}$ was calculated as -3.70
$\text{k}_{\text{B}}\text{T}$, assuming that the guest ligand in the fully
decoupled state could sample the entire simulation box, as shown by testing.
Consequently, the total binding free energy was estimated to be -4.73 $\pm
0.22$ $\text{k}_{\text{B}}\text{T}$.

Lastly, from the 2D simulation, we additionally recovered the 2D free energy
surface of the system, with the dependence between the alchemical variable and
the configurational CV removed by shifting the whole free energy surface by
its 1D projection along the alchemical direction (see Figure 8). This free
energy surface exhibits two high-probability regions in the sampling space,
which include the narrow ditch approximately at $N=1$ to $N=1.5$, and a deep
basin in higher alchemical states with large $N$ values. The extended ditch
clearly corresponds to the bound state, while the deep basin in the weakly
coupled states reflects the preference of hydrated configurations over the
other end. In the fully coupled state, the free energy barrier peaks around
$N=6$, which separates the bound (approximately $N<1.5$) and unbound
(approximately $8<N<9$). Notably, in the fully coupled state, crossing the
free energy barrier from the bound state to the unbound state is harder than
in the opposite direction, which caused the extended time scale of water
exclusion. On the other hand, water exclusion is the slowest degree of freedom
in weakly coupled states that can be reached in alchemical sampling. By
flattening out the free energy surface, 2D alchemical metadynamics shortened
both time scales. Importantly, the free energy basins shown in Figure 8
categorized System 3 as Scenario B in 1, so the escape from such free energy
basins demonstrates that alchemical metadynamics that biases configurational
CV can address the two most common challenges shown in 1 that fail most
traditional alchemical free energy methods.

Figure 8: (A) The 2D free energy surface and (B) contour plot of System 3. The
2D free energy surface was computed from the 2D alchemical metadynamics but
the marginalized free energy as a function of $\lambda$ was removed to make it
possible to see the free energy in the water occupancy direction as a function
of lambda.

##  5 Conclusion

In this study, we proposed alchemical metadynamics, which expanded the
configurationally defined sampling space allowed in traditional metadynamics
with an additional sampling direction. Alchemical metadynamics is most useful
when the CV space is multidimensional. With the configurational bias, it
encourages the system to escape from configurational metastable subspace that
could have easily trap the system. It retains the advantages of traditional
alchemical free energy methods, but also enables higher flexibility in
sampling rough free energy surfaces.

With different test systems and alchemical processes, we showed that 1D
alchemical metadynamics had at least comparable performance as expanded
ensemble, and was able to accurately calculate the solvation free energy of an
argon atom. We also showed that 2D alchemical metadynamics could eliminate the
dependence of free energy calculations on the starting metastable state due to
restricted configurational sampling in System 2. In System 3, we demonstrated
that 2D alchemical metadynamics could accelerate the water exclusion by
shortening its time scale approximately by a factor of 10. More importantly,
the success in both Systems 2 and 3 manifests the usage of alchemical
metadynamics in overcoming challenges that can frustrate traditional
alchemical free energy methods that do not bias configurational CVs. We
conclude that alchemical metadynamics is promising in enhancing sampling in
challenging systems, such as highly flexible protein-peptide binding
complexes, or protein-nucleic acid binding complexes. The method can be
trivially combined with more sophisticated algorithms, such as path collective
variable 49 tICA 50, SGOOP 51, RAVE 52, or other similar machine learning
methods 53, 54, 55, 56, 57, 58.

## Author Contributions

W.-T.H., M.R.S. primarily conceptualized the project and designed the
methodology, with additional contributions from P.T.M and G.B.; P.T.M.
implemented the sampling methods in PLUMED and GROMACS, with contributions
from G.B. Experiments were performed by W.-T.H.; experiments were analyzed and
validated by W.-T.H. with contributions from M.R.S.; W.-T.H wrote the original
manuscript draft; editing and review of the manuscript was done by M.R.S.,
G.B. and P.T.M. M.R.S. supervised the project and obtained resources.

This study was supported by the grant from National Science Foundation
(OAC-1835720) (W.-T.S and M.R.S) and an MolSSI Software Fellowship (P.T.M.).
MolSSI is funded by NSF CHE-2136142. The computational work done in this
publication used resources provided from the Extreme Science and Engineering
Discovery Environment (XSEDE), which is supported by National Science
Foundation grant number ACI-1548562. Specifically, it used the Bridges-2
system, which is supported by NSF ACI-1928147, located at the Pittsburgh
Supercomputing Center (PSC).

## References

  * Hénin et al. 2022 Hénin, J.; Lelièvre, T.; Shirts, M. R.; Valsson, O.; Delemotte, L. Enhanced sampling methods for molecular dynamics simulations. _arXiv preprint arXiv:2202.04164_ 2022,
  * Torrie and Valleau 1977 Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. _Journal of Computational Physics_ 1977, _23_ , 187–199
  * Darve et al. 2008 Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. _The Journal of chemical physics_ 2008, _128_ , 144120
  * Laio and Parrinello 2002 Laio, A.; Parrinello, M. Escaping free-energy minima. _Proceedings of the National Academy of Sciences_ 2002, _99_ , 12562–12566
  * Bussi et al. 2006 Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M. Free-energy landscape for $\beta$ hairpin folding from combined parallel tempering and metadynamics. _Journal of the American Chemical Society_ 2006, _128_ , 13435–13441
  * Bonomi and Parrinello 2010 Bonomi, M.; Parrinello, M. Enhanced sampling in the well-tempered ensemble. _Physical review letters_ 2010, _104_ , 190601
  * Piana and Laio 2007 Piana, S.; Laio, A. A bias-exchange approach to protein folding. _The journal of physical chemistry B_ 2007, _111_ , 4553–4559
  * Marinari and Parisi 1992 Marinari, E.; Parisi, G. Simulated tempering: a new Monte Carlo scheme. _EPL (Europhysics Letters)_ 1992, _19_ , 451
  * Sugita and Okamoto 1999 Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. 1999
  * Sugita et al. 2000 Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional replica-exchange method for free-energy calculations. _The Journal of Chemical Physics_ 2000, _113_ , 6042–6051
  * Lyubartsev et al. 1992 Lyubartsev, A.; Martsinovski, A.; Shevkunov, S.; Vorontsov-Velyaminov, P. New approach to Monte Carlo calculation of the free energy: Method of expanded ensembles. _The Journal of chemical physics_ 1992, _96_ , 1776–1783
  * Knight and Brooks III 2009 Knight, J. L.; Brooks III, C. L. $\lambda$-Dynamics free energy simulation methods. _Journal of computational chemistry_ 2009, _30_ , 1692–1700
  * Oshima et al. 2019 Oshima, H.; Re, S.; Sugita, Y. Replica-exchange umbrella sampling combined with Gaussian accelerated molecular dynamics for free-energy calculation of biomolecules. _Journal of Chemical Theory and Computation_ 2019, _15_ , 5199–5208
  * Knight and Brooks III 2011 Knight, J. L.; Brooks III, C. L. Multisite $\lambda$ dynamics for simulated structure–activity relationship studies. _Journal of chemical theory and computation_ 2011, _7_ , 2728–2739
  * Marques et al. 2019 Marques, S. M.; Bednar, D.; Damborsky, J. Computational study of protein-ligand unbinding for enzyme engineering. _Frontiers in chemistry_ 2019, _6_ , 650
  * Wang et al. 2017 Wang, Y.; Martins, J. M.; Lindorff-Larsen, K. Biomolecular conformational changes and ligand binding: from kinetics to thermodynamics. _Chemical science_ 2017, _8_ , 6466–6473
  * Trapl et al. 2020 Trapl, D.; del Río, C. C.; Kříž, P.; Spiwok, V. Prediction of pKa in a system with high orthogonal barriers: Alchemical flying Gaussian method. _Chemical Physics Letters_ 2020, _760_ , 138012
  * Hahn et al. 2020 Hahn, D. F.; König, G.; Hünenberger, P. H. Overcoming orthogonal barriers in alchemical free energy calculations: On the relative merits of $\lambda$-variations, $\lambda$-extrapolations, and biasing. _Journal of Chemical Theory and Computation_ 2020, _16_ , 1630–1645
  * Bieler et al. 2014 Bieler, N. S.; Häuselmann, R.; Hünenberger, P. H. Local elevation umbrella sampling applied to the calculation of alchemical free-energy changes via $\lambda$-dynamics: The $\lambda$-LEUS scheme. _Journal of Chemical Theory and Computation_ 2014, _10_ , 3006–3022
  * Bieler et al. 2015 Bieler, N. S.; Tschopp, J. P.; Hünenberger, P. H. Multistate $\lambda$-local-elevation umbrella-sampling (MS-$\lambda$-LEUS): Method and application to the complexation of cations by crown ethers. _Journal of Chemical Theory and Computation_ 2015, _11_ , 2575–2588
  * Zheng et al. 2008 Zheng, L.; Chen, M.; Yang, W. Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems. _Proceedings of the National Academy of Sciences_ 2008, _105_ , 20227–20232
  * Zheng and Yang 2012 Zheng, L.; Yang, W. Practically efficient and robust free energy calculations: double-integration orthogonal space tempering. _Journal of Chemical Theory and Computation_ 2012, _8_ , 810–823
  * Hayes et al. 2017 Hayes, R. L.; Armacost, K. A.; Vilseck, J. Z.; Brooks III, C. L. Adaptive landscape flattening accelerates sampling of alchemical space in multisite $\lambda$ dynamics. _The Journal of Physical Chemistry B_ 2017, _121_ , 3626–3635
  * Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Computer Physics Communications_ 2014, _185_ , 604–613
  * Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. _Physical review letters_ 2008, _100_ , 020603
  * Hastings 1970 Hastings, W. K. Monte Carlo sampling methods using Markov chains and their applications. 1970,
  * Barker 1965 Barker, A. A. Monte carlo calculations of the radial distribution functions for a proton? electron plasma. _Australian Journal of Physics_ 1965, _18_ , 119–134
  * Geman and Geman 1984 Geman, S.; Geman, D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. _IEEE Transactions on pattern analysis and machine intelligence_ 1984, 721–741
  * Liu and Liu 2001 Liu, J. S.; Liu, J. S. _Monte Carlo strategies in scientific computing_ ; Springer, 2001; Vol. 10
  * Liu 1996 Liu, J. S. Peskun’s theorem and a modified discrete-state Gibbs sampler. _Biometrika_ 1996, _83_
  * Chodera 2016 Chodera, J. D. A simple method for automated equilibration detection in molecular simulations. _Journal of chemical theory and computation_ 2016, _12_ , 1799–1805
  * Desgranges and Delhommelle 2012 Desgranges, C.; Delhommelle, J. Evaluation of the grand-canonical partition function using expanded Wang-Landau simulations. I. Thermodynamic properties in the bulk and at the liquid-vapor phase boundary. _The Journal of Chemical Physics_ 2012, _136_ , 184107
  * Wang and Landau 2001 Wang, F.; Landau, D. P. Efficient, multiple-range random walk algorithm to calculate the density of states. _Physical review letters_ 2001, _86_ , 2050
  * Belardinelli and Pereyra 2007 Belardinelli, R.; Pereyra, V. Fast algorithm to calculate density of states. _Physical Review E_ 2007, _75_ , 046701
  * Dama et al. 2014 Dama, J. F.; Rotskoff, G.; Parrinello, M.; Voth, G. A. Transition-tempered metadynamics: robust, convergent metadynamics via on-the-fly transition barrier estimation. _Journal of Chemical Theory and Computation_ 2014, _10_ , 3626–3633
  * Branduardi et al. 2012 Branduardi, D.; Bussi, G.; Parrinello, M. Metadynamics with adaptive Gaussians. _Journal of chemical theory and computation_ 2012, _8_ , 2247–2254
  * Micheletti et al. 2004 Micheletti, C.; Laio, A.; Parrinello, M. Reconstructing the density of states by history-dependent metadynamics. _Physical Review Letters_ 2004, _92_ , 170601
  * Bussi and Laio 2020 Bussi, G.; Laio, A. Using metadynamics to explore complex free-energy landscapes. _Nature Reviews Physics_ 2020, _2_ , 200–212
  * Bussi et al. 2007 Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. _The Journal of chemical physics_ 2007, _126_ , 014101
  * Berendsen et al. 1984 Berendsen, H. J.; Postma, J. v.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. _The Journal of chemical physics_ 1984, _81_ , 3684–3690
  * Parrinello and Rahman 1980 Parrinello, M.; Rahman, A. Crystal structure and pair potentials: A molecular-dynamics study. _Physical review letters_ 1980, _45_ , 1196
  * Parrinello and Rahman 1981 Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. _Journal of Applied physics_ 1981, _52_ , 7182–7190
  * Essmann et al. 1995 Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. _The Journal of chemical physics_ 1995, _103_ , 8577–8593
  * Hess et al. 1997 Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: a linear constraint solver for molecular simulations. _Journal of computational chemistry_ 1997, _18_ , 1463–1472
  * Shirts and Chodera 2008 Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. _The Journal of chemical physics_ 2008, _129_ , 124105
  * Rizzi et al. 2018 Rizzi, A.; Murkli, S.; McNeill, J. N.; Yao, W.; Sullivan, M.; Gilson, M. K.; Chiu, M. W.; Isaacs, L.; Gibb, B. C.; Mobley, D. L., et al. Overview of the SAMPL6 host–guest binding affinity prediction challenge. _Journal of computer-aided molecular design_ 2018, _32_ , 937–963
  * Rizzi et al. 2020 Rizzi, A.; Jensen, T.; Slochower, D. R.; Aldeghi, M.; Gapsys, V.; Ntekoumes, D.; Bosisio, S.; Papadourakis, M.; Henriksen, N. M.; De Groot, B. L., et al. The SAMPL6 SAMPLing challenge: assessing the reliability and efficiency of binding free energy calculations. _Journal of computer-aided molecular design_ 2020, _34_ , 601–633
  * Gilson et al. 1997 Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. The statistical-thermodynamic basis for computation of binding affinities: a critical review. _Biophysical journal_ 1997, _72_ , 1047–1069
  * Branduardi et al. 2007 Branduardi, D.; Gervasio, F. L.; Parrinello, M. From A to B in free energy space. _The Journal of chemical physics_ 2007, _126_ , 054103
  * M. Sultan and Pande 2017 M. Sultan, M.; Pande, V. S. tICA-metadynamics: accelerating metadynamics by using kinetically selected collective variables. _Journal of chemical theory and computation_ 2017, _13_ , 2440–2447
  * Tiwary and Berne 2016 Tiwary, P.; Berne, B. Spectral gap optimization of order parameters for sampling complex molecular systems. _Proceedings of the National Academy of Sciences_ 2016, _113_ , 2839–2844
  * Ribeiro et al. 2018 Ribeiro, J. M. L.; Bravo, P.; Wang, Y.; Tiwary, P. Reweighted autoencoded variational Bayes for enhanced sampling (RAVE). _The Journal of chemical physics_ 2018, _149_ , 072301
  * Rohrdanz et al. 2011 Rohrdanz, M. A.; Zheng, W.; Maggioni, M.; Clementi, C. Determination of reaction coordinates via locally scaled diffusion map. _The Journal of chemical physics_ 2011, _134_ , 03B624
  * Sultan and Pande 2018 Sultan, M. M.; Pande, V. S. Automated design of collective variables using supervised machine learning. _The Journal of chemical physics_ 2018, _149_ , 094106
  * McCarty and Parrinello 2017 McCarty, J.; Parrinello, M. A variational conformational dynamics approach to the selection of collective variables in metadynamics. _The Journal of chemical physics_ 2017, _147_ , 204109
  * Chen and Ferguson 2018 Chen, W.; Ferguson, A. L. Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration. _Journal of computational chemistry_ 2018, _39_ , 2079–2102
  * Wang et al. 2019 Wang, Y.; Ribeiro, J. M. L.; Tiwary, P. Past–future information bottleneck for sampling molecular reaction coordinate simultaneously with thermodynamics and kinetics. _Nature communications_ 2019, _10_ , 1–8
  * Wehmeyer and Noé 2018 Wehmeyer, C.; Noé, F. Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics. _The Journal of chemical physics_ 2018, _148_ , 241703


# FDNet: Feature Decoupled Segmentation Network for Tooth CBCT Image

###### Abstract

Precise Tooth Cone Beam Computed Tomography (CBCT) image segmentation is
crucial for orthodontic treatment planning. In this paper, we propose FDNet, a
Feature Decoupled Segmentation Network, to excel in the face of the variable
dental conditions encountered in CBCT scans, such as complex artifacts and
indistinct tooth boundaries. The Low-Frequency Wavelet Transform (LF-Wavelet)
is employed to enrich the semantic content by emphasizing the global
structural integrity of the teeth, while the SAM encoder is leveraged to
refine the boundary delineation, thus improving the contrast between adjacent
dental structures. By integrating these dual aspects, FDNet adeptly addresses
the semantic gap, providing a detailed and accurate segmentation. The
framework’s effectiveness is validated through rigorous benchmarks, achieving
the top Dice and IoU scores of 85.28% and 75.23%, respectively. This
innovative decoupling of semantic and boundary features capitalizes on the
unique strengths of each element to significantly elevate the quality of
segmentation performance.

Index Terms—  CBCT, Tooth Segmentation, Wavelet Transform, SAM Encoder,
Feature Decoupling

##  1 Introduction

Tooth segmentation from Cone Beam Computed Tomography (CBCT) scans forms the
crux of dental diagnostics, orthodontic treatment planning, and dental
restoration procedures undertaken by stomatologists and dentists. The creation
of precise tooth models is incumbent upon accurate segmentation, which has
been traditionally a manual and tedious endeavor. The prevalence of low-
quality images in CBCT scans, characterized by complex artifacts and blurred
boundaries, complicates the task of segmentation, often leading to unclear
demarcations and challenges in extracting accurate semantic information. The
exigency for a robust, automated segmentation framework that can adeptly
handle these prevalent image quality issues while ensuring precise tooth
delineation is palpable.

Existing segmentation methodologies, notably U-Net [1] and its derivatives [2,
3], have shown effectiveness across various medical imaging domains. However,
when applied to CBCT tooth segmentation, especially in 2D scenarios, a
significant semantic gap between the encoder and decoder stages emerges,
detrimentally impacting the segmentation performance [4, 5]. Furthermore, the
inadequate fusion of semantic and boundary information, essential for
precisely distinguishing adjacent dental structures, compounds the challenge
citemea. This highlights the pressing need for a bespoke solution designed to
navigate the distinct hurdles associated with CBCT tooth imaging.

Addressing the outlined challenges, we introduce FDNet, a pioneering triple-
branch network that effectively decouples and processes semantic and boundary
features to overcome the challenges of CBCT tooth imaging. At its core are the
LF-Wavelet and the SAM encoder [6], tailored for semantic enrichment and
boundary delineation in images, respectively. LF-Wavelet operates within the
Low-frequency and Image Fusion (LIF) module, focusing on semantic information
amplification, while the SAM encoder, a pre-trained Vision Transformer, excels
at boundary feature extraction. Their synergistic operation mitigates the
semantic gap between encoder and decoder stages and enhances resilience
against common image defects. Through FDNet, we aim to mend the semantic gap
and bolster the framework against typical image defects, advancing
segmentation accuracy in CBCT tooth imaging.

The core contributions of our study are as follows:

  * •

We introduce FDNet, a distinctive triple-branch network that improves
segmentation outcome by effectively decoupling and integrating semantic and
boundary information.

  * •

FDNet utilizes LIF to lessen the semantic gap between the encoder and decoder
stages, which is a common drawback in preceding models.

  * •

FDNet achieves superior performance on DICE and IoU over existing state-of-
the-art models, providing a potential new baseline for tooth CBCT image
segmentation.

Fig. 1: The architecture of the proposed FDNet segmentation model.
Specifically, the LIF module fuses features from the $U_{w}$ and $U_{o}$
encoders before passing them to the CTB module.

##  2 METHODS

###  2.1 Architecture Overview

Our FDNet is a triple-branch model engineered for precision in tooth CBCT
image segmentation, as shown in Fig. 1, epitomizes our innovative approach to
feature decoupling, specifically targeting the nuanced intricacies of tooth
CBCT images. It harnesses the LF-Wavelet and original image data, channeling
these into the encoder before proceeding through the LIF’s Channel-wise Cross-
attention Transformer Block (CTB) module for enhanced feature integration,
capturing varied scale information effectively $(R_{j},j\in\\{1,2,3,4\\})$.
Subsequently, the refined features are seamlessly combined with the SAM
encoder’s output, feeding into the RDS decoder to accomplish superior
segmentation accuracy and boundary clarity.

###  2.2 LF-Wavelet Transform

In the initial phase of our FDNet framework, we adopt an LF-Wavelet to address
the semantic segmentation needs in CBCT tooth imaging. Accurate segmentation
hinges on the inclusion of low-frequency (LF) semantics, such as overall shape
and color consistency, which are crucial for understanding the global context
of dental structures. To remedy this, we apply a wavelet decomposition to
isolate and fuse the low-frequency details back into the original image,
enriching the data fed into our deep-learning model. This process imbues the
model with a fused representation of intricate semantic and structural
information. Figure 2 displays this process, allowing for a visual
appreciation of how LF-Wavelet Transform prioritizes.

Fig. 2: Comparison of low-frequency transformed and original images, with each
low-frequency counterpart on the top and its original image at the bottom.
Fig. 3: The results of segmentation using different methods. It shows our
model’s adept performance in achieving precise segmentation across four dental
conditions, including scenarios with complex artifacts and blurred tooth
boundaries.

The mathematical expression for this LF enhancement is articulated as follows:

$\displaystyle
X_{i}^{\prime}=W^{-1}\left(LF\left(W\left(X_{i}\right)\right)\right).$ (1)

Here, $X_{i}$ is the original input image, $W$ represents the wavelet
transform, $W^{-1}$ is its inverse, and $LF$ is the extraction of low-
frequency components. The enhanced image $X_{i}^{\prime}$ is thus a composite,
fortified with low-frequency information, enhancing the segmentation process.
By concatenating these low-frequency embeddings with the original image data,
we construct a composite input that capitalizes on the strengths of both
spectral domains in the first stage. Through this methodology, we mitigate the
potential loss of contextual semantics and augment the model’s proficiency in
translating global features into precise local patterns, which is pivotal for
refined and robust segmentation outcomes.

Moreover, the uniformity offered by the low-frequency details helps address
segmentation challenges, such as small data interruptions - gaps or holes
caused by high-frequency noise. Without this focus on low frequencies, such
disruption could lead to fragmented segmentation, adversely affecting dental
diagnostics and treatment planning.

###  2.3 SAM Encoder

After the low-frequency wavelet transformation, our FDNet framework
incorporates the SAM Encoder to enhance feature representation for precise
CBCT tooth segmentation. Pre-trained on the vast SA-1B segmentation dataset,
the SAM Encoder is adept at general feature extraction, a capability essential
for diverse imaging conditions.

The SAM Encoder’s edge lies in its ability to highlight boundary features,
which is critical for the clear delineation of adjacent dental structures
within CBCT images. While the SAM Encoder is a powerful tool for feature
extraction, its usage is optimized when guided by specific prompts, as it does
not inherently understand the segmentation goals or the semantics involved.

We address this by employing a channel-wise cross attention (CCA) module [7]
that synergistically combines the SAM Encoder’s boundary-focused embeddings
with the feature maps $R_{j}$ from the LIS module. The revised formula that
depicts this integration is:

$\displaystyle\left\\{\begin{array}[]{l}R_{j}{}^{\prime}=CCA\left(\operatorname{Concat}\left(\operatorname{Sam}\left(X_{i}\right),R_{j}\right),D_{j}\right),\\\
D_{j-1}=\operatorname{Conv}\left(\operatorname{Concat}\left(R_{j}^{\prime},D_{j}\right)\right).\quad
j=2,3,4\end{array}\right.$ (4)

Here, $R_{j}{}^{\prime}$ is the enhanced feature map produced by the CCA
module, which harmonizes the SAM Encoder’s boundary embeddings with CTB output
features $R_{j}$ and the decoder features $D_{j}$. The decoder feature map
$D_{j-1}$, resulting from this concatenation at levels $j=2,3,4$, embodies a
holistic representation of both global and local contextual details, enabling
an enriched feature landscape for the decoder’s operation.

The redefined concatenation process thus formulates a composite feature map
$S$ that integrates the SAM Encoder’s boundary delineation capabilities with
the semantic and structural insights derived from the LF-Wavelet and the
decoder features. Such an integrated feature map is instrumental in capturing
the nuances of CBCT tooth imaging, providing the RDS module, which serves as
the decoder in our FDNet framework, with a rich tapestry of features to
effectively discern and segment dental anatomy.

##  3 EXPERIMENT

###  3.1 Dataset

Given the advancements in three-dimensional dental CBCT technology and the
escalating importance of dentistry driven by economic upturn and an aging
population, the task of tooth segmentation has gained prominence. In light of
this, we leveraged a dataset provided by Hangzhou Dental Hospital,
encompassing 9,000 labeled CBCT images, to accomplish our segmentation
objectives. From this dataset, a subset of 288 images was randomly selected to
constitute our training set, while a distinct set of 400 images was earmarked
for testing purposes.

Table 1: Evaluation comparison on multiple tooth segmentation models.

| Dice (%) | IoU (%) | Recall | Precision
---|---|---|---|---
U-Net [1] | $81.69\pm 10.66$ | $69.45\pm 13.32$ | $82.14\pm 10.06$ | $82.20\pm 14.20$
U-Net++ [8] | $81.44\pm 10.75$ | $69.93\pm 13.69$ | $82.30\pm 10.41$ | $82.47\pm 14.31$
Attention U-Net [9] | $81.45\pm 10.66$ | $69.91\pm 13.58$ | $83.10\pm 9.49$ | $81.79\pm 14.78$
MedT [10] | $72.37\pm 15.97$ | $58.85\pm 17.20$ | $61.66\pm 17.98$ | $\mathbf{92.19\pm 13.51}$
Swin-Unet [11] | $75.23\pm 14.51$ | $62.19\pm 16.43$ | $75.90\pm 16.39$ | $77.44\pm 14.59$
UCTransNet [7] | $\underline{83.46\pm 9.60}$ | $\underline{73.88\pm 13.21}$ | $\mathbf{90.15\pm 10.86}$ | $81.02\pm 13.03$
FDNet (Ours) | $\mathbf{85.28\pm 9.78}$ | $\mathbf{75.23\pm 13.08}$ | $\underline{89.77\pm 10.99}$ | $\underline{82.67\pm 12.47}$

###  3.2 Main Results

To evaluate the performance of our proposed FDNet framework, we conducted a
comparative analysis against several advanced models, including U-Net [1],
U-Net++ [8], Attention U-Net [9], MedT [10], Swin-Uet [11], and UCTransNet [7]
on our tooth CBCT image dataset. All models were meticulously trained and
tested based on their official source code and preset parameters. We used a
test set comprising 400 patient tooth samples. The comparative visual results
are illustrated in Fig. 3, while the quantitative evaluation is in Table 1.

To highlight our model’s consistent performance across various patient cases
and imaging conditions, we computed the average and standard deviation from
multiple experiments. Our analysis revealed that FDNet achieved superior
segmentation accuracy, surpassing the second-ranked model, UCTransNet, in both
Dice and IoU metrics—indicators of segmentation correctness and the overlap
between predicted and actual segmentation, respectively. While FDNet’s average
Recall is marginally lower than that of UCTransNet, it showcases superior
Precision, with only MedT having a slightly higher value. However, it is
essential to note that MedT’s higher Precision is countered by its lower
Recall, suggesting a trade-off between these metrics. Our model’s balanced
performance in both metrics means it maintains a harmonious balance between
identifying relevant tooth structures (Recall) and minimizing false positives
(Precision). This balance is critical in dental CBCT segmentation, where both
over-segmentation and under-segmentation can have detrimental implications for
subsequent dental treatment planning.

###  3.3 Ablation Study

An ablation study was orchestrated to dissect the individual contributions of
our core modules: SAM encoder, LIF and LW-Wavelet, in the realm of tooth
segmentation. As illustrated in Table 2, this study illuminated the pivotal
roles of each module in augmenting the segmentation quality. It underscored
that the confluence of semantic information through LW-Wavelet, the robust
feature extraction by the SAM encoder, and the fusion functionality of the LIF
module synergistically drive superior segmentation outcomes.

Table 2: Ablation Study Outlining the Impact of SAM Encoder, LIF, and LF-
Wavelet on Segmentation Performance.

Methods | Dice (%) | IoU (%) | Recall | Precision
---|---|---|---|---
FDNet (w/o SAM) | 85.24 | 75.14 | 87.56 | $\boldsymbol{84.77}$
FDNet (w/o LIF) | 83.94 | 74.01 | 87.43 | 82.31
FDNet (w/o LW) | 84.33 | 74.12 | 86.47 | 83.98
FDNet (Ours) | $\boldsymbol{85.64}$ | $\boldsymbol{75.72}$ | $\boldsymbol{88.86}$ | $84.38$

##  4 CONCLUSION

In this work, we introduce FDNet, an innovative approach for tooth
segmentation in CBCT imaging, addressing the inherent challenges presented by
diverse dental conditions, indistinct tooth boundaries, and semantic paucity
in CBCT images. FDNet’s architectural ingenuity lies in its unique feature
decoupling strategy. By incorporating the Low-Frequency Wavelet Transform
within our framework, we enhance the semantic richness by capturing the global
structural details from low-frequency components. Concurrently, the SAM
encoder extracts boundary features critical for distinguishing adjacent dental
structures. This strategic fusion of components facilitates a semantic bridge
between the encoding and decoding phases, markedly elevating the segmentation
accuracy, particularly in the face of complex dental scenarios. Our extensive
testing underscores FDNet’s superior performance over leading-edge models,
highlighting its capacity to revolutionize automated tooth segmentation across
various patient dental conditions and positioning it as a versatile and
reliable framework in the realm of dental diagnostics and orthodontic
treatment planning.

##  5 COMPLIANCE WITH ETHICAL STANDARDS

This study was performed in line with the principles of the Declaration of
Helsinki. Approval was granted by the Ethics Committee of Lishui University
Medical School (N0.2022YR014).

##  6 ACKNOWLEDGMENTS

This research is supported in part by the National Natural Science Foundation
of China [(No.62206242) and (No.62201323)].

## References

  * [1] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18. Springer, 2015, pp. 234–241.
  * [2] Sixiao Zheng, Jiachen Lu, Hengshuang Zhao, Xiatian Zhu, Zekun Luo, Yabiao Wang, Yanwei Fu, Jianfeng Feng, Tao Xiang, Philip HS Torr, et al.,  “Rethinking semantic segmentation from a sequence-to-sequence perspective with transformers,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2021, pp. 6881–6890.
  * [3] Yunxiang Li, Shuai Wang, Jun Wang, Guodong Zeng, Wenjun Liu, Qianni Zhang, Qun Jin, and Yaqi Wang,  “Gt u-net: A u-net like group transformer network for tooth root segmentation,”  in Machine Learning in Medical Imaging: 12th International Workshop, MLMI 2021, Held in Conjunction with MICCAI 2021, Strasbourg, France, September 27, 2021, Proceedings 12. Springer, 2021, pp. 386–395.
  * [4] Nabil Ibtehaz and M Sohel Rahman,  “Multiresunet: Rethinking the u-net architecture for multimodal biomedical image segmentation,”  Neural networks, vol. 121, pp. 74–87, 2020.
  * [5] Yanwei Pang, Yazhao Li, Jianbing Shen, and Ling Shao,  “Towards bridging semantic gap to improve semantic segmentation,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 4230–4239.
  * [6] Alexander Kirillov, Eric Mintun, Nikhila Ravi, Hanzi Mao, Chloe Rolland, Laura Gustafson, Tete Xiao, Spencer Whitehead, Alexander C. Berg, Wan-Yen Lo, Piotr Dollár, and Ross Girshick,  “Segment anything,”  arXiv:2304.02643, 2023.
  * [7] Haonan Wang, Peng Cao, Jiaqi Wang, and Osmar R. Zaiane,  “Uctransnet: Rethinking the skip connections in u-net from a channel-wise perspective with transformer,”  Proceedings of the AAAI Conference on Artificial Intelligence, vol. 36, no. 3, pp. 2441–2449, Jun. 2022.
  * [8] Zongwei Zhou, Md Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang,  “Unet++: A nested u-net architecture for medical image segmentation,”  in Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support: 4th International Workshop, DLMIA 2018, and 8th International Workshop, ML-CDS 2018, Held in Conjunction with MICCAI 2018, Granada, Spain, September 20, 2018, Proceedings 4. Springer, 2018, pp. 3–11.
  * [9] Ozan Oktay, Jo Schlemper, Loic Le Folgoc, Matthew Lee, Mattias Heinrich, Kazunari Misawa, Kensaku Mori, Steven McDonagh, Nils Y Hammerla, Bernhard Kainz, et al.,  “Attention u-net: Learning where to look for the pancreas,”  arXiv preprint arXiv:1804.03999, 2018.
  * [10] Jeya Maria Jose Valanarasu, Poojan Oza, Ilker Hacihaliloglu, and Vishal M. Patel,  “Medical transformer: Gated axial-attention for medical image segmentation,”  in Medical Image Computing and Computer Assisted Intervention – MICCAI 2021, Cham, 2021, pp. 36–46, Springer International Publishing.
  * [11] Hu Cao, Yueyue Wang, Joy Chen, Dongsheng Jiang, Xiaopeng Zhang, Qi Tian, and Manning Wang,  “Swin-unet: Unet-like pure transformer for medical image segmentation,”  in European conference on computer vision. Springer, 2022, pp. 205–218.