11institutetext: NVIDIA
11email<EMAIL_ADDRESS>
# Automated 3D Segmentation of Kidneys and Tumors in MICCAI KiTS 2023
Challenge

Andriy Myronenko  Dong Yang  Yufan He  Daguang Xu

###### Abstract

Kidney and Kidney Tumor Segmentation Challenge (KiTS) 2023 [6] offers a
platform for researchers to compare their solutions to segmentation from 3D
CT. In this work, we describe our submission to the challenge using automated
segmentation of Auto3DSeg111https://monai.io/apps/auto3dseg available in
MONAI222https://github.com/Project-MONAI/MONAI. Our solution achieves the
average dice of 0.835 and surface dice of 0.723, which ranks first and wins
the KiTS 2023 challenge333https://kits-challenge.org/kits23/#kits23-official-
results.

###### Keywords:

Auto3DSeg MONAI Segmentation.

##  1 Introduction

Almost half a million people are diagnosed with kidney cancer annually. Each
year, a larger number of kidney tumors are detected, and currently, it is
difficult to determine whether a tumor is malignant or benign using
radiographic methods. The risk of metastatic progression remains a serious
concern, highlighting the need for reliable systems to objectively
characterize kidney tumor images and predict treatment outcomes.

For almost five years, the KiTS [5] initiative has maintained and expanded a
publicly available collection of hundreds of segmented CT scans featuring
kidney tumors. This year’s KiTS’23 [6] competition includes an expanded
training set consisting of 489 cases. The goal of the challenge is to develop
an automated method to segment kidneys, tumors and cysts.

##  2 Methods

We implemented our approach with MONAI [1] using Auto3DSeg open source
project. Auto3DSeg is an automated solution for 3D medical image segmentation,
utilizing open source components in MONAI, offering both beginner and advanced
researchers the means to effectively develop and deploy high-performing
segmentation algorithms.

The minimal user input to run Auto3DSeg for KiTS’23 datasets, is

⬇

1#!/bin/bash

2python -m monai.apps.auto3dseg AutoRunner run \

3 \--input="./input.yaml"

where a user provided input config (input.yaml) includes only a few lines:

⬇

1# This is the YAML file "input.yaml"

2modality: CT

3datalist: "./dataset.json"

4dataroot: "/data/kits23"

5

6class_names:

7- { "name": "kidney_and_mass", "index": [1,2,3] }

8- { "name": "mass", "index": [2,3] }

9- { "name": "tumor", "index": [2] }

10sigmoid : true

When running this command, Auto3DSeg will analyze the dataset, generate
hyperparameter configurations for several supported algorithms, train them,
and produce inference and ensembling. The system will automatically scale to
all available GPUs and also supports multi-node training.

The 3 minimum user options (in input.yaml) are data modality (CT in this
case), location of the downloaded KiTS’23 dataset (dataroot), and the list of
input filenames with an associated fold number (dataset.json). We generate the
5-fold split assignments randomly. Since KiTS defines its specific label
mapping (from integer class labels to 3 subregions, see Fig. 1), we have to
define it in the config, and since these subregions are overlapping, we use
”sigmoid: true” to indicate multi-label segmentation, where the final
activation is sigmoid (instead of the default softmax).

Currently, the default Auto3DSeg setting trains three 3D segmentation
algorithms: SegResNet [8], DiNTS [4] and SwinUNETR [3, 9] with their unique
training recipes. SegResNet and DiNTS are convolutional neural network (CNN)
based architectures, whereas SwinUNETR is based on transformers. Each is
trained using 5-fold cross validation.

Figure 1: KiTS’23 data example of an axial slice with the provided annotations
of kidneys (red), tumors (green) and cysts (blue). The classes of interest
that KiTS tasks to segment are: a) all foreground combined b) tumors + cysts
(green and blue) c) tumors only (green).

For model inference, a sliding-window scheme is used to create probability
maps, which are re-sampled back to its original spacing. This allows
ensembling prediction from different algorithms even if there were trained at
different resolutions

The simplicity of Auto3DSeg is a very minimal user input, which allows even
non-expert users to achieve a great baseline performance. The system will take
care of most of the heavy lifting to analyze, configure and optimally utilize
the available GPU resources. And for expert users, there are many
configuration options that can be manually provided to override the automatic
values, for better performance tuning.

In the final prediction, we ensemble the best model checkpoints only from
SegResNet and DiNTS algorithms, since they performed better during cross-
validation. We also applied a few small customizations to the baseline
Auto3DSeg workflow. We describe the baseline Auto3DSeg method and the
customization below.

###  2.1 Training and Validation Data

Our submission made use of the official KiTS’23 training set alone.

###  2.2 SegResNet

SegResNet444https://docs.monai.io/en/stable/networks.html is an encode-decoder
based semantic segmentation network based on [8]. It is a U-net based
convolutional neural network with deep supervision (see Figure 2).

Figure 2: SegResNet network configuration. The network uses repeated ResNet
blocks with batch normalization and deep supervision

The default Auto3DSeg SegResNet configuration was used, which includes 5
levels of 1, 2, 2, 4, 4 blocks. It follows a common CNN approach to downsize
image dimensions by 2 progressively (down to $16\times$ smaller) and
simultaneously increase feature size by. All convolutions are $3\times 3\times
3$ with an initial number of filters equal to 32. The encoder is trained with
a $256\times 256\times 256$ input region. The decoder structure is similar to
the encoder one, but with a single block per each spatial level. Each decoder
level begins with upsizing with transposed convolution: reducing the number of
features by a factor of 2 and doubling the spatial dimension, followed by the
addition of encoder output of the equivalent spatial level. The number of
levels and the region size is automatically configured. We use spatial
augmentation including random affine and flip in all axes,random intensity
scale, shift, noise and blurring. We use the dice loss, and sum it over all
deep-supervision sublevels:

$Loss=\sum_{i=0}^{4}\frac{1}{2^{i}}Loss(pred,target^{\downarrow})$ (1)

where the weight $\frac{1}{2^{i}}$ is smaller for each sublevel (smaller image
size) $i$. The target labels are downsized (if necessary) to match the
corresponding output size using nearest neighbor interpolation

We use the AdamW optimizer with an initial learning rate of $2e^{-4}$ and
decrease it to zero at the end of the final epoch using the Cosine annealing
scheduler. We use batch size of 1 per GPU, We use weight decay regularization
of $1e^{-5}$. Input images were re-scaled from $\left[-54,242\right]$ to the
$\left[-1,1\right]$ CT interval, followed by a sigmoid function. The range was
determined automatically by the data analysis step to include the intensity
pattern variations within the foreground regions.

###  2.3 DiNTS

DiNTS stands for Differentiable Network Topology Search (DiNTS) scheme, an
advanced methodology that fosters more dynamic topologies and an integrated
two-level search. DiNTS has demonstrated superior performance, achieving top-
tier results in the Medical Segmentation Decathlon (MSD) challenge [2].

The DiNTS algorithm utilizes a densely-connected lattice-based network,
training with a $96\times 96\times 96$ model input for both training and
inference. It leverages automatic mixed precision (AMP) and the SGD optimizer,
with an initial learning rate of 0.2 and a loss defined by the Dice plus focal
Loss. For the data processing, we utilize intensity normalization and random
cropping, as well as random rotation, zoom, Gaussian smoothing, intensity
scaling and shifting, Gaussian noising, and flipping.

In the quest for enhanced model performance, we fine-tuned the checkpoints
that were initially trained using the default training recipe, adjusting them
with various patch sizes for a span of 25 epochs. Our observations suggest
that a larger patch size often leads to improved model performance. Taking
into account our computational budget, we have selected a patch size from a
range between $192^{3}$ and $192\times 192\times 288$ for each fold of the
model (based on validation Dice scores) as the configuration in our final
model inference.

###  2.4 Metrics

The output of the network has 3 channels followed by a sigmoid, to segment
each of the 3 KiTS’23 expected classes: a) Kidney + Tumor + Cyst; b) Tumor +
Cyst; c) Tumor only. This creates a multi-label segmentation, where each voxel
can belong to more than one label. We use an average dice metric of these 3
classes to select the best validation checkpoints (without considering surface
distances).

###  2.5 Auto3DSeg customizations for KiTS’23

Even though the default Auto3DSeg configuration achieved a good baseline
cross-validation performance automatically, we did a few customization
including cropping to kidneys region and post-processing.

Training on the full size 3D CT images can be time consuming, so we pre-
cropped the images around the kidneys region. A simple rectangular box was
used, based on the ground truth labels. Since we used only 1 cropping per
image, the cropped region included not only kidneys but everything in between
including the spine. Training on such cropped images has 2 advantages:
firstly, it allowed for faster training, since smaller images can be cached in
RAM, and secondly, it simplified the task for the network. The disadvantage of
such an approach is that it requires finding the bounding box of the kidneys
region first.

We trained a separate segmentation network to find the foreground and
calculate the kidneys bounding box coordinates. For all the tasks we used the
same exact network architecture, trained all at $0.78\times 0.78\times
0.78mm^{3}$ CT resolution. Arguably, bounding box detection could have been
done faster, using a simpler detection network and at a lower CT resolution,
but here we saved on coding time, by reusing the framework. We trained the
first round of models fast (using a smaller number of epochs), to be used as a
bounding box detector. And after that trained longer, using only the cropped
(around kidneys) regions. This approach is somewhat similar in spirit to the
KiTS 2021 champion solution of coarse-to-fine training [10] (based on nnU-net
[7]), but here we do not re-use or concatenate masks detected at a coarse
level, we simply use it to detect the bounding box for faster training.

We also added binary post-processing on the final segmentation masks. Firstly,
we remove small connected components (smaller then 100 voxels total) based on
the foreground (merged labels). Secondly, we correct for ”outline” of some
tumor region. During a training stage, we noticed that on a small set of
images, network predictions of the tumor label have a small rim (1-2 voxels)
of cyst label. This happened mostly because the network was trained as a
multi-label task, where each voxel can be assigned to several classes. Since,
it’s not possible for a tumor to have a ”cystic” outline (even if it looks
like a cyst image pattern), by definition, we decided to correct such cases
with a simple binary post-processing. In our cross-validation tests, this
final post-processing did not actually affect the accuracy metrics, but we
still decided to include it.

Finally. we increased the size of the network input patch during training to
$256\times 256\times 256$ for SegResNet and to $192\times 192\times 288$ for
DiNTS, which allowed for faster training and also slightly increased the
cross-validation performance.

###  2.6 Optimization

