# Learned Mappings for Targeted Free Energy Perturbation between Peptide
Conformations

Soohaeng Yoo Willow  Department of Chemistry, Illinois Institute of
Technology, Chicago, Illinois, 60616  Lulu Kang  Department of Applied
Mathematics, Illinois Institute of Technology, Chicago, Illinois, 60616  David
D. L. Minh<EMAIL_ADDRESS>Department of Chemistry, Department of Biology, and
Center for Interdisciplinary Scientific Computation, Illinois Institute of
Technology, Chicago, Illinois, 60616

###### Abstract

Targeted free energy perturbation uses an invertible mapping to promote
configuration space overlap and the convergence of free energy estimates.
However, developing suitable mappings can be challenging. Wirnsberger et al.
(2020) demonstrated the use of machine learning to train deep neural networks
that map between Boltzmann distributions for different thermodynamic states.
Here, we adapt their approach to free energy differences of a flexible bonded
molecule, deca-alanine, with harmonic biases with different spring centers.
When the neural network is trained until “early stopping” - when the loss
value of the test set increases - we calculate accurate free energy
differences between thermodynamic states with spring centers separated by 1 Å
and sometimes 2 Å. For more distant thermodynamic states, the mapping does not
produce structures representative of the target state and the method does not
reproduce reference calculations.

##  I Introduction

Free energy calculations are a powerful tool that are increasingly used to
design materials [1] and drugs [2]. Accurately calculating free energy
differences between a pair of thermodynamic states often entails performing
simulations of a multiple intermediate states along a thermodynamic process
connecting the states [3]. Simulating multiple intermediates promotes overlap
between the most important configuration space in neighboring thermodynamic
states, a requirement for the convergence of free energy estimates [4, 5, 6].

Circumventing the simulation of these intermediates could improve simulation
accuracy and reduce resource requirements. Simulations of intermediates can
have problems that degrade the accuracy of free energy estimates. For example,
alchemical processes in which particles are created or destroyed are known to
suffer from pathologies such as the end-point catastrophe - there can be poor
configuration space overlap between states where particles nearly appear or
disappear, leading to spurious discontinuities in free energy estimates along
a thermodynamic process - and the presence of artificial energy minima [7].
Moreover, intermediate states are usually of no particular scientific
interest. Avoiding them could lower resource requirements, allowing more
scientists to access the tool and reducing the consumption of computing and
energy resources.

Targeted free energy perturbation (TFEP) [8] is an approach that can bypass
the simulation of intermediate states. Rather than simulating many
intermediates, TFEP achieves configuration space overlap via an invertible
mapping between the end states. In 2002, Jarzynski [8] formulated TFEP as an
extension of the classic free energy perturbation (FEP) identity [9]. In FEP,
configurations are sampled from a single state. The free energy is an
exponential average of the difference between the potential energy of each
sampled configuration in the target and sampled state (Equation 2). On the
other hand, instead of using the same configuration for both potential energy
evaluations, TFEP uses the potential energy of the _mapped_ configuration in
the target state. Along with the Jacobian of the mapping, the potential energy
of the mapped configuration is incorporated into a generalized work (Equation
6) in the exponential average (Equation 5).

There were several precedents to the formalism of Jarzynski [8]. In 1985,
Voter [10] performed free energy calculations with volume-preserving
translations. In 1995, Severance  _et al._ [11] considered transformations of
bond lengths and angles based on harmonic force field parameters as well as
the rotation of a dihedral angle. Although they did not explicitly mention
Jacobians, all of their example mappings had a Jacobian determinant of unity;
the term was not needed to reproduce analytical free energies. TFEP can be
considered as a generalization of the celebrated Jarzynski identity [12, 13]
that was published in 1997; in the Jarzynski identity, mapping is performed by
a nonequilibrium switching process [14]. In 2000, Miller and Reinhardt [15]
proposed that such switching processes can incorporate scaling
transformations.

There have been several extensions and generalizations to TFEP. Hahn and Then
[16] implemented the bidirectional estimator that Jarzynski [8] suggested in
the original TFEP paper. Soon after Jarzynski [8], Meng and Schilling [17]
described a multistate generalization of TFEP to arbitrary statistical
distributions, which Paliwal and Shirts [18] applied to molecular systems. In
addition to free energies, the multistate generalization is capable of
estimating arbitrary expectation values [17, 18]. Subsequently, TFEP and its
extensions have been applied to spin systems and other statistical
distributions, but our focus here will be on molecular systems.

Although TFEP is a compelling concept, it has been difficult to apply because
of the lack of suitable mapping functions for molecular systems. Mappings have
been developed for expanding a cavity in a fluid [8], inserting a particle
into a fluid [16], the interconversion of water models [18], and between
crystals at different temperature [19] and volume [20]. However, the mappings
are not broadly generalizable; human intuition and creativity have been
required for every distinct application.

While human intelligence is still important, we are now in the age of
artificial intelligence. Apropos to the emerging era, a few groups have used
deep learning to train maps between molecular systems for TFEP. Wirnsberger
_et al._ [21], from the Google subsidiary DeepMind, who also developed
AlphaFold [22], learned mappings to ensembles of fluids containing a solute
with different radii. Ding and Zhang [23] trained mappings between a Boltzmann
and tractable distribution - for which the normalized density is known and
independent and identically distributed samples can be easily generated - to
compute absolute free energies for different conformations of di-alanine and
temperatures of deca-alanine. Subsequently, they applied the same approach to
a host-guest system to compute binding free energies [24]. Rizzi  _et al._
[25] trained neural networks for mapping between two levels of quantum theory
for a simple chemical reaction in the gas phase. Wirnsberger  _et al._ [26]
learned mappings between a lattice model with random perturbations to
different phases of Lennard-Jones systems.

In recent years, the use of learned mappings in molecular simulation beyond
TFEP has been pioneered by Frank Noé and collaborators. In 2019, Noé  _et al._
[27] introduced Boltzmann generators, which use a reference state and a
learned mapping to generate samples that may be reweighed to the Boltzmann
distribution. They demonstrated the method on generating samples of a bistable
dimer in fluid (initially described in an paper coauthored by the
corresponding author [28]) and a small protein. While their publication
primarily focused on sampling opposed to free energies, they did report
computing free energy differences between independent Boltzmann generators
using the average (opposed to the exponential average from TFEP [8]) of the
generalized work. Subsequently, Sbailò  _et al._ [29] described the use of
learned mappings for Monte Carlo moves, as demonstrated in the bistable dimer
system. Finally, Invernizzi  _et al._ [30] used learned mappings for replica
exchange, performing simulations of di-alanine and tetra-alanine. Besides Noé
and collaborators, Mahmoud  _et al._ [31] have developed a hierarchical
sampling procedure to generate samples of a medium-sized (106 residue)
protein.

Here, we apply learned mappings to free energy differences between different
conformations of deca-alanine. Our approach is similar to Wirnsberger  _et
al._ [21] but the system is qualitatively different. Wirnsberger  _et al._
[21] modeled a cavity in a fluid of neutral monatomic molecules, which has a
simpler potential energy function due to the lack of bonded (bond length, bond
angle, and torsion) and Coulomb interactions. Our system is similar to Ding
and Zhang [23], but the approach is qualitatively different. While Ding and
Zhang [23] learned mappings between Boltzmann and tractable distributions, we
learn mappings between Boltzmann distributions corresponding to different
thermodynamic states (different harmonic biases). Mapping to a tractable
distribution has an advantage in that the tractable distribution can be
selected to have significant configuration space overlap with the unmapped
molecular distribution. Moreover, it is generally faster to evaluate the
probability density and its derivatives for a tractable distribution than the
potential energy of a molecular system. On the other hand, directly mapping
between pairs of molecular distributions can permit a subset of coordinates to
be perturbed. For processes in which changes are limited to a subset of
coordinates, such as those in the focused confinement method [32, 33], such
mappings may be particularly fruitful. Hence both approaches to mapping may
find separate niches.

The remainder of the paper is as follows. In Section II, we review the theory
of TFEP and the loss function. We then describe computational methods (Section
III), including the deca-alanine model system, the neural network
architecture, training procedure, and free energy estimators. We then report
results along with discussion of their implications (Section IV), and lastly
describe our conclusions (Section V).

##  II Theory

###  II.1 Targeted free energy perturbation

Let us first define notation. For a thermodynamic state $A$, a configuration
$\bm{x}$ has the equilibrium probability density
$\rho_{A}(\bm{x})=\exp\left[-\beta U_{A}(\bm{x})\right]/Z_{A}$. In this
expression, $\beta$ is the inverse of the temperature and Boltzmann’s constant
$k_{B}$, $\beta=\left(k_{B}T\right)^{-1}$, $U_{A}(\bm{x})$ is the potential
energy of the configuration in state $A$, and $Z_{A}=\int\exp\left[-\beta
U_{A}(\bm{x})\right]d\bm{x}$ is the configurational integral of the state, an
integral over all space. Analogous definitions apply to a second thermodynamic
state $B$. The objective of our calculation is the Helmholtz free energy
difference between states $A$ and $B$, defined as,

$\Delta F=F_{B}-F_{A}=-\beta^{-1}\ln\left(\frac{Z_{B}}{Z_{A}}\right).$ (1)

In principle, a free energy difference $\Delta F$ may be calculated based on
the classic free energy perturbation identity [9],

$\displaystyle\mathbb{E}_{A}\left[e^{-\beta\Delta U}\right]=e^{-\beta\Delta
F},$ (2)

where $\mathbb{E}_{A}[\cdot]$ represents an expectation value in the state $A$
and $\Delta U(\bm{x})=U_{B}(\bm{x})-U_{A}(\bm{x})$ is the potential energy
difference of the _same_ configuration between the two states. (For notational
simplicity, the configuration dependence is implicit within the expectation
value). Specifically, a free energy difference may be estimated using the
sample mean of the observable $\exp\left(-\beta\Delta U(\bm{x})\right)$ over
configurations drawn from state $A$. As exponential averages suffer from
finite-sample bias [34], better estimation performance may be achieved using
samples from both $A$ and $B$ [35] or from multiple states (including $A$ and
$B$) [36]. Regardless of the estimator, free energy calculations do not
reliably converge to the true values unless there is overlap between the
highest-probability regions of configuration space for the pair (or series) of
states [4, 5, 6].

Configuration space overlap may be improved by a bijective (one-to-one) and
invertible mapping (Figure 1). A mapping is simply a function that changes a
configuration $\bm{x}$ into a new configuration $\bm{x}^{\prime}=M(\bm{x})$.
For an invertible mapping, the inverse function transforms $\bm{x}^{\prime}$
into $\bm{x}$, $\bm{x}=M^{-1}(\bm{x}^{\prime})$. Applying a mapping to state
$A$ creates a new state $A^{\prime}$. The probability of the new state is
given by the change of variables formula,

$\displaystyle\rho_{A^{\prime}}(M(\bm{x}))$ $\displaystyle=$
$\displaystyle\frac{\rho_{A}(\bm{x})}{|J_{M}(\bm{x})|},$ (3)

where $|J_{M}(\bm{x})|$ is the determinant of the Jacobian of the mapping.
Similarly, the inverse mapping can be applied to the state $B$ to create a new
state $B^{\prime}$ with density,

$\displaystyle\rho_{B^{\prime}}(M^{-1}(\bm{x}))$ $\displaystyle=$
$\displaystyle\frac{\rho(\bm{x})}{|J_{M^{-1}}(\bm{x})|}.$ (4)

A mapping can be defined such that $A^{\prime}$ is close to $B$ and conversely
$B^{\prime}$ is close to $A$ (Figure 1).

Figure 1: Schematic showing how mapping can improve configuration space
overlap. Each shape represents the highest-probability regions of the
respective states.

Mappings may be incorporated into free energy calculations via a
generalization of Equation 2, TFEP [8],

$\displaystyle\mathbb{E}_{A}\left[e^{-\beta\Phi_{F}}\right]=e^{-\beta\Delta
F}.$ (5)

Equation 5 replaces the $\Delta U$ in Equation 2 with a generalized work,

$\displaystyle\Phi_{F}(\bm{x})$ $\displaystyle=$ $\displaystyle
U_{B}(M(\bm{x}))-U_{A}(\bm{x})-\beta^{-1}\log|J_{M}(\bm{x})|,$ (6)

in which the configurations in state $A$ and $B$ may be different. We refer to
this term as work because it is a generalization [14] of the nonequilibrium
work from Jarzynski’s identity [12, 13]. The work includes a subscript $F$ to
denote a “forward” mapping from $A$ to $A^{\prime}$. For the “reverse” mapping
from $B$ to $B^{\prime}$, the work includes the subscript $R$,

$\displaystyle\Phi_{R}(\bm{x})$ $\displaystyle=$ $\displaystyle
U_{A}(M^{-1}(\bm{x}))-U_{B}(\bm{x})-\beta^{-1}\log|J_{M^{-1}}(\bm{x})|.$ (7)

The designation of a direction of forward or reverse is arbitrary. In the case
of an identity mapping, the work is the conventional potential energy
difference and Equation 5 reduces to Equation 2. However, a mapping that
increases configuration space overlap can improve the convergence of free
energy estimates. Indeed, if a mapping perfectly transforms between $A$ and
$B$,

$\displaystyle\frac{\rho_{A^{\prime}}(M(\bm{x}))}{\rho_{B}(\bm{x})}=\frac{\rho_{B^{\prime}}(M^{-1}(\bm{x}))}{\rho_{A}(\bm{x})}=1,$
(8)

then only a single sample is needed to estimate the free energy difference [8,
25]. In addition to the sample mean estimator based on Equation 5, mappings
have also been incorporated into bidirectional [16] and multistate [17, 18]
estimators for free energies and expectation values.

###  II.2 Loss functions

While it can be challenging to _design_ mappings that improve configuration
space overlap and accelerate the convergence of free energy calculations,
Wirnsberger  _et al._ [21] demonstrated that machine learning may be used to
_train_ mappings. Inspired by configuration space overlap, they developed loss
functions for their machine learning models based on a statistical distance
between mapped and targeted distributions, the Kullback-Leibler (KL)
divergence [37]. The quality of the forward mapping is quantified by the KL
divergence between the mapped density $\rho_{A^{\prime}}(M(\bm{x}))$ and the
target density $\rho_{B}(\bm{x})$,

