# Condensed-phase molecular representation to link structure and
thermodynamics in molecular dynamics

Bernadette Mohr<EMAIL_ADDRESS>Diego van der Mast  Tristan Bereau  Max
Planck Institute for Polymer Research, 55128 Mainz, Germany.  Van’t Hoff
Institute for Molecular Sciences and Informatics Institute, University of
Amsterdam, Amsterdam 1098 XH, The Netherlands.

###### Abstract

Materials design requires systematic and broadly applicable methods to extract
structure–property relationships. In the context of molecules, physics-based
representations are key ingredients for the subsequent machine learning of
various properties. We focus here on thermodynamic properties of (complex)
molecular liquids: establishing links between condensed-phase intermolecular
structures and thermodynamic properties. Compared to electronic properties,
they pose two challenges: ($i$) a Boltzmann-ensemble representation and ($ii$)
the relevance of the environment. We address these aspects by adapting the
Spectrum of London and Axilrod-Teller-Muto representation (SLATM). Atomic
representations describe both covalent and non-covalent interactions, which we
sum over all constituents of a molecule. We further ensemble-average the
representation to relate it to thermodynamic properties. As an application, we
focus here on a challenging biomolecular system: lipid selectivity of small
molecules in mitochondrial membranes. We rely on coarse-grained molecular
dynamics simulations. Dimensionality reduction of the ensemble SLATM using
principal components reveals simple, interpretable relationships between
structure and free energy. Our methodology offers a systematic route to
connect higher-order intermolecular interactions with thermodynamics.

##  I Introduction

Computational molecular design is rapidly becoming one of the most exciting
fields of our time thanks to its impressive developments and broad
applicability[1, 2, 3, 4, 5, 6]. The idea is simple: look for molecules or
materials with desirable properties. In practice, solving the underlying
inverse problem remains challenging, requiring ample computational resources,
and an approach that makes large use of the underlying physics and chemistry
at hand. Electronic properties have spearheaded the movement: quantum-
mechanical (QM) calculations (e.g., density-functional theory) over large
numbers of molecules have been successfully used in the context of machine
learning to predict various properties with increasing accuracy and
generalization [7]. In no small part is this success due to the development of
molecular representations: they exploit physically meaningful degrees of
freedom (e.g., $r^{-1}$ scaling for Coulombic interactions) and account for
symmetries via invariances [8]. The picture appears less optimistic when
considering thermodynamic properties, in particular in the context of
condensed-phase liquids.

While electronic properties typically focus on single molecules in the gas,
our consideration of thermodynamic properties brings two specificities:

  1. 1.

Thermodynamics underlines the role of conformational entropy. Beyond a static
structure, the diversity of conformations heavily impact the energetics;

  2. 2.

The condensed phase involves a molecule embedded in a dense environment,
highlighting the balance of covalent and non-covalent interactions.

Traditional approaches often rely on the radial distribution function (RDF) to
analyze the pairwise structure of a liquid. RDFs remain one of the hallmarks
of analyzing molecular dynamics (MD) and Monte Carlo simulations alike,
parametrizing coarse-grained (CG) force fields, and more recently, the
aforementioned materials design[9, 10, 11]. RDFs and machine learning-type
molecular representations share their emphasis on body-order interactions—the
latter having generalized the RDF’s pairwise focus to include higher-order
many-body interactions. Thus, it is to be expected that molecular
representations developed for electronic properties ought to be a solid
starting point for thermodynamics.

In this work, we adapt molecular representations to thermodynamic properties
of molecular liquids via the two criteria we identified above, the
conformational ensemble and the complex balance of covalent and non-covalent
interactions. We focus on the spectrum of London Axilrod-Teller-Muto (SLATM)
molecular representation, used predominantly in machine learning for QM
calculations[12, 13]. The SLATM representation accounts for interatomic
interactions expanded along the 1-, 2-, and 3-body contributions. Our first
requirement pertains to the encoding of an ensemble average. Previous studies
have proposed summing over representations of individual conformations to
reproduce the appropriate Boltzmann weight, but neglected the role of the
environment[14, 15]. SLATM, however, does not distinguish between covalent and
non-covalent interactions, it thus naturally fulfills our second requirement
for a condensed-phase description.

We demonstrate the utility of the SLATM representation for linking molecular
structure to a thermodynamic property in a complex condensed-phase liquid. We
base our investigation on coarse-grained molecular dynamics simulations of a
challenging biomolecular system: Lipid selectivity of small molecules in
mitochondrial membranes. The problem involves the subtle identification of
preferential interactions between two similar lipids: cardiolipin (CL) and
phosphatidylglycerol (PG)[16, 17, 18, 19, 20, 21, 22, 23]. In a previous
study, we devised a rigorous discovery pipeline combining CG simulations,
free-energy calculations, and active learning, which, taken together, led to
the identification of design rules[24]. We further showed that such CG
simulations can be used to propose small-molecule probes for experimental
validation, with exciting results both _in vitro_ and _in vivo_[25].