We train the method on an 8-GPU 48GB NVIDIA A40 machine, with a batch size of
1 per GPU, which is equivalent to batch size of 8 single GPU training.
Auto3DSeg caches on-the-fly all the resampled data in RAM during the first
training epoch, when sufficient amount of RAM is available (otherwise a
fraction of the data is cached). This way only the first epoch suffers a slow-
down due to disk i/o and resampling, and the rest of the training process is
fast.

##  3 Results

Based on our random 5-fold split, the average dice scores per fold are shown
in Table 1. For the final submission we used an ensemble of 15 models: 10
models of SegResNet (5 folds trained twice), and 5 models of DiNTS. SegResNet
A and B training runs in Table 1 had the same configurations.

| Fold 1 | Fold 2 | Fold 3 | Fold 4 | Fold 5 | Average
---|---|---|---|---|---|---
SegResNet A | 0.8997 | 0.8739 | 0.8923 | 0.8911 | 0.8892 | 0.88924
SegResNet B | 0.8995 | 0.8773 | 0.8913 | 0.889 | 0.8865 | 0.88872
DiNTS | 0.8810 | 0.8647 | 0.8806 | 0.8752 | 0.8822 | 0.8767
Table 1: Avgerage Dice results of the 15 trained models based on our 5-fold
data split.

On the final hidden challenge dataset, our submission achieved an average Dice
score of 0.835, which ranked first among other submissions.

##  4 Conclusion

We described our winning solution to KiTS 2023 challenge using Auto3DSeg from
MONAI. Our final submission is en ensemble of 15 CNN models, 10 of SegResNet
and 5 of DiNTS. We hope that open source tools in MONAI will help more
researchers to achieve good baseline 3D segmentation results on their
particular task. Our solution achieves the average dice of 0.835 and surface
dice of 0.723, which ranks first on the KiTS 2023 leaderboard555https://kits-
challenge.org/kits23/#kits23-official-results.

## References

  * [1] Project-monai/monai, https://doi.org/10.5281/zenodo.5083813
  * [2] Antonelli, M., Reinke, A., Bakas, S., Farahani, K., Kopp-Schneider, A., Landman, B.A., Litjens, G., Menze, B., Ronneberger, O., Summers, R.M., et al.: The medical segmentation decathlon. Nature communications 13(1), 4128 (2022)
  * [3] Hatamizadeh, A., Nath, V., Tang, Y., Yang, D., Roth, H.R., Xu, D.: Swin unetr: Swin transformers for semantic segmentation of brain tumors in mri images. In: International MICCAI Brainlesion Workshop. pp. 272–284. Springer (2021)
  * [4] He, Y., Yang, D., Roth, H., Zhao, C., Xu, D.: Dints: Differentiable neural network topology search for 3d medical image segmentation. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 5841–5850 (2021)
  * [5] Heller, N., Isensee, F., Maier-Hein, K.H., Hou, X., Xie, C., Li, F., Nan, Y., Mu, G., Lin, Z., Han, M., et al.: The state of the art in kidney and kidney tumor segmentation in contrast-enhanced ct imaging: Results of the kits19 challenge. Medical Image Analysis 67, 101821 (2021)
  * [6] Heller, N., Wood, A., Isensee, F., Rädsch, T., Tejpaul, R., Papanikolopoulos, N., Weight, C.: The 2023 kidney and kidney tumor segmentation challenge, https://kits-challenge.org/kits23/
  * [7] Isensee, F., Jaeger, P.F., Kohl, S.A.A., Petersen, J., Maier-Hein, K.H.: nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation. Nature Methods 18, 203–211 (2021)
  * [8] Myronenko, A.: 3D MRI brain tumor segmentation using autoencoder regularization. In: International MICCAI Brainlesion Workshop. pp. 311–320. Springer (2018)
  * [9] Tang, Y., Yang, D., Li, W., Roth, H.R., Landman, B., Xu, D., Nath, V., Hatamizadeh, A.: Self-supervised pre-training of swin transformers for 3d medical image analysis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 20730–20740 (2022)
  * [10] Zhao, Z., Chen, H., Wang, L.: A coarse-to-fine framework for the 2021 kidney and kidney tumor segmentation challenge. In: Heller, N., Isensee, F., Trofimova, D., Tejpaul, R., Papanikolopoulos, N., Weight, C. (eds.) Kidney and Kidney Tumor Segmentation. pp. 53–58 (2022)


# Solvent Effects on Extractant Conformational Energetics in Liquid-Liquid
Extraction: A Simulation Study of Molecular Solvents and Ionic Liquids

Xiaoyu Wang  Chemical Sciences and Engineering Division, Argonne National
Laboratory, 9700 S Cass Ave, Lemont, IL 60439<EMAIL_ADDRESS>Srikanth
Nayak  Chemical Sciences and Engineering Division, Argonne National
Laboratory, 9700 S Cass Ave, Lemont, IL 60439  Richard E. Wilson  Chemical
Sciences and Engineering Division, Argonne National Laboratory, 9700 S Cass
Ave, Lemont, IL 60439  L. Soderholm  Chemical Sciences and Engineering
Division, Argonne National Laboratory, 9700 S Cass Ave, Lemont, IL 60439
Michael J. Servis  Chemical Sciences and Engineering Division, Argonne
National Laboratory, 9700 S Cass Ave, Lemont, IL 60439<EMAIL_ADDRESS>
###### Abstract

Extractant design in liquid-liquid extraction (LLE) is a research frontier of
metal ion separations that typically focuses on the direct extractant-metal
interactions. However, a more detailed understanding of energetic drivers of
separations beyond primary metal coordination is often lacking, including the
role of solvent in the extractant phase. In this work, we propose a new
mechanism for enhancing metal-complexant energetics with nanostructured
solvents. Using molecular dynamics simulations with umbrella sampling, we find
that the organic solvent can reshape the energetics of the extractant’s
intramolecular conformational landscape. We calculate free energy profiles of
different conformations of a representative bidentate extractant,
n-octyl(phenyl)-N,N-diisobutyl carbamoyl methyl phosphinoxide (CMPO), in four
different solvents: dodecane, tributyl phosphate (TBP), and dry and wet ionic
liquid (IL) 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide
([EMIM][Tf_2N]). By promoting reorganization of the extractant molecule into
its binding conformation, our findings reveal how particular solvents can
ameliorate this unfavorable step of the metal separation process. In
particular, the charge alternating nanodomains formed in ILs substantially
reduce the free energy penalty associated with extractant reorganization.
Importantly, using alchemical free energy calculations, we find that this
stabilization persists even when we explicitly include the extracted cation.
These findings provide insight into the energic drivers of metal ion
separations and potentially suggest a new approach to designing effective
separations using a molecular-level understanding of solvent effects.

###### keywords:

American Chemical Society, LaTeX

[supp-]supplementary Uranyl,CMPO,Ionic Liquid,Molecular Dynamics,IR

##  1 Introduction

Liquid-liquid extraction (LLE) leverages the relative solubilities of chemical
species between two immiscible phases, usually an aqueous and an organic
phase, to selectively recover a wide range of critical materials, including
rare earths, actinides, and platinum group elements. 1, 2, 3 In
hydrometallurgical processes, the extractant, an amphiphilic metal complexant,
selectively binds the targeted ions and solubilizes them in the nonpolar
organic phase. LLE is a free energy driven process, relying on small
differences in metal solvation free energies between the two phases. For many
separation processes, these differences are often only a few kJ/mol: many
large energetic factors contribute to the extraction and cancel out, resulting
in a total free energy of extraction that can be much smaller than each
contribution. For example, enthalpic metal ion binding with the complexant
generally favors phase transfer, whereas extractant-solvent interactions,
extractant conformational reorganization4 and self-assembly5 into molecular
complexes inhibit the process.

Approaches to rational process design are typically simplified by focusing on
the most significant enthalpic interactions. For example, recent studies have
focused on improving the chelating functionality of the complexant for
targeted metal recovery. 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 Such works
typically focus on understanding and improving local metal coordination
binding and extractant conformational energetics. Although not all
interactions are explicitly and properly considered when following this
approach, the role of organic phase hierarchical structuring has received
limited recent attention.16, 17, 18, 19 In general, complexant-solvent
interactions are accounted for with simple continuum models, despite the well-
documented solvent effects on metal separation energetics. 7, 20, 5, 21, 22
Herein we discuss one such structural factor—multidentate-complexant
conformation stability—and its potential impact on the energetics driving a
separation.

Multidentate amphiphilic molecules are a common class of extractants that
utilize chelation effects to enhance the energetics of separation and to
improve the selectivity of metal ions. 23, 24, 15, 25 However, for many
multidentate extractants, the binding conformation is not the most favorable
for the free, unbound molecule. 23 Generally, multidentate complexants with
chelating oxygen atoms bind metal ions with oxygen atoms aligned in the cis
conformation. In contrast, free (unbound) extractant molecules typically favor
the trans conformation, where the repulsion between electronegative oxygen
atoms is minimized. For example, Kapoor et al. studied single-crystal
structures of unbound CMPO, as well as its chelate complex with uranyl
nitrate, and found that the CMPO in the urnayl chelate complex exhibits a cis
conformation, while the unbound CMPO remains at a trans conformation. 26 To
quantify the energetic differences between the two molecular conformations—cis
and trans—we consider a representative bidentate extractant with well-
characterized separation behavior, CMPO. The molecular structure and the cis,
gauche and trans conformations are shown in Figure 1. Our gas phase density
functional theory (DFT) calculations (see below) find that the potential
energy difference between cis and trans CMPO is significant (up to 50 kJ/mol)
compared to the overall separation free energy. Such a high energy barrier
strongly suggests that, though the cis conformation is relevant while
bidentate chelating the metal ion, the trans conformation dominates for the
free complexant in typical nonpolar organic solvents such as dodecane.27, 28

Figure 1: The cis, gauche and trans conformations of CMPO are shown, with
corresponding oxygen-oxygen distances. Color scheme: red – oxygen, blue –
nitrogen, gray – carbon, orange – phosphorus.