$\displaystyle D_{\mathrm{KL}}[\rho_{A^{\prime}}||\rho_{B}]$ $\displaystyle=$
$\displaystyle\beta(\mathbb{E}_{A}[\Phi_{F}]-\Delta F).$ (9)

Conversely, the quality of the reverse mapping is described by the KL
divergence between $\rho_{B^{\prime}}(M^{-1}(\bm{x}))$ and $\rho_{A}(\bm{x})$,

$D_{\mathrm{KL}}[\rho_{B^{\prime}}||\rho_{A}]=\beta(\mathbb{E}_{B}[\Phi_{R}]+\Delta
F).$ (10)

As $\beta$ and $\Delta F$ are constants with respect to mapping, they do not
need to be included in loss functions. Hence, a loss function suitable for
training only the forward mapping is,

$\mathcal{L}_{F}=\mathbb{\bar{E}}_{A}[\Phi_{F}],$ (11)

where $\mathbb{\bar{E}}_{A}[\cdot]$ is a sample mean estimate of an
expectation value in state $A$. Wirnsberger  _et al._ [21] also introduced and
recommended a bidirectional loss function,

$\mathcal{L}=\mathbb{\bar{E}}_{A}[\Phi_{F}]+\mathbb{\bar{E}}_{B}[\Phi_{R}],$
(12)

which we use in the present work.

##  III Computational Methods

###  III.1 Model system

Alanine deca-peptide (ACE-(ALA)9-NME) was selected as a model system. The
system was modeled in gas phase at room temperature, where it is a stable
$\alpha$ helix. The initial structure was modified from Nguyen and Minh [38]
by hand to include the end caps ACE and NME. The system has 102 atoms.
Parameters from the AMBER ff14SB force field [39] were used. Harmonic
restraints were added to the carbon of the first peptide bond and the nitrogen
atom of the last peptide bond, leading to the total potential energy,

$U_{\lambda}(x)=U_{0}(x)+k||\bm{x}_{C}-\bm{x}_{C}^{0}||^{2}+k||\bm{x}_{N}-\bm{x}_{N}^{\lambda}||^{2},$
(13)

where $U_{0}(x)$ is the standard gas-phase AMBER potential energy. In the
harmonic potential, $k=50$ kJ/mol/Å2 is the force constant, $\bm{x}_{C}$ and
$\bm{x}_{N}$ are coordinates of the restrained carbon and nitrogen atoms,
respectively, and $||\cdot||$ represents the Euclidean norm, the length of the
vector. The spring center $\bm{x}_{C}^{0}$ is at the origin and
$\bm{x}_{N}^{\lambda}=(0,0,\lambda)$, where
$\lambda\in\left\\{14,14.2,14.4,...,20.6,20.8,21\right\\}$ Å.

For each $\lambda$, an ensemble of structures was generated with a molecular
dynamics simulation. Isothermal molecular dynamics simulations were performed
using the Langevin integrator (LangevinIntegrator) in OpenMM 7.7 [40] at a
temperature of $T=300$ K with a time step of 1 fs for a duration of 20 ns. The
last 10000 configurations (every 1 ps for the last 10 ns) were saved.

###  III.2 Neural network architecture

Mappings were performed using real-valued non-volume preserving (real NVP)
transformations [41] between Cartesian coordinates. Real NVP uses a deep
neural network composed of a stack of simple bijections that include scaling
and translation. It is especially well-suited for TFEP because it is stably
invertible and because the Jacobian of the transformation can be efficiently
computed.

To adapt real NVP to molecular systems in Cartesian coordinates, we used
separate coupling layers for $x$, $y$, and $z$ dimensions (Figure 2). The
first coupling layer $M_{x}$ transforms the $x$ coordinates while leaving the
$y$ and $z$ coordinates invariant,

$\displaystyle
M_{x}(\\{x,y,z\\})=\left\\{x^{\prime}=xe^{s(y,z)}+t(y,z),y,z\right\\}$ (14)

where $s$ and $t$ are functions encoded in the neural network. The inverse of
this mapping is
$M_{x}^{-1}(\\{x^{\prime},y,z\\})=\left\\{x=\left[x^{\prime}-t(y,z)\right]e^{-s(y,z)},y,z\right\\}$.
Jacobian determinants of the mappings $M_{x}$ and $M_{x}^{-1}$ are,

$\displaystyle\log|J_{M_{x}}|=\sum_{i=1}^{N_{\mathrm{atoms}}}s_{i}(y,z)$ (15)
$\displaystyle\log|J_{M_{x}^{-1}}|=-\sum_{i=1}^{N_{\mathrm{atoms}}}s_{i}(y,z),$
(16)

where the sum is over the atoms and the limit the total number of atoms
$N_{\mathrm{atoms}}$. Similarly, the second coupling layer $M_{y}$ and the
third coupling layer $M_{z}$ transform the $y$ and $z$ coordinates with
$(x,z)$ and $(x,y)$ invariant, respectively. Each coupling layer
$M_{k}:\mathbb{T}^{3N}\rightarrow\mathbb{T}^{3N}$ used the same architecture
but different learnable parameters $\theta_{k}$. The coupling layers consisted
of 10 hidden layers with 128 dimensions (compared to 102 atoms).

(a) |
---|---
|
(b) |
|
Figure 2: Schematic of the neural network architecture. (a) In the first
coupling layer $M_{x}$, $y$ and $z$ are left unchanged and the learnable
parameters $\theta$ for $s$ and $t$ functions are trained to transform the $x$
coordinates. (b) The entire real NVP network consists of $M=M_{z}\circ
M_{y}\circ M_{x}$.

The overall mapping $M$ of the real NVP transformation is given as

$\displaystyle M=M_{z}\circ M_{y}\circ M_{x},$ (17)

where $\circ$ means that the output of the right side is input for the left
side. Its inverse is,

$\displaystyle M^{-1}=M_{x}^{-1}\circ M_{y}^{-1}\circ M_{z}^{-1}.$ (18)

The Jacobian determinant of the mapping $M$ and the inverse mapping $M^{-1}$
are,

$\log|J_{M}|=\sum_{\nu\in\\{x,y,z\\}}\log|J_{M_{\nu}}|,\\\
\log|J_{M^{-1}}|=\sum_{\nu\in\\{x,y,z\\}}\log|J_{M_{\nu}^{-1}}|$ (19)

The mapping was implemented in JAX [42].

###  III.3 Training

The real NVP map was initialized as an identity map. This was achieved by
initializing the transformation $t$ and scaling $s$ factors at zero, such that
$x^{\prime}=x\exp(0)+0=x$. Randomness in initial conditions was introduced by
drawing all other initial parameters from the standard normal distribution.

The real NVP neural network was trained by minimizing the value of the loss
function (Equation 12). As the loss function (Equation 12) is based on the
molecular mechanics energies of mapped configurations, we implemented the
AMBER force field in JAX. This implementation is freely available at
https://github.com/swillow/jax_amber. To ensure that mapped $\bm{x}_{N}$ are
similar to $\bm{x}_{N}^{\lambda}$, we used $k=200$ kJ/mol/Å2 instead of $k=50$
kJ/mol/Å2 in the loss function. Data were divided into a training set (80%)
and a test set (20%). Out of the last 10 ns of simulation, 2000 configurations
(from 10 ns to 12 ns) were designated as test sets while the remaining 8000
configurations (from 12 ns to 20 ns) were designated as training sets. The
training loss was minimized using the adaptive moment estimation (Adam)
optimizer [43] with a learning rate of 1.0e-4.

We performed different amounts of training in initial and later calculations.
Initially, training was performed for 150000 steps. We compared results from
mappings obtained after different amounts of training:

  1. 1.

Early stopping, with the minimum loss of the test set.

  2. 2.

L $\sim$ 0, where the loss of the training set was closest to zero.

  3. 3.

Complete, after all 150000 training steps.

Subsequently, after observing superior free energy estimation from early
stopping, we stopped training once the loss of the test set in the last 200
steps was over $3.0~{}k_{B}T$ greater than the minimum.

###  III.4 Free energy estimation

Reference free energies as a function of $\lambda$ were calculated the
multistate Bennett acceptance ratio (MBAR) estimator [36] without mapping.
They made us of the entire ensemble of equilibrated structures and energies.
Mapped free energies were estimated based on data from pairs of $\lambda$
using the generalization of the Bennett Acceptance Ratio [35] described by
Hahn and Then [16]. Thus, the mapped calculations were based on a small subset
of the structures and energies used in the reference calculations.

##  IV Results and Discussion

###  IV.1 The free energy landscape of deca-alanine has multiple barriers

The free energy of deca-alanine as as function of the spring position
generally increases between $\lambda=14$ and $\lambda=21$ Å (Figure 3). There
are local minima near 14, 16, and 18 Å that are separated by barriers. The
barriers correspond to local peaks in the potential energy, suggesting that
crossing them corresponds to breaking intramolecular interactions such as
hydrogen bonds. As the potential energy barriers are larger than the free
energy barriers, they are partially compensated for by entropic increases.

Figure 3: The reference free energy surface and average potential energy
difference compared to the state with $\lambda=14$ Å.

###  IV.2 Overlap decreases as $\lambda$ separation increases

While there may be some overlap between thermodynamic states with $\lambda$
separated by 1 Å, there does not appear to be any overlap between states where
it is separated by 2 Å or more ($\Delta\lambda\geq 2$ Å). The overlap between
thermodynamic states may be evaluated based on the histograms of $\Phi_{F}$
and $-\Phi_{R}$. Without mapping (or an identity map), the work is simply the
potential energy difference of the same configuration in the pair of states.
For both separations, distributions of the potential energy difference are
unimodal. For $\Delta\lambda=1$ Å, there is some overlap between the tails of
the distributions of $\Phi_{F}$ and $-\Phi_{R}$. For $\Delta\lambda=2$ Å, the
distributions of potential energy differences are broader and further
separated such that there is no evident overlap. Due to the lack of overlap
between these distributions, the free energy difference cannot be accurately
calculated (without a mapping that improves overlap).

Figure 4:  Normalized histograms of work without mapping for $\lambda_{A}=20$
Å and $\lambda_{B}=21$ Å (solid line) and for $\lambda_{A}=19$ Å and
$\lambda_{B}=21$ Å (dashed line). For both pairs, $\beta\Phi_{F}$ are colored
green with forward diagonal hatches and $-\beta\Phi_{R}$ are colored red with
backwards diagonal hatches.

Even for $\Delta\lambda\leq 2$ Å, our unmapped work distributions are more
distinct from one another than those reported in Wirnsberger  _et al._ [21].
Wirnsberger  _et al._ [21] focused on a solute in a Lennard-Jones fluid, in
which the states $A$ and $B$ were generated by changing the solute radii from
$R_{A}=2.5974~{}\sigma$ to $R_{B}=2.8444~{}\sigma$. For $\sigma=3.15$ Å, the
radius of the solute in state $A$ is 8.18 Å and that of the solute in state
$B$ is 8.96 Å, for a difference of 0.78 Å. The distance between unmapped work
distributions was about
$\beta(\mathbb{E}_{A}[\Phi_{F}]+\mathbb{E}_{B}[\Phi_{R}])\sim 40$. In
contrast, we observe
$\beta(\mathbb{E}_{A}[\Phi_{F}]+\mathbb{E}_{B}[\Phi_{R}])\sim 200$ and 600 for
1 Å and 2 Å, respectively.

###  IV.3 Loss values for the training and test set diverge after overfitting

In all of the learned mappings, the loss values for the training and test set
exhibit different behavior from other another (Figure 5). At the start of
training, both loss values quickly decrease. With additional training steps,
the training loss continues to trend downward and even becomes negative. On
the other hand, the test loss approaches a minimum and sharply increases after
around 2000 steps.

Figure 5: Loss values for the mapping between states with $\lambda_{A}=20$ Å
and $\lambda_{B}=21$ Å for training (blue) and test (red) data. The inset is a
subset of the range emphasizing the minimum of the test loss.

A deviation between the training and test loss is common in machine learning,
including the learned mapping in Wirnsberger  _et al._ [21] and Rizzi  _et
al._ [25] \- and is an indicator of overfitting to the training set. For this
reason, Wirnsberger  _et al._ [21] adopted an early stopping criterion. Rizzi
_et al._ [25] suggested that overtraining leads to a systematic bias in the
free energy estimate and thus they estimated free energies based on evaluation
sets separate from their training sets.

###  IV.4 Mapping increases the overlap between forward and reverse work
distributions

With early stopping, mapping improves overlap between the work distributions
and the feasibility of free energy estimation (Figure 6a). At this mapping,
which provides a minimum loss for the test set, there is considerable overlap
between the distributions of $\Phi_{F}$ and $-\Phi_{R}$. Trends in the average
work and free energy difference are consistent with applying Jensen’s
inequality to the expectation values in the TFEP expression, Equation 5. We
anticipate that
$\mathbb{E}_{A}[\Phi_{F}]>-\beta^{-1}\ln\mathbb{E}_{A}[e^{-\beta\Phi_{F}}]=\Delta
F$. Similarly, $\mathbb{E}_{B}[\Phi_{R}]>-\Delta F$. Hence, expectation values
of the work bound the free energy according to
$-\mathbb{E}_{B}[\Phi_{R}]<\Delta F<\mathbb{E}_{A}[\Phi_{F}]$. While these
inequalities are true for expectation values, they do not necessarily hold for
estimates based on finite sample sizes, which are not only imprecise but
systematically biased [34]. With early stopping, however, the order of
expectation values and the free energy follows the anticipated trend.
Consequently, the free energy difference, which is at the value of $\Phi$ in
which the densities are equal [16], is straightforward to identify. However,
further training leads to unexpected trends in the distributions of
generalized energies.