Our application introduced here consists of extracting structure–property
relationships for small-molecule probes that are selective to CL (against PG).
The ensemble SLATM representation will encode the structural aspects of small
molecules in their complex environments via the body-order decomposition. CL
selectivity is measured via a difference in free energy between the two
environments. Because the resulting molecular representation is high-
dimensional, we apply dimensionality reduction using principal component
analysis (PCA). PCA linearly relates the encoded features, offering good
performance in the low-data regime and thereby increased interpretability.

We propose the SLATM representation in combination with dimensionality
reduction to systematically investigate the connection between body-order
interactions, the molecular structure, and thermodynamic properties from MD
trajectories.

##  II Methods

###  II.1 Coarse-grained molecular dynamics trajectories

We generated the trajectory data used in this work during a previous study on
data-driven high throughput screening[24]. We measured the binding selectivity
of small, drug-like compounds to CL by calculating each candidate compound’s
CG partitioning free energies from a water phase into a CL membrane. We used
coarse-grained representations for all components of the systems. The
comparison with partitioning free energies into a membrane comprising the
phospholipid phosphatidylglycerol (PG) allowed us to monitor CL-selectivity
accurately since CL and PG are structurally very similar[19, 16, 26, 27]. The
all-atom and CG structures of two lipids highlight the major difference in the
CG representations, namely the Nda particle in the CL headgroup versus the P4
particle in the PG headgroup (Figure 1).

Figure 1: The structural formulae of cardiolipin (CL) and phosphatidylglycerol
(PG) next to their coarse-grained (CG) representation, highlighting the
position of the CG beads in the lipid structures and the physicochemical
properties represented by the corresponding bead type. The noteworthy
difference between the CG representations of CL and PG is the Nda bead
representing the CL headgroup vs. the P4 bead comprising the PG headgroup. The
structures were rendered using ChemSketch, and the images of the CG
representations were created with VMD.[28, 29]

####  II.1.1 Physics-based coarse-grained representation

Computational compound screening still suffers from insufficient experimental
reference measurements and the immense computational cost of probing the
chemical space of small organic molecules[6, 30]. An elegant method to remedy
this is to reduce the chemical space’s vastness by translating it into a
lower-dimensional set of coarse-grained components. At the same time, coarse-
graining directs special attention to the predominant intermolecular driving
forces that are an important aspect of molecular design applications[31, 32,
33, 34, 35, 36].
The lipids, water, and ion particles were represented using the Martini 2 CG
force field[37, 38, 39]. In our previous study, we constructed the candidate
compounds as graph structures from random combinations of different CG beads
types as nodes, connected by edges representing bonds, to ensure a
sufficiently even coverage of chemical compound space[24]. The bead types used
for the candidate compounds are defined by the reduced resolution _5+1_ bead-
type force field developed for the previous study, {T1, T2, T3, T4, T5,
Q0$\pm$}—a Martini alternative aimed at screening small molecules[40, 24]. We
limited the maximum number of beads per compound to five in order to stay
within the molecular weight limit defined in Lipinsky’s rule of five for drug-
likeness of small molecules[41]. An illustration of all bead types
parameterized in the applied force fields is included in Figure S1.
Combining the _5+1_ bead types with the Martini 2 bead types constituting the
CL and the PG model (Figure 1), and the central beads of the water and ion
models {POL, PQd}, leads to a set of 14 different beads in the trajectory
data: {T1, T2, T3, T4, T5, $\text{{Q0}}\pm$, Nda, P4, Qa, Na, C1, C3, POL,
PQd}.

####  II.1.2 Measuring target-selectivity as free-energy differences

As a measure for CL selectivity, we used the difference between the
partitioning free energies of a candidate compound into a CL and a PG membrane
from previous studies,

$\Delta\Delta G_{\text{CL}\rightarrow\text{PG}}=\Delta
G^{\text{CL}}_{\text{W}\rightarrow\text{M}}-\Delta
G^{\text{PG}}_{\text{W}\rightarrow\text{M}}.\text{\cite[cite]{[\@@bibref{Number}{mohr2022data,
kleinwachter2022clib}{}{}]}}$ (1)

We calculated the partitioning free energies by successively coupling the Van
der Waals and electrostatic interactions between each compound and a pure
water environment, as well as two systems containing a CL and a PG membrane,
respectively[42, 43, 44]. We applied 40 intermediate coupling steps for each
interaction type to ensure adequate sampling. Subsequently, we estimated the
free energies of solvation via the MBAR method and calculated the free energy
differences following Equation (1)[45, 46].

####  II.1.3 Preparation of the trajectories