One ligand design strategy used to avoid the energetic cost of the
extractant’s conformation is to fix the extractant in the chelating
conformation by adding a rigid backbone 6, which comes at the cost of
flexibility in the binding pocket and leads to more complicated synthetic
routes for the extractant molecule. Here, we suggest that the energetic cost
of such conformational reorganization can be reduced instead by tuning
solvent-extractant interactions. To broaden our understanding of the potential
role of the organic solvent, we first consider two conventional solvents,
dodecane and TBP, specifically chosen to cover a range of polarities, and in
addition they are well-characterized and commonly employed with CMPO in LLE
processes. 29 Furthermore, we consider two model IL systems – dry and “wet”
(water saturated) [EMIM][Tf_2N]. ILs are a class of solvents that have
received substantial attention in recent years due in part to their unique
liquid nanostructure:30, 31, 32, 33 they form polar/apolar and
cationic/anionic domains, 30, 34, 31, 19 resulting in drastically different
nanostructuring compared to traditional organic molecular solvents. ILs have a
broad appeal, finding their use in a wide range of applications through highly
customizable organic cations and anions, which enable tailoring of
physicochemical properties based on specific application needs. They also
serve as low vapor pressure replacements for the volatile organic solvents
typically employed in separations.35, 36 ILs have also shown potential to
enhance metal ion separations. 37, 38, 39, 40, 41, 42 For example, with CMPO,
the distribution ratio for uranyl extraction increases from $\sim$1 to
$\sim$20 (for an initial aqueous of 0.01 M HNO3) when the extractant phase
solvent is changed from dodecane to IL. 37, 43 The increase in this case is
due, at least in part, to the different extraction mechanisms enabled by the
IL, such as cation exchange using the IL cation, which makes direct comparison
to conventional solvents challenging. Here, we compare solvents in a manner
that is agnostic to the particular extraction mechanism.

In the following, we present MD simulations with a free energy estimation
technique, umbrella sampling, to show how, compared to typical solvents such
as dodecane, TBP, and, to a greater extent, IL, can significantly bias the
conformational free energy profile of the extractant, stabilizing the metal-
binding cis-CMPO conformation. To verify that this effect persists even when
the extractant binds with the metal ion, we also carry out IL-phase solvation
free-energy calculations in which an urnayl ion is annihilated while bound to
both cis-restrained and free CMPO. By comparing the two end states, restrained
cis and free CMPO, we are able to separate the stabilization effects
contributed by the metal-extractant complexation and the IL-phase solvent
effect, finding that the solvent effect still contributes substantially to
favoring the binding conformation. Overall, in contrast to traditional
research approaches focused on designing task-specific extractants by
tailoring their molecular structures and improving their functionalities, our
findings suggest how structure-directing solvents could be used to control
extractant conformational energetics, thereby suggesting a new route for
enhancing LLE separation efficiency.

##  2 Experiments

###  2.1 Materials

Caution! U-238 is an alpha-emitting isotope. All experiments described were
performed in specially designed laboratories with negative pressure fume hoods
and gloveboxes, using strict radiological controls.

CMPO was obtained from the Argonne National Laboratory stock. The
[EMIM][Tf_2N], n-dodecane and TBP were obtained from Sigma-Aldrich. NaNO3 and
HNO3 were obtained from Fisher Chemicals. UO2(NO3)2.6H2O was obtained from the
Argonne National Laboratory stock. All the chemicals were used as received.

###  2.2 FT-IR Spectra by Density Functional Theory (DFT) Calculations

In the DFT calculations, we use a complex consisting of one uranyl, with one
bidentate CMPO and two bidentate nitrate anions to fill the equatorial
coordination vacancies and neutralize the complex. The geometry of the uranyl-
CMPO complex is first optimized with the hybrid functional theory M06-2X 44 in
Gaussian 16 45; the Stuttgart/Dresden relativistic large core potential (RLC)
basis set (ECP), obtained from Basis Set Exchange 46, is used for uranium,
while the basis set 6-311++G(d,p) is used for other elements in the complex.
Followed by geometry optimization, single-point energy with frequency
calculation is deployed. The C=O stretching frequencies are summarized in
Table S1.

###  2.3 FT-IR Spectra

Samples for ATR-FTIR were prepared as follows. CMPO was dissolved in IL at
0.1M and contacted with an aqueous solution of 1M NaNO3 and 1mM HNO3 at equal
phase volumes. This “pre-equilibrated” organic phase was collected and used as
the extractant solution for an aqueous solution containing 20mM UO2(NO3)2, 1M
NaNO3 and 1mM HNO3. The volumes of the two phases for the extraction were
equal. The phases were contacted for 30 minutes and centrifuged. The separated
organic phase with uranyl ion and the uranyl-free, pre-equilibrated organic
phases were then analyzed by FT-IR.

ATR-FTIR spectra were collected with a Nicolet iS50 ATR spectrometer equipped
with a 2.8 mm round, type IIa diamond crystal, KBr beamsplitter, and a DTGS
detector. Spectra were collected with a resolution of 4 cm-1, a zero-filling
factor of 2, and 32 scans per spectrum in the wavenumber range of 4000 to 400
cm-1. Ambient background was subtracted from the sample spectra. The ATR
crystal was carefully cleaned with methanol between the samples.

###  2.4 MD Simulations

The General AMBER Force Field (GAFF2) 47, 48 is used to describe intra- and
inter-molecular interactions for CMPO, TBP, and [EMIM][Tf_2N]. Atomic partial
charges for these molecules are derived from Austin Model 1 with bond charge
corrections (AM1-BCC) method. 49, 50 The parameter set for dodecane is taken
from an optimized GAFF model to prevent freezing at room temperature. 51 The
SPC/E water model 52 is used to describe water. 53 A uranyl model by Kerisit
and Liu has shown structural, kinetic, and thermodynamic properties that agree
with experimental data, so we use this model for the uranyl ion in the
alchemical simulations.

Geometric mixing rules are used for all cross-term Lennard-Jones interactions.
An energy cutoff of 12 Å is set for nonbonded Lennard-Jones and electrostatic
interactions. Long-range electrostatic interactions are accounted for by a
particle-particle particle-mesh solver. 55 The leap-frog stochastic dynamics
integrator at 350 K along with a timestep of 1 femtosecond is used for the
thermostat 56; the Parrinello-Rahman barostat 57, 58 is applied at 1
atmosphere. All simulations are performed in the GROMACS MD software package.
59

####  2.4.1 Umbrella sampling

A cis-CMPO is packed with 100 solvent molecules using PACKMOL 60: 100
dodecane; 100 TBP; 100 [EMIM][Tf_2N] for dry IL simulation; 70 [EMIM][Tf_2N]
and 30 water moleucles for wet IL simulation. The umbrella sampling starts
with an oxygen-oxygen distance of 2.8 Å and reaches 5.1 Å, with an increment
of 0.1 Å for each window using a harominic bias potential of 100 kcal/mol⋅Å.
In each window, a 20 nanosecond simulation is performed, whereas the oxygen-
oxygen distance is sampled every 0.05 picosecond. All free energy profiles are
analyzed using the Weighted Histogram Analysis Method (WHAM). Histograms are
shown in Figures S2-S6. Sample GROMACS input scripts and force field files
(.itp) are provided at Zenodo: doi.org/10.5281/zenodo.8378117.

####  2.4.2 Alchemical free energy estimations

In order to isolate the impact of metal complexation on the solvent effect, we
conduct alchemical transformation calculations on a uranyl nitrate ion pair
bound to the CMPO molecule. A UO_2(NO_3)^+ cation complex is annihilated from
the IL phase to mimic the cation exchange mechanism in the uranyl extraction.
43 The annihilation of a charged particle is a well-known size-effect problem
in alchemical free energy calculation. 61 In order to minimize this size
effect, simulations are performed for one CMPO in a box with 300 ion pairs of
[EMIM][Tf_2N], which sizes up to roughly 50 Å on each dimension. Such a box
size is shown to minimize the size effect after application of a uniform
background charge, which we implement in the long-range electrostatics. 62 The
free energy difference between the start and the end states is measured using
the TI-CUBIC method by fitting the curve $\langle(\frac{\partial
U(\lambda)}{\partial\lambda})\rangle_{\lambda_{i}}$ vs $\lambda_{i}$ to a
natural cubic spline and then analytically integrating through the cubic
equation. 63 This methodology has been implemented in the Alchemical Analysis
code. 64 $\langle(\frac{\partial
U(\lambda)}{\partial\lambda})\rangle_{\lambda_{i}}$ and $\Delta G$ in each
window of perturbing $\lambda_{i}$ are shown in Figures S7-S11.

##  3 Results and Discussions

###  3.1 FT-IR

First, to demonstrate that CMPO binds uranyl in a bidentate manner with the
cis-conformation even in the relatively polar IL environment, we measure FT-IR
spectra after pre-equilibration with the metal-free aqueous phase and again
after contact with the uranyl-containing aqueous phase. Although P=O
stretching is in the region around 1200 cm-1, which is difficult to
distinguish from features associated with pure IL 65, EXFAS studies have
clearly shown P=O binds strongly with uranyl. 43 So, here we will only focus
on the C=O region: if C=O binds with uranyl, we will have strong evidence that
CMPO coordinates the uranyl in a bidentate manner. Figure 2 shows the region
of 1500–1700 cm-1, which includes the C=O stretching region of interest, with
the difference between the IR spectra with and without uranyl shown in the
right panel. (FTIR spectra for the entire wavenumber range are shown in Figure
S1.)

Although C=O stretching is also overlapped by the strong signal of IL
imidazolium C=C stretching 65, the difference in IR after uranyl extraction
can provide direct evidence of the interaction between the C=O group and the
uranyl cation. As shown in Figure 2, the C=O stretching mode shifts upon
uptake of uranyl into the organic phase, with $\Delta\nu_{\ce{C=O}}=0.42$
cm-1. This is consistent with the shift reported in the literature for the
extraction of CMPO of other metal types that is attributed to bindentate
complexation. 66 The peak shift in the C=O vibrational mode also corresponds
directly to our DFT-calculated peak shift upon uranyl binding (see Table S1).
To ensure that the magnitude of the shifted peaks is not small compared to the
amount of uranyl in the organic phase, we also show the change in IR spectra
after extraction of uranyl in the U=O region (900–1000 cm-1) in Figure 2. This
peak associated with the uranyl stretching mode is similar in intensity to the
shifted C=O peak. Although these intensities cannot be directly and
quantitatively compared, this still indicates that the approximate amount of
shifted C=O bonds is of the same order as the total amount of extracted uranyl
molecules. Overall, this suggests that CMPO forms bidentate-extracted
complexes with uranyl.

Figure 2: Left panel: IR spectra of C=O region of interest for 0.1 M CMPO in
[EMIM][Tf_2N] for the pre-equilibrated organic phase (solid line) and uranyl-
contacted organic phase (dashed line). Middle panel: subtracted IR spectra of
the uranyl-contacted organic phase by the pre-equibrated organic phase. Right
panel: subtracted IR spectra of U=O region of interest.

###  3.2 Free energy calculations and reaction coordinate