Figure 6: Effect of stopping criterion on work distributions for training data
between $\lambda_{A}=20$ Å and $\lambda_{B}=21$ Å. Normalized histograms are
shown for mapping functions from (a) early stopping, (b) when the loss is
closest to zero, and (c) after 150000 training steps. $\beta\Phi_{F}$ are
colored green with forward diagonal hatches and $-\beta\Phi_{R}$ are colored
red with backwards diagonal hatches. The reference free energy difference
$\Delta F$ obtained via the MBAR estimator is represented by the vertical
dashed line. A similar figure for the testing data is available as Figure S1
in the Supplementary Material.

When the loss is close to zero, we observe that both estimated averages are
less than the free energy difference,
$\mathbb{\hat{E}}_{B}[-\Phi_{R}]\sim\mathbb{\hat{E}}_{A}[\Phi_{F}]<\Delta F$
(Figure 6b). With a perfect mapping, we would anticipate that
$\mathbb{E}_{B}[-\Phi_{R}]=\mathbb{E}_{A}[\Phi_{F}]=\Delta F$. The deviation
of the expectation values from the free energy difference demonstrates that
perfect mapping is not achieved. The potential energy distributions of the
mapped configurations mostly overlap with the targeted (unmapped) potential
energy distribution (Figure S2 of the Supplementary Material). This overlap
between potential energy distributions shows that the source of the observed
inequality is from the Jacobian; $\log|J_{M}|\neq 0$ and $\log|J_{M^{-1}}|\neq
0$.

Finally, at the end of training, when the training loss values are negative
($\mathcal{L}<0$), the forward work values are actually much smaller than the
reverse work values. As shown in Figure 6 (c), the order of averages and free
energies is flipped such that $\mathbb{\hat{E}}_{A}[\Phi_{F}]\ll\Delta
F\ll\mathbb{\hat{E}}_{B}[-\Phi_{R}]$. This flipping occurs because mapped
energies $U_{B}(M(\bm{x}_{A}))$ and $U_{A}(M^{-1}(\bm{x}_{B}))$ are lower than
the targeted energies $U_{B}(\bm{x}_{B})$ and $U_{A}(\bm{x}_{A})$ (Figure S2
of the Supplementary Material). Molecular mechanics potential energies,
especially Lennard-Jones repulsion terms, are very sensitive to small changes
in coordinates; there are many opportunities for energy optimization relative
to samples from the Boltzmann distribution. While these low energies
demonstrate that the optimizer is working well, the generated structures are
not representative of the target distribution. Moreover, the overlap between
work distributions is reduced compared to when the loss is close to zero.

Except with early stopping, trends in the work distributions of the test set
do not match that of the training set (Figure S1 in the Supplementary
Material.) With early stopping, the training and test set have comparable work
distributions. When $\mathcal{L}\sim 0$, the work distributions for the test
set are similar to the distributions after early stopping;
$\rho(\beta\Phi_{F})$ and $\rho(-\beta\Phi_{R})$ of the test set do not have
nearly complete overlap. For complete training, the anticipated order of
averages and free energies $\mathbb{\hat{E}}_{B}[-\Phi_{R}]<\Delta
F<\mathbb{\hat{E}}_{A}[\Phi_{F}]$ is preserved, but the distributions are very
broad. The corresponding large test loss indicates that there is overfitting
to the training set.

###  IV.5 Mapping enables accurate free estimation for smaller $\Delta\lambda$

For most pairs of states where $\Delta\lambda=1$ Å, TFEP is able to accurately
reproduce reference free energy differences. Mapped estimates are most
accurate between pairs of states with $\Delta F$ less than 2 kJ/mol, less
accurate when $\Delta F$ is between 2 and 4 kJ/mol, and least accurate for
$\Delta F$ greater than 4 kJ/mol (Figure 7b).

Figure 7: Comparison of free energy estimates for a separation of 1 Å. Mapped
free energy differences were computed between $\lambda_{A}$ and
$\lambda_{B}=\lambda_{A}+1$ Å for $\lambda_{A}\in\\{14,15,...,20\\}$ Å. (a)
Free energy surfaces were computed as a cumulative sum. (b) Individual mapped
$\Delta\hat{F}$ were compared to reference values. In both panels, estimates
were computed with maps from early stopping (magenta squares) and
$\mathcal{L}\sim 0$ (green circles).

Among stopping criterion, TFEP is more reliable when using a learned mapping
with early stopping. When the loss value is near zero, TFEP performance is
sometimes improved over early stopping, but in several cases with large free
energy differences, the estimate significantly deviates from the reference
value. When $-\mathbb{E}_{B}[\Phi_{R}]=\mathbb{E}_{A}[\Phi_{F}]=\Delta F$, we
would expect that the free energies are easy to predict. However, we have
already shown (Figure 6b), that the learned mappings do not provide perfect
mappings and that there is overfitting to the training set.

Differences between stopping criterion are more evident for free energy
surface reconstruction (Figure 7a). With early stopping, the estimated free
energy surface reproduces the local minimum near 16 Å and the barrier between
that state and at 18 Å. Finer details in the free energy surface are not
evident because the mapped estimate was performed at intervals of 1 Å opposed
to 0.2 Å. In contrast to the accurate surface reconstruction using the mapping
from early stopping, when the training loss is near zero, the estimated free
energy surfaces has few features.

These results demonstrate the ability to achieve configuration space overlap
and accurate free energy estimation through learned mapping, opposed to
through sampling of additional intermediate states. Comparable free energy
differences are obtained using a fifth of the number of samples between
adjacent $\lambda$ values. While the mapped estimation procedure incurs
additional costs in training and is more complex, it may lead to an overall
reduction in computational cost.

Using early stopping, reasonable mapped free energy estimates were also
obtained for select pairs of states for which $\Delta\lambda=2$. For
$\lambda_{A}=18$ Å and $\lambda_{B}=20$ Å, the mapped estimate based on early
stopping (3.1 kJ/mol) is better than the $\mathcal{L}\sim 0$ estimate (7.4
kJ/mol) at reproducing the reference free energy (3.6 kJ/mol). For
$\lambda_{A}=19$ Å and $\lambda_{B}=21$ Å, the early stopping estimate (5.5
kJ/mol) is also closer to the reference free energy (4.8 kJ/mol) than the
$\mathcal{L}\sim 0$ estimate (6.6 kJ/mol). These pairs of states are not
separated by potential energy barriers (Figure 3).

###  IV.6 There is poor mapping for more distant states

In contrast to the accurate performance for smaller $\Delta\lambda$, TFEP does
not provide accurate free energy differences for pairs of states in which the
spring center is separated by 4 Å (Table 1). Not only are TFEP results
significantly different from the reference calculation, the two independent
training replicates diverge from one another. This suggests that the
optimization leads to maps that are quite different from one another. Thus,
these maps are at local opposed to global minima of the loss function. To
pinpoint the causes and effects of poor mapping we further analyzed the work,
potential energy, and structural distributions of the mapping between
$\lambda_{A}=14$ Å and $\lambda_{B}=18$ Å .

Table 1:  Comparison of TFEP and MBAR free energy estimates for 4 Å separations. TFEP estimates were based on learned mappings with two independent trainings. $(\lambda_{A},\lambda_{B})$ | $\Delta F_{\mathrm{ref}}$ |  $\Delta F_{\mathrm{TFEP}}$ (Training 1) |  $\Delta F_{\mathrm{TFEP}}$ (Training 2)
---|---|---|---
(14, 18) | 8.461 | -20.419 | -35.181
(15, 19) | 8.691 | -25.454 | -11.418
(16, 20) | 9.269 | 11.294 | 7.835
(17, 21) | 8.154 | 4.757 | 17.162

For this pair of states, TFEP fails to reproduce reference free energies
because the work distributions are still not overlapping (Figure 8). Without
mapping, there is a gap of over 2000 $k_{B}T$ between the peaks of the forward
and reverse work distributions. Mapping significantly reduces this gap, but
there is not significant probability density at the free energy difference,
making it difficult to estimate accurately.

Figure 8:  Normalized histograms of work (a) without and (b) with mapping
between $\lambda_{A}=14$ Å and $\lambda_{B}=18$ Å. $\beta\Phi_{F}$ are colored
green with forward diagonal hatches and $-\beta\Phi_{R}$ are colored red with
backwards diagonal hatches.

The failure to achieve work distribution overlap occurs because mapping does
not reproduce the targeted potential energy distribution (Figure 9). Overall,
the mapped structures have a higher potential energy than the target
distribution. For the forward mapping, the potential energy distribution is
shifted higher and broadened. There is significant overlap with the target
potential energy distribution. For the reverse mapping, there is a more
pronounced shifting and broadening and the high-energy tail of the
distribution is heavier. These suggest that the forward mapping is more
successful than the reverse mapping.

Figure 9:  Normalized histograms of potential energies $U$ for (a)
$\lambda_{A}=14$ and (b) $\lambda_{B}=18$ Å, based on samples from the
distribution (line) or mapped samples (dashed line).

From a structural perspective, the failure of reverse mapping is related to
the formation of hydrogen bonds. While state $A$ is characterized by a stable
$\alpha$-helix, state $B$ is comparatively unfolded. Near the ends of the
helix, structures mapped from state $B$ to state $A$ are similar to structures
from state $A$ (Figure 10). However, the middle of the helix is compressed and
lacks hydrogen bonds. The mapped states are in a metastable state opposed to
the most stable structure of the peptide.

Figure 10:  Comparison of the targeted structures (red color) with the mapped
structures (green color) in the reverse process from $\lambda_{B}=18$ Å to
$\lambda_{A}=14$ Å.

Given the shortcomings of our more distant mappings, the development of
improved training algorithms to map between distant thermodynamic states
across rugged potential energy surfaces is an intriguing area for future
research. One possible approach to better training of mappings may be to
borrow ideas from sampling. Sbailò  _et al._ [29] demonstrated that learned
mappings may be used to propose Monte Carlo moves between metastable states of
a dimer. Mappings capable of generating reasonable structures in the target
distribution will probably work well for TFEP.

As mentioned in the introduction, an alternative to mapping between distant
molecular distributions is to map each molecular distribution to a comparable
tractable distribution. Using this strategy, Ding and Zhang [23] were able to
reproduce free energy differences between deca-alanine at very different
temperatures, 300 K and 500 K. Subsequently, Ding and Zhang [24] applied the
same strategy to obtaining binding free energies for a host-guest complex.
Reference calculations were performed with an attach-pull-release procedure
with multiple intermediate states, suggesting that the end states occupy quite
different configuration space. However, it is unclear whether their methods
will work well on larger and more flexible protein-ligand systems. Moreover,
as mentioned in the introduction, mapping between molecular distributions may
be the only way to map subsets of coordinates.

##  V Conclusions

We have successfully used learned mappings in TFEP to reproduce reference free
energy differences between conformations of a flexible bonded molecule with a
rugged free energy surface, deca-alanine. The learned mappings were most
successful when the unmapped states were similar. On the other hand, our
procedure was unable to learn a mapping from an unfolded to alpha helical
structure of the model peptide.

##  VI Data Availability Statement

Code to perform mapped free energy calculations and sample data are available
at https://github.com/swillow/jax_amber.

##  VII Acknowledgement

We thank Yindong Chen and Joseph DePaolo-Boisvert for helpful discussions and
comments on the manuscript. This research was financially supported in part by
National Institutes of Health grant R01GM127712 to DDLM and National Science
Foundation grants DMS 1916467 and DMS 2153029 to LK. The content is solely the
responsibility of the authors and does not necessarily represent the official
views of the National Institutes of Health or the National Science Foundation.