In this work, we only analyze the trajectories generated in the last, fully
coupled state of the solvation free-energy calculations in each of the two
lipid environments. To prevent the candidate compounds from collapsing and the
coordinates of the CG beads subsequently overlapping, we now introduced
appropriate angles and constraints to the compound structures according to
their geometry (Figure S2). Finally, we generated new trajectories in both
lipid environments for all 439 compounds with a length of 20 nanoseconds after
minimizing and equilibrating the systems using Gromacs 2020[47]. The standard
simulation settings for the Martini force field were applied according to the
literature[48].
We corrected the resulting trajectories for periodic boundary conditions and
centered the entire systems around the compounds. To facilitate statistical
analysis, we selected 200 frames from each trajectory according to the highest
density of the compound’s center of mass position over the entire trajectory.
We preprocessed the trajectories using the MDAnalysis package[49, 50]. The
solutes were extracted together with all membrane, solvent and ion particles
surrounding their center of mass within the long-range interaction cutoff
distance of 1.1 nm from each of the 200 previously selected frames per
trajectory, again with tools from MDAnalysis.

###  II.2 Analysis

####  II.2.1 The Spectrum of London Axilrod-Teller-Muto representation

The Spectrum of London Axilrod-Teller-Muto (SLATM) representation describes a
molecule as a vector of atomic environments encoding all 1-, 2-, and 3-body
interactions occurring within the scope of the simulated long-range
interactions (Figure 2)[12, 13]. Akin to other molecular representations,
SLATM ignores the notion of covalent bonding. The representation is
translationally, rotationally, and permutationally invariant. The 1-body
interactions uniquely identify atom types (or, in this work, CG bead types) of
the force field. 2- and 3-body interactions are inspired by London dispersion
forces and the Axilrod-Teller-Muto potential[51, 52, 53]. For interactions not
found in a particular system, the corresponding section of the representation
is padded with zeros.

Figure 2: Schematic of a SLATM histogram representation and cartoon
representation of a CG molecule and two associated body-order interactions.
The SLATM vector is decomposed in 1-, 2-, and 3- body contributions: Particle
counts (purple), pairwise interactions (blue), and triplets (red),
respectively. Above: cartoon representation of a molecule and exemplary body-
order interactions around the T5 particle.

We worked with the "atomic" representation (aSLATM) that is centered around
one particle in a molecule (Figure 2), following an existing
implementation[54]. We adjusted the aSLATM parameters to scale from the
original settings for atoms to representing CG beads: sigmas = [0.3, 0.2]
(width of the Gaussian kernel), dgrids = [0.2, 0.2], (bin widths of the binned
2-body and 3-body interaction distances) rcut = 8.0 (maximum interaction
distance) and rpower = 6 (exponent of the leading order term in the London
potential)[54]. We ensemble average the aSLATM contribution of each CG in a
compound bead over 200 snapshots of the corresponding MD simulation. The
ensemble representation, $\langle\text{aSLATM}\rangle$, is summed over all
particles of the molecule to yield a single representation, which we denote
$\langle\mathcal{X}\rangle$.

####  II.2.2 Post-processing of the SLATM representations

Many-body interactions are evaluated over a distance r = (0, rcut], divided
into equal steps with specified step sizes. Consequently, at this stage we
obtained ensemble-averaged interaction spectra with multiple counts for each
many-body interaction at increasing distances from the evaluated particle
(Figure 2). In order to obtain a single value for each many-body interaction,
we calculated the weighted arithmetic mean over the corresponding spectrum
$\bar{x}=\sum_{i=1}^{n}w^{\prime}_{i}x_{i}$ (Equation S1).

The averaged observations of the individual interactions differ by several
orders of magnitude (Figure S3 (a) and (c)). We therefore scaled the
representations by taking the logarithm (Figure S3 (b) and (d)); further
details are included in the SI, section 2.2. Given an environment, say CL, the
representation consists of the logarithm of the ensemble-averaged SLATM with a
mean value for each observed many-body interaction,
$\ln\langle\mathcal{X}\rangle_{\text{CL}}$.

Focusing on the _difference_ in observed interactions between CL and PG
environments, our quantity of interest is the difference between the two log-
SLATMs (Figure S4), leading to

$\Delta\langle\mathcal{X}\rangle=\ln\frac{\langle\mathcal{X}\rangle_{\text{CL}}}{\langle\mathcal{X}\rangle_{\text{PG}}}.$
(2)

Our dataset consists of 439 compounds, for which we compute a
$\Delta\langle\mathcal{X}\rangle$, each of which is a vector of size 1480, to
account for all 14 CG bead types present.

####  II.2.3 Principal component analysis