Previously, we applied metadynamics to sample the whole free energy landscape
of extractant conformation. 67 From there, we have seen significant free
energy penalties for the bound-state cis conformation in dodecane (up to 30
kJ/mol). This large difference in energies between conformations and high
energetic barriers between conformations necessitate advanced sampling
techniques, such as umbrella sampling, 68 to adequately sample the entire
conformational space. Compared with the metadynamics technique in our previous
study, umbrella sampling, with harmonic bias potentials, can better control
collective variables along a reaction coordinate. However, mapping this entire
2D Ramachandran plot using umbrella samling is not computationally feasible,
so we instead construct a 1D reaction coordinate that differentiates the three
molecular conformations of interest: the phosphine oxide to carbonyl oxygen
distance. This distance is a meaningful reaction coordinate because it readily
distinguishes cis, gauche and trans conformations. By this construction of the
reaction coordinate, we can reduce the 2D Ramachandran plot to a 1D map. By
performing umbrella sampling along this reaction coordinate, we can calculate
the free energy differences between each conformation from MD simulations
using explicitly modeled solvent.

###  3.3 Gas phase, dodecane and TBP solvents

The free energy profile in the gas phase provides information on the intrinsic
energetics of the molecule in the absence of any solvent interactions. This
free energy profile along the reaction coordinate is shown in Figure 3 (purple
line). The cis-CMPO conformation corresponds to the left end of each profile
at 2.8 Å. The free energy levels off around 3.8 Å, corresponding to the gauche
conformation, while the global free energy minimum at the trans conformation
is around 4.4 Å. At higher distances, the free energy continues to increase as
the molecule is distorted. Since the umbrella sampling provides the relative
free energy rather than the absolute one, we anchor all profiles at $\Delta
G_{trans}$ = 0 kJ/mol for cross-comparisons. Compared to the gas phase, the
dodecane solvent phase (green line) has a similar profile, with a primary
difference of the dodecane profile being slightly lower energies of the gauche
and cis regions compared to the gas phase. Overall, both in the gas phase and
in the dodecane, CMPO shows a significant free energy penalty for the cis
conformation relative to trans (roughly ΔG _cis-trans=30 kJ/mol). Such a large
energy penalty works against the formation of metal-extractant complexes where
the free, unbound trans-CMPO rearranges into the cis conformation to bind the
metal.

Figure 3: Free energy profiles in each solvent are plotted along the oxygen-
oxygen distance reaction coordinate, with distances corresponding to each of
the CMPO conformations highlighted.

While dodecane might be expected to act similarly to the gas phase, TBP is a
partially polar solvent and may affect the conformational energetics of CMPO.
A snapshot of CMPO in TBP is shown in Figure 4a. For TBP, in Figure 3 (blue
line) we find that the free energy profile has been comparatively reduced,
with ΔG _cis-trans and ΔG _gauche-trans decreasing by 8 kJ/mol and 10 kJ/mol,
respectively. The difference between trans and cis states is roughly twice as
large between TBP and dodecane than dodecane and the gas phase. A shallow
local minimum is also observed around 3.4 Å, indicating that the gauche
conformation is stabilized by the TBP solvent. In contrast to dodecane, TBP is
partially polar and, as we find, can partially stabilize cis and gauche
conformations through dipole alignment with the CMPO head group dipole; such
solvent arrangement is shown in Figure 4b. To quantify this effect, we
calculate spatial distribution functions (SDFs) to show the ensemble average
alignment of TBP with CMPO for each conformation. SDFs, often used to
understand solvation in molecular simulations, show the spatial distribution
of a particular atom or molecule of interest surrounding a reference atom or
molecule. Figure 4c-e show the SDFs of cis, gauche and trans CMPO,
respectively. For the cis and gauche conformations, it is clear that TBP’s
$sp^{2}$ oxygen, phosphorus, and $sp^{3}$ oxygen have structured distributions
above and below the CMPO head group, indicating solvent TBP molecules align
with the CMPO dipole. The trans-CMPO has a weaker intrinsic dipole because the
polar O=P and O=C bonds point in opposite directions. This prevents solvent
organization around the CMPO, as evidenced by the relative lack of structure
in the SDF shown in Figure 4e. Overall, these results show how organization of
the partially polar TBP solvent molecules around the CMPO stabilizes the more
polar cis, gauche-CMPO conformations.

Figure 4: CMPO in TBP. a) shows the overall simulation box; b) shows a zoom-in
inset of TBP alignment with cis-CMPO, and the dashed arrow represents the
overall dipole orientation of the TBP molecule; c)-e) show the spatial
distribution functions (SDF) of TBP double-bonded $sp^{2}$ oxygen (red),
phosphorus (orange) and single-bonded $sp^{3}$ oxygen (yellow) centered by
cis, gauche and trans CMPO respectively.

###  3.4 IL solvents

Next, we consider IL solvents, which are structurally unique compared to
traditional molecular solvents. In this case, the right panel in Figure 3
shows the cis and gauche free energies are dramatically reduced. Instead of
the global minimum at the trans, the gauche becomes the most stable
conformation. The ΔG _gauche-trans is –5 kJ/mol—a reduction of 8 kJ/mol
compared to dodecane; the ΔG _cis-trans is lowered from 30 kJ/mol in dodecane
to 7 kJ/mol in the dry IL. We attribute this significant stabilization of the
cis and gauche conformations to the alignment of the CMPO dipole with the
cation/anion nanodomains characteristic of ILs. 31, 34 Figure 5a-b shows a
representative snapshot in which CMPO is surrounded by cationic/anionic
domains. More importantly, the SDFs in Figure 5 show that, for both O=P and
O=C groups, strong negative partial charges of the $sp^{2}$ oxygen atoms lead
to cationic [EMIM] domains surrounding those oxygen atoms. Conversely, the
other partially positive end of the O=P and O=C dipoles associate with the
anionic [Tf_2N] domains. The cationic and anionic layers therefore induce
dipole alignment of the embedded extractant molecule, favoring the cis and
gauche conformations. Considering that multidentate extractants, such as CMPO,
may need to overcome a high energy barrier to reach their metal-binding
conformation, such significant stabilization induced by the solvent could
potentially improve separations performance. This observation implies new
opportunities for further control of separation energetics, such as the
customization of nanodomains through IL design.

Figure 5: CMPO in dry [EMIM][Tf_2N]. a) shows the cationic (blue) and anionic
(red) domains in the simulation box; b) shows a zoom-in inset of cis-CMPO
surrounded by cationic/anionic domains; c)-e) show the spatial distribution
functions (SDF) of cationic EMIM (blue) and anionic Tf_2N (red) centered by
cis, gauche and trans CMPO respectively.

The solubility of water in ILs can be significant and, in LLE applications
where the IL is in direct contact with an aqueous phase, nonnegligible amounts
of water transfer to the IL phase. For example, the saturated concentration of
water in [EMIM][Tf_2N] is 0.2982 mole fraction, or 2.89% volume fraction
(vol%). 69 Therefore, to account for the realistic solution environment
encountered in LLE, we also consider the water-saturated IL. Despite hydrogen
bonding between water and CMPO’s polar O=P and O=C sites, the relative free
energies of the CMPO conformations are nearly the same as the dry IL (Figure
3). A zoom-in inset in Figure 6b shows two water molecules are forming
hydrogen bonds with the $sp^{2}$ oxygen atoms on CMPO. SDFs in Figure 6c-e
further illustrate such molecular events: we can see that water distributes
(cyan) around the two $sp^{2}$ oxygen atoms for all three cis, gauche and
trans CMPO conformations. This finding suggests that any impact of water
hydrogen bonding favoring cis or gauche conformations in the CMPO is either
not significant or roughly cancelled out by the water interfering with the
dipole alignment induced by the IL nanodomains.

Figure 6: CMPO in wet [EMIM][Tf_2N]. a) shows the cationic (blue), anionic
(red), and water (cyan) domains in the simulation box; b) shows a zoom-in
inset of cis-CMPO surrounded by cationic/anionic domains while forming
hydrogen bond with water molecules; c)-e) show the spatial distribution
functions (SDF) of cationic EMIM (blue), anionic Tf_2N (red) and water (cyan)
centered by cis, gauche and trans CMPO respectively.

###  3.5 Alchemical free energy calculations with uranyl ion

The above analysis considers the impact of the solvent in the absence of metal
ions. In an actual extraction process, the metal ion will participate in the
coordination with the extractant molecule. Such participation of the metal ion
will potentially diminish the solvent preorganization effect onto the
extrantant energetics, especially when one considers that part of the
extractant will be binding to the metal ion, potentially suppressing
interactions with the surrounding solvent. To directly test the solvent effect
with the presence of the metal ion, we perform a sequence of alchemical
transformation calculations. By calculating the solvation free energy of a
uranyl nitrate bound to CMPO with and without a bias potential on the
extractant conformation, we can ascertain the free energy cost of constraining
the extractant to the binding conformation and compare it to the same cost in
the absence of the metal. The difference, then, reveals how much the cation
has diminished the effect of the solvent. By comparing this to the free energy
profiles above obtained by umbrella sampling, we can identify which fraction
of the solvent reorganization free energy resulting from the solvent persists
in the presence of the metal. A schematic flow chart of how alchemical
transformations are carried out in this work is shown in Figure 7. Details
associated with each alchemical transformation are described below.

Figure 7: Schematic flow chart of alchemical transformations calculated. Step
(a1): Annihilating UO_2(NO_3)^+ with a restraint placed to keep CMPO stay at
cis-conformation. Step (a2): Annihilating UO_2(NO_3)^+ without any restraint
on CMPO. Step (b): Annihilating the restraint on CMPO to let it relax from
cis-conformation to free.

  * (a1)

The annihilation UO_2(NO_3)^+ with a cis-CMPO: One UO_2^2+ and one NO_3^- are
contained in the simulation box. In the initial configuration of the metal
complexant UO_2[NO_3][CMPO], CMPO and NO_3^- chelate with uranyl in a
bidentate fashion. The UO_2^2+ and NO_3^- are decoupled from the simulation
box by mutating $\lambda_{vdw}$ and $\lambda_{coul}$ in the soft-core version
of the 12-6-1 potential. $\lambda_{coul}$ is decoupled first followed by
$\lambda_{vdw}$, in order to prevent the “end-point catastrophes”. 12-point
scheme is chosen for $\lambda_{coul}$ (0.00922, 0.04794, 0.11505, 0.20634,
0.31608, 0.43738, 0.56262, 0.68392, 0.79366, 0.88495, 0.95206, 0.99078) to
minimize quadrature errors. For $\lambda_{vdw}$, an increment of 0.1 is chosen
between each window. To keep CMPO at cis-comformation, a harmonic restraint of
100 kcal/mol⋅Å is placed between phosphine oxide and carbonyl oxygen to
maintain the distance at 2.8 Å.

  * (a2)