## References

  * Tolborg  _et al._ [2022] K. Tolborg, J. Klarbring, A. M. Ganose, and A. Walsh, Free energy predictions for crystal stability and synthesisability, Digital Discovery 1, 586 (2022).
  * Armacost  _et al._ [2020] K. A. Armacost, S. Riniker, and Z. Cournia, Novel Directions in Free Energy Methods and Applications, Journal of Chemical Information and Modeling 60, 1 (2020).
  * Pohorille  _et al._ [2010] A. Pohorille, C. Jarzynski, and C. Chipot, Good Practices in Free-Energy Calculations, The Journal of Physical Chemistry B 114, 10235 (2010).
  * Pearlman and Kollman [1989] D. A. Pearlman and P. A. Kollman, The lag between the Hamiltonian and the system configuration in free energy perturbation calculations, The Journal of Chemical Physics 91, 7831 (1989).
  * Wu and Kofke [2005a] D. Wu and D. A. Kofke, Phase-space overlap measures. II. Design and implementation of staging methods for free-energy calculations, The Journal of Chemical Physics 123, 084109 (2005a).
  * Wu and Kofke [2005b] D. Wu and D. A. Kofke, Phase-space overlap measures. I. Fail-safe bias detection in free energies calculated by molecular simulation, The Journal of Chemical Physics 123, 054103 (2005b).
  * Lee  _et al._ [2020] T.-S. Lee, B. K. Allen, T. J. Giese, Z. Guo, P. Li, C. Lin, T. D. McGee, D. A. Pearlman, B. K. Radak, Y. Tao, H.-C. Tsai, H. Xu, W. Sherman, and D. M. York, Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery, Journal of Chemical Information and Modeling 60, 5595 (2020).
  * Jarzynski [2002] C. Jarzynski, Targeted free energy perturbation, Physical Review E 65, 046122 (2002).
  * Zwanzig [1954] R. Zwanzig, High-temperature equation of state by a perturbation method. I. Nonpolar gases, Journal of Chemical Physics 22, 1420 (1954).
  * Voter [1985] A. F. Voter, A Monte Carlo method for determining free-energy differences and transition state theory rate constants, The Journal of Chemical Physics 82, 1890 (1985).
  * Severance  _et al._ [1995] D. L. Severance, J. W. Essex, and W. L. Jorgensen, Generalized alteration of structure and parameters: A new method for free-energy perturbations in systems containing flexible degrees of freedom, Journal of Computational Chemistry 16, 311 (1995), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.540160306 .
  * Jarzynski [1997a] C. Jarzynski, Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach, Physical Review E 56, 5018 (1997a).
  * Jarzynski [1997b] C. Jarzynski, Nonequilibrium equality for free energy differences, Physical Review Letters 78, 2690 (1997b).
  * Oberhofer  _et al._ [2005] H. Oberhofer, C. Dellago, and P. L. Geissler, Biased Sampling of Nonequilibrium Trajectories: Can Fast Switching Simulations Outperform Conventional Free Energy Calculation Methods?, The Journal of Physical Chemistry B 109, 6902 (2005).
  * Miller and Reinhardt [2000] M. A. Miller and W. P. Reinhardt, Efficient free energy calculations by variationally optimized metric scaling: Concepts and applications to the volume dependence of cluster free energies and to solid–solid phase transitions, The Journal of Chemical Physics 113, 7035 (2000).
  * Hahn and Then [2009] A. M. Hahn and H. Then, Using bijective maps to improve free-energy estimates, Physical Review E 79, 011113 (2009).
  * Meng and Schilling [2002] X.-L. Meng and S. Schilling, Warp Bridge Sampling, Journal of Computational and Graphical Statistics 11, 552 (2002).
  * Paliwal and Shirts [2013] H. Paliwal and M. R. Shirts, Multistate reweighting and configuration mapping together accelerate the efficiency of thermodynamic calculations as a function of molecular geometry by orders of magnitude, The Journal of Chemical Physics 138, 154108 (2013).
  * Tan  _et al._ [2010] T. B. Tan, A. J. Schultz, and D. A. Kofke, Efficient calculation of temperature dependence of solid-phase free energies by overlap sampling coupled with harmonically targeted perturbation, The Journal of Chemical Physics 133, 134104 (2010).
  * Moustafa  _et al._ [2015] S. G. Moustafa, A. J. Schultz, and D. A. Kofke, Very fast averaging of thermal properties of crystals by molecular simulation, Physical Review E 92, 043303 (2015).
  * Wirnsberger  _et al._ [2020] P. Wirnsberger, A. J. Ballard, G. Papamakarios, S. Abercrombie, S. Racanière, A. Pritzel, D. Jimenez Rezende, and C. Blundell, Targeted free energy estimation via learned mappings, Journal of Chemical Physics 153, 144112 (2020).
  * Jumper  _et al._ [2021] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, A. Bridgland, C. Meyer, S. A. A. Kohl, A. J. Ballard, A. Cowie, B. Romera-Paredes, S. Nikolov, R. Jain, J. Adler, T. Back, S. Petersen, D. Reiman, E. Clancy, M. Zielinski, M. Steinegger, M. Pacholska, T. Berghammer, S. Bodenstein, D. Silver, O. Vinyals, A. W. Senior, K. Kavukcuoglu, P. Kohli, and D. Hassabis, Highly accurate protein structure prediction with AlphaFold, Nature 596, 583 (2021).
  * Ding and Zhang [2020] X. Ding and B. Zhang, Computing Absolute Free Energy with Deep Generative Models, Journal of Physical Chemistry B 124, 10166 (2020).
  * Ding and Zhang [2021] X. Ding and B. Zhang, DeepBAR: A Fast and Exact Method for Binding Free Energy Computation, Journal of Physical Chemistry Letters 12, 2509 (2021).
  * Rizzi  _et al._ [2021] A. Rizzi, P. Carloni, and M. Parrinello, Targeted free energy perturbation revisited: Accurate free energies from mapped reference potentials, Journal of Physical Chemistry Letters 12, 9449 (2021), arxiv:2106.09295 .
  * Wirnsberger  _et al._ [2022] P. Wirnsberger, G. Papamakarios, B. Ibarz, S. Racanière, A. J. Ballard, A. Pritzel, and C. Blundell, Normalizing flows for atomic solids, Machine Learning: Science and Technology 3, 025009 (2022).
  * Noé  _et al._ [2019] F. Noé, S. Olsson, J. Köhler, and H. Wu, Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning, Science 365, eaaw1147 (2019).
  * Nilmeier  _et al._ [2011] J. P. Nilmeier, G. E. Crooks, D. D. L. Minh, and J. D. Chodera, PNAS Plus: Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation, Proceedings of the National Academy of Sciences of the United States of America 108, E1009 (2011), arxiv:1105.2278 .
  * Sbailò  _et al._ [2021] L. Sbailò, M. Dibak, and F. Noé, Neural mode jump Monte Carlo, Journal of Chemical Physics 154, 074101 (2021).
  * Invernizzi  _et al._ [2022] M. Invernizzi, A. Krämer, C. Clementi, and F. Noé, Skipping the Replica Exchange Ladder with Normalizing Flows, The Journal of Physical Chemistry Letters 13, 11643 (2022).
  * Mahmoud  _et al._ [2022] A. H. Mahmoud, M. Masters, S. J. Lee, and M. A. Lill, Accurate Sampling of Macromolecular Conformations Using Adaptive Deep Learning and Coarse-Grained Representation, Journal of Chemical Information and Modeling 62, 1602 (2022).
  * van der Vaart  _et al._ [2019] A. van der Vaart, P. B. Orndorff, and S. T. Le Phan, Calculation of Conformational Free Energies with the Focused Confinement Method, Journal of Chemical Theory and Computation 15, 6760 (2019).
  * Orndorff  _et al._ [2020] P. B. Orndorff, S. T. Le Phan, K. H. Li, and A. van der Vaart, Conformational Free-Energy Differences of Large Solvated Systems with the Focused Confinement Method, Journal of Chemical Theory and Computation 16, 5163 (2020).
  * Zuckerman and Woolf [2004] D. M. Zuckerman and T. B. Woolf, Systematic finite-sampling inaccuracy in free energy differences and other nonlinear quantities, Journal of Statistical Physics 114, 1303 (2004).
  * Bennett [1976] C. H. Bennett, Efficient Estimation of Free-Energy Differences from Monte Carlo Data, Journal of Computational Physics 22, 245 (1976).
  * Shirts and Chodera [2008] M. R. Shirts and J. D. Chodera, Statistically optimal analysis of samples from multiple equilibrium states, Journal of Chemical Physics 129, 124105 (2008).
  * Kullback [1959] S. Kullback,  _Information Theory and Statistics_ (John Wiley, 1959).
  * Nguyen and Minh [2016] T. H. Nguyen and D. D. L. Minh, Intermediate Thermodynamic States Contribute Equally to Free Energy Convergence: A Demonstration with Replica Exchange, Journal of Chemical Theory and Computation 12, 2154 (2016).
  * Maier  _et al._ [2015] J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser, and C. Simmerling, ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB, Journal of Chemical Theory and Computation 11, 3696 (2015).
  * Eastman  _et al._ [2017] P. Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A. Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, R. P. Wiewiora, B. R. Brooks, and V. S. Pande, OpenMM 7: Rapid development of high performance algorithms for molecular dynamics, PLOS Computational Biology 13, e1005659 (2017).
  * Dinh  _et al._ [2017] L. Dinh, J. Sohl-Dickstein, and S. Bengio, Density estimation using Real NVP, arXiv:1605.08803 [cs, stat] (2017), arxiv:1605.08803 [cs, stat] .
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: Composable transformations of Python+NumPy programs (2018).
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization (2017), arxiv:1412.6980 [cs] .

## Supplementary Material

Figure S1: Effect of stopping criterion on estimated probability densities of
$\beta\Phi$ for test data between $\lambda_{A}=20$ Å and $\lambda_{B}=21$ Å.
Work distributions for mapping functions from (a) early stopping, (b) when the
loss is closest to zero, and (c) after 150000 training steps. The reference
free energy difference $\Delta F$ obtained via the MBAR estimator is
represented by the vertical dashed line. A similar figure for the training
data is available as Fig 6. Figure S2:  Estimated probability densities of (a)
and (c) potential energies and (b) and (d) terminal nitrogen positions for
$\lambda_{A}=20$ Å and $\lambda_{B}=21$ Å. Normalized histograms are based on
simulated (blue line) and mapped samples based on early stopping (green dashed
line) and complete training (red dotted line).


$\tilde{\rho}\to\infty$. Doing this we get

$c_{1}=5,\qquad c_{2}=-{5\pi\over 4},\qquad
c_{3}=-\frac{25+13\pi}{2\sqrt{2}},\qquad c_{4}=\frac{13}{4\sqrt{2}},\qquad
c_{5}=c_{6}=0\,\,.$ (B.16)

Let us now write the final form of our ${\cal O}(\epsilon)$ solution. The
dilaton $\phi(\tilde{\rho})$, the warp factor $h(\tilde{\rho})$, and the
functions $\alpha(\tilde{\rho})$ and $G(\tilde{\rho})$ are written in (3) in
terms of the function $\Omega(\tilde{\rho})$ defined in (3.6). The internal
deformation functions $F(\tilde{\rho})$ and $S(\tilde{\rho})$ are written in
(3.8) in terms of the first-order functions $F_{1}(\tilde{\rho})$ and
$S_{1}(\tilde{\rho})$. In order to present compact expressions for these
functions, let us define the integrals ${\cal I}_{P}(x)$ and ${\cal I}_{Q(x)}$
as

${\cal
I}_{P}^{(1)}(x)\,=\,\int_{1}^{x}z^{2}P(z)\,{\rm{d}}z\,\,,\qquad\qquad{\cal
I}_{Q}^{(1)}(x)\,=\,\int_{x}^{\infty}z^{2}Q(z)\,{\rm{d}}z\,\,.$ (B.17)

This procedure finally yields

$\displaystyle F_{1}(\tilde{\rho})$ $\displaystyle=$
$\displaystyle{3\sqrt{2}\,\pi\over 8}\,\Big{(}P(\tilde{\rho})\,{\cal
I}_{Q}^{(1)}(\tilde{\rho})+Q(\tilde{\rho})\,{\cal
I}_{P}^{(1)}(\tilde{\rho})\Big{)}+{1\over 2}\,\tilde{\rho}^{3}-{1\over
8}\,+\,{1\over 10}\,(\tilde{\rho}^{4}+2)\,\Omega(\tilde{\rho})$ $\displaystyle
S_{1}(\tilde{\rho})$ $\displaystyle=$ $\displaystyle-{3\sqrt{2}\,\pi\over
32}\,\Big{(}P(\tilde{\rho})\,{\cal
I}_{Q}^{(1)}(\tilde{\rho})+Q(\tilde{\rho})\,{\cal
I}_{P}^{(1)}(\tilde{\rho})\Big{)}+{1\over 2}\,\tilde{\rho}^{3}-{1\over
8}\,+\,{1\over 10}\,(\tilde{\rho}^{4}+2)\,\Omega(\tilde{\rho})\ .$ (B.18)

##  Appendix C Solution with non-zero chemical potential

Let us now add non-zero chemical potential to the black hole solution of
appendix B. We will first expand the dilaton in powers of $\epsilon$ and
$\delta$ as:

$\phi=\epsilon\,\phi_{1}(\tilde{\rho})\,+\,\epsilon\,\delta^{2}\,\phi_{2}(\tilde{\rho})\,\,,$
(C.1)

where $\phi_{1}(\tilde{\rho})$ is the solution found in appendix B,in (B.7).
It turns out that the equation for $\phi_{2}$ is completely decoupled and
given by

$\phi_{2}^{\prime\prime}+\frac{\left(5\tilde{\rho}^{4}-1\right)}{\tilde{\rho}\left(\tilde{\rho}^{4}-1\right)}\phi_{2}^{\prime}-\frac{Q_{\textrm{\scriptsize
c}}^{2}}{120\,\rho_{\textrm{\scriptsize
h}}^{6}\,\tilde{\rho}^{3}\left(\tilde{\rho}^{4}-1\right)}=0\,\,.$ (C.2)

The solution for this equation is:

$\phi_{2}=d_{2}+d_{1}\log\left(1-\frac{1}{\tilde{\rho}^{4}}\right)-\frac{Q_{\textrm{\scriptsize
c}}^{2}}{480\rho_{\textrm{\scriptsize
h}}^{6}\tilde{\rho}}\left(\tilde{\rho}\left(\log\left(\frac{\tilde{\rho}-1}{\tilde{\rho}+1}\right)+2\arctan(\tilde{\rho})\right)+4\right)\
,$ (C.3)

where $d_{1}$ and $d_{2}$ are integration constants. Similarly, we expand
$\alpha(\tilde{\rho})$ and $b(\tilde{\rho})$ as:

$\alpha(\tilde{\rho})\,=\,\alpha_{0}\,+\,\epsilon\,\alpha_{1}(\tilde{\rho})\,+\,\epsilon\,\delta^{2}\,\alpha_{2}(\tilde{\rho})\,\,,\qquad\qquad
b(\tilde{\rho})\,=\,b_{0}(\tilde{\rho})\,+\,\epsilon\,\delta^{2}\,b_{2}(\tilde{\rho})\,\,,$
(C.4)

where $\alpha_{0}$, $\alpha_{1}(\tilde{\rho})$, and $b_{0}(\tilde{\rho})$ are
the zeroth and first order terms in (3) per subscript labeling. Notice that
$b(\tilde{\rho})$ does not contain ${\mathcal{O}}(\epsilon))$ terms, in
agreement with our result (3) for the black hole with vanishing chemical
potential. The equations for $\alpha_{2}(\tilde{\rho})$ and
$b_{2}(\tilde{\rho})$ are:

$\displaystyle\alpha_{2}^{\prime\prime}+\frac{5\tilde{\rho}^{4}-1}{\tilde{\rho}\left(\tilde{\rho}^{4}-1\right)}\alpha_{2}^{\prime}-\frac{Q_{\textrm{\scriptsize
c}}^{2}}{120\rho_{\textrm{\scriptsize
h}}^{6}\,\tilde{\rho}^{3}\left(\tilde{\rho}^{4}-1\right)}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle
b_{2}^{\prime\prime}-\frac{5\tilde{\rho}^{4}-9}{\tilde{\rho}\left(\tilde{\rho}^{4}-1\right)}b_{2}^{\prime}+\frac{16}{\tilde{\rho}^{2}\left(\tilde{\rho}^{4}-1\right)^{2}}b_{2}-\frac{Q_{\textrm{\scriptsize
c}}^{2}}{30\,\rho_{\textrm{\scriptsize h}}^{6}\,\tilde{\rho}^{7}}$
$\displaystyle=$ $\displaystyle 0\ ,$ (C.5)