Principal component analysis (PCA) aims at finding the subspace in a set of
multivariate samples with the highest overall variance, extending its
applicability from dimensionality reduction to feature selection and feature
extraction[55, 56, 57, 58]. The original variables are linearly transformed
into uncorrelated new variables corresponding to the highest variance in the
original samples, the principal components (PC). The lower-dimensional PCs can
be correlated to the features responsible for the corresponding fraction of
variance. In this manner, the interpretability of the information can be
increased while information loss is minimized[59, 2, 60]. Apart from the PCs,
PCA yields the corresponding eigenvectors and eigenvalues of the covariance
matrix of the original variables. The eigenvectors indicate the directions of
high variance in a set of samples, and the eigenvalues represent the
proportion of the variance explained by the PC corresponding to the respective
eigenvector. To highlight the contribution of each original variable on a PC,
we focus on the loadings

$l_{i}=\sqrt{\lambda_{i}}\Psi_{i},$ (3)

where $\lambda_{i}$ and $\Psi_{i}$ denote the $i^{\textrm{th}}$ eigenvalue and
eigenvector, respectively. The square root of the eigenvalue corresponds to a
scaling by the standard deviation of the underlying data.

We used the PCA implementation of the scikit-learn package with the random
seed set to a constant value for reproducibility[61, 62, 63]. Default settings
were used for the remaining hyperparameters.

##  III Results and Discussion

We employ PCA to obtain an unsupervised (low-dimensional) embedding that
facilitates analysis and interpretation of the originally high-dimensional MD
trajectories.

###  III.1 Low-dimensional embedding

The magnitude of the eigenvalues, $\\{\lambda_{i}\\}$, show little clear
variation but instead a virtually continuous decay[60]. The first six PCs
represent $\sim 77\%$ of the overall variance, which we focus on below. The
amount of variance explained by each of the first 10 PCs is shown in Figure
S5. Because all small molecules probe the same two lipid environments, the
1-body contributions cancel out from $\Delta\langle\mathcal{X}\rangle$. For
this reason, we limit our evaluation to the 2- and 3-body interactions. To
interpret the first six PCs, we cross-correlate them with different
physicochemical descriptors described in the following paragraphs. We use
linear regression to measure the correlation, quantified by the coefficient of
determination $R^{2}$[64].

The difference in ensemble-SLATM environments engineered in
$\Delta\langle\mathcal{X}\rangle$ only retains the difference in observed
interactions between the two lipid environments and is therefore likely to
tease out selectivity effects. Indeed we find PC3 to correlate significantly
with selectivity, as measured from the difference in free energy,
$\Delta\Delta G$ of Equation (1) (Figure 3 (a) and Figure S6). Simultaneously,
PC3 also most strongly correlates with hydrophobicity, as monitored through
the water–octanol partitioning free energy $\Delta
G_{\text{W}\rightarrow\text{Ol}}$ used in the parameterization of the _5+1_
bead-type force field (Figure S7)[40]. PC3 also largely correlates to the
ratio of polar beads in each solute and the number of charges per solute
(Figures S8 and S9). While we rationalize the physicochemical nature of PCs 3
to 6 below, it is not clear to us what PC1 and PC2 primarily represent.

Figure 3: Cross-correlation of (a) the third principal component (PC3) to CL
selectivity $\Delta\Delta G$, and (b) PC5 to the ratio of hydrogen-bonding
beads in each solute. The color gradients further visualize the respective
physicochemical descriptor represented on the vertical axis.

We find the strongest correlation between the ratio of hydrogen-bonding beads
per solute and PC5 (Figure 3 (b), Figure S10). The magnitude of the distance
vector $\Delta\langle\mathcal{X}\rangle$ ($L_{2}$-norm) quantifies the
difference between the observed interactions in the CL and the PG environment.
Here, we find the strongest correlation to PC6 (Figure S11). We also observe
additional, weaker correlations between the descriptors mentioned above and
the other PCs considered. Hydrophobicity also correlates to PC4 and PC5
(Figure S7). PC1 is correlated to the ratio of polar beads per compound
(Figure S8), and PC4 and PC5 to the ratio of charged beads per compound
(Figure S9). PC2 is correlated to the ratio of hydrogen bonding beads (Figure
S10), and PC5 is also correlated to the $L_{2}$-norm, Figure S11. The
correlation between PCs likely results from the low dimensionality of the CG
model itself, which heavily compresses the diversity of physicochemical
interactions in a few bead types.

###  III.2 From many-body interactions to molecular structures

We next investigate the set of interactions that critically impact the
resulting thermodynamics of binding to a specific lipid environment. We
construct a graph consisting of nodes and edges, where nodes are bead types,
and edges are the 2- and (projected) 3-body interactions (Figure 4). The
thickness of the edges is to represent the frequency of sampling of the
interaction. We compute this frequency by summing the occurrence of
interactions between the same CG bead types. Interactions with the 50 largest
absolute loading values, $|l_{i}|$, (Equation (3)) are included in the
analysis (Figures S12 and S13).