The annihilation of UO_2(NO_3)^+ with a free CMPO: This alchemical
transformation is the same as in the previous step, except for the restraint
on the CMPO. We do not apply any restraint on the CMPO molecules, thus, it can
freely move in solution, while the UO_2^2+ and NO_3^- are phased out.

  * (b)

The annihilation of the restraint on cis-CMPO: In this step, in order to test
how much energy can be recovered from relaxing the CMPO conformation, we phase
out the restraint from the end state in step 1. The restraint is decoupled by
linearly scaling down the restraint to 0 kcal/mol⋅Å with 21 windows.

From the alchemical transformations outlined above, we define the relevant
free energies given in Equations 1 and 2.

$\Delta\Delta G^{\text{total}}_{cis-\text{free}}=\Delta G_{\text{(a1)}}-\Delta
G_{\text{(a2)}}$ (1)

$\Delta\Delta G_{\text{UO}_{2}}=\Delta\Delta
G^{\text{total}}_{cis-\text{free}}-\Delta\Delta G_{\text{(b)}}$ (2)

$\Delta\Delta G^{\text{total}}_{cis-\text{free}}$ is the free energy
difference of solvating the uranyl between the cis-constrained and
unconstrained CMPO molecule. $\Delta\Delta G_{\text{UO}_{2}}$ is the
contribution to $\Delta\Delta G^{\text{total}}_{cis-\text{free}}$ that results
from the presence of the uranyl, i.e., the contribution that does not come
from just the solvent in the absence of the uranyl. ΔG_(a1), ΔG_(a2), and
ΔΔG_(b) are from the alchemical free energy calculations in steps (a1), (a2),
and (b) as we describe above. We visually show the relationship of these free
energy terms in Figure 8. From the same initial configuration, steps (a1)
(annihilation of uranyl with CMPO restrained) and (a2) (annihilation of uranyl
with CMPO free) result in different end states — cis and free CMPO in IL,
respectively. These two alchemical transformations have also resulted in a
notable difference in free energies, with a ΔΔG^total_cis-free of 20.92
kJ/mol. Noting that annihilation of the uranyl nitrate with CMPO restrained to
the cis conformation requires a higher energy (ΔG_(a2) $<$ ΔG_(a1)), this
observation indicates that the uranyl ion is better solvated by cis-CMPO. To
prove that these numbers are not affected by the error when the curves of
$\langle(\frac{\partial U(\lambda)}{\partial\lambda})\rangle_{\lambda_{i}}$
vs. $\lambda_{i}$ are fitted to cubic splines, in addition to TI-CUBIC, we
have carried out extra free energy estimations using BAR 70 and MBAR 71. The
comparisons have been shown in Figure S8, S10; differences between these
methodologies are minimal.

The difference in free energy between two end states, ΔΔG^total_cis-free,
includes the stabilization effects of the uranyl metal ion and the IL
cationic/anionic nanodomain. Although these two end states do not contain the
uranyl metal ion, its influence comes from the steps (a1) and (a2) during the
annihilation of metal ion. By then phasing out the restraint that keeps the
CMPO at the cis conformation, 10.70 kJ/mol of free energy (ΔΔG_(b)) is
recovered in step (b). It is worth mentioning that, in this case, TI-CUBIC
could give sizable uncertainty when fitting the curve of
$\langle(\frac{\partial U(\lambda)}{\partial\lambda})\rangle_{\lambda_{i}}$
vs. $\lambda_{i}$ using a cubic spline, which is shown in Figure S11. This
uncertainty becomes more relevant in the last window, where an abrupt decrease
of $\langle(\frac{\partial U(\lambda)}{\partial\lambda})\rangle_{\lambda_{i}}$
vs. $\lambda_{i}$ occurs. As a result, fitting $\langle(\frac{\partial
U(\lambda)}{\partial\lambda})\rangle_{\lambda_{i}}$ in this region using a
cubic spline could potentially lead to error in the estimated free energy.
Therefore, to avoid this issue, we use MBAR, which has shown to provide
minimal variance by making use of equilibrium data collected from multiple
states, 71 to estimate the free energy for this step, finding a value of
ΔΔG_(b) = 10.70 kJ/mol. To show the reliability of MBAR over TI-CUBIC, we also
show the overlap matrix in Figure 12S, which has satisfied overlapping for the
last several $\lambda_{i}$s. This ΔΔG_(b) from the alchemical calculations is
slightly smaller than ΔG_cis-gauche in Figure 3 from umbrella sampling (showed
by the orange line), which is presumably due in part to how, when the bias
potential is removed from the cis-CMPO, we are sampling an ensemble average of
all lower energetic configurations for free the CMPO, rather than the lowest
gauche configuration in Figure 3. Although alchemical and umbrella sampling
are two types of free energy methodologies based on different theories, the
observation that ΔΔG_(b) and ΔG_cis-gauche agree with each other in a
reasonable way provides extra confidence in both sets of free energy
calculations.

To verify that ΔΔG_(b) comes from the solvent effects, we also perform step
(b) in the gas phase. With the change in solvent background, ΔΔG_(b) increases
to 27.68 kJ/mol, again agreeing the umbrella sampling result for ΔG_cis-trans
in Figure 3 (shown by the purple line). After cross-checking the values for
ΔΔG_(b) from the solvent effect, we can subtract this contribution from
ΔΔG^total_cis-free, resulting in 10.22 kJ/mol (ΔΔG_UO2) from stabilization of
the cis conformation through metal-extractant binding. From breaking down each
contribution and comparing ΔΔG^total_cis-free, ΔΔG_(b) and ΔΔG_UO2, we
therefore conclude that roughly half of the stabilization of the cis
conformation from the IL solvent is retained in the presence of the extracted
metal ion. If the stabilization imparted by the anion domain is unaffected by
the replacement of the IL cation domain with the metal cation upon binding,
perhaps approximately half of the total energetic stabilization of the cis
conformation would be expected.

Figure 8: Free energies of each step in the alchemcial transformation
calculations. Steps (a1) and (a2) involving annihilating the uranyl ion have
been colored by blues, while step (b) involving the extractant pre-
organization without the uranyl ion is colored by red.

##  4 Conclusions

Reversible LLE processes are typically driven by small free energy differences
in the solubility of the targeted solute between immiscible phases; however,
individual steps in the metal complexation and extraction process are often
characterized by relatively larger energetic contributions that mostly cancel
out. These include the electrostatic interactions between extractant and
metal, as well as the assembly of multiple extractants into the final
extracted complex. In this study, we find that the choice of organic solvent
induces substantial changes in the energetic cost associated with organizing
the extractant into the binding conformation. For example, the roughly 20
kJ/mol free energy difference between binding and non-binding conformations we
find between dodecane and IL is large compared to the typical total extraction
free energies. Then, by directly introducing the metal ion in the simulations
and calculating the free energy solvation from the alchemical transformation
of annihilating the metal ion with both free and cis-restrained CMPO, we are
able to separate stabilization effects of the metal from those of the IL
solvent. These alchemical free energy calculations show that roughly half of
the extractant reorganizational free energy results from just the solvent
effects measured with umbrella sampling even in the presence of the metal. The
magnitude of this effect, with a calculated value of 10.79 kJ/mol, is
sufficiently large to have a substantial impact on practical separations
performance.

Our molecular dynamics simulations reveal that the mechanism for stabilizing
the binding conformation in the IL is alignment with the anion/cation
nanodomains, which create layers of positive and negative charge that promote
dipole alignment of the extractant molecule, thereby stabilizing the more
strongly dipolar cis conformation. Overall, these results suggest that
considering solvent effects on extractant conformation directly and explicitly
is essential to broadly understand and manipulate the delicate balance of
extraction energetics. Such a mechanism by which the solvent induces a
significant change in the extractant’s conformational free energy landscape
indicates broad possibilities for tailoring solvents to manipulate separation
energetics.

This work was supported by U.S. Department of Energy (DOE), Office of Science,
Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and
Biosciences Division, Separation Science Program, under Contract DE-
AC02-06CH11357 to UChicago Argonne, LLC, Operator of Argonne National
Laboratory. We gratefully acknowledge the computing resources provided on
Bebop, a high-performance computing cluster operated by the Laboratory
Computing Resource Center at Argonne National Laboratory. We thank Dr. Peter
Tkac and Dr. M. Alex Brown for generously providing CMPO.

IR spectra experiments; histograms from umbrella sampling along the reaction
coordinate in dodecane, TBP, dry and wet IL; plots of $\langle(\frac{\partial
U(\lambda)}{\partial\lambda})\rangle_{\lambda_{i}}$ and $\Delta G$ vs
$\lambda$ in the alchemical free energy calculations.