and, remarkably, can be solved in terms of $\phi_{2}(\tilde{\rho})$ as

$\alpha_{2}(\tilde{\rho})\,=\,\phi_{2}(\tilde{\rho})\,\,,\qquad\qquad
b_{2}(\tilde{\rho})\,=\,4\Big{(}1-{1\over\tilde{\rho}^{4}}\Big{)}\,\phi_{2}(\tilde{\rho})\,\,.$
(C.6)

Similarly, we can expand the warp factor $h(\tilde{\rho})$ as:

$h(\tilde{\rho})\,=\,h_{0}(\tilde{\rho})\,+\,\epsilon\,h_{1}(\tilde{\rho})\,+\,\epsilon\,\delta^{2}\,h_{2}(\tilde{\rho})\,\,,$
(C.7)

where $h_{0}(\tilde{\rho})$ and $h_{1}(\tilde{\rho})$ can be extracted from
(3). The second order contribution $h_{2}(\tilde{\rho})$ satisfies the
equation

$\displaystyle
h_{2}^{\prime\prime}+\frac{13\tilde{\rho}^{4}-9}{\tilde{\rho}\left(\tilde{\rho}^{4}-1\right)}h_{2}^{\prime}-\frac{16}{\tilde{\rho}^{2}\left(\tilde{\rho}^{4}-1\right)}h_{2}$
$\displaystyle\quad\qquad+\frac{4Q_{\textrm{\scriptsize
c}}\,\tilde{\rho}^{2}}{\rho_{\textrm{\scriptsize
h}}^{4}\left(\tilde{\rho}^{4}-1\right)^{2}}b_{2}+\frac{8Q_{\textrm{\scriptsize
c}}}{\rho_{\textrm{\scriptsize
h}}^{4}\,\tilde{\rho}^{2}\left(\tilde{\rho}^{4}-1\right)}\alpha_{2}-\frac{Q_{\textrm{\scriptsize
c}}^{3}}{160\rho_{\textrm{\scriptsize
h}}^{10}\,\tilde{\rho}^{7}\left(\tilde{\rho}^{4}-1\right)}=0\,\,,$ (C.8)

which can also be solved in terms of the dilaton function
$\phi_{2}(\tilde{\rho})$ as:

$h_{2}(\tilde{\rho})=\frac{3Q_{\textrm{\scriptsize
c}}}{4\rho_{\textrm{\scriptsize
h}}^{4}\tilde{\rho}^{4}}\,\phi_{2}(\tilde{\rho})\,\,.$ (C.9)

Let us next represent $F(\tilde{\rho})$ and $S(\tilde{\rho})$ as in (4.7). The
differential equations for $F_{2}(\tilde{\rho})$ and $S_{2}(\tilde{\rho})$ are
coupled and contain $\alpha_{2}(\tilde{\rho})$ and $b_{2}(\tilde{\rho})$ as
sources. These equations can be decoupled by performing the same linear
combinations as in (B.8), i.e., by introducing new functions
$g_{1}(\tilde{\rho})$ and $g_{2}(\tilde{\rho})$ as

$\displaystyle
F_{2}(\tilde{\rho})={1\over\sqrt{2}}g_{1}(\tilde{\rho})-{4\over\sqrt{17}}g_{2}(\tilde{\rho})\qquad\qquad
S_{2}(\tilde{\rho})={1\over\sqrt{2}}g_{1}(\tilde{\rho})+{1\over\sqrt{17}}g_{2}(\tilde{\rho})\,\,.$
(C.10)

The equations satisfied by $g_{1}$ and $g_{2}$ are:

$\displaystyle
g_{1}^{\prime\prime}+\frac{1-5\tilde{\rho}^{4}}{\tilde{\rho}-\tilde{\rho}^{5}}\,g_{1}^{\prime}-\frac{32\tilde{\rho}^{2}}{\tilde{\rho}^{4}-1}\,g_{1}-{24\,\sqrt{2}\,\tilde{\rho}^{2}\over\tilde{\rho}^{4}-1}\Big{(}d_{1}\log\left(1-\frac{1}{\tilde{\rho}^{4}}\right)+d_{2}\Big{)}$
$\displaystyle\qquad+{Q_{\textrm{\scriptsize c}}^{2}\,\tilde{\rho}^{3}\over
10\sqrt{2}\,\rho_{\textrm{\scriptsize
h}}^{6}\,(\tilde{\rho}^{4}-1)}\Big{(}\log\left(\frac{\tilde{\rho}-1}{\tilde{\rho}+1}\right)+2\arctan(\tilde{\rho})\Big{)}+{Q_{\textrm{\scriptsize
c}}^{2}\,(80\tilde{\rho}^{4}+3)\over 200\sqrt{2}\,\rho_{\textrm{\scriptsize
h}}^{6}\,\tilde{\rho}^{3}(\tilde{\rho}^{4}-1)}=0$ $\displaystyle
g_{2}^{\prime\prime}+\frac{1-5\tilde{\rho}^{4}}{\tilde{\rho}-\tilde{\rho}^{5}}\,g_{2}^{\prime}-\frac{12\tilde{\rho}^{2}}{\tilde{\rho}^{4}-1}\,g_{2}+\frac{\sqrt{17}\,Q_{\textrm{\scriptsize
c}}^{2}}{1200\,\rho_{\textrm{\scriptsize
h}}^{6}\,\tilde{\rho}^{3}\left(\tilde{\rho}^{4}-1\right)}=0\,\,.$ (C.11)

The equation for $g_{1}$ can be explicitly integrated as:

$\displaystyle g_{1}(\tilde{\rho})$ $\displaystyle=$
$\displaystyle-\frac{1}{2400\sqrt{2}\rho_{\textrm{\scriptsize
h}}^{6}\tilde{\rho}}\Bigg{(}7200\rho_{\textrm{\scriptsize
h}}^{6}\tilde{\rho}^{5}\left(d_{1}\log\left(1-\frac{1}{\tilde{\rho}^{4}}\right)+2d_{1}+d_{2}\right)+4Q_{\textrm{\scriptsize
c}}^{2}\left(5\tilde{\rho}^{4}-9\right)$ (C.12)
$\displaystyle+5Q_{\textrm{\scriptsize
c}}^{2}\tilde{\rho}\left(\tilde{\rho}^{4}-2\right)\left(\log\left(\frac{\tilde{\rho}-1}{\tilde{\rho}+1}\Bigg{)}+2\arctan(\tilde{\rho})\right)\right)+d_{3}\left(1-2\tilde{\rho}^{4}\right)$
$\displaystyle+\frac{1}{2}d_{4}\left(2-\left(1-2\tilde{\rho}^{4}\right)\log\left(1-\frac{1}{\tilde{\rho}^{4}}\right)\right)\,\,,$

where $d_{3}$ and $d_{4}$ are new integration constants. To integrate the
equation for $g_{2}$ in (C.11) we again make use of the Wronskian method. The
homogeneous part of the equation is the same as the one in the equation for
$f_{2}$ in (B.9) and, therefore, is solved by the functions $P$ and $Q$
defined in (B.11). To write a special solution of the second-order
inhomogeneous equation in (C.11), we define new functions ${\cal
I}_{P}^{(2)}(x)$ and ${\cal I}_{Q}^{(2)}(x)$ as

${\cal I}_{P}^{(2)}(x)\,=\,\int_{1}^{x}{P(z)\over
z^{2}}\,{\rm{d}}z\,\,,\qquad\qquad{\cal
I}_{Q}^{(2)}(x)\,=\,\int_{x}^{\infty}{Q(z)\over z^{2}}\,{\rm{d}}z\,\,.$ (C.13)

Then, the general solution for $g_{2}(\tilde{\rho})$ can be written as:

$g_{2}(\tilde{\rho})=d_{5}\,P(\tilde{\rho})+d_{6}\,Q(\tilde{\rho})+{\sqrt{17}\,\pi\,Q_{\textrm{\scriptsize
c}}^{2}\over 9600\,\sqrt{2}\,\rho_{\textrm{\scriptsize
h}}^{6}}\Big{(}P(\tilde{\rho})\,{\cal
I}_{Q}^{(2)}(\tilde{\rho})\,+\,Q(\tilde{\rho})\,{\cal
I}_{P}^{(2)}(\tilde{\rho})\Big{)}\,\,,$ (C.14)

where $d_{5}$ and $d_{6}$ are two new integration constants. By imposing
regularity of $\phi_{2}$, $g_{1}$, and $g_{2}$, both at the horizon
$\tilde{\rho}=1$ and at the UV $\tilde{\rho}\to\infty$, we can determine the
constants of integration, which can be shown to be given by:

$\displaystyle d_{1}\,=\,{d_{2}\over\pi}\,=\,{Q_{\textrm{\scriptsize
c}}^{2}\over 480\,\rho_{\textrm{\scriptsize h}}^{6}}\ \ ,\qquad\qquad
d_{3}\,=\,-{3+2\pi\over 480\,\sqrt{2}}\,{Q_{\textrm{\scriptsize
c}}^{2}\over\rho_{\textrm{\scriptsize h}}^{5}}$ $\displaystyle
d_{4}\,=\,{Q_{\textrm{\scriptsize c}}^{2}\over
120\,\sqrt{2}\,\rho_{\textrm{\scriptsize h}}^{5}}\ \ ,\qquad\qquad
d_{5}\,=\,d_{6}\,=\,0\,\,.$ (C.15)

Using these values for the $d_{i}$’s we get the expressions of
$\phi(\tilde{\rho})$, $\alpha(\tilde{\rho})$, $b(\tilde{\rho})$, and
$h(\tilde{\rho})$ displayed in eqs. (4.5) and (4.6) of the main text. The
expression for $G(\tilde{\rho})$ in 4.6) can be obtained from the other
functions by using (2.8). The functions $F(\tilde{\rho})$ and
$G(\tilde{\rho})$ can be written as in (4.7), with $F_{2}(\tilde{\rho})$ and
$S_{2}(\tilde{\rho})$ given by

$\displaystyle F_{2}(\tilde{\rho})$ $\displaystyle=$ $\displaystyle-{\pi\over
4\sqrt{2}}\,\Big{(}P(\tilde{\rho})\,{\cal
I}_{Q}^{(2)}(\tilde{\rho})+Q(\tilde{\rho})\,{\cal
I}_{P}^{(2)}(\tilde{\rho})\Big{)}\,-\,{5\over 2}\,\tilde{\rho}^{3}+{5\over
8}\,+\,{9\over 2\,\tilde{\rho}}+{1\over
2}(2-\tilde{\rho}^{4})\,\Omega(\tilde{\rho})$ $\displaystyle
S_{2}(\tilde{\rho})$ $\displaystyle=$ $\displaystyle{\pi\over
16\sqrt{2}}\,\Big{(}P(\tilde{\rho})\,{\cal
I}_{Q}^{(2)}(\tilde{\rho})+Q(\tilde{\rho})\,{\cal
I}_{P}^{(2)}(\tilde{\rho})\Big{)}\,-\,{5\over 2}\,\tilde{\rho}^{3}+{5\over
8}\,+\,{9\over 2\,\tilde{\rho}}+{1\over
2}(2-\tilde{\rho}^{4})\,\Omega(\tilde{\rho})\ .$ (C.16)

To complete the solution it remains to find $j(\tilde{\rho})$. Its equation at
leading order can be easily found from (A.8) and is given by

$j^{\prime\prime}+\frac{\tilde{\rho}^{4}+3}{\tilde{\rho}\left(\tilde{\rho}^{4}-1\right)}j^{\prime}-\frac{9\,\tilde{\rho}^{2}}{\tilde{\rho}^{4}-1}j-\frac{Q_{\textrm{\scriptsize
c}}\,\tilde{\rho}}{2\rho_{\textrm{\scriptsize
h}}\left(\tilde{\rho}^{4}-1\right)}=0\ .$ (C.17)

The homogeneous version of this equation has the following two independent
solutions:

$\displaystyle
J_{1}(\tilde{\rho})=F\left(-\frac{3}{4},\frac{3}{4};1;1-\tilde{\rho}^{4}\right)\,\,,\qquad
J_{2}(\tilde{\rho})=\left(\tilde{\rho}^{4}-1\right)^{-3/4}\,\,F\left(\frac{3}{4},\frac{3}{4};\frac{5}{2};\frac{1}{1-\tilde{\rho}^{4}}\right)\,\,.\qquad\qquad$
(C.18)

The general solution for $j(\tilde{\rho})$ can be found by adding a linear
combination of $J_{1}(\tilde{\rho})$ and $J_{2}(\tilde{\rho})$ and a special
solution of the complete equation obtained again by the Wronskian method.
Notice that $J_{1}(\tilde{\rho})$ is regular at $\tilde{\rho}=1$ and blows up
at $\tilde{\rho}\to\infty$, whereas $J_{2}(\tilde{\rho})$ has the opposite
behavior and diverges at the horizon and is regular at the UV. By taking these
facts into account, it is not difficult to find the regular solution for
$j(\tilde{\rho})$. If we define ${\cal I}_{1}(\tilde{\rho})$ and ${\cal
I}_{2}(\tilde{\rho})$ as

${\cal I}_{1}(\tilde{\rho})\,=\,\int_{1}^{\tilde{\rho}}{J_{1}(z)\over
z^{2}}{\rm{d}}z\,\,,\qquad\qquad{\cal
I}_{2}(\tilde{\rho})\,=\,\int_{\tilde{\rho}}^{\infty}{J_{2}(z)\over
z^{2}}{\rm{d}}z$ (C.19)

then we have

$\displaystyle j(\tilde{\rho})\,=\,-\frac{Q_{\textrm{\scriptsize
c}}\big{[}\Gamma\left(\frac{3}{4}\right)\big{]}^{2}}{8\sqrt{\pi}\,\rho_{\textrm{\scriptsize
h}}}\left(J_{1}(\tilde{\rho})\,{\cal
I}_{2}(\tilde{\rho})+J_{2}(\tilde{\rho})\,{\cal I}_{1}(\tilde{\rho})\right)\
.$ (C.20)