An illustration of the interactions analyzed to construct the graphs in Figure
4 is shown in Figure S14. All interactions analyzed involved either the Nda or
the P4 bead, emphasizing their relevance for selectivity engineering. Similar
graphs for the remaining PCs are shown in Figures S15 to S20.

The information contained in the graphs is also illustrated in Figure 5, where
we emulate a three-dimensional positioning of the solute beads around the
lipids. Here, we manually position the solute bead types around the lipid
beads they predominantly interact with. The arrow widths reflect the
interaction frequency in the same way we use to weigh the edges in 4.

Figure 4: Visualization of the dominant interactions of PC3 with negative
loadings. The graph in (a) shows interactions involving the CL headgroup bead
Nda, whereas the graph in (b) includes all interactions involving the PG
headgroup bead P4. The width of the edges scales with the occurrence of an
interaction between two beads in this subset. Edges highlighted in orange show
interactions between solute beads, and dashed edges represent interactions
with beads used to model water or sodium ions. Figure 5: Illustration of the
three-dimensional information extracted from the largest negative loadings in
PC3. The hypothetic solutes in (a) and (b) illustrate the dominant
interactions involving CL and PG, respectively (Figure 4).

###  III.3 Structure–property relationships

We visualize all pairwise combinations of the six main PCs and the maximum
coefficients of their respective eigenvectors as biplots to gain insight into
the low-dimensional embedding of$\Delta\langle\mathcal{X}\rangle$ (Figure 6
and Figures S21 to S32)[65, 66]. We find the same correlations to features as
in section III.1 (Figure 3 and Figures S6 to S11). Color gradients
representing the physicochemical descriptors used for cross-correlation are
most distinctly visible in plots where both PCs are correlated to the applied
descriptor. We find the clearest effect in the biplots of PC3 and PC5 as a
distinct separation between two groups of the 439 compounds. A clear
difference is visible between neutral and non-neutral compounds (Figure 6
(a)). The link to hydrophobicity is less visible (Figures S23 and S24). While
the separation by selectivity $\Delta\Delta G$ is much less pronounced, it is
still apparent that increasing CL selectivity is predominantly found from
positive to negative scores along PC3 (Figure 6 (b)). In contrast, they
cluster around zero along PC5. The Q0-Nda and Q0-P4 interactions point in
opposite directions of increased CL selectivity in the case of Q0-Nda and
reduced CL selectivity for Q0-P4. The solid diagonal lines along the gradients
of the two descriptors (i.e., the loadings) are both nearly parallel to the
dashed diagonal line. This geometric alignment leads to an appealing feature:
a linear relation between PC3 and PC5 acts as a proxy for selectivity.

This apparent linear relationship between PC3 and PC5 against CL selectivity
gives us the means to build a _proxy_ for further prediction. Figure 6 (c)
includes three test compounds, for which we _predicted_ CL selectivity, i.e.,
without expensive calculation of reference free-energy differences
$\Delta\Delta G$. We constructed the new compounds A and B to follow the
three-dimensional structural aspects we identify as indicative of CL
selectivity in section III.2 (Figures 4 (a) and 5 (a)). Compound C was
originally eliminated from the previous study because it was deemed unlikely
to properly insert at the membrane interface[24]. Here, we introduce it as a
negative test. We now generate one new trajectory for each of the three test
compounds in each lipid environment and translate these trajectories into
SLATM representations $\Delta\langle\mathcal{X}\rangle$. We transform the
three new $\Delta\langle\mathcal{X}\rangle$ using the PCA model of our
original set of 439 compounds. Compound C is featured in the bottom-right part
of the biplot, indicative of poor selectivity. Compounds A and B, on the other
hand, are on the other side of the boundary line tangent to the dashed
diagonal line, suggesting high CL selectivity. Indeed both compounds were
created to follow the structural properties identified in section III.2. The
positions of compound A and compound B along PC5 reflect its correlation to
the number of hydrogen bonding sites in the compounds. For both A and B, the
positions on the biplot fall outside the range of the training set. This
effect may be due to the weak correlation of PC3 to the ratio of hydrogen
bonding beads. Further, the positioning may also suggest undersampling of the
chemical space considered[67, 68].

Figure 6: Biplots of PC3 and PC5, colored by the ratio of charged beads per
solute (a), and the selectivity $\Delta\Delta G$ (b) and (c). The values of
the principal components are scaled to the interval $[-1,1]$. In (a) and (b),
we show the six highest eigenvector coefficients of 2-body interactions. In
(c), three example compounds are classified for CL selectivity by PCA, without
calculating their respective partitioning free-energy difference $\Delta\Delta
G$.