## References

  * Sholl and Lively 2016 Sholl, D. S.; Lively, R. P. Seven chemical separations to change the world. _Nature_ 2016, _532_ , 435–437
  * Sun et al. 2012 Sun, X.; Luo, H.; Dai, S. Ionic liquids-based extraction: a promising strategy for the advanced nuclear fuel cycle. _Chemical Reviews_ 2012, _112_ , 2100–2128
  * De Jesus et al. 2021 De Jesus, K.; Rodriguez, R.; Baek, D.; Fox, R.; Pashikanti, S.; Sharma, K. Extraction of lanthanides and actinides present in spent nuclear fuel and in electronic waste. _Journal of Molecular Liquids_ 2021, _336_ , 116006
  * Lumetta et al. 2002 Lumetta, G. J.; Rapko, B. M.; Garza, P. A.; Hay, B. P.; Gilbertson, R. D.; Weakley, T. J.; Hutchison, J. E. Deliberate design of ligand architecture yields dramatic enhancement of metal ion affinity. _Journal of the American Chemical Society_ 2002, _124_ , 5644–5645
  * Špadina et al. 2021 Špadina, M.; Dufrêche, J.-F.; Pellet-Rostaing, S.; Marčelja, S.; Zemb, T. Molecular forces in liquid–liquid extraction. _Langmuir_ 2021, _37_ , 10637–10656
  * Lewis et al. 2011 Lewis, F. W.; Harwood, L. M.; Hudson, M. J.; Drew, M. G.; Desreux, J. F.; Vidick, G.; Bouslimani, N.; Modolo, G.; Wilden, A.; Sypula, M., et al. Highly efficient separation of actinides from lanthanides by a phenanthroline-derived bis-triazine ligand. _Journal of the American Chemical Society_ 2011, _133_ , 13093–13102
  * Ansari et al. 2012 Ansari, S. A.; Pathak, P.; Mohapatra, P. K.; Manchanda, V. K. Chemistry of diglycolamides: promising extractants for actinide partitioning. _Chemical Reviews_ 2012, _112_ , 1751–1772
  * McCann et al. 2016 McCann, B. W.; Silva, N. D.; Windus, T. L.; Gordon, M. S.; Moyer, B. A.; Bryantsev, V. S.; Hay, B. P. Computer-aided molecular design of bis-phosphine oxide lanthanide extractants. _Inorganic Chemistry_ 2016, _55_ , 5787–5803
  * Ivanov and Bryantsev 2016 Ivanov, A. S.; Bryantsev, V. S. A Computational Approach to Predicting Ligand Selectivity for the Size-Based Separation of Trivalent Lanthanides. _European Journal of Inorganic Chemistry_ 2016, _2016_ , 3474–3479
  * Brigham et al. 2017 Brigham, D. M.; Ivanov, A. S.; Moyer, B. A.; Delmau, L. H.; Bryantsev, V. S.; Ellis, R. J. Trefoil-shaped outer-sphere ion clusters mediate lanthanide (III) ion transport with diglycolamide ligands. _Journal of the American Chemical Society_ 2017, _139_ , 17350–17358
  * Morgenstern et al. 2020 Morgenstern, A.; Lilley, L. M.; Stein, B. W.; Kozimor, S. A.; Batista, E. R.; Yang, P. Computer-assisted design of macrocyclic chelators for actinium-225 radiotherapeutics. _Inorganic Chemistry_ 2020, _60_ , 623–632
  * Bessen et al. 2021 Bessen, N. P.; Popov, I. A.; Heathman, C. R.; Grimes, T. S.; Zalupski, P. R.; Moreau, L. M.; Smith, K. F.; Booth, C. H.; Abergel, R. J.; Batista, E. R., et al. Complexation of lanthanides and heavy actinides with aqueous sulfur-donating ligands. _Inorganic Chemistry_ 2021, _60_ , 6125–6134
  * Liu et al. 2022 Liu, T.; Johnson, K. R.; Jansone-Popova, S.; Jiang, D.-e. Advancing Rare-Earth Separation by Machine Learning. _JACS Au_ 2022, _2_ , 1428–1434
  * Driscoll et al. 2022 Driscoll, D. M.; Liu, H.; Reinhart, B.; Popovs, I.; Bocharova, V.; Jansone-Popova, S.; Jiang, D.-e.; Ivanov, A. S. Noncoordinating Secondary Sphere Ion Modulates Supramolecular Clustering of Lanthanides. _The Journal of Physical Chemistry Letters_ 2022, _13_ , 12076–12081
  * Zhang et al. 2022 Zhang, X.; Adelman, S. L.; Arko, B. T.; De Silva, C. R.; Su, J.; Kozimor, S. A.; Mocko, V.; Shafer, J. C.; Stein, B. W.; Schreckenbach, G., et al. Advancing the Am Extractant Design through the Interplay among Planarity, Preorganization, and Substitution Effects. _Inorganic Chemistry_ 2022, _61_ , 11556–11570
  * Yamamuro et al. 2011 Yamamuro, O.; Yamada, T.; Kofu, M.; Nakakoshi, M.; Nagao, M. Hierarchical structure and dynamics of an ionic liquid 1-octyl-3-methylimidazolium chloride. _The Journal of Chemical Physics_ 2011, _135_ , 054508
  * Baldwin et al. 2018 Baldwin, A. G.; Ivanov, A. S.; Williams, N. J.; Ellis, R. J.; Moyer, B. A.; Bryantsev, V. S.; Shafer, J. C. Outer-sphere water clusters tune the lanthanide selectivity of diglycolamides. _ACS central science_ 2018, _4_ , 739–747
  * Servis et al. 2020 Servis, M. J.; Martinez-Baez, E.; Clark, A. E. Hierarchical phenomena in multicomponent liquids: simulation methods, analysis, chemistry. _Physical Chemistry Chemical Physics_ 2020, _22_ , 9850–9874
  * Wang et al. 2020 Wang, Y.-L.; Li, B.; Sarman, S.; Mocci, F.; Lu, Z.-Y.; Yuan, J.; Laaksonen, A.; Fayer, M. D. Microstructural and dynamical heterogeneities in ionic liquids. _Chemical Reviews_ 2020, _120_ , 5798–5877
  * Sprakel and Schuur 2019 Sprakel, L. M.; Schuur, B. Improving understanding of solvent effects on intermolecular interactions in reactive liquid–liquid extraction with Isothermal Titration Calorimetry and molecular modeling. _Journal of Industrial and Engineering Chemistry_ 2019, _72_ , 364–373
  * Guilbaud and Zemb 2015 Guilbaud, P.; Zemb, T. Depletion of water-in-oil aggregates from poor solvents: Transition from weak aggregates towards reverse micelles. _Current Opinion in Colloid & Interface Science_ 2015, _20_ , 71–77
  * Li et al. 2019 Li, Z.; Zhang, Z.; Smolders, S.; Li, X.; Raiguel, S.; Nies, E.; De Vos, D. E.; Binnemans, K. Enhancing metal separations by liquid–liquid extraction using polar solvents. _Chemistry–A European Journal_ 2019, _25_ , 9197–9201
  * Hay et al. 2001 Hay, B. P.; Dixon, D. A.; Vargas, R.; Garza, J.; Raymond, K. N. Structural Criteria for the Rational Design of Selective Ligands. 3. Quantitative Structure- Stability Relationship for Iron (III) Complexation by Tris-Catecholamide Siderophores. _Inorganic Chemistry_ 2001, _40_ , 3922–3935
  * Boehme and Wipff 2002 Boehme, C.; Wipff, G. Carbamoylphosphine oxide complexes of trivalent lanthanide cations: role of counterions, ligand binding mode, and protonation investigated by quantum mechanical calculations. _Inorganic Chemistry_ 2002, _41_ , 727–737
  * Ellis et al. 2017 Ellis, R. J.; Brigham, D. M.; Delmau, L.; Ivanov, A. S.; Williams, N. J.; Vo, M. N.; Reinhart, B.; Moyer, B. A.; Bryantsev, V. S. “Straining” to separate the rare earths: how the lanthanide contraction impacts chelation by diglycolamide ligands. _Inorganic Chemistry_ 2017, _56_ , 1152–1160
  * Kapoor et al. 2005 Kapoor, R. N.; Guillory, P.; Schulte, L.; Cervantes-Lee, F.; Haiduc, I.; Parkanyi, L.; Pannell, K. H. Di (p-tert-butylphenyl)-N, N-di-(iso-butyl) carbamoylmethylphosphine oxide and its organotin and uranyl adducts: structural and spectroscopic characterization. _Applied Organometallic Chemistry_ 2005, _19_ , 510–517
  * Chaumont and Wipff 2004 Chaumont, A.; Wipff, G. Solvation of uranyl (II) and europium (III) cations and their chloro complexes in a room-temperature ionic liquid. a theoretical study of the effect of solvent “humidity”. _Inorganic Chemistry_ 2004, _43_ , 5891–5901
  * Chaumont and Wipff 2006 Chaumont, A.; Wipff, G. Solvation of uranyl–CMPO complexes in dry vs. humid forms of the [BMI][PF_6] ionic liquid. A molecular dynamics study. _Physical Chemistry Chemical Physics_ 2006, _8_ , 494–502
  * Ellis et al. 2012 Ellis, R. J.; Audras, M.; Antonio, M. R. Mesoscopic aspects of phase transitions in a solvent extraction system. _Langmuir_ 2012, _28_ , 15498–15504
  * Canongia Lopes and Pádua 2006 Canongia Lopes, J. N.; Pádua, A. A. Nanostructural organization in ionic liquids. _The Journal of Physical Chemistry B_ 2006, _110_ , 3330–3335
  * Hayes et al. 2015 Hayes, R.; Warr, G. G.; Atkin, R. Structure and nanostructure in ionic liquids. _Chemical Reviews_ 2015, _115_ , 6357–6426
  * Yan et al. 2019 Yan, X.; Anguille, S.; Bendahan, M.; Moulin, P. Ionic liquids combined with membrane separation processes: A review. _Separation and Purification Technology_ 2019, _222_ , 230–253
  * Lyu et al. 2020 Lyu, Y.; Brennecke, J. F.; Stadtherr, M. A. Review of recent aromatic–aliphatic–ionic liquid ternary liquid–liquid equilibria and their modeling by COSMO-RS. _Industrial & Engineering Chemistry Research_ 2020, _59_ , 8871–8893
  * Kashyap et al. 2013 Kashyap, H. K.; Santos, C. S.; Murthy, N. S.; Hettige, J. J.; Kerr, K.; Ramati, S.; Gwon, J.; Gohdo, M.; Lall-Ramnarine, S. I.; Wishart, J. F.; Margulis, C. J.; Castner, E. W. J. Structure of 1-alkyl-1-methylpyrrolidinium bis (trifluoromethylsulfonyl) amide ionic liquids with linear, branched, and cyclic alkyl groups. _The Journal of Physical Chemistry B_ 2013, _117_ , 15328–15337
  * Dietz and Stepinski 2008 Dietz, M. L.; Stepinski, D. C. Anion concentration-dependent partitioning mechanism in the extraction of uranium into room-temperature ionic liquids. _Talanta_ 2008, _75_ , 598–603
  * Billard 2013 Billard, I. _Handbook on the Physics and Chemistry of Rare Earths_ ; Elsevier, 2013; Vol. 43; pp 213–273
  * Visser and Rogers 2003 Visser, A. E.; Rogers, R. D. Room-temperature ionic liquids: new solvents for f-element separations and associated solution chemistry. _Journal of Solid State Chemistry_ 2003, _171_ , 109–113
  * Dai et al. 1999 Dai, S.; Ju, Y.; Barnes, C. Solvent extraction of strontium nitrate by a crown ether using room-temperature ionic liquids. _Journal of the Chemical Society, Dalton Transactions_ 1999, 1201–1202
  * Kubota et al. 2011 Kubota, F.; Shimobori, Y.; Baba, Y.; Koyanagi, Y.; Shimojo, K.; Kamiya, N.; Goto, M. Application of ionic liquids to extraction separation of rare earth metals with an effective diglycol amic acid extractant. _Journal of Chemical Engineering of Japan_ 2011, _44_ , 307–312
  * Rout et al. 2009 Rout, A.; Venkatesan, K.; Srinivasan, T.; Rao, P. V. Extraction of americium (III) from nitric acid medium by CMPO-TBP extractants in ionic liquid diluent. _Radiochimica Acta_ 2009, _97_ , 719–725
  * Shen et al. 2011 Shen, Y.; Wang, S.; Zhu, L.; Wang, J.; Wu, W. Extraction of Th (IV) from an HNO3 solution by diglycolamide in ionic liquids. _Industrial & Engineering Chemistry Research_ 2011, _50_ , 13990–13996
  * Turanov et al. 2020 Turanov, A.; Karandashev, V.; Boltoeva, M. Solvent extraction of intra-lanthanides using a mixture of TBP and TODGA in ionic liquid. _Hydrometallurgy_ 2020, _195_ , 105367
  * Visser et al. 2003 Visser, A. E.; Jensen, M. P.; Laszak, I.; Nash, K. L.; Choppin, G. R.; Rogers, R. D. Uranyl coordination environment in hydrophobic ionic liquids: an in situ investigation. _Inorganic Chemistry_ 2003, _42_ , 2197–2199
  * Zhao and Truhlar 2008 Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. _Theoretical Chemistry Accounts_ 2008, _120_ , 215–241
  * Frisch et al. 2016 Frisch, M. J. et al. Gaussian˜16 Revision C.01. 2016; Gaussian Inc. Wallingford CT
  * Pritchard et al. 2019 Pritchard, B. P.; Altarawy, D.; Didier, B.; Gibson, T. D.; Windus, T. L. New basis set exchange: An open, up-to-date resource for the molecular sciences community. _Journal of Chemical Information and Modeling_ 2019, _59_ , 4814–4820
  * Wang et al. 2004 Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. _Journal of Computational Chemistry_ 2004, _25_ , 1157–1174
  * He et al. 2020 He, X.; Man, V. H.; Yang, W.; Lee, T.-S.; Wang, J. A fast and high-quality charge model for the next generation general AMBER force field. _The Journal of Chemical Physics_ 2020, _153_
  * Jakalian et al. 2000 Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. _Journal of Computational Chemistry_ 2000, _21_ , 132–146
  * Jakalian et al. 2002 Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. _Journal of Computational Chemistry_ 2002, _23_ , 1623–1641
  * Vo et al. 2015 Vo, Q. N.; Hawkins, C. A.; Dang, L. X.; Nilsson, M.; Nguyen, H. D. Computational study of molecular structure and self-association of tri-n-butyl phosphates in n-dodecane. _The Journal of Physical Chemistry B_ 2015, _119_ , 1588–1597
  * Berendsen et al. 1987 Berendsen, H.; Grigera, J.; Straatsma, T. The missing term in effective pair potentials. _Journal of Physical Chemistry_ 1987, _91_ , 6269–6271
  * Izadi and Onufriev 2016 Izadi, S.; Onufriev, A. V. Accuracy limit of rigid 3-point water models. _The Journal of Chemical Physics_ 2016, _145_ , 074501
  * Kerisit and Liu 2013 Kerisit, S.; Liu, C. Structure, kinetics, and thermodynamics of the aqueous uranyl (VI) cation. _The Journal of Physical Chemistry A_ 2013, _117_ , 6421–6432
  * Hockney and Eastwood 1989 Hockney, R. W.; Eastwood, J. W. _Computer simulation using particles_ ; Adam Hilger: New York, 1989
  * Goga et al. 2012 Goga, N.; Rzepiela, A.; De Vries, A.; Marrink, S.; Berendsen, H. Efficient algorithms for Langevin and DPD dynamics. _Journal of Chemical Theory and Computation_ 2012, _8_ , 3637–3649
  * Parrinello and Rahman 1981 Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. _Journal of Applied Physics_ 1981, _52_ , 7182–7190
  * Nosé and Klein 1983 Nosé, S.; Klein, M. Constant pressure molecular dynamics for molecular systems. _Molecular Physics_ 1983, _50_ , 1055–1076
  * Abraham et al. 2015 Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. _SoftwareX_ 2015, _1_ , 19–25
  * Martínez et al. 2009 Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A package for building initial configurations for molecular dynamics simulations. _Journal of Computational Chemistry_ 2009, _30_ , 2157–2164
  * Rocklin et al. 2013 Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hünenberger, P. H. Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects. _The Journal of Chemical Physics_ 2013, _139_
  * Bogusz et al. 1998 Bogusz, S.; Cheatham III, T. E.; Brooks, B. R. Removal of pressure and free energy artifacts in charged periodic systems via net charge corrections to the Ewald potential. _The Journal of Chemical Physics_ 1998, _108_ , 7070–7084
  * Bruckner and Boresch 2011 Bruckner, S.; Boresch, S. Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration. _Journal of Computational Chemistry_ 2011, _32_ , 1320–1333
  * Klimovich et al. 2015 Klimovich, P. V.; Shirts, M. R.; Mobley, D. L. Guidelines for the analysis of free energy calculations. _Journal of Computer-Aided Molecular Design_ 2015, _29_ , 397–411
  * Hofft et al. 2008 Hofft, O.; Bahr, S.; Kempter, V. Investigations with infrared spectroscopy on films of the ionic liquid [EMIM] Tf2N. _Langmuir_ 2008, _24_ , 11562–11566
  * Tkac et al. 2012 Tkac, P.; Vandegrift, G. F.; Lumetta, G. J.; Gelis, A. V. Study of the interaction between HDEHP and CMPO and its effect on the extraction of selected lanthanides. _Industrial & Engineering Chemistry Research_ 2012, _51_ , 10433–10444
  * Wang and Servis 2023 Wang, X.; Servis, M. J. Using Metadynamics to Reveal Extractant Conformational Free Energy Landscapes. 2023, arXiv:2309.06400 [physics.chem-ph]
  * Torrie and Valleau 1977 Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. _Journal of Computational Physics_ 1977, _23_ , 187–199
  * Freire et al. 2008 Freire, M. G.; Carvalho, P. J.; Gardas, R. L.; Marrucho, I. M.; Santos, L. M.; Coutinho, J. A. Mutual solubilities of water and the [C_n mim][Tf_2N] hydrophobic ionic liquids. _The Journal of Physical Chemistry B_ 2008, _112_ , 1604–1610
  * Bennett 1976 Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo data. _Journal of Computational Physics_ 1976, _22_ , 245–268
  * Shirts and Chodera 2008 Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. _The Journal of chemical physics_ 2008, _129_