On the other hand, we can use these results to find the worldvolume gauge
field $A_{t}^{\prime}$. Let us write:

$\displaystyle{1\over 2\pi\alpha^{\prime}}A_{t}^{\prime}(\tilde{\rho})={1\over
2\pi\alpha^{\prime}}\Big{(}A_{t}^{\prime\,(0)}(\tilde{\rho})+\epsilon\,A_{t}^{\prime\,(1)}(\tilde{\rho})\Big{)}=\tilde{\delta}\,\left(a_{0}(\tilde{\rho})+\epsilon\,a_{1}(\tilde{\rho})\right)\
,$ (C.21)

where

$\displaystyle a_{0}(\tilde{\rho})$ $\displaystyle=$
$\displaystyle-\frac{\rho_{\textrm{\scriptsize
h}}}{\sqrt{6}\pi\alpha^{\prime}}{1\over\tilde{\rho}^{2}}$ $\displaystyle
a_{1}(\tilde{\rho})$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{3}\,\rho_{\textrm{\scriptsize
h}}}{16\,\alpha^{\prime}\tilde{\rho}^{2}}\Big{(}P(\tilde{\rho})\,I_{Q}(\tilde{\rho})+Q(\tilde{\rho})\,I_{P}(\tilde{\rho})-\frac{\sqrt{2}\left(3\left(1-4\tilde{\rho}^{3}\right)-\left(3\tilde{\rho}^{4}+1\right){4\over
5}\Omega(\tilde{\rho})\right)}{3\pi}$ (C.22)
$\displaystyle+\frac{80\sqrt{2}\,\rho_{\textrm{\scriptsize
h}}}{\pi\,Q_{\textrm{\scriptsize c}}}j(\tilde{\rho})\Big{)}\ .$

The UV value of $A_{t}/(2\pi\alpha^{\prime})$ is related to the chemical
potential $\mu$. This value can be obtained by integrating (C.21) with respect
to $\tilde{\rho}$. Doing this integration with the condition that
$A_{t}(\tilde{\rho}=1)=0$, we get

${A_{t,UV}\over 2\pi\alpha^{\prime}}\,=\,{A_{t}(\tilde{\rho}\to\infty)\over
2\pi\alpha^{\prime}}\,=\,{1\over
2\pi\alpha^{\prime}}\int_{1}^{\infty}\,A_{t}^{\prime}(\tilde{\rho})\,{\rm{d}}\tilde{\rho}\,=\,\tilde{\delta}\,\left(a_{0,UV}\,+\,\epsilon\,a_{1,UV}\right)\,\,,$
(C.23)

where $a_{0,\text{UV}}$ and $a_{1,\text{UV}}$ are given by

$\displaystyle a_{0,\text{UV}}$ $\displaystyle=$
$\displaystyle\int_{1}^{\infty}\,a_{0}(\tilde{\rho})\,{\rm{d}}\tilde{\rho}\,=\,-{\rho_{\textrm{\scriptsize
h}}\over\sqrt{6}\,\pi\alpha^{\prime}}$ $\displaystyle a_{1,\text{UV}}$
$\displaystyle=$
$\displaystyle\int_{1}^{\infty}\,a_{1}(\tilde{\rho})\,{\rm{d}}\tilde{\rho}\,\approx\,-1.06\,a_{0,\text{UV}}\
,$ (C.24)

The integration of $a_{0,\text{UV}}$ in (C.24) can be done analytically,
whereas the one needed to get $a_{1,\text{UV}}$ has to be done numerically. By
comparing these results with the chemical potential calculated in (5.40), we
notice that the ${\cal O}(\epsilon\tilde{\delta})$ terms are different (see
(5.42)). This difference can be attributed to the backreaction of the flavors.

##  Appendix D Dimensional reduction

The ten-dimensional supergravity action can be dimensionally reduced both to
four and five dimensions. The details of these constructions were given in
[51] and will be briefly reviewed in this appendix. Since this reduced theory
corresponds to vanishing chemical potential, in this appendix we restrict
ourselves to this $\mu=0$ case. In the reduction to five dimensions we adopt
the following ansatz for the metric

${\rm{d}}s^{2}_{10}=e^{{10\over
3}\gamma}g_{pq}{\rm{d}}z^{p}{\rm{d}}z^{q}+e^{-2(\gamma+\lambda)}{\rm{d}}s^{2}_{\rm{KE}}+e^{2(4\lambda-\gamma)}\left({\rm{d}}\tau+{\cal
A}\right)^{2}\,\,,$ (D.1)

where $g_{pq}$ is a metric for the coordinates
$z^{p}=(t,x^{1},x^{2},x^{3},\rho)$ and $(\gamma,\lambda)$ are scalar fields
depending on the 5d coordinates $z^{p}$. In our ansatz $\gamma$ and $\lambda$
depend only on the radial coordinate $\rho$ and are related to the functions
of the ten-dimensional metric as

$\lambda={1\over 5}\left(\log F-\log S\right)\,\,,\qquad\qquad\gamma={1\over
10}\left(-2\log F-8\log S-{5\over 2}\log h\right)\,\,.$ (D.2)

Moreover, the reduced 5d metric takes the form

${\rm{d}}s^{2}_{5}=-c_{1}^{2}\,{\rm{d}}t^{2}+c_{2}^{2}\left(({\rm{d}}x^{1})^{2}+({\rm{d}}x^{2})^{2}\right)+c_{3}^{2}\,({\rm{d}}x^{3})^{2}+c_{\rho}^{2}\,{\rm{d}}\rho^{2}\,\,,$
(D.3)

where the different coefficients are:

$\displaystyle c_{1}=F^{1/3}S^{4/3}b^{1/2}h^{1/6}\,\,,\qquad\qquad
c_{2}=F^{1/3}S^{4/3}h^{1/6}$ $\displaystyle
c_{3}=\alpha\,F^{1/3}S^{4/3}h^{1/6}\,\,,\qquad\qquad
c_{\rho}=\frac{\alpha\,F^{4/3}h^{2/3}S^{16/3}}{\rho^{5}b^{1/2}}\,\,.$ (D.4)

We will use the reduced theory to compute the VEV of the stress-energy tensor.
This calculation makes use of the holographic renormalization formalism to
regulate the on-shell action. The VEV of the stress-energy tensor of the dual
theory is obtained by taking the functional derivative of the renormalized
theory with respect to the boundary metric. The result is [51]:

$\displaystyle\langle T^{\mu}{}_{\nu}\rangle={V_{5}\over
2\kappa_{10}^{2}}\sqrt{\gamma}\left(-2K^{\mu}{}_{\nu}+\delta^{\mu}{}_{\nu}\left(2K+W_{5d}\right)\right)\Big{|}_{\rho\to\infty}+Q_{\textrm{\scriptsize
f}}{V_{5}\over\kappa_{10}^{2}}\sqrt{\hat{\gamma}}e^{2\gamma+2\lambda+{\phi\over
2}}\delta^{\mu}{}_{\nu}\Big{|}_{\rho\to\infty}\,\,,$ (D.5)

where $\gamma$ is the determinant of the metric induced on constant $\rho$
slices, $K^{\mu}{}_{\nu}$ is the extrinsic curvature of these constant $\rho$
surfaces ($K=K^{\mu}{}_{\mu}$), and $W_{5d}$ is the superpotential, given by

$\displaystyle W_{5d}=-6e^{{8\over 3}\gamma-4\lambda}-4e^{{8\over
3}\gamma+6\lambda}+Q_{\textrm{\scriptsize c}}e^{{20\over 3}\gamma}\,\,.$ (D.6)

In the last term in (D.5), $\hat{\gamma}$ is the determinant of the metric
induced on constant $\rho$ and $x^{3}$ slices. It is not difficult to show
that, at first order in $\epsilon$, that

$\langle
T^{\mu}{}_{\nu}\rangle=\text{diag}(-\mathcal{E},p_{xy},p_{xy},p_{z})=\mathcal{E}\,\text{diag}\Bigg{(}-1\,,\,{1\over
3}\,+\,{10\over 9}\epsilon\,,\,{1\over 3}\,+\,{10\over 9}\epsilon\,,\,{1\over
3}\,-\,{20\over 9}\epsilon\Bigg{)}\,\,,$ (D.7)

where $\mathcal{E}$ is the energy density:

$\displaystyle\mathcal{E}={3\,{\rm Vol}\big{(}{\cal
M}_{5}\big{)}\over(2\pi)^{7}\,\alpha^{\prime 4}\,g_{\textrm{\scriptsize
s}}^{2}}\,\rho_{\textrm{\scriptsize
h}}^{4}\left(1\,-\,\frac{5}{6}\epsilon\right)\ .$ (D.8)

This result matches with the one found by thermodynamic methods in section 5.

###  D.1 Hydrodynamics in the shear channel

The five-dimensional supergravity theory can be further reduced along the
$x^{3}$-direction. The $5d\to 4d$ ansatz for the reduction takes the form[51]:

${\rm{d}}s^{2}_{5}\,=\,e^{-\beta}\,{\rm{d}}s^{2}_{4}\,+\,e^{2\beta}\,({\rm{d}}x^{3})^{2}\
,$ (D.9)

where $\beta$ is a new scalar field of the $4d$ supergravity. For our ansatz
$e^{\beta}\,=\,c_{3}$, where $c_{3}$ is the function displayed in (D.4). This
$4d$ theory can be used to compute the VEV of the stress-energy tensor of the
dual theory along the $(t,x^{1},x^{2})$-directions. The result of this
calculation coincides with the one displayed in (D.7) and will not be repeated
here. Moreover, we can use the $4d$ theory to obtain the transport
coefficients for perturbations propagating along the ($x^{1},x^{2}$)-plane. In
this appendix we will analyze the so-called shear channel. With this purpose
we fluctuate the $4d$ metric about the background solution via the
substitution

$g_{mn}\,\to\,g_{mn}\,+\,h_{mn}\,\,,$ (D.10)

where $g_{mn}$ is the reduced metric of the first-order solution at zero
chemical potential. We will work in the radial gauge for the metric, in which

$h_{m\rho}=0\,\,,\qquad\qquad(m=t,x^{1},x^{2},\rho)\,\,.$ (D.11)

Without loss of generality we consider a perturbation propagating along the
$x^{2}$-direction. In the shear channel only the metric fluctuations
$h_{tx^{1}}$ and $h_{x^{1}\,x^{2}}$ are excited. Let us assume that these
fluctuations (in the Fourier space) have frequency $\omega$ and momentum $q$
and, accordingly, let us parametrize them as

$h_{tx^{1}}=e^{-i\left(\omega
t-qx^{2}\right)}c_{2}^{2}(\rho)\,H_{tx^{1}}(\rho)\,\,,\qquad\qquad
h_{x^{1}\,x^{2}}=e^{-i\left(\omega
t-qx^{2}\right)}c_{2}^{2}(\rho)\,H_{x^{1}\,x^{2}}(\rho)\,\,.$ (D.12)

It turns out that the equations for the fluctuations can be reduced to a
single second-order differential equation for a gauge invariant combination
$X=X_{0}+\epsilon\,X_{1}$, defined as

$X(\rho)=q\,H_{tx^{1}}(\rho)+\omega H_{x^{1}\,x^{2}}(\rho)\,\,.$ (D.13)

In order to write the fluctuation equation in a simple form, let us perform
the following change of the radial variable

$x\,=\,\sqrt{b(\rho})\,=\,\sqrt{1-{\rho_{\textrm{\scriptsize
h}}^{4}\over\rho^{4}}}\,\,.$ (D.14)

In this new variable, the horizon is at $x=0$, while the boundary resides at
$x=1$. Moreover, let us work with the dimensionless frequency and momentum
$\hat{q}$ and $\hat{\omega}$ defined in (6.2). Then, if now the prime denotes
derivative with respect to $x$, the fluctuation equation becomes:

$\displaystyle
X^{\prime\prime}\,-\,\frac{\hat{q}^{2}x^{2}+\hat{\omega}^{2}}{x\left(\hat{q}^{2}x^{2}-\hat{\omega}^{2}\right)}\,X^{\prime}\,+\,\Gamma\,\frac{\hat{q}^{2}x^{2}-\hat{\omega}^{2}}{x^{2}\sqrt{1-x^{2}}\,k(x)}\,X=0\,\,,$
(D.15)

where $\Gamma$ is the constant:

$\displaystyle\Gamma=\frac{\pi^{2}T^{2}Q_{\textrm{\scriptsize
c}}}{\rho_{\textrm{\scriptsize h}}^{2}}=4\,-\,15\,\epsilon\,\,,$ (D.16)

and $k(x)$ is the following function

$\displaystyle
k(x)=4\left(x^{2}-1\right)+5\epsilon\left(x^{2}-1+4\left(1-x^{2}\right)^{1/4}+{4\over
5}\,x^{2}\,\Omega\left({1\over\left(1-x^{2}\right)^{1/4}}\right)\right)\,\,.$
(D.17)

Notice that $\Gamma$ and $k(x)$ contain terms which are zeroth-order in the
flavor parameter, together with the ones that are first-order in $\epsilon$.
We want to find solutions of (D.15) which satisfy infalling boundary
conditions at the horizon $x=0$, as well as Dirichlet boundary conditions at
the UV boundary $x=1$. These solutions only exist if certain dispersion
relation $\hat{\omega}=\hat{\omega}(\hat{q})$ is satisfied. Our objective in
this section is to find such dispersion relation in a power series of
$\hat{q}$. Actually, we will only keep quadratic and quartic terms of
$\hat{q}$, as in (6.4), which contains two transport coefficients
$\hat{D}_{\eta}$ and $\hat{\tau}_{s}$.

In order to impose infalling boundary conditions at the horizon $x=0$, we
adopt an ansatz for $X(x)$

$X(x)\,=\,x^{-i\hat{\omega}}\,S(x)\,\,,$ (D.18)