The predictive approach significantly reduces the required computational time
compared to the calculation of partitioning free energy differences
$\Delta\Delta G$, making it an exciting candidate for computational high-
throughput screening. The previously reported $\sim 24\text{ GPU h}$ for a
neutral compound and $\sim 48\text{ GPU h}$ for a charged compound on a Tesla
GPU reduced to $\sim 0.3\text{ GPU h}$ and $\sim 0.7\text{ GPU h}$ on an
Nvidia GTX 980, respectively, a decrease of $\sim$ 98.5%[24]. The selectivity
prediction through PCA also mostly eliminates the considerable difference in
computational cost between charged and neutral compounds. We can identify the
2-body and 3-body interactions that strongly link to CL selectivity and
predict new candidates through the direct linear interpretability of PCA based
on coarse-grained MD trajectories and the SLATM representations retaining the
relevant thermodynamical and structural information.

##  IV Conclusion

In the present work, we evaluate the Spectrum of London Axilrod-Teller-Muto
(SLATM) representation for solving the relationship between molecular
structure and thermodynamics. The SLATM representations provide a systematic
link between physicochemical interactions present in molecular dynamics (MD)
trajectories and body-order decomposed structural features. Our adapted
ensemble-SLATM representation addresses both aspects of non-covalent
interactions in condensed-phase environments and Boltzmann-averages for
thermodynamic properties. The methodology was applied to coarse-grained
molecular dynamics (MD) simulations of small molecules targeting phospholipid
cardiolipin (CL) membrane environments. Based on these trajectories, we
investigate the link between (ensemble) structure and target (thermodynamic)
selectivity across a subset of chemical space.

We translate snapshots of these trajectories into SLATM representations to
capture intermolecular interactions and average over the resulting conformers
to account for the statistical ensemble. Because direct analysis of the SLATM
representations is complex due to their high dimensionality, we apply
principal component analysis (PCA) for dimensionality reduction. Through PCA,
we are able to connect high variance in the SLATM representations to features
we had identified in the previous study. The lower-dimensional coordinates
correlated to these features allow us to extract selectivity-relevant 2-body
and 3-body intermolecular interactions, which we use to generate CG candidate
structures with potentially improved CL selectivity. By visualizing the main
principal coordinates pairwise, we observe a linear separation of points in
the lower-dimensional space between the coordinates correlated to CL
selectivity and the fraction of hydrogen-bonding sites per small molecule. We
use this linear relationship as the starting point for a predictive model to
propose candidate compounds with improved selectivity. Predicting novel
compounds based on a pre-trained PCA model provides a reduction of the
computational time by $\approx$98.5% compared to free energy calculations,
suggesting a potential use for virtual high-throughput screening endeavors.

Our results demonstrate a novel approach to solving the structure–property
relationship for thermodynamic properties from MD trajectories using the SLATM
representation and dimensionality reduction.

###### Acknowledgements.

We sincerely thank Ioana Ilie, Joseph Rudzinski and Jocelyne Vreede for
critical reading of the manuscript. We acknowledge support from the Sectorplan
Bèta & Techniek of the Dutch Government. This work was completed in part with
resources provided by the Dutch national e-infrastructure with the support of
SURF Cooperative.