# The layer-wise L1 Loss Landscape of Neural Nets is more complex around local
minima

Peter Hinz  ETH Zürich

###### Abstract

For fixed training data and network parameters in the other layers the L1 loss
of a ReLU neural network as a function of the first layer’s parameters is a
piece-wise affine function. We use the Deep ReLU Simplex algorithm to
iteratively minimize the loss monotonically on adjacent vertices and analyze
the trajectory of these vertex positions. We empirically observe that in a
neighbourhood around a local minimum, the iterations behave differently such
that conclusions on loss level and proximity of the local minimum can be made
before it has been found: Firstly the loss seems to decay exponentially slow
at iterated adjacent vertices such that the loss level at the local minimum
can be estimated from the loss levels of subsequently iterated vertices, and
secondly we observe a strong increase of the vertex density around local
minima. This could have far-reaching consequences for the design of new
gradient-descent algorithms that might improve convergence rate by exploiting
these facts.

Key words: Deep ReLU programming, Loss surface, local minima, ReLU neural
network, vertex structure

##  I INTRODUCTION AND MOTIVATION

Due to its relevance for the design of efficient neural network training
algorithms, a deeper understanding of the loss landscape of neural networks is
of great interest and a topic of active research. However, being a non-convex
function of a typically very high number of parameters which is usually
explored only point-wise or path-wise by training algorithms, insight in its
structure requires sophisticated ideas and techniques.

Related literature covers both, theoretical results for specific settings and
experimental studies. In [1] it is shown that for squared error loss, specific
networks and limited training samples the training loss behaves nicely in the
sense that almost all local minima are global. The work [2] studies loss
visualization techniques and qualitatively relates geometrical properties,
generalization error and trainability. The authors of [3] and [4]
experimentally analyze the loss landscape, its properties such as curvature
and Hessian matrix along the trajectories of various optimization algorithms
and compare their resulting converged local minima. In this work we perform a
novel analysis of the vertex structure of the layer-wise L1 loss of feed-
forward ReLU neural networks.

Feed-forward ReLU neural networks with $L\in\mathbb{N}$ hidden layers and
layer widths $(n_{0},\ldots,n_{L+1})\in\mathbb{N}^{L+2}$ are functions
$f_{\theta}:\mathbb{R}^{n_{0}}\to\mathbb{R}^{n_{L+1}}$ of the form

$f_{\theta}(x)=W_{L+1}g_{\theta_{L}}\circ\cdots\circ
g_{\theta_{1}}(x)+b_{L+1}\quad\textnormal{ for }x\in\mathbb{R}^{n_{0}}$ (1)

with layer transition functions
$g_{\theta_{l}}:\mathbb{R}^{n_{l-1}}\to\mathbb{R}^{n_{l}},x\mapsto\textnormal{ReLU}(W_{l}x+b_{l})$,
$l\in\left\\{1,\cdots,L\right\\}$ that use a coordinate-wise application of
the activation function $\textnormal{ReLU}:t\mapsto\max(0,t)$. The parameter
$\theta$ determines the weight and bias matrices, i.e.
$\theta=(\theta_{1},\ldots,\theta_{L+1})$ with
$\theta_{l}=(W_{l},b_{l})\in\mathbb{R}^{n_{l}\times
n_{l-1}}\times\mathbb{R}^{n_{l}}$, $l\in\left\\{1,\ldots,L+1\right\\}$. These
functions are piece-wise affine on convex regions which have the form of a
feasible region of linear inequalities.

We recently published the _Deep ReLU Simplex (DRLSimplex)_ method [5] which is
a novel algorithm for the iterative minimization of such functions. It
iterates on vertices of feed-forward ReLU neural networks and is closely
related to the simplex algorithm from linear programming. The main extension
of the DRLSimplex over the standard simplex algorithm is the ability to
efficiently change the feasible region at a vertex. The path of visited
vertices has the important property that the corresponding function values are
non-increasing. Figure 1 qualitatively depicts the iteration behaviour: In the
initial phase, a set of active constraints is constructed while the initial
affine region determined by the starting position is not left and takes the
role of the feasible region in linear programming (steps depicted in green).
After a vertex has been found, the feasible region is changed to some of the
adjacent regions and the iteration continues by finding an adjacent vertex and
changing the feasible region (steps depicted in red). The next vertex is
chosen such that the objective function is decreased and the iteration
terminates when the objective function value is larger at all adjacent
vertices. This way, after a finite number of iterations the DRLSimplex
algorithm converges at a local minimum.