where $S(x)$ must be regular at $x=0$. Let us expand $S(x)$ in powers of
$\hat{q}$ as:

$S(x)\,=\,S_{0}(x)\,+\,\hat{q}^{2}\,S_{2}(x)\,+\,\ldots\,\,.$ (D.19)

Plugging (D.18) and the expansion (D.19) into (D.15) we get the following two
differential equations for $S_{0}(x)$ and $S_{2}(x)$:

$\displaystyle S_{0}^{\prime\prime}-{1\over x}\,S_{0}^{\prime}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
S_{2}^{\prime\prime}-\frac{1}{x}S_{2}^{\prime}$ $\displaystyle=$
$\displaystyle-\Bigg{[}{2\hat{D}_{\eta}\over
x^{2}}+{\Gamma\over\sqrt{1-x^{2}}\,k(x)}\Bigg{]}S_{0}-\frac{2\hat{D}_{\eta}\,\left(\hat{D}_{\eta}-x^{2}\right)}{x^{3}}\,S_{0}^{\prime}\,\,.$
(D.20)

We can immediately integrate the equation for $S_{0}$ in (D.20), yielding

$\displaystyle S_{0}(x)\,=\,c_{1}{x^{2}\over 2}\,+\,c_{2}\ ,$ (D.21)

where $c_{1}$ and $c_{2}$ are integration constants. The differential equation
for $S_{2}$ in (D.20) contains terms which depend on the flavor expansion
parameter $\epsilon$. Since we are working at first order in $\epsilon$, it is
convenient to expand the transport coefficients in powers of $\epsilon$ as:

$\hat{D}_{\eta}\,=\,\hat{D}_{0}\,+\,\epsilon\,\hat{D}_{1}\,\,,\qquad\qquad\hat{\tau}_{s}\,=\,\hat{\tau}_{0}\,+\,\epsilon\,\hat{\tau}_{1}\,\,,$
(D.22)

as well as for the function $S_{2}$:

$S_{2}(x)\,=\,s(x)+\epsilon\,\tilde{s}(x)\,\,.$ (D.23)

Proceeding in this way we generate two second-order differential equations for
$s(x)$ and $\tilde{s}(x)$, which can be analytically integrated (these
functions have rather long expressions and will not be written here). This
integration generates two new integration constants $c_{3}$ and $c_{4}$ for
$s(x)$ and another two $c_{5}$ and $c_{6}$ for $\tilde{s}(x)$. Expanding our
analytic result for $s(x)$ near $x=0$, up to second order in $x$, we get

$\displaystyle s(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{(}2\hat{D}_{0}(c_{1}\hat{D}_{0}+c_{2})\log
x-c_{1}-c_{2}+2c_{4}\Big{)}$ (D.24)
$\displaystyle+\frac{1}{4}x^{2}\left(c_{1}\hat{D}_{0}(2\log
x-1)+c_{1}+2c_{2}\log\left(\frac{x}{2}\right)+c_{2}+2c_{3}\right)+\ldots\,\,.$

Similarly, $\tilde{s}(x)$ near $x=0$ takes the form

$\displaystyle\tilde{s}(x)$ $\displaystyle=$
$\displaystyle-\Big{[}(6\hat{D}_{1}+5i\pi-5-10\log
2)c_{1}+\big{(}(15+25i)\pi-70+40\log 2\big{)}c_{2}-12c_{5}\Big{]}{x^{2}\over
24}$ (D.25) $\displaystyle+\Big{[}{c_{1}\over
2}x^{2}+2c_{1}\hat{D}_{0}+c_{2}\Big{]}\,\hat{D}_{1}\,\log x\,+\,{1\over
12}\big{(}10i\pi+35-20\log 2)c_{1}$
$\displaystyle\qquad\qquad\qquad\qquad+{1\over 12}\big{(}5i\pi-5-10\log
2\big{)}c_{2}\,+\,c_{6}+\ldots\ .$

Combining all these equations, we get
$S(x)=S_{0}(x)+\hat{q}^{2}(s(x)+\epsilon\tilde{s}(x))$ for small $\hat{q}$ and
small $x$. Let us now obtain an expression by expanding $S(x)$ in opposite
order and demand that it should be consistent with the above. We first expand
near $x=0$ and write

$\displaystyle S(x)=1+\sigma_{2}\,x^{2}+\sigma_{4}\,x^{4}+\ldots\ .$ (D.26)

Plugging the expansion (D.26) into (D.15) we get $\sigma_{2}$ and $s_{4}$ in
terms of $\hat{\omega}$ and $\hat{q}$. Expanding this result for small
$\hat{q}$ and $\epsilon$, $\sigma_{2}$, we obtain

$\displaystyle\sigma_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\hat{D}_{0}}+\frac{1}{4}\hat{q}^{2}(2\hat{\tau}_{0}-1)+\epsilon\left(\frac{\hat{D}_{1}}{2\hat{D}_{0}^{2}}+\frac{\hat{q}^{2}\hat{\tau}_{1}}{2}\right)+\ldots$
$\displaystyle\sigma_{4}$ $\displaystyle=$
$\displaystyle-\frac{(4-12\hat{D}_{0})\hat{q}^{2}}{64\hat{D}_{0}}-\epsilon\frac{\hat{q}^{2}\left(5\pi\hat{D}_{0}^{2}-30\hat{D}_{0}^{2}+30\hat{D}_{0}^{2}\log
2-4\hat{D}_{1}\right)}{64\hat{D}_{0}^{2}}+\ldots\ .$ (D.27)

Thus, an expansion for small $x$ and small $\hat{q}$ (D.26) leads to the
following expressions for $S_{0}(x)$ and $S_{2}(x)$:

$\displaystyle S_{0}(x)$ $\displaystyle=$ $\displaystyle
1-\frac{x^{2}}{2\hat{D}_{0}}+\epsilon\,\frac{\hat{D}_{1}x^{2}}{2\hat{D}_{0}^{2}}$
$\displaystyle S_{2}(x)$ $\displaystyle=$
$\displaystyle{12x^{4}\hat{D}_{0}^{2}\,+\,16(2\hat{\tau}_{0}-1)\hat{D}_{0}^{2}\,x^{2}\,-\,4\hat{D}_{0}\,x^{4}\over
64\hat{D}_{0}^{2}}$ (D.28) $\displaystyle+{\epsilon\over
64\,\hat{D}_{0}^{2}}\,\Big{[}(30-5\pi-30\log
2)\hat{D}_{0}^{2}\,x^{4}\,+\,4\hat{D}_{1}\,x^{4}\,+\,32\hat{\tau}_{1}\,\hat{D}_{0}^{2}\,x^{2}\Big{]}\
.$

Let us now compare (D.28) and (D.21). Since the right-hand side of (D.21) does
not depend on $\epsilon$, we conclude that

$\hat{D}_{1}=0\,\,.$ (D.29)

Moreover, by comparing the expression for $S_{2}$ in (D.28) and the one that
follows from (D.24) and (D.25), we can relate the different integration
constants and the transport coefficients as

$\displaystyle c_{1}=-{1\over\hat{D}_{0}}\,\,,\qquad c_{2}=1\qquad
c_{3}=\frac{\hat{D}_{0}(2\hat{\tau}_{0}-3+\log 4)+1}{2\hat{D}_{0}}\,\,,\qquad
c_{4}={1\over 2}-{1\over 2\hat{D}_{0}}$ $\displaystyle
c_{5}=\frac{\hat{D}_{0}(12\hat{\tau}_{1}+\pi(15+25i)-70+20\log
4)+5(1-i\pi+\log 4)}{12\hat{D}_{0}}$ $\displaystyle
c_{6}=\frac{5(\hat{D}_{0}(2-2i\pi+\log(16))+14+4i\pi-\log(256))}{24\hat{D}_{0}}\,\,.$
(D.30)

With these constants, the function $S_{0}$ becomes

$S_{0}(x)=1-{1\over 2\,\hat{D}_{0}}x^{2}\,\,.$ (D.31)

Imposing the Dirichlet boundary condition $S_{0}(x=1)=0$ at $x=1$ leads to the
result for $\hat{D}_{0}$,

$\hat{D}_{0}={1\over 2}\,\,.$ (D.32)

Using this new condition and (D.1) we arrive at the following expressions for
$s(x)$ and $\tilde{s}(x)$