## References

  * Curtarolo  _et al._ [2013] S. Curtarolo, G. L. Hart, M. B. Nardelli, N. Mingo, S. Sanvito, and O. Levy, Nature materials 12, 191 (2013).
  * Ferguson [2017] A. L. Ferguson, Journal of Physics: Condensed Matter 30, 043002 (2017).
  * Sidky  _et al._ [2020] H. Sidky, W. Chen, and A. L. Ferguson, Molecular Physics 118, e1737742 (2020).
  * Gkeka  _et al._ [2020] P. Gkeka, G. Stoltz, A. Barati Farimani, Z. Belkacemi, M. Ceriotti, J. D. Chodera, A. R. Dinner, A. L. Ferguson, J.-B. Maillet, H. Minoux, _et al._ , Journal of chemical theory and computation 16, 4757 (2020).
  * Dijkstra and Luijten [2021] M. Dijkstra and E. Luijten, Nature materials 20, 762 (2021).
  * Bereau [2021] T. Bereau, Modelling and Simulation in Materials Science and Engineering 29, 023001 (2021).
  * von Lilienfeld  _et al._ [2020] O. A. von Lilienfeld, K.-R. Müller, and A. Tkatchenko, Nature Reviews Chemistry 4, 347 (2020).
  * Musil  _et al._ [2021] F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Chemical Reviews 121, 9759 (2021).
  * Terban and Billinge [2021] M. W. Terban and S. J. Billinge, Chemical Reviews 122, 1208 (2021).
  * Peter and Kremer [2009] C. Peter and K. Kremer, Soft Matter 5, 4357 (2009).
  * Zhu  _et al._ [2021] H. Zhu, Y. Huang, J. Ren, B. Zhang, Y. Ke, A. K.-Y. Jen, Q. Zhang, X.-L. Wang, and Q. Liu, Advanced Science 8, 2003534 (2021).
  * Huang  _et al._ [2018] B. Huang, N. O. Symonds, and O. A. von Lilienfeld, arXiv preprint arXiv:1807.04259 (2018).
  * Huang  _et al._ [2020] B. Huang, N. O. Symonds, and O. A. von Lilienfeld, Handbook of Materials Modeling: Methods: Theory and Modeling , 1883 (2020).
  * Rauer and Bereau [2020] C. Rauer and T. Bereau, The Journal of chemical physics 153, 014101 (2020).
  * Weinreich  _et al._ [2021] J. Weinreich, N. J. Browning, and O. A. von Lilienfeld, The Journal of Chemical Physics 154, 134113 (2021).
  * Dudek [2017] J. Dudek, Frontiers in cell and developmental biology 5, 90 (2017).
  * Paradies  _et al._ [2014] G. Paradies, V. Paradies, V. De Benedictis, F. M. Ruggiero, and G. Petrosillo, Biochimica et Biophysica Acta (BBA)-Bioenergetics 1837, 408 (2014).
  * Elías-Wolff  _et al._ [2019] F. Elías-Wolff, M. Lindén, A. P. Lyubartsev, and E. G. Brandt, Soft Matter 15, 792 (2019).
  * Houtkooper and Vaz [2008] R. Houtkooper and F. Vaz, Cellular and Molecular Life Sciences 65, 2493 (2008).
  * Pennington  _et al._ [2019] E. R. Pennington, K. Funai, D. A. Brown, and S. R. Shaikh, Biochimica et Biophysica Acta (BBA)-Molecular and Cell Biology of Lipids 1864, 1039 (2019).
  * Paradies  _et al._ [2019] G. Paradies, V. Paradies, F. M. Ruggiero, and G. Petrosillo, Cells 8, 728 (2019).
  * Gonzalvez  _et al._ [2013] F. Gonzalvez, M. DÁurelio, M. Boutant, A. Moustapha, J.-P. Puech, T. Landes, L. Arnauné-Pelloquin, G. Vial, N. Taleux, C. Slomianny, _et al._ , Biochimica et Biophysica Acta (BBA)-Molecular Basis of Disease 1832, 1194 (2013).
  * Yi  _et al._ [2022] Q. Yi, S. Yao, B. Ma, and X. Cang, Journal of Lipid Research , 100227 (2022).
  * Mohr  _et al._ [2022] B. Mohr, K. Shmilovich, I. S. Kleinwächter, D. Schneider, A. L. Ferguson, and T. Bereau, Chemical Science 13, 4498 (2022).
  * Kleinwächter  _et al._ [2022] I. Kleinwächter, B. Mohr, A. Joppe, N. Hellmann, T. Bereau, H. D. Osiewacz, and D. Schneider, RSC Chemical Biology 3, 941 (2022).
  * Chen  _et al._ [2018] W.-W. Chen, Y.-J. Chao, W.-H. Chang, J.-F. Chan, and Y.-H. H. Hsu, Scientific reports 8, 1 (2018).
  * Pokorná  _et al._ [2016] L. Pokorná, P. Čermáková, A. Horváth, M. G. Baile, S. M. Claypool, P. Griač, J. Malínskỳ, and M. Balážová, Biochimica et Biophysica Acta (BBA)-Bioenergetics 1857, 34 (2016).
  * Advanced Chemistry Development Inc, Toronto [ON] C. Advanced Chemistry Development Inc, Toronto (ON), Acd/chemsketch freeware.
  * Humphrey  _et al._ [1996] W. Humphrey, A. Dalke, and K. Schulten, Journal of Molecular Graphics 14, 33 (1996).
  * Dobson  _et al._ [2004] C. M. Dobson _et al._ , Nature 432, 824 (2004).
  * Menichetti  _et al._ [2017] R. Menichetti, K. H. Kanekal, K. Kremer, and T. Bereau, The Journal of Chemical Physics 147, 125101 (2017).
  * Menichetti  _et al._ [2019] R. Menichetti, K. H. Kanekal, and T. Bereau, ACS central science 5, 290 (2019).
  * Voth [2008] G. A. Voth,  _Coarse-graining of condensed phase and biomolecular systems_ (CRC press, 2008).
  * Noid [2013] W. G. Noid, The Journal of chemical physics 139, 09B201_1 (2013).
  * Alessandri  _et al._ [2021] R. Alessandri, F. Grünewald, and S. J. Marrink, Advanced Materials 33, 2008635 (2021).
  * Kanekal  _et al._ [2022] K. H. Kanekal, J. F. Rudzinski, and T. Bereau, arXiv preprint arXiv:2203.07487 (2022).
  * Marrink  _et al._ [2007] S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. De Vries, The Journal of Physical Chemistry B 111, 7812 (2007).
  * Michalowsky  _et al._ [2017] J. Michalowsky, L. V. Schäfer, C. Holm, and J. Smiatek, The Journal of Chemical Physics 146, 054501 (2017).
  * Michalowsky  _et al._ [2018] J. Michalowsky, J. Zeman, C. Holm, and J. Smiatek, The Journal of Chemical Physics 149, 163319 (2018).
  * Kanekal and Bereau [2019] K. H. Kanekal and T. Bereau, The Journal of chemical physics 151, 164106 (2019).
  * Lipinski [2004] C. A. Lipinski, Drug Discovery Today: Technologies 1, 337 (2004).
  * Torrie and Valleau [1977] G. M. Torrie and J. P. Valleau, Journal of Computational Physics 23, 187 (1977).
  * Chipot and Pohorille [2007] C. Chipot and A. Pohorille,  _Free Energy Calculations_ (Springer, 2007).
  * Mey  _et al._ [2020] A. S. Mey, B. Allen, H. E. B. Macdonald, J. D. Chodera, M. Kuhn, J. Michel, D. L. Mobley, L. N. Naden, S. Prasad, A. Rizzi, _et al._ , arXiv preprint arXiv:2008.03067 (2020).
  * Shirts and Chodera [2008] M. R. Shirts and J. D. Chodera, The Journal of Chemical Physics 129, 124105 (2008).
  * [46] Beauchamp, Kyle A and Chodera, John D and Naden, Levi N and Shirts, Michael R, Python implementation of the multistate bennett acceptance ratio (mbar), published under the MIT license.
  * Abraham  _et al._ [2015] M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1, 19 (2015).
  * De Jong  _et al._ [2016] D. H. De Jong, S. Baoukina, H. I. Ingólfsson, and S. J. Marrink, Computer Physics Communications 199, 1 (2016).
  * Richard J. Gowers  _et al._ [2016] Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan DomaÅ„ski, David L. Dotson, SÃ©bastien Buchoux, Ian M. Kenney, and Oliver Beckstein, in  _Proceedings of the 15th Python in Science Conference_, edited by Sebastian Benthall and Scott Rostrup (2016) pp. 98 – 105\.
  * Michaud-Agrawal  _et al._ [2011] N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein, Journal of computational chemistry 32, 2319 (2011).
  * London [1930] F. London, Zeitschrift für Physik 63, 245 (1930).
  * Axilrod and Teller [1943] B. Axilrod and E. Teller, The Journal of Chemical Physics 11, 299 (1943).
  * Muto [1943] Y. Muto, J. Phys. Math. Soc. Jpn 17, 629 (1943).
  * Christensen  _et al._ [2017] A. Christensen, F. Faber, B. Huang, L. Bratholm, A. Tkatchenko, K. Muller, and O. von Lilienfeld, URL https://github. com/qmlcode/qml (2017).
  * Pearson [1901] K. Pearson, The London, Edinburgh, and Dublin philosophical magazine and journal of science 2, 559 (1901).
  * Hotelling [1933] H. Hotelling, Journal of educational psychology 24, 417 (1933).
  * Wold  _et al._ [1987] S. Wold, K. Esbensen, and P. Geladi, Chemometrics and intelligent laboratory systems 2, 37 (1987).
  * Jolliffe [2002] I. Jolliffe,  _Principal Component Analysis_ , Springer Series in Statistics (Springer, 2002).
  * Jolliffe and Cadima [2016] I. T. Jolliffe and J. Cadima, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 20150202 (2016).
  * Glielmo  _et al._ [2021] A. Glielmo, B. E. Husic, A. Rodriguez, C. Clementi, F. Noé, and A. Laio, Chemical Reviews 121, 9722 (2021).
  * Tipping and Bishop [1999] M. E. Tipping and C. M. Bishop, Neural computation 11, 443 (1999).
  * Bishop [2016] C. Bishop,  _Pattern Recognition and Machine Learning_ , Information Science and Statistics (Springer New York, 2016).
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Journal of Machine Learning Research 12, 2825 (2011).
  * Camacho  _et al._ [2010] J. Camacho, J. Picó, and A. Ferrer, Chemometrics and Intelligent Laboratory Systems 100, 48 (2010).
  * Gower and Hand [1995] J. C. Gower and D. J. Hand,  _Biplots_ , Vol. 54 (CRC Press, 1995).
  * Blasius  _et al._ [2009] J. Blasius, P. H. Eilers, and J. Gower, Computational Statistics & Data Analysis 53, 3145 (2009).
  * Hughes [1968] G. Hughes, IEEE transactions on information theory 14, 55 (1968).
  * Chandrasekaran and Jain [1974] B. Chandrasekaran and A. K. Jain, IEEE Transactions on Computers 100, 102 (1974).