Figure 1: Qualitative illustration of the DRLSimplex algorithm iterating on
vertices of a continuous, piece-wise affine function given in the form of a
feed-forward ReLU neural network. In the first phase (green) a vertex is found
similarly to standard linear programming and in subsequent iterations then the
feasible regions is changed.

The L1 training loss

$L(\theta)=\sum_{i=1}^{N}\sum_{j=1}^{n_{L+1}}|(y_{i})_{j}-f_{\theta}(x_{i})_{j}|$
(2)

of $N\in\mathbb{N}$ samples
$(x_{1},y_{1}),\ldots,(x_{N},y_{N})\in\mathbb{R}^{n_{0}}\times\mathbb{R}^{n_{L+1}}$
itself can be written as a ReLU feed-forward neural network when only one
layer shall be optimized while keeping the other layer parameters fixed.
Therefore, we can apply the DRLSimplex algorithm to generate a sequence of
network parameters (weight matrix and bias vector entries) corresponding to a
specific layer which have non-increasing L1 loss value. We analyze the
behaviour of this sequence of vertex positions and loss values as we
approximate a local minimum. The DRLSimplex algorithm is well suited for this
kind of analysis due to the following reasons:

  1. 1.

Compared to gradient descent, where the step size is usually a function of
past step sizes or decaying at a predetermined rate, the step size of the
DRLSimplex algorithm is automatically determined in every iteration by the
proximity of the next vertex in the selected direction. This way we directly
get access to local properties of the underlying objective function, i.e. the
loss landscape.

  2. 2.

Gradient descent algorithms do not guarantee that the training loss values
decay in every iteration. Without proper step size control, they do not
necessarily converge. In contrast, the DRLSimplex algorithm is guaranteed to
terminate in a local minimum without any adjustments.

From these two facts it is clear that the DRLSimplex algorithm provides a
finite loss-decreasing trajectory of adjacent vertices which has a direct
relation to the converged local minimum and is _not_ affected by an arbitrary
externally defined factors such as step size control rules. This makes this
algorithm perfectly suited for an analysis on how local vertex density and
corresponding quantities behave as we approximate local minima.

##  II SIMULATION AND RESULTS

For simplicity we only consider the case of first layer L1 optimization, i.e.
we search a local minimum of the function

$L_{1,\theta^{*}}:\theta_{1}\mapsto
L\left((\theta_{1},\theta^{*}_{2},\ldots,\theta^{*}_{L+1})\right).$ (3)

with fixed weight matrices and bias vectors
$\theta_{2}^{*},\ldots,\theta^{*}_{L+1}$ for the other layers. Note that this
is not a strong restriction since by equations (2) and (1) for every layer the
L1 loss optimization of this layer’s parameters can be interpreted as first-
layer parameter training of the smaller network starting with that layer and
using different training predictors transformed by the lower layers instead of
the original predictors.

As we have shown in [5] our DRLSimplex can be implemented efficiently by
making use of the insight and best-practices concerning performance and
numerical stability that evolved over decades in the field of linear
programming. However, such an implementation still has to be written and
tested such that at the present state we rely on complete matrix inversion
instead of simplex tableau modifications to favour correctness over
performance. Therefore we were are only able perform our analysis on a small
scale.

Our simulation was carried out as follows: We picked the network architecture
$L=4$ and $(n_{0},\ldots,n_{5})=(4,5,4,3,2,1)$ such that there are
$n_{1}(n_{0}+1)=25$ first layer parameters consisting of the $5\times 4$
weight matrix and the $5$-element bias vector that are applied to the network
input. For the fixed other layers’ parameters
$\theta^{*}_{2},\ldots,\theta^{*}_{L+1}$ we sampled all entries from a uniform
distribution on $[-1,1]$. We then generated $N=500$ random training samples
$(x_{1},y_{1}),\ldots,(x_{500},y_{500})\in\mathbb{R}^{4}\times\mathbb{R}$ all
with entries from independent uniform distributions on $[-3,3]$. Finally we
sampled a random initial position with uniformly distributed entries on
$[-20,20]$ in the $25$-dimensional parameter space of first layer parameters
to start the iteration there using our DRLSimplex algorithm. We performed our
analysis several times also with a different number of layers and other layer
widths and the results were always similar. We will illustrate our findings in
a representative simulated case of such a optimization trajectory of vertices
ending in a local minimum $\hat{\theta}_{1}$ generated by our DRLSimplex
algorithm.

In Figures 2, 3 and 4 we depict the path of the training loss, a smoothed
curve of the L2 distance of consecutive iterated parameter estimates and the
distance to the local minimum vertex that our algorithm has finally converged
to. At the beginning of the training process, our DRLSimplex algorithm
iterates in the boundary of the same initial feasible region while
successively adding hyperplane constraints until a vertex has been found. This
process generates the first 25 iteration positions. In consequence, all
iterated parameter estimates starting with index 25 are vertices of the layer-
wise loss (3).

Figure 2: First layer training loss (3) at iterated first layer parameters
returned by our DRLSimplex algorithm against iteration number in logarithmic
scale. The loss sequence is non-increasing by construction. Figure 3: Running
mean of 40 consecutive L2-distances of iterated first layer parameters (the
first 39 mean values are based on fewer observations).  Figure 4: L2-distance
of iterated first-layer parameters to the final converged local minimum
$\hat{\theta}_{1}$ (last element of sequence of iterated adjacent vertices)
returned by the DRLSimplex algorithm.

After a number of initial training iterations we observe that the remaining
training process can be partitioned in two phases:

  1. 1.

In the _exponential decay phase_ we observe that the training loss, the L2
distance between successive iterated parameter estimates and the distance to
the finally converged minimum all decrease roughly exponentially against the
iteration index as can be seen in Figures 2, 3 and 4 respectively. This is a
typical observation we made in each of our repeated simulations. In this phase
we approximate the local minimum in a systematic and stable way. The fact that
our considered iterations are adjacent vertices of $L_{1,\theta^{*}}$ allows
us to conclude that the vertex density increases exponentially during this
phase as we approximate the local minimum $\hat{\theta}_{1}$.

  2. 2.

In the _fine tuning phase_ we observe a slower loss decrease than in the
exponential decay tuning phase. The step size of consecutive vertices
iterations and the distance to the final local minimum do not show a
systematic behaviour across our simulations. We conclude that in contrast to
the exponential decay phase, this phase is highly individual and case-
specific.

##  III INTERPRETATION

The decreasing distance observed in the trajectory of adjacent vertices
approaching a local minimum corresponds to a higher vertex density, i.e. more
and smaller affine regions around local minima. If the number of regions is
considered as a natural measure of function complexity, this observation
indicates that the function complexity of the considered layer-wise L1 loss
function is greater around local minima.

Furthermore, the two phases we observed in the training path of iterated
parameter estimates indicate that the neighbourhood of local minima of the
layer-wise L1 loss is also two fold as shown in Figure 5:

Figure 5: Schematic interpretation of neighbourhood of layer-wise L1 loss
local minima $\hat{\theta}_{1}$ based on our results.

  1. 1.

There is an _outer exponential neighbourhood_ that behaves essentially
similarly for all local minima we observed in the sense that the vertex
distance, the distance to the local minimum and the loss level decay
exponentially on a sequence of adjacent vertices approaching the local
minimum.

  2. 2.

Furthermore there is an _inner specific neighbourhood_ where the loss decay of
such a sequence is slower and both, the vertex distance and the distance to
the local minimum are individual and do not have a typical behaviour.

The general behaviour repeatedly observed in the outer exponential
neighbourhood suggests that the vertex structure of most local minima behaves
consistently in a similar way in their neighbourhood. This fact could be used
to extract information about a local minimum based on the iterated parameter
values. For example, we can already give a good estimate on the loss value at
a local minimum $\hat{\theta}_{1}$ that the DRLSimplex algorithm finally
converges to before that local minimum has been found, based on an exponential
decay fit as soon as this decay regime is recognized. In other words, without
having found a local minimum the trajectory of vertices returned by our
DRLSimplex algorithm already reveals insight in its loss level.

In future research one could analyze this _information leak_ of local minima
onto their neighbourhood. Furthermore, in the light of neural network training
one can speculate on the relation of over-fitting and the inner specific
neighbourhood, and a deeper analysis would be a second interesting research
direction.

##  IV CONCLUSION

In this work we considered the L1 training loss of ReLU feed-forward neural
network as a function of the first layer’s parameters and for fixed parameters
in the other layers. This function is piece-wise affine on convex regions and
induces a vertex structure. We used the DRLSimplex algorithm to obtain a loss-
decreasing sequence of adjacent vertices ending in a local minimum. Our
numerical experiments results suggest that locally around local minima the
vertex density of the layer-wise L1 loss is larger, which indicates that there
are more and smaller regions located around local minima and this can be
interpreted as higher function complexity.

A further insight is the fact that that layer-wise L1 loss local minima seem
to leak information to their neighbourhood and that this information can be
retrieved from the steps of an iterating optimization procedure. We used our
DRLSimplex algorithm to get access to local properties of the loss function
such as vertex density and in future research it might be interesting to
analyze to what extent also the trajectory of iteration parameters generated
by popular training procedures carries information on proximity, direction or
loss level of nearby local minima that could be used to adjust step size
control in gradient-descent like minimization procedures and ultimately
improve convergence rate of neural network training algorithms.

## References

  * [1] Q. Nguyen and M. Hein, “The Loss Surface of Deep and Wide Neural Networks,” vol. 70, pp. 2603–2612, 06–11 Aug 2017.
  * [2] H. Li, Z. Xu, G. Taylor, C. Studer, and T. Goldstein, “Visualizing the loss landscape of neural nets,” arXiv preprint arXiv:1712.09913, 2017.
  * [3] D. J. Im, M. Tao, and K. Branson, “An empirical analysis of the optimization of deep network loss surfaces,” arXiv preprint arXiv:1612.04010, 2017\.
  * [4] Q. Yuan and N. Xiao, “Experimental exploration on loss surface of deep neural network,” International Journal of Imaging Systems and Technology, vol. 30, no. 4, pp. 860–873, 2020.
  * [5] P. Hinz and S. van de Geer, “Deep ReLU Programming,” arXiv preprint arXiv:2011.14895, 2021.