$\displaystyle s(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{(}x^{2}\left(\tau_{0}-\log\big(x/2\big{missing})\right)-1+\sqrt{1-x^{2}}+{1\over
2}\log\left(\frac{1-\sqrt{1-x^{2}}}{\sqrt{1-x^{2}}+1}\right)$
$\displaystyle+\left(1-x^{2}\right){\rm
arctanh}\left(\sqrt{1-x^{2}}\right)\Bigg{)}$ $\displaystyle\tilde{s}(x)$
$\displaystyle=$
$\displaystyle-\frac{5\left(1-x^{2}\right)^{3/4}}{x^{2}}\left(\frac{1}{3}\,F\left(1,1;\frac{7}{4};1-\frac{1}{x^{2}}\right)-\frac{1}{21}2\left(1-x^{2}\right)\,F\left(1,2;\frac{11}{4};1-\frac{1}{x^{2}}\right)+\frac{2x^{2}}{3}\right)$
(D.33)
$\displaystyle+\frac{1}{4}x^{2}(2\tau_{1}+5(\log(4)-2))-\sqrt{1-x^{2}}\left(\Omega\left(\frac{1}{(1-x^{2})^{{1\over
4}}}\right)+\frac{15}{4}\right)$
$\displaystyle+\frac{5}{4}\left(x^{2}+2\right){\rm
arctanh}\left(\sqrt{1-x^{2}}\right)-\frac{5}{4}(\pi-5+\log(4))\ .$

Requiring that $s(x=1)=0$ determines the value of $\hat{\tau}_{0}$, which
turns out to be

$\hat{\tau}_{0}=1-\log 2\ .$ (D.34)

Moreover, the Dirichlet boundary condition $\tilde{s}(x=1)=0$ is satisfied if
$\hat{\tau}_{1}$ is given by

$\hat{\tau}_{1}=\frac{5}{2}(\pi-3)\ .$ (D.35)

## References

  * [1] J. M. Maldacena, _The large N limit of superconformal field theories and supergravity_ , _Adv. Theor. Math. Phys._ 2 (1998) 231 [9711200].
  * [2] J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal and U. A. Wiedemann, _Gauge/String Duality, Hot QCD and Heavy Ion Collisions_. Cambridge University Press, 2014, 10.1017/CBO9781139136747, [1101.0618].
  * [3] J. McGreevy, _Holographic duality with a view toward many-body physics_ , _Adv. High Energy Phys._ 2010 (2010) 723105 [0909.0518].
  * [4] A. V. Ramallo, _Introduction to the AdS/CFT correspondence_ , _Springer Proc. Phys._ 161 (2015) 411 [1310.4319].
  * [5] E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisenstein and A. P. MacKenzie, _Nematic Fermi Fluids in Condensed Matter Physics_ , _Annual Review of Condensed Matter Physics_ 1 (2010) 153 [0910.4166].
  * [6] K. Yagi and N. Yunes, _Relating follicly-challenged compact stars to bald black holes: A link between two no-hair properties_ , _Phys. Rev._ D91 (2015) 103003 [1502.04131].
  * [7] K. Yagi and N. Yunes, _I-Love-Q Relations: From Compact Stars to Black Holes_ , _Class. Quant. Grav._ 33 (2016) 095005 [1601.02171].
  * [8] R. A. Janik and P. Witaszczyk, _Towards the description of anisotropic plasma at strong coupling_ , _JHEP_ 09 (2008) 026 [0806.2141].
  * [9] P. M. Chesler and L. G. Yaffe, _Horizon formation and far-from-equilibrium isotropization in supersymmetric Yang-Mills plasma_ , _Phys. Rev. Lett._ 102 (2009) 211601 [0812.2053].
  * [10] T. Azeyanagi, W. Li and T. Takayanagi, _On String Theory Duals of Lifshitz-like Fixed Points_ , _JHEP_ 06 (2009) 084 [0905.0688].
  * [11] D. Mateos and D. Trancanelli, _The anisotropic N=4 super Yang-Mills plasma and its instabilities_ , _Phys. Rev. Lett._ 107 (2011) 101601 [1105.3472].
  * [12] D. Mateos and D. Trancanelli, _Thermodynamics and Instabilities of a Strongly Coupled Anisotropic Plasma_ , _JHEP_ 07 (2011) 054 [1106.1637].
  * [13] J.-i. Koga, K. Maeda and K. Tomoda, _Holographic superconductor model in a spatially anisotropic background_ , _Phys. Rev._ D89 (2014) 104024 [1401.6501].
  * [14] S. Jain, N. Kundu, K. Sen, A. Sinha and S. P. Trivedi, _A Strongly Coupled Anisotropic Fluid From Dilaton Driven Holography_ , _JHEP_ 01 (2015) 005 [1406.4874].
  * [15] L. Cheng, X.-H. Ge and S.-J. Sin, _Anisotropic plasma at finite $U(1)$ chemical potential_, _JHEP_ 07 (2014) 083 [1404.5027].
  * [16] E. Banks and J. P. Gauntlett, _A new phase for the anisotropic N=4 super Yang-Mills plasma_ , _JHEP_ 09 (2015) 126 [1506.07176].
  * [17] D. Roychowdhury, _On anisotropic black branes with Lifshitz scaling_ , _Phys. Lett._ B759 (2016) 410 [1509.05229].
  * [18] D. Roychowdhury, _Holography for anisotropic branes with hyperscaling violation_ , _JHEP_ 01 (2016) 105 [1511.06842].
  * [19] V. Jahnke and A. S. Misobuchi, _Probing strongly coupled anisotropic plasmas from higher curvature gravity_ , _Eur. Phys. J. C_ 76 (2016) 309 [1510.03774].
  * [20] E. Banks, _Phase transitions of an anisotropic N=4 super Yang-Mills plasma via holography_ , _JHEP_ 07 (2016) 085 [1604.03552].
  * [21] D. Ávila, D. Fernández, L. Patiño and D. Trancanelli, _Thermodynamics of anisotropic branes_ , _JHEP_ 11 (2016) 132 [1609.02167].
  * [22] A. Donos, J. P. Gauntlett and O. Sosa-Rodriguez, _Anisotropic plasmas from axion and dilaton deformations_ , _JHEP_ 11 (2016) 002 [1608.02970].
  * [23] D. Giataganas, U. Gürsoy and J. F. Pedraza, _Strongly-coupled anisotropic gauge theories and holography_ , _Phys. Rev. Lett._ 121 (2018) 121601 [1708.05691].
  * [24] G. Itsios, N. Jokela, J. Järvelä and A. V. Ramallo, _Low-energy modes in anisotropic holographic fluids_ , _Nucl. Phys._ B940 (2019) 264 [1808.07035].
  * [25] I. Aref’eva and K. Rannu, _Holographic Anisotropic Background with Confinement-Deconfinement Phase Transition_ , _JHEP_ 05 (2018) 206 [1802.05652].
  * [26] A. Karch and A. O’Bannon, _Metallic AdS/CFT_ , _JHEP_ 09 (2007) 024 [0705.3870].
  * [27] T. Albash, V. G. Filev, C. V. Johnson and A. Kundu, _Finite temperature large N gauge theory with quarks in an external magnetic field_ , _JHEP_ 07 (2008) 080 [0709.1547].
  * [28] T. Albash, V. G. Filev, C. V. Johnson and A. Kundu, _Quarks in an external electric field in finite temperature large N gauge theory_ , _JHEP_ 08 (2008) 092 [0709.1554].
  * [29] J. Erdmenger, R. Meyer and J. P. Shock, _AdS/CFT with flavour in electric and magnetic Kalb-Ramond fields_ , _JHEP_ 12 (2007) 091 [0709.1551].
  * [30] E. D’Hoker and P. Kraus, _Magnetic Brane Solutions in AdS_ , _JHEP_ 10 (2009) 088 [0908.3875].
  * [31] E. D’Hoker and P. Kraus, _Charged Magnetic Brane Solutions in AdS (5) and the fate of the third law of thermodynamics_ , _JHEP_ 03 (2010) 095 [0911.4518].
  * [32] K. Jensen, A. Karch and E. G. Thompson, _A Holographic Quantum Critical Point at Finite Magnetic Field and Finite Density_ , _JHEP_ 05 (2010) 015 [1002.2447].
  * [33] N. Evans, A. Gebauer, K.-Y. Kim and M. Magou, _Phase diagram of the D3/D5 system in a magnetic field and a BKT transition_ , _Phys. Lett._ B698 (2011) 91 [1003.2694].
  * [34] E. D’Hoker and P. Kraus, _Holographic Metamagnetism, Quantum Criticality, and Crossover Behavior_ , _JHEP_ 05 (2010) 083 [1003.1302].
  * [35] K.-Y. Kim, B. Sahoo and H.-U. Yee, _Holographic chiral magnetic spiral_ , _JHEP_ 10 (2010) 005 [1007.1985].
  * [36] C. Hoyos, T. Nishioka and A. O’Bannon, _A Chiral Magnetic Effect from AdS/CFT with Flavor_ , _JHEP_ 10 (2011) 084 [1106.4030].
  * [37] M. Ammon, V. G. Filev, J. Tarrío and D. Zoakos, _D3/D7 Quark-Gluon Plasma with Magnetically Induced Anisotropy_ , _JHEP_ 09 (2012) 039 [1207.1047].
  * [38] U. Gürsoy, M. Järvinen, G. Nijs and J. F. Pedraza, _Inverse Anisotropic Catalysis in Holographic QCD_ , _JHEP_ 04 (2019) 071 [1811.11724].
  * [39] N. Evans, K.-Y. Kim, J. P. Shock and J. P. Shock, _Chiral phase transitions and quantum critical points of the D3/D7(D5) system with mutually perpendicular E and B fields at finite temperature and density_ , _JHEP_ 09 (2011) 021 [1107.5053].
  * [40] D. E. Kharzeev and H.-U. Yee, _Chiral helix in AdS/CFT with flavor_ , _Phys. Rev._ D84 (2011) 125011 [1109.0533].
  * [41] A. Donos, J. P. Gauntlett and C. Pantelidou, _Magnetic and Electric AdS Solutions in String- and M-Theory_ , _Class. Quant. Grav._ 29 (2012) 194006 [1112.4195].
  * [42] N. Jokela and A. V. Ramallo, _Universal properties of cold holographic matter_ , _Phys. Rev._ D92 (2015) 026004 [1503.04327].
  * [43] G. Itsios, N. Jokela and A. V. Ramallo, _Collective excitations of massive flavor branes_ , _Nucl. Phys._ B909 (2016) 677 [1602.06106].
  * [44] S. S. Gubser and S. S. Pufu, _The Gravity dual of a p-wave superconductor_ , _JHEP_ 11 (2008) 033 [0805.2960].
  * [45] M. Ammon, J. Erdmenger, M. Kaminski and P. Kerner, _Superconductivity from gauge/gravity duality with flavor_ , _Phys. Lett._ B680 (2009) 516 [0810.2316].
  * [46] P. Basu, J. He, A. Mukherjee and H.-H. Shieh, _Superconductivity from D3/D7: Holographic Pion Superfluid_ , _JHEP_ 11 (2009) 070 [0810.3970].
  * [47] N. Iizuka, S. Kachru, N. Kundu, P. Narayan, N. Sircar and S. P. Trivedi, _Bianchi Attractors: A Classification of Extremal Black Brane Geometries_ , _JHEP_ 07 (2012) 193 [1201.4861].
  * [48] A. Donos and J. P. Gauntlett, _Helical superconducting black holes_ , _Phys. Rev. Lett._ 108 (2012) 211601 [1203.0533].
  * [49] N. Iizuka and K. Maeda, _Study of Anisotropic Black Branes in Asymptotically anti-de Sitter_ , _JHEP_ 07 (2012) 129 [1204.3008].
  * [50] E. Conde, H. Lin, J. M. Penín, A. V. Ramallo and D. Zoakos, _D3/D5 theories with unquenched flavors_ , _Nucl. Phys._ B914 (2017) 599 [1607.04998].
  * [51] J. M. Penín, A. V. Ramallo and D. Zoakos, _Anisotropic D3-D5 black holes with unquenched flavors_ , _JHEP_ 02 (2018) 139 [1710.00548].
  * [52] N. Jokela, J. M. Penín, A. V. Ramallo and D. Zoakos, _Gravity dual of a multilayer system_ , _JHEP_ 03 (2019) 064 [1901.02020].
  * [53] U. Gran, N. Jokela, D. Musso, A. V. Ramallo and M. Tornsö, _Holographic fundamental matter in multilayered media_ , _JHEP_ 12 (2019) 38 [1909.01864].
  * [54] C. Hoyos, N. Jokela, J. M. Penín and A. V. Ramallo, _Holographic spontaneous anisotropy_ , _JHEP_ 04 (2020) 062 [2001.08218].
  * [55] C. Hoyos, N. Jokela, J. M. Penín, A. V. Ramallo and J. Tarrío, _Risking your NEC_ , _JHEP_ 10 (2021) 112 [2104.11749].
  * [56] A. Karch and E. Katz, _Adding flavor to AdS / CFT_ , _JHEP_ 06 (2002) 043 [hep-th/0205236].
  * [57] O. DeWolfe, D. Z. Freedman and H. Ooguri, _Holography and defect conformal field theories_ , _Phys. Rev._ D66 (2002) 025009 [hep-th/0111135].
  * [58] J. Erdmenger, Z. Guralnik and I. Kirsch, _Four-dimensional superconformal theories with interacting boundaries or defects_ , _Phys. Rev. D_ 66 (2002) 025020 [hep-th/0203020].
  * [59] S. K. Domokos and A. B. Royston, _Supersymmetry of the D3/D5 Defect Field Theory_ , 2207.07671.
  * [60] C. Nunez, A. Paredes and A. V. Ramallo, _Unquenched Flavor in the Gauge/Gravity Correspondence_ , _Adv. High Energy Phys._ 2010 (2010) 196714 [1002.1088].
  * [61] S. Kobayashi, D. Mateos, S. Matsuura, R. C. Myers and R. M. Thomson, _Holographic phase transitions at finite baryon density_ , _JHEP_ 02 (2007) 016 [hep-th/0611099].
  * [62] D. Mateos, S. Matsuura, R. C. Myers and R. M. Thomson, _Holographic phase transitions at finite chemical potential_ , _JHEP_ 11 (2007) 085 [0709.1225].
  * [63] F. Bigazzi, A. L. Cotrone, J. Mas, A. Paredes, A. V. Ramallo and J. Tarrío, _D3-D7 Quark-Gluon Plasmas_ , _JHEP_ 11 (2009) 117 [0909.2865].
  * [64] F. Bigazzi, A. L. Cotrone, J. Mas, D. Mayerson and J. Tarrío, _D3-D7 Quark-Gluon Plasmas at Finite Baryon Density_ , _JHEP_ 04 (2011) 060 [1101.3560].
  * [65] M. M. Caldarelli, R. Emparan and B. Van Pol, _Higher-dimensional Rotating Charged Black Holes_ , _JHEP_ 04 (2011) 013 [1012.4517].
  * [66] S. S. Gubser, _Einstein manifolds and conformal field theories_ , _Phys. Rev. D_ 59 (1999) 025006 [hep-th/9807164].
  * [67] V. Balasubramanian and P. Kraus, _A Stress tensor for Anti-de Sitter gravity_ , _Commun. Math. Phys._ 208 (1999) 413 [hep-th/9902121].
  * [68] C. Hoyos, N. Jokela, D. Rodríguez Fernández and A. Vuorinen, _Breaking the sound barrier in AdS/CFT_ , _Phys. Rev. D_ 94 (2016) 106008 [1609.03480].
  * [69] C. Ecker, C. Hoyos, N. Jokela, D. Rodríguez Fernández and A. Vuorinen, _Stiff phases in strongly coupled gauge theories with holographic duals_ , _JHEP_ 11 (2017) 031 [1707.00521].
  * [70] C. Hoyos, N. Jokela and A. Vuorinen, _Holographic approach to compact stars and their binary mergers_ , _Prog. Part. Nucl. Phys._ 126 (2022) 103972 [2112.08422].
  * [71] M. Natsuume, _Causal hydrodynamics and the membrane paradigm_ , _Phys. Rev. D_ 78 (2008) 066010 [0807.1392].
  * [72] J. M. Maldacena, _Wilson loops in large N field theories_ , _Phys. Rev. Lett._ 80 (1998) 4859 [hep-th/9803002].
  * [73] S.-J. Rey and J.-T. Yee, _Macroscopic strings as heavy quarks in large N gauge theory and anti-de Sitter supergravity_ , _Eur. Phys. J. C_ 22 (2001) 379 [hep-th/9803001].
  * [74] S. Ryu and T. Takayanagi, _Holographic derivation of entanglement entropy from AdS/CFT_ , _Phys. Rev. Lett._ 96 (2006) 181602 [hep-th/0603001].
  * [75] S. Ryu and T. Takayanagi, _Aspects of Holographic Entanglement Entropy_ , _JHEP_ 08 (2006) 045 [hep-th/0605073].
  * [76] W. Fischler, A. Kundu and S. Kundu, _Holographic Mutual Information at Finite Temperature_ , _Phys. Rev. D_ 87 (2013) 126012 [1212.4764].
  * [77] J. Erdmenger and N. Miekley, _Non-local observables at finite temperature in AdS/CFT_ , _JHEP_ 03 (2018) 034 [1709.07016].
  * [78] N. Jokela and A. Pönni, _Notes on entanglement wedge cross sections_ , _JHEP_ 07 (2019) 087 [1904.09582].
  * [79] N. Jokela and A. Pönni, _Towards precision holography_ , _Phys. Rev. D_ 103 (2021) 026010 [2007.00010].
  * [80] M. M. Wolf, F. Verstraete, M. B. Hastings and J. I. Cirac, _Area Laws in Quantum Systems: Mutual Information and Correlations_ , _Phys. Rev. Lett._ 100 (2008) 070502 [0704.3906].
  * [81] M. Headrick, _Entanglement Renyi entropies in holographic theories_ , _Phys. Rev. D_ 82 (2010) 126010 [1006.0047].
  * [82] S. A. Hartnoll and P. Petrov, _Electron star birth: A continuous phase transition at nonzero density_ , _Phys. Rev. Lett._ 106 (2011) 121601 [1011.6469].
  * [83] V. G. M. Puletti, S. Nowling, L. Thorlacius and T. Zingg, _Holographic metals at finite temperature_ , _JHEP_ 01 (2011) 117 [1011.6261].
  * [84] O. Bergman, N. Jokela, G. Lifschytz and M. Lippert, _Quantum Hall Effect in a Holographic Model_ , _JHEP_ 10 (2010) 063 [1003.4965].
  * [85] N. Jokela, M. Järvinen and M. Lippert, _A holographic quantum Hall model at integer filling_ , _JHEP_ 05 (2011) 101 [1101.3329].
  * [86] Y. Bea, N. Jokela, M. Lippert, A. V. Ramallo and D. Zoakos, _Flux and Hall states in ABJM with dynamical flavors_ , _JHEP_ 03 (2015) 009 [1411.3335].
  * [87] B. S. DiNunno, N. Jokela, J. F. Pedraza and A. Pönni, _Quantum information probes of charge fractionalization in large- $N$ gauge theories_, _JHEP_ 05 (2021) 149 [2101.11636].
  * [88] A. F. Faedo, D. Mateos and J. Tarrío, _Three-dimensional super Yang-Mills with unquenched flavor_ , _JHEP_ 07 (2015) 056 [1505.00210].
  * [89] N. Jokela, J. M. Penín and K. C. S. Rigatos, _Gravity duals for defect quivers in the Veneziano limit_ , 2112.14677.
  * [90] F. Benini, _A Chiral cascade via backreacting D7-branes with flux_ , _JHEP_ 10 (2008) 051 [0710.0374].