11institutetext: Technical University of Darmstadt, Karolinenpl. 5, 64289
Darmstadt, Germany
11email<EMAIL_ADDRESS>
# How Reliable Are Out-of-Distribution Generalization Methods for Medical
Image Segmentation? ††thanks:  Supported by the Bundesministerium für
Gesundheit (BMG) with grant [ZMVI1-2520DAT03A]. The final authenticated
version of this manuscript will be published in Lecture Notes in Pattern
recognition in the life and natural sciences - DAGM GCPR 2021 (doi will
follow).

Antoine Sanner  0000-0002-4917-9529 Camila González  0000-0002-4510-7309
Anirban Mukhopadhyay  0000-0003-0669-4018

###### Abstract

The recent achievements of Deep Learning rely on the test data being similar
in distribution to the training data. In an ideal case, Deep Learning models
would achieve Out-of-Distribution (OoD) Generalization, i.e. reliably make
predictions on out-of-distribution data. Yet in practice, models usually fail
to generalize well when facing a shift in distribution. Several methods were
thereby designed to improve the robustness of the features learned by a model
through Regularization- or Domain-Prediction-based schemes. Segmenting medical
images such as MRIs of the hippocampus is essential for the diagnosis and
treatment of neuropsychiatric disorders. But these brain images often suffer
from distribution shift due to the patient’s age and various pathologies
affecting the shape of the organ. In this work, we evaluate OoD Generalization
solutions for the problem of hippocampus segmentation in MR data using both
fully- and semi-supervised training. We find that no method performs reliably
in all experiments. Only the V-REx loss stands out as it remains easy to tune,
while it outperforms a standard U-Net in most cases.

###### Keywords:

Semantic segmentation Medical Images Out-of-Distribution Generalization

##  1 Introduction

Semantic segmentation of medical images is an important step in many clinical
procedures. In particular, the segmentation of the hippocampus from MRI scans
is essential for the diagnosis and treatment of neuropsychiatric disorders.
Automated segmentation methods have improved vastly in the past years and now
yield promising results in many medical imaging applications [11]. These
methods can technically exploit the information contained in large datasets.
However, no matter how large a training dataset is or how good the results on
the in-distribution data are, methods may fail on Out-of-Distribution (OoD)
data. OoD Generalization remains crucial for the reliability of deep neural
networks, as insufficient generalization may vastly limit their implementation
in practical applications.

Distribution shifts occur when the data at test time is different in
distribution from the training data. In the context of hippocampus
segmentation, the age of the patient [13] and various pathologies [4] can
affect the shape of the organ. Using a different scanner or simply using
different acquisition parameters can also cause a distribution shift [5].

OoD Generalization alleviates this issue by training a model such that it
generalizes well to new distributions at test time without requiring any
further training. Several strategies exist to approach this. Regularization-
based approaches enforce the learning of robust features across training
datasets [2, 9], which can reliably be used to produce an accurate prediction
regardless of context. On the other hand, Domain-Prediction-based methods
focus on harmonizing between domains [6, 7] by including a domain predictor in
the architecture.

In this work, we perform a thorough evaluation of several state-of-the-art OoD
Generalization methods for segmenting the hippocampus. We find that no method
performs reliably across all experiments. Only the Regularization-based V-REx
loss stands out as it remains easy to tune, while its worst results remains
relatively good.

##  2 Methods

The setting for Out-of-Distribution Generalization has been defined [2] as
follows. Data is collected from multiple environments and the source of each
data point is known. An environment describes a set of conditions under which
the data has been measured. One can for instance obtain a different
environment by using a different scanner or studying a different group of
patients. These environments contain spurious correlations due for instance to
dataset biases.

Only fully labeled datasets from a limited set of environments is available
during training, and the goal is to learn a predictor, which performs well on
all environments. As such, a model trained using this method will
theoretically perform well on unseen but semantically related data.

Figure 1:  Network architecture for the segmentation task with Domain-
prediction mechanism. The U-Net generates features, which are then used by the
segmentor for the main task and by the auxiliary head for domain prediction.
This schema contains some simplifications: the segmentor and domain predictor
are implemented respectively using a convolutional layer and a CNN. The
architecture displayed is 2D, while the implementation is 3D.

###  2.1 Regularization-based methods

Arjovsky et al. [2] formally define the problem of OoD Generalization as
follows. Consider the datasets
$D_{e}:=\\{(x_{i}^{e},y_{i}^{e})\\}^{n_{e}}_{i=1}$ collected under multiple
training environments $e\in\mathcal{E}_{tr}$. These environments describe the
same pair of random variables measured under different conditions. The dataset
$D_{e}$, from environment $e$, contains examples identically and independently
distributed according to some probability distribution $P(X_{e},Y_{e})$. The
goal is then to learn a predictor $Y\approx f(X)$, which performs well across
all unseen related environments $\mathcal{E}_{all}\supset\mathcal{E}_{tr}$.
The objective is to minimize:

$R^{OoD}(f)=\max_{e\in\mathcal{E}_{all}}\mathcal{R}^{e}(f)$

where $\mathcal{R}^{e}(f):=\mathbb{E}_{X^{e},Y^{e}}[\ell(f(X^{e}),Y^{e})]$ is
the risk under environment $e$.

IRMv1 [2] and Risk Extrapolation (V-REx and MM-Rex) [9] add regularization to
the training loss to enforce strict equality between training risks by finding
a data representation common to all environments. IRM Games [1] poses the
problem as finding the Nash equilibrium of an ensemble game.

More precisely, Krueger et al. [9] propose to solve this problem by minimizing
the following loss (V-REx):

$\min_{\mathcal{X}\rightarrow\mathcal{Y}}\lambda\cdot\mathrm{Var}\\{R^{1}(\Phi),...,R^{n_{e}}(\Phi)\\}+\frac{1}{|\mathcal{E}_{tr}|}\sum_{e\in\mathcal{E}_{tr}}R^{e}(\Phi)$

where $\Phi$ is the invariant predictor, the scalar $\lambda\geq 0$ controls
the balance between reducing average risk and enforcing the equality of risks,
and $Var$ stands for the variance between the risks across training
environments. To the best of our knowledge, these approaches have not been
evaluated yet on image segmentation tasks.

###  2.2 Domain-Prediction-based methods

Given data from multiple sources, Domain-Prediction-based methods find a
harmonized data representation such that all information relating to the
source domain of the image is removed [6, 7, 14]. This goal can be achieved by
appending another head to the network, which acts as a domain classifier.
During training, the ability of the domain classifier to predict the domain is
minimized to random chance, thus reducing the domain-specific information in
the data representation. Domain-Prediction-based methods differ from
Regularization-based methods in that they rather remove the need to annotate
all images from a new target dataset to be able to train a model on new
unlabeled data, and so allow leveraging non-annotated data.

As described by Dinsdale et al. [6], the network shown in Fig. 1 is a modified
U-Net with a second head which acts as a domain predictor. The final goal is
to find a representation that maximizes the performance on a segmentation task
with input images $\mathbf{X_{d}}\in\mathbb{R}^{B_{d}\times W\times H\times
D}$ and task labels $\mathbf{Y_{d}}\in\mathbb{R}^{B_{d}\times W\times H\times
D\times C}$ while minimizing the performance of the domain predictor. This
network is composed of an encoder, a segmentor, and a domain predictor with
respective weights $\Theta_{repr}$, $\Theta_{seg}$, and $\Theta_{dp}$.

The network is trained iteratively by minimizing the loss on the segmentation
task $L_{s}$, the domain loss $L_{d}$, and the confusion loss $L_{conf}$. The
confusion loss penalizes a divergence of the domain predictor’s prediction
from a uniform distribution and is used to remove source-related information
from $\Theta_{repr}$. It is also important that the segmentation loss $L_{s}$
be evaluated separately for the data from each dataset to prevent the
performance being driven by only one dataset, if their sizes vary
significantly. In our work, the segmentation loss $L_{s}$ takes the form of
the sum of a Sorensen-Dice loss and a binary cross-entropy loss. The domain
loss $L_{d}$ is used to assess how much information remains in about the
domains.

The Domain-Prediction method thus minimizes the total loss function:

$\displaystyle
L(\textbf{X},\textbf{Y},\textbf{D},\Theta_{repr},\Theta_{seg},\Theta_{dp})=$
$\displaystyle\sum_{d\in\mathcal{E}_{tr}}L_{s}(\mathbf{X_{d}},\mathbf{Y_{d}};\Theta_{repr},\Theta_{seg})$
$\displaystyle+\alpha\cdot
L_{d}(\textbf{X},\textbf{D},\Theta_{repr};\Theta_{dp})$
$\displaystyle+\beta\cdot
L_{conf}(\textbf{X},\textbf{D},\Theta_{dp};\Theta_{repr})$

where $\alpha$ and $\beta$ represent weights of the relative contributions for
the different loss functions. $(\mathbf{X_{d}},\mathbf{Y_{d}})$ corresponds to
the labeled data available for the domain $d$, while $(\mathbf{X},\mathbf{D})$
corresponds to all the images (labeled and unlabeled) available and their
corresponding domain. This is interesting as training the domain predictor
does not require labeled data and semi-supervised training can be used to find
a harmonized data representation.

###  2.3 A Combined Method for OoD Generalization

The losses introduced in the context of Regularization-based OoD
Generalization can be combined to the Domain-Prediction mechanism. We choose
to apply the V-REx loss to the segmentor only. The intuition is that if the
segmentor learns robust features through regularization, then there will be
less source-related information in the extracted features and can make the
Domain-Prediction part of the method easier. Besides, Domain Prediction
schemes already require to compute the segmentation loss for each domain
separately, which is also a requirement for Regularization-based methods. This
method aims to minimize the following loss function:

$\displaystyle\min_{\mathcal{X}\rightarrow\mathcal{Y}}\lambda$
$\displaystyle\cdot\mathrm{Var}\\{L_{s}(\mathbf{X_{1}},\mathbf{Y_{1}};\Theta_{repr},\Theta_{seg}),...,L_{s}(\mathbf{X_{n_{d}}},\mathbf{Y_{n_{d}}};\Theta_{repr},\Theta_{seg})\\}$
$\displaystyle+\sum_{d=1}^{n_{d}}L_{s}(\mathbf{X_{d}},\mathbf{Y_{d}};\Theta_{repr},\Theta_{seg})$
$\displaystyle+\alpha\cdot
L_{d}(\mathbf{X},\textbf{D},\Theta_{repr};\Theta_{dp})$
$\displaystyle+\beta\cdot
L_{conf}(\mathbf{X},\textbf{D},\Theta_{dp};\Theta_{repr})$

##  3 Experimental Setup

In the following, we present the datasets used in this work and detail the
evaluation strategy used.

###  3.1 Datasets

We use a corpus of three datasets for the task of hippocampus segmentation
from various studies. The Decath dataset [12] contains MR images of both
healthy adults and schizophrenia patients. The HarP dataset [3] is the product
of an effort from the European Alzheimer’s Disease Consortium and Alzheimer’s
Disease Neuroimaging Initiative to harmonize the available protocols and
contains MR images of senior healthy subjects and Alzheimer’s disease
patients. Lastly, the Dryad dataset [10] contains MR images of healthy young
adults.

Since both HarP and Dryad datasets provide whole-head MR images at the same
resolution, we crop each scan at two fixed positions, respectively for the
right and left hippocampus. The shape of the resulting crops (64, 64, 48) fits
every hippocampus in both datasets. Since segmentation instances from the
Decath dataset are smaller, they are centered and zero-padded to the same
shape. Similarly, as all the datasets do not provide the same number of
classes in their annotation, we restrict the classes to “background” and
“hippocampus”.

###  3.2 Evaluation

We first train three instances of a 3D U-Net, using 3 layers for encoding, and
decoding with no dropout, respectively on each of the datasets to assess how
well the features learned on a dataset generalize to the other datasets. We
then evaluate the proposed method with five-fold cross-validation. In the
context of fully supervised training, a model is trained on two of the
datasets and evaluated on the third dataset. 10% of the data in each fold is
used for validation and guiding the training schedule. The dataset that is not
used to train the model is considered entirely as test data. Due to the nature
of the task, we need to find a hyper-parameter/data augmentation configuration
that allows the model to generalize well on the test dataset no matter which
pair of training dataset is used.

Figure 2:  Schema of testing/validation/training splits in the context of
fully supervised and semi-supervised training. The test splits remain equal
while the segmentor’s training and validation splits are subsets of the splits
for the domain predictor. “D.P.” refers to the domain predictor and “Seg.” to
the segmentor.

For semi-supervised training, we train our models on all 3 datasets. The
splits for the domain predictor are computed the same way as during fully
supervised training. As for the segmentor, its training and validation splits
are subsets of the domain predictor’s respective splits, as shown in Fig. 2.
The test splits for both heads remain equal. The splits are computed only once
and reused for each method to reduce the variance during testing.

The network is implemented using Python 3.8 and PyTorch 1.6.0. Some light data
augmentation is also used: RandomAffine, RandomFlip, RandomMotion. The batch
size for each dataset is selected so that the training duration is the same
with all datasets.

##  4 Results

In this section, we first evaluate the generality of the features learned by a
standard U-Net on each dataset. We then compare the different methods in a
fully supervised settings, and we outline their limitations. Finally, we
assess whether Domain-Prediction-based methods can leverage their ability of
training on unlabeled data.

###  4.1 U-Net results on each dataset

We train three instances of a standard 3D U-Net respectively on each of the
datasets, in order to assess how well features learned by a model on a given
dataset generalize well to the other two datasets. The results can be seen in
Table 1. Each of three trained models achieves a mean Dice score in the range
of $84$ to $90$ with a low standard deviation below 2. As a comparison,
Isensee et al. [8], Carmo et al. [4], and Zhu and al. [15] achieve
respectively on Decath, HarP, and Dryad a $90$, $86$, and $89$ Dice scores.
So, these results are close to state-of-the-art, which is not our goal per se,
but they should give a good insight on how features learned in one dataset
generalize to the other ones.

Training dataset | Decath | HarP | Dryad
---|---|---|---
Decath only | 88.5 ± 1.4 | 28.7 ± 12.0 | 17.1 ± 14.3
HarP only | 76.5 ± 2.4 | 84.6 ± 1.4 | 80.0 ± 1.7
Dryad only | 59.0 ± 4.8 | 57.9 ± 1.9 | 85.2 ± 1.9
Table 1:  Dice coefficient, comparing generalization capability of the
features learned by the U-Net in each dataset using a standard no OoD
Generalization mechanism. The results on the diagonal are for the respective
in-distribution dataset of each model and otherwise the results correspond to
OoD testing.

The performance of the three models on OoD Generalization differ greatly.
While the model trained on HarP already achieves great results, the other two
models seem to struggle much more to generalize and attain a much lower
average and higher standard deviation. So, training each of the three pairs of
datasets provides a different setting and allows testing multiple scenarios.

###  4.2 Fully supervised Training

For all dataset configuration, we use two datasets for training and the third
one for OoD testing. The training datasets are also referred to as in-
distribution datasets.

####  4.2.1 Training on Decath and HarP:

In this first setup, we train all models on Decath and HarP. The results can
be seen in Table 2. Training a U-Net only on HarP already learns features that
generalize well, so we are interested in seeing how adding “lesser quality”
data (Section 4.1) influences the results on Dryad. The results are shown in
Table 2. If we first focus on the results on the in-distribution datasets, we
observe that all methods achieve state-of-the-art results. All methods achieve
satisfying results on the out-of-distribution dataset with the combined method
having a lead. The domain prediction accuracy reaches near random results for
all Domain-Prediction-based methods, which indicates that the domain-
identifying information is removed.

Method | Decath | HarP | Dryad | Dom. Pred. Acc.
---|---|---|---|---
U-Net | 89.5 ± 0.4 | 85.7 ± 0.7 | 81.2 ± 1.5 | -
U-Net + V-REx | 89.4 ± 0.6 | 85.3 ± 1.2 | 81.7 ± 0.5 | -
Domain-Prediction | 90.0 ± 0.4 | 86.8 ± 1.3 | 82.8 ± 2.5 | 48.6 ± 6.2
Combined | 88.9 ± 0.7 | 85.2 ± 0.8 | 84.6 ± 0.4 | 56.1 ± 10.9
Table 2:  Dice coefficients and domain prediction accuracy (when applicable).
A bold dataset name denotes that the dataset is used for training.

####  4.2.2 Training on Decath and Dryad:

This training setup is perhaps the most interesting one of the three, as we
have seen that training a model on either of Decath and Dryad does not yield
good generalization results. As shown in Table 3, all methods yield state-of-
the-art results on both in-distribution datasets. Surprisingly, all methods
perform well on HarP, although the regular U-Net to a lesser degree. Regarding
domain prediction accuracy, the Domain-Prediction method reaches once more
near random chance results. The domain prediction accuracy is particularly
high for the combined method, since these results stem from an intermediary
training stage where the segmentor and domain predictor are trained jointly.

Method | Decath | HarP | Dryad | Dom. Pred. Acc.
---|---|---|---|---
U-Net | 88.7 ± 0.7 | 69.8 ± 2.4 | 89.2 ± 0.5 | -
U-Net + V-REx | 89.3 ± 0.2 | 73.7 ± 1.9 | 89.9 ± 0.3 | -
Domain-Prediction | 89.6 ± 0.4 | 72.9 ± 2.4 | 90.5 ± 0.3 | 56.5 ± 11.5
Combined | 89.6 ± 0.2 | 73.2 ± 1.5 | 90.5 ± 0.3 | 90.0 ± 6.3
Table 3:  Dice coefficients and domain prediction accuracy (when applicable).
A bold dataset name denotes that the dataset is used for training.

####  4.2.3 Training on HarP and Dryad:

In this last setup, we train the models on HarP and Dryad and test them on
Decath. The results can be seen in Table 4. All methods perform almost equally
well on in-distribution datasets. However, they differ vastly on the first
column, with the U-Net + V-REx method being the clear winner. The regular
U-Net and the Domain-Prediction method perform similarly, with the latter
having a smaller variance. Removing domain-related information from the U-Net
features does not seem to be enough as both Domain-Prediction methods achieve
a near random domain prediction accuracy, which indicates that domain-
identifying information has been removed, but the combined method vastly
underperforms.

Method | Decath | HarP | Dryad | Dom. Pred. Acc.
---|---|---|---|---
U-Net | 65.6 ± 18.1 | 84.9 ± 0.8 | 89.8 ± 0.7 | -
U-Net + V-REx | 78.1 ± 2.6 | 85.3 ± 0.8 | 90.1 ± 0.3 | -
Domain-Prediction | 69.4 ± 6.4 | 86.3 ± 0.8 | 90.6 ± 0.6 | 54.4 ± 7.0
Combined | 46.2 ± 10.6 | 84.2 ± 1.1 | 89.7 ± 0.9 | 52.6 ± 9.3
Table 4:  Dice coefficients and domain prediction accuracy (when applicable).
The usual data augmentation scheme is used. A bold dataset name denotes that
the dataset is used for training.

If we observe the segmentation masks predicted by models from different folds
of the same cross-validation, we can see in Fig. 3 that while the first
prediction matches the ground truth nicely, the second prediction fails to
fully recognize the hippocampus and is very noisy. What is troublesome is that
both models achieve Dice scores of $86\%$ and $90\%$ respectively for the test
splits of in-distribution datasets. The domain prediction accuracy is reduced
to random chance without causing any significant drop in performance on the
segmentation task. As such, the Domain-Prediction method reduces domain
information in the model features as expected. Reducing domain prediction
accuracy does not seem to be enough however, as the combined method also
achieves low domain detection results, but fares worse than a regular U-Net on
the test dataset. Overall, we have no way of predicting whether a model is
usable on new data by looking at the metrics on the validation data.

Figure 3:  Segmentation mask prediction for an instance from Decath. Both
models are trained using the Domain-Prediction-based method but with different
training data folds and achieve respectively $78.9\%$ and $59.9\%$ Dice score
on the out-of-distribution dataset Decath. Dice scores on the test splits of
the in-distribution datasets are within $1\%$ of each other.

It is worth noting that by adding some RandomBiasField and RandomNoise to the
data augmentation scheme, both Domain-Prediction-based methods perform
significantly better, as can be seen in Table 5. The results on the training
datasets remain stable with a decrease in standard deviation. The Dice scores
on Decath, however, significantly improve. However, this change cannot be
justified by just looking at the metrics on validation data. Furthermore, this
scheme yields significantly worse results when training on any other dataset
pair, which is why it is not used anywhere else.

Method | Decath | HarP | Dryad | Dom. Pred. Acc.
---|---|---|---|---
Domain-Prediction | 75.1 ± 1.7 | 84.2 ± 0.9 | 90.0 ± 0.3 | 57.8 ± 7.6
Combined | 74.1 ± 4.3 | 82.3 ± 1.2 | 89.0 ± 0.6 | 61.2 ± 3.8
Table 5:  Dice coefficients and domain prediction accuracy (when applicable)
for a different data augmentation scheme. A bold dataset name denotes that the
dataset is used for training.

####  4.2.4 Discussion

OoD Generalization methods seem to be quite limiting in the sense that they
require the user to find a hyper-parameter configuration that works well in
all setups. The results on the validation data are too good to allow us to
fine tune the hyper-parameters for each method in each of the settings. As
such, no method tested gives reliable results, although the U-Net + V-REx
method always performs at least on par, and often better than, the reference
U-Net.

We also observe that we have no way of predicting whether a model is usable on
new data by looking at the metrics on the validation data. The Dice scores on
the validation data are not representative of how well a model generalizes on
out-of-distribution data. A near random guessing domain prediction accuracy is
also not sufficient, as seen in Section 4.2.3.

###  4.3 Semi-supervised Training

During the fully supervised training, we noted that Domain-Prediction-based
methods trained using HarP and Dryad can fail to fully recognize the
hippocampus and produce very noisy segmentations. The main goal of these
experiments is to test whether adding a few labeled data points from Decath to
the training set can alleviate this issue.

Being able to train on unlabeled data is also a very interesting feature for a
method, considering that image annotation is a time-consuming task. In the
next titles, if present, the numbers between square brackets next to a
dataset’s name mean refer respectively to the number of labeled data points
used for training and validation from this dataset in addition to the
unlabeled data points.

####  4.3.1 Training on Decath[5,5], HarP and Dryad:

In this experiment, we train the models on all datasets but only HarP and
Dryad are fully annotated. Decath is almost fully unlabeled. Only 10 labeled
points are considered and are split between training and validation. The
results are shown in Table 6.

Method | Decath[5, 5] | HarP | Dryad | Dom. Pred. Acc.
---|---|---|---|---
U-Net | 88.1 ± 0.6 | 85.9 ± 1.0 | 90.5 ± 0.5 | -
U-Net + V-REx | 87.8 ± 0.8 | 85.7 ± 1.0 | 90.5 ± 0.4 | -
Domain-Prediction | 87.4 ± 0.7 | 85.4 ± 1.4 | 90.0 ± 0.5 | 56.6 ± 5.2
Combined | 86.9 ± 1.7 | 85.8 ± 1.1 | 90.3 ± 0.6 | 57.2 ± 1.5
Table 6:  Dice coefficients and domain prediction accuracy (when applicable).
A bold dataset name denotes that the dataset is used for training. If present,
the numbers between brackets next to an in-distribution dataset’s name refer
respectively to the number of data points used for training and validation
from this dataset.

All methods perform similarly on the segmentation task and reach state-of-the-
art results on all datasets. The Domain-Prediction-based methods decrease the
domain prediction accuracy, but it remains far better than random chance
($55\%$ vs $33\%$). During training, sudden variations in metrics for the
segmentation tasks were observed even at learning rates as low as $5\cdot
10^{-5}$. However, the results achieved here are only slightly worse than the
one reached during fully-supervised training. As such, we expect that refining
this tuning only leads to marginal gains on the segmentation task. Overall, no
method here is able to outperform the others and yield significant improvement
over using a standard U-Net.

####  4.3.2 Training on Decath[6,4], HarP[6,4] and Dryad[6,4]:

From what we have seen in Section 4.3.1, training models on two fully labeled
datasets and one partially labeled one does not allow discriminating between
methods. Instead, we can consider that we have multiple dataset available, but
only a few labeled instances for each of them. With only 18 instances being
used for training, we want to see whether one method outperforms the others.

Method | Decath[6, 4] | HarP[6, 4] | Dryad[6, 4] | Dom. Pred. Acc.
---|---|---|---|---
U-Net | 86.7 ± 0.7 | 80.2 ± 1.4 | 87.3 ± 1.2 | -
U-Net + V-REx | 86.9 ± 0.5 | 80.7 ± 1.0 | 87.5 ± 1.0 | -
Domain-Prediction | 82.9 ± 3.3 | 75.7 ± 4.8 | 85.6 ± 1.5 | 56.1 ± 4.4
Combined | 85.6 ± 0.5 | 78.4 ± 2.0 | 86.4 ± 0.7 | 54.7 ± 3.7
Table 7:  Dice coefficients label and domain prediction accuracy (when
applicable). A bold dataset name denotes that the dataset is used for
training. If present, the numbers between parentheses next to an in-
distribution dataset’s name refer respectively to the number of data points
used for training and validation from this dataset.

As shown in Table 7, all methods perform almost equally well. The results for
the U-Net and U-Net + V-REx methods are impressive, considering that there is
only a slight drop in performance now that we are training these models using
far fewer labeled data points. The results for the Domain-Prediction-based
methods are slightly worse, and the domain prediction accuracy remains
relatively high for both methods. These results are coherent with the ones
obtained during the previous experiment as the same hyper-parameter
combination is used, although the combined method achieves better results than
the standard Domain-Prediction method on the second column.

####  4.3.3 Discussion

Overall, no method reliably outperforms the standard U-Net. Domain-Prediction-
based methods actually performs worse than the reference in the last
experiment. Either removing domain related information does not scale well
with the number of training domains, or a much finer hyper-parameter tuning is
required to allow these methods to leverage their ability to train on
unlabeled data. The feasibility of this tuning can also be put in question due
to the increase in instability observed during training.

##  5 Conclusion

In this work, we evaluate a variety of methods for Out-of-Distribution
Generalization on the task of hippocampus segmentation. In particular, we
evaluate Regularization-based methods which regularize the performance of a
model across training environments. We also consider a Domain-Prediction-based
method, which adds a domain classifier to the architecture, its goal being to
reduce the prediction accuracy of this classifier. Lastly, we explore how well
a method uniting both approaches performs.

We compare these methods in a fully supervised setting, and we observe the
limitations of OoD Generalization methods. No method performs reliably in all
experiments. Only the V-REx loss stands out as it remains easy to tune, while
its worst results remains close to the reference U-Net.

To gauge the ability of Domain-Prediction-based methods to train on unlabeled
data, we subsequently evaluate all methods in a semi-supervised settings,
using the minimum number of labeled images required to retain a stable
training trajectory. The model trained with a V-REx loss maintains impressive
results despite the lack of training data, while Domain-Prediction-based
methods show their limits when training on an increased number of domains.

The combined method achieves good results in some settings, but suffers from
instability in others. Domain-Prediction-based methods have a lot of
potential, but require a lot of fine-tuning and can cause the training to slow
down.

In the future, we wish to evaluate these methods on another corpus of
datasets. This different setting will allow us to evaluate whether the results
observed here still hold, especially whether Domain-Prediction-based methods
can make use of training on unlabeled data during semi-supervised training. We
also wish to confirm whether using the V-REx loss remains a good option in
another setting, as it achieves similar or better results than a standard loss
on this corpus of datasets.

## References

  * [1] Ahuja, K., Shanmugam, K., Varshney, K.R., Dhurandhar, A.: Invariant risk minimization games. CoRR (2020), http://arxiv.org/abs/2002.04692
  * [2] Arjovsky, M., Bottou, L., Gulrajani, I., Lopez-Paz, D.: Invariant risk minimization. ArXiv (2020), http://arxiv.org/abs/1907.02893
  * [3] Boccardi, M., Bocchetta, M., Morency, F.C., Collins, D.L., Nishikawa, M., Ganzola, R., Grothe, M.J., Wolf, D., Redolfi, A., Pievani, M., Antelmi, L., Fellgiebel, A., Matsuda, H., Teipel, S., Duchesne, S., Jack, C.R., Frisoni, G.B.: Training labels for hippocampal segmentation based on the EADC-ADNI harmonized hippocampal protocol. Alzheimer’s & Dementia 11(2), 175–183 (2015). https://doi.org/10.1016/j.jalz.2014.12.002, https://linkinghub.elsevier.com/retrieve/pii/S155252601402891X
  * [4] Carmo, D., Silva, B., Yasuda, C., Rittner, L., Lotufo, R.: Hippocampus segmentation on epilepsy and alzheimer’s disease studies with multiple convolutional neural networks. Heliyon (2020), http://arxiv.org/abs/2001.05058
  * [5] Castro, D.C., Walker, I., Glocker, B.: Causality matters in medical imaging. Nature Communications 11(1) (Jul 2020). https://doi.org/10.1038/s41467-020-17478-w, https://doi.org/10.1038/s41467-020-17478-w
  * [6] Dinsdale, N.K., Jenkinson, M., Namburete, A.I.L.: Deep learning-based unlearning of dataset bias for mri harmonisation and confound removal. bioRxiv (2020). https://doi.org/10.1101/2020.10.09.332973, https://www.biorxiv.org/content/early/2020/12/14/2020.10.09.332973
  * [7] Ganin, Y., Lempitsky, V.: Unsupervised domain adaptation by backpropagation. ArXiv (2015), http://arxiv.org/abs/1409.7495
  * [8] Isensee, F., Petersen, J., Klein, A., Zimmerer, D., Jaeger, P.F., Kohl, S., Wasserthal, J., Koehler, G., Norajitra, T., Wirkert, S., Maier-Hein, K.H.: nnU-net: Self-adapting framework for u-net-based medical image segmentation. ArXiv (2018), http://arxiv.org/abs/1809.10486
  * [9] Krueger, D., Caballero, E., Jacobsen, J.H., Zhang, A., Binas, J., Priol, R.L., Courville, A.: Out-of-distribution generalization via risk extrapolation (REx). CoRR (2020), http://arxiv.org/abs/2003.00688
  * [10] Kulaga-Yoskovitz, J., Bernhardt, B.C., Hong, S.J., Mansi, T., Liang, K.E., van der Kouwe, A.J., Smallwood, J., Bernasconi, A., Bernasconi, N.: Multi-contrast submillimetric 3 tesla hippocampal subfield segmentation protocol and dataset. Scientific Data 2(1), 150059 (2015). https://doi.org/10.1038/sdata.2015.59, http://www.nature.com/articles/sdata201559
  * [11] Litjens, G., Kooi, T., Bejnordi, B.E., Setio, A.A.A., Ciompi, F., Ghafoorian, M., van der Laak, J.A., van Ginneken, B., Sánchez, C.I.: A survey on deep learning in medical image analysis. Medical Image Analysis 42 (Dec 2017). https://doi.org/10.1016/j.media.2017.07.005, https://doi.org/10.1016/j.media.2017.07.005
  * [12] Simpson, A.L., Antonelli, M., Bakas, S., Bilello, M., Farahani, K., van Ginneken, B., Kopp-Schneider, A., Landman, B.A., Litjens, G., Menze, B., Ronneberger, O., Summers, R.M., Bilic, P., Christ, P.F., Do, R.K.G., Gollub, M., Golia-Pernicka, J., Heckers, S.H., Jarnagin, W.R., McHugo, M.K., Napel, S., Vorontsov, E., Maier-Hein, L., Cardoso, M.J.: A large annotated medical image dataset for the development and evaluation of segmentation algorithms. CoRR (2019), http://arxiv.org/abs/1902.09063
  * [13] Xu, Y., Valentino, D.J., Scher, A.I., Dinov, I., White, L.R., Thompson, P.M., Launer, L.J., Toga, A.W.: Age effects on hippocampal structural changes in old men: The haas. NeuroImage 40(3), 1003–1015 (2008). https://doi.org/https://doi.org/10.1016/j.neuroimage.2007.12.034, https://www.sciencedirect.com/science/article/pii/S105381190701141X
  * [14] Xue, Y., Feng, S., Zhang, Y., Zhang, X., Wang, Y.: Dual-task self-supervision for cross-modality domain adaptation. In: Martel, A.L., Abolmaesumi, P., Stoyanov, D., Mateus, D., Zuluaga, M.A., Zhou, S.K., Racoceanu, D., Joskowicz, L. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2020, vol. 12261, pp. 408–417. Springer International Publishing (2020), http://link.springer.com/10.1007/978-3-030-59710-8˙40, series Title: Lecture Notes in Computer Science
  * [15] Zhu, H., Shi, F., Wang, L., Hung, S.C., Chen, M.H., Wang, S., Lin, W., Shen, D.: Dilated dense u-net for infant hippocampus subfield segmentation. Frontiers in Neuroinformatics 13, 30 (2019). https://doi.org/10.3389/fninf.2019.00030, https://www.frontiersin.org/article/10.3389/fninf.2019.00030/full


# Time Reversal of the Overdamped Langevin Equation and Fixman’s Law

Oliver T. Dyer  Robin C. Ball<EMAIL_ADDRESS>Dept of Physics,
University of Warwick, Coventry CV4 7AL, UK.

###### Abstract

We discuss how the first order Langevin equation for the overdamped dynamics
of an interacting system has a natural time reversal of simple but surprising
form, with consequences for correlation functions. This leads to the
correlation of interactions as a strictly restraining term in the time-
dependent diffusion tensor of the system, deriving the relation first
suggested by Fixman. Applying this to the time-dependent diffusion of dilute
polymer coils leads to the quantitative calibration of Kirkwood’s
approximation for their hydrodynamic radius. We find the generalized “Fixman
Law” for dissipation with a memory kernel, which has revealing causal
structure, and we also discuss the case of the second order Langevin Equation.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

The first order Langevin equation Langevin (1908)

$\frac{dx_{i}}{dt}=M_{ij}f_{j}(x(t))+k_{\textrm{B}}T\partial_{k}M_{ki}+u_{i}(t)$
(1)

describes the dynamics of a subsystem of coordinates $x_{i}(t)$ with
interaction forces $f_{i}(x)$, embedded in (or part of) a larger background
system in thermal equilibrium at temperature $T$ with short timescale
response. Being rooted in equilibrium it should have a time reversal symmetry,
and we show how to find this and non-trivial fluctuation theorem results which
follow. In later sections we show a notable application and discuss the
generalisation of our results to both inertial memory and Langevin dynamics in
memory media.

The background system being in thermal equilibrium contributes random velocity
contributions $u_{i}(t)$ with correlation $\langle
u_{i}(t)u_{j}(t^{\prime})\rangle=2k_{\textrm{B}}TM_{ij}(x)\delta(t-t^{\prime})$
related to the mobility $M_{ij}$ by the Fluctuation-Dissipation Theorem Kubo
(1957). We adopt Itô calculus in taking the $u_{i}(t)$ to be uncorrelated with
present as well as past positions, requiring us to include the explicit
Brownian drift term $k_{\textrm{B}}T\partial_{k}M_{ki}$ Lau and Lubensky
(2007).

We now seek a time reversal of Eq. (1) with $t^{R}=-t$ and $x^{R}(t^{R})=x(t)$
obeying

$\frac{dx_{i}^{R}}{dt^{R}}=M_{ij}f_{j}(x^{R}(t^{R}))+k_{\textrm{B}}T\partial_{k}M_{ki}+u_{i}^{R}(t^{R})$
(2)

which requires that the time reversal of the random velocity contribution be
given by

$u_{i}^{R}(t^{R})=-u_{i}(t)-2v_{i}(t),$ (3)

where $v_{i}(t)=M_{ij}f_{j}(x(t))+k_{\textrm{B}}T\partial_{k}M_{ki}$ is the
deterministic contribution to the forwards motion of the system.

At first sight the deterministic contribution to the time reverse random
velocity is puzzling because the random velocity should be unbiased, but this
is conditional on the direction of time: $u(t)$ is uncorrelated with past and
present $x(t^{\prime}\leq t)$. Correspondingly $u^{R}(t^{R})$ should be
uncorrelated with future and present $x(t^{\prime}\geq t)$, and the difference
of condition matters because $x(t^{\prime}\geq t)$ clearly cumulates influence
from the earlier random velocity terms $u(t)$.

Useful correlation identities follow from Eq. (3) by considering that for
$t^{R}\geq t^{\prime R}$ we have $\langle u_{i}^{R}(t^{R})v_{j}^{R}(t^{\prime
R})\rangle=0$ where $v_{i}^{R}(t^{R})=v_{i}(t)$. Substituting back in terms of
unreversed quantities then leads to

$\langle u_{i}(t)v_{j}(t^{\prime})\rangle=-2\langle
v_{i}(t)v_{j}(t^{\prime})\rangle,\quad t<t^{\prime}$ (4)

whereas $\langle u_{i}(t)v_{j}(t^{\prime})\rangle=0$ for $t>t^{\prime}$. From
this the full velocity autocorrelation of the interacting system can be
expressed as

$\langle\dot{x_{i}}(t)\dot{x_{j}}(t^{\prime})\rangle=2M_{ij}/(k_{\textrm{B}}T)\,\delta(t-t^{\prime})-\langle
v_{i}(t)v_{j}(t^{\prime})\rangle.$ (5)

The power of Eq. (5) is that it expresses the diffusion of the interacting
system in terms of the bare non-interacting value _minus_ a restraining
contribution from the deterministic motion. This result was first proposed by
Fixman by Diffusion Equation arguments Fixman (1981), in the polymer context
(without Brownian drift) which we discuss below, but in a manner which did not
convince later authors albeit they found some numerical evidence to support it
when $t\simeq t^{\prime}$ Liu and Dünweg (2003). As it is now proved, we will
refer to it as _Fixman’s Law_.

## Application to Polymer Diffusion

The archetypal application is to the motion of colloid and polymer systems in
a Newtonian solvent, where the $x_{i}(t)$ are the $3N$ coordinates of
$N$particle position vectors $\vec{r}_{n}(t)$ and radius $a_{n}$, with joint
mobility tensor given by $O_{nn}=I/\left(6\pi\eta a_{n}\right)$ for diagonal
blocks (with $I$ the three dimensional identity matrix) and essentially that
of Oseen with $M\rightarrow
O_{mn}=\left(I+\hat{r}\hat{r}\right)/\left(8\pi\eta|\vec{r}_{m}-\vec{r}_{n}|\right)$
for off-diagonal blocks $m\neq n$, but matching onto the diagonal case at
close approach Rotne and Prager (1969); Dyer and Ball (2017). Being based on
divergence-free Stokes flow, these mobility tensors have zero Brownian drift.
For an isolated polymer coil of $N$ beads Eq. (5) then becomes

$\langle\dot{\vec{r}}_{n}(t)\dot{\vec{r}}_{m}(t^{\prime})\rangle=2k_{\textrm{B}}T\,\langle
O_{mn}\rangle\delta(t-t^{\prime})-\langle\vec{v}_{n}(t)\vec{v}_{m}(t^{\prime})\rangle$
(6)

where the leading term is Kirkwood’s approximationKirkwood (1954) and all the
time dependent memory is captured in the counter terms correlating the
velocity contribution $\vec{v}_{m}(t)=\sum_{n}O_{mn}\cdot\vec{F}_{n}(t)$
driven by the conservative forces holding the polymer chain together.

Focussing for simplicity on the coil centre of mass
$\vec{R}(t)=\sum_{n}\vec{r}_{n}(t)/N$, the corresponding velocity
autocorrelation function is then given by

$\langle\vec{V}(t)\vec{V}(t^{\prime})\rangle=2D_{\textrm{K}}\delta(t-t^{\prime})I-A_{\textrm{F}}(t-t^{\prime})I$
(7)

where $D_{\textrm{K}}=k_{\textrm{B}}TN^{-2}\sum_{mn}\langle O_{mn}\rangle$ is
the classical Kirkwood diffusivityKirkwood (1954) and all the memory in the
polymer centre of mass motion comes from the Fixman counter-term

$A_{\textrm{F}}(t-t^{\prime})I=N^{-2}\\!\\!\\!\\!\sum_{mnm^{\prime}n^{\prime}}\\!\langle
O_{mm^{\prime}}(t)\cdot\vec{F}_{m^{\prime}}(t)\vec{F}_{n^{\prime}}(t^{\prime})\cdot
O_{n^{\prime}n}(t^{\prime})\rangle.$ (8)

The long time polymer coil diffusivity is then given by

$D_{L}=D_{\textrm{K}}-\Delta=D_{\textrm{K}}-\int_{0}^{\infty}A_{\textrm{F}}(\tau)d\tau$
(9)

and given that $D_{L}$ is the natural experimental measurement whilst
$D_{\textrm{K}}$ is a direct configurational moment, it is important to
understand their difference $\Delta$, all arising from the Fixman term.

Figure 1: Autocorrelation of the force driven Centre of Mass velocities
$A_{\textrm{F}}(t)$ (as defined in Eq. 8) _vs_ time $t$ for swollen chains
(lower curves) and for Gaussian chains (upper curves with $A_{\textrm{F}}$
shifted up by a factor of 10). With units naturally scaled in terms of
measured Zimm time $\tau_{\textrm{Z}}$ and Kirkwood (short time) diffusion
coefficient $D_{\textrm{K}}$, there is good superposition for the wide range
of chain length $N$ shown. The plots exclude measurements at times short in
absolute (rather than scaled) units, for which each curve breaks away from the
superposition (see text).

Figure 2: Upper curves: fractional decrease in diffusivity, found by
numerically integrating the data in Fig. 1 for the autocorrelation function
$A_{\textrm{F}}(t)$, plotted against the upper time limit of integration.
Chain lengths are as per legend in Fig. 1. Data has been truncated when noise
begins to dominate additional contributions leading to incomplete curves,
especially for swollen chains. The lower curves show the same data shifted
vertically to a common master curve, with plateau value zero and shown in
gray, the shift for each curve then giving an estimate of the eventual long
time plateau of the relative decrease in diffusivity. The sections of each
curve below the gray arcs in the upper plots are omitted to make the master
curves clearer. Fig. 5 in the Supplementary Material gives a different
presentation of these master curves.

In Fig. 1 we show the time-dependence of $A_{\textrm{F}}(t)$ in Gaussian and
swollen chains measured in Wavelet Monte Carlo dynamics simulations Dyer and
Ball (2017); Dyer (2019). In each case the longer time data is tolerably
consistent with the natural scaling form
$A_{\textrm{F}}(t)=D_{\textrm{K}}/\tau_{\textrm{Z}}~{}\mathfrak{h}(t/\tau_{\textrm{Z}})$
as plotted where we have used the relaxation time $\tau_{\textrm{Z}}$ of the
longest Rouse mode as the natural scale of time, and this behaviour on its own
would lead upon time integration of $A_{\textrm{F}}$ to a contribution to
$\Delta$ proportional to $D_{\textrm{K}}$. Outside the range of data shown in
Fig. 1 we also see short time structure
$A_{\textrm{F}}(t)=\mathfrak{g}(N)\mathfrak{f}(t)$ on monomeric timescales
which on its own would lead to a contribution to $\Delta$ proportional to
$\mathfrak{g}(N)$, and for the particular case $t=0$ earlier work Liu and
Dünweg (2003) showed that $\mathfrak{g}(N)$ does appear to approach $1/N$
albeit slowly.

In Fig. 2 we show the corresponding cumulative integrals for
$\Delta(t)=\int_{0}^{t}A_{\textrm{F}}(t^{\prime})dt^{\prime}$ in scaled
variables, as far as we can measure them above noise. Assuming these plots do
go far enough to capture all the short time parts, the total values of
$\Delta_{N}/D_{\textrm{K}}$ inclusive of the scaling parts are estimated by
finding the vertical shifts to align to a master curve in the scaled variables
with plateau value zero.

The measured estimates of $\Delta_{N}/D_{\textrm{K}}$ are shown in Fig. 3
plotted against the anticipated correction to scaling which is
$N^{-1}/R_{H}^{-1}$. These enable us to give the first quantitative estimates
of the asymptotic values for true long time diffusion coefficient compared to
the Kirkwood short time formulaKirkwood (1954),

$\lim_{N\to\infty}\frac{D_{L}-D_{\textrm{K}}}{D_{\textrm{K}}}=\begin{cases}(-3.4\pm
0.5)\%&\textrm{swollen chains}\\\ (-8.9\pm 0.5)\%&\textrm{Gaussian
chains}.\end{cases}$ (10)

These values rmk  considerably firm up earlier ones, where for phantom chains
Fixman Fixman (1986) concluded the deficit was _around_ 8% whilst for swollen
chains Liu & Dünweg Liu and Dünweg (2003) concluded it was _at least_ 3.5%.

These corrections matter wherever a measured chain diffusivity or mobility,
such as might be observed by dynamic light scattering or sedimentation
respectively, is compared with the Kirkwood formula given in terms of
configurational statistics. Direct determination of these corrections to our
accuracy from simulation without use of the Fixman law would be very hard, as
this would require direct measurement of $D$ to better than $1\%$ which would
entail simulation (with full hydrodynamics included) in excess of $10^{4}$
chain relaxation times for each chain length. Direct determination by
experiment would face the added difficulties of determining the
configurational statistics to better than $1\%$ and controlling any influence
of polydispersity.

Figure 3:  Plateau values of $\Delta/D_{\textrm{K}}$ for different polymer
chain lengths, plotted against the correction to scaling variable
$1/(ND_{\textrm{K}})$. Data for long chains were obtained using the master
curves in Fig. 2 to extrapolate out to the plateau. The error bars on these
data show the difference between this master curve extrapolation and the
estimated upper bound from a linear extrapolation. For the data that did reach
$t/\tau_{\textrm{Z}}=1$, the error bars show twice the standard deviation of
$\Delta/D_{\textrm{K}}$ measured across all simulations contributing to each
data point. The dashed lines to the left of the data indicate fits of the 4
leftmost data points with the function $m/(ND_{\textrm{K}})+c$, and the data
at $1/(ND_{\textrm{K}})=0$ show the extrapolated infinite chain value with
error bars calculated as $2\sqrt{\sigma_{cc}^{2}}$, with $\sigma^{2}$ the
covariance matrix for $m$ and $c$.

##  I Second Order Langevin Equation

The basis for the second order Langevin equation is slightly different, but it
leads to a matching time reversal result. Working in terms of forces we have

$m_{ij}\frac{d^{2}x_{j}}{dt^{2}}+Z_{ij}\frac{dx_{j}}{dt}=f_{i}(x(t))+\phi_{i}(t).$
(11)

with random forces $\phi_{i}(t)=Z_{ij}u_{j}$ where $Z=M^{-1}$ is the impedance
matrix and $m$ an inertial tensor.

We now pose the analogous time reversed equation

$m_{ij}\frac{d^{2}x^{R}_{j}}{(dt^{R})^{2}}+Z_{ij}\frac{dx^{R}_{j}}{dt^{R}}=f_{i}(x^{R}(t^{R}))+\phi^{R}_{i}(t^{R}),$
(12)

with $x^{R}(t^{R}=-t)=x(t)$ as before, and by inspection this requires

$\phi^{R}_{i}(t^{R})=\phi_{i}(t)-2Z_{ij}\frac{dx_{j}}{dt}.$ (13)

This looks quite different from the first order case, but it turns out to
exactly agree in the limit of inertia being negligible. Multiplying through by
the mobility matrix to obtain the random currents, eliminating $Z\cdot dx/dt$
using the second order Langevin equation (11), and then noting that
$M_{ij}\ddot{x}_{j}\rightarrow-k_{\textrm{B}}T\partial_{j}M_{ji}$ in the limit
for the first order Langevin equation leads back to the time reversal of the
random currents as in Eq. (3).

##  II Generalisation to Memory Media

Finally we have considered the generalisation to the first order Langevin
equation for motion in a time dependent (or memory) medium, seeking insight
into the nature of time reversal, and status of Fixman’s law in relation to
other fluctuation-dissipation type results. The generalisation turns out to be
quite distinctive, confirming that Fixman’s law (5) is self-standing result
and not a rewriting of prior fluctuation theorems. We start from

$\frac{dx}{dt}=\int G(t-t^{\prime})f(x(t^{\prime}))dt^{\prime}+u(t)=G*f+u$
(14)

where $G(t-t^{\prime})$ is the response function of the medium which we now
take to depend only on explicit time difference and not on $x(t)$, so there is
no Brownian drift. Causality dictates that $G(t-t^{\prime})=0$ for
$t<t^{\prime}$, and for the forces we will use the shorthand notation
$f(t)=f(x(t))$. We leave indices implicit now but for the multivariate case
$G$ is a matrix, in terms of which the Fluctuation-Dissipation Theorem Kubo
(1957) gives

$\left\langle
u(t)u^{T}(t^{\prime})\right\rangle=kT\left(G(t-t^{\prime})+G^{T}(t^{\prime}-t)\right).$
(15)

We will allow the memory medium to be non-reversible, in particular in terms
of its off-diagonal elements of $\left\langle
u(t)u^{T}(t^{\prime})\right\rangle$. This then means that whilst we seek a
time reversal transformation preserving the form of the Langevin Equation as a
class, it will no longer in general be a symmetry for each case. The
transformation is found as before by negating all time arguments to $t^{R}=-t$
in our equation of motion, leading to $-dx^{R}/dt^{R}=G^{R}*f^{R}+u(-t^{R})$,
where $x^{R}(t^{R})=x(t)=x(-t^{R})$ and similarly for $f^{R}$ and $G^{R}\equiv
G^{R}(t^{R})=G(t)=G(-t^{R})$. This matches the natural time reverse equation

$dx^{R}/dt^{R}={G^{T}}*f^{R}+u^{R}(t^{R})$ (16)

by taking

$u^{R}(t^{R})=-u(-t^{R})-(G^{T}+G^{R})*f^{R}.$ (17)

This time reversal of the random velocities (17) features terms $G^{T}*f^{R}$
which is causal in reversed time and $G^{R}*f^{R}$ which is anti-causal in
their propagation of the conservative forces, complicating derivation of
results from it. We show in the supplementary material that it follows at some
length that

$\langle uf^{T}\rangle=-(G+G^{TR})*\langle ff^{T}\rangle_{-}$ (18)

where $\langle f(t)f^{T}(t^{\prime})\rangle_{-,+}$ are respectively the
$t<t^{\prime}$ and $t>t^{\prime}$ parts of $\langle
f(t)f^{T}(t^{\prime})\rangle$. It is then a matter of straightforward
substitution to find the generalisation of Eq. (8) as

$\begin{split}A_{\textrm{F}}&=\left\langle\dot{x}\dot{x}^{T}\right\rangle_{\textrm{free}}-\left\langle\dot{x}\dot{x}^{T}\right\rangle\\\
&=G*\left\langle ff^{T}\right\rangle_{+}*G+G^{TR}*\left\langle
ff^{T}\right\rangle_{-}*G^{TR}.\end{split}$ (19)

Figure 4:  Diagrammatic representation of Fixman Law for a memory medium. Full
lines denote the total correlation whilst dashed lines denote their
correlations without interaction, as per the background medium. All the
correlations are consistently directed in time as shown by the arrows.

This memory medium version of Fixman’s law (19) is richer in structure than
conventional Fluctuation-Dissipation Theorem type results. Specialising to
$t>t^{\prime}$ for clarity, $A_{\textrm{F}}(t-t^{\prime})$ can be expressed as

$\int\displaylimits_{t>u>u^{\prime}>t^{\prime}}\mkern-20.0mudu\,du^{\prime}\left\langle
G(t-u)\\!\cdot\\!f(u)\;(G^{R}(t^{\prime}-u^{\prime})\\!\cdot\\!f(u^{\prime}))^{T}\right\rangle$

which is the correlation of the causal response to the applied causal forces
with the corresponding anti-causal response, with the added restriction that
only causally ordered factors of forces contribute. Substituting the $G$ in
terms of correlation functions the full result is represented diagrammatically
in Fig. 4, where the time ordering of its terms can be seen to be natural.

Further specialisations, generalisations and applications await attention.
Without memory media Langevin equations describe a Markov process so
everything can be worked through in terms of Fokker-Planck equations. For the
dynamics in memory media we have only considered their response depending on
the time lag, whereas it is natural to consider complex fluid media with
response dependent on the full configurational history. The application to
colloid dynamics is interesting because Fixman’s law separates the
contribution of hydrodynamic and conservative interactions, whilst for
polymers we would like to generalise to encompass rheological response
functions.

###### Acknowledgements.

This work has been supported by the Engineering and Physical Sciences Research
Council (EPSRC), Grant No.s EP/L505110/1 and EP/M508184/1 (PhD studentship for
O.T.D.). Computing facilities were provided by the Scientific Computing
Research Technology Platform of the University of Warwick.

## References

  * Langevin (1908) P. Langevin, C. R. Acad. Sci. (Paris) 146, 530 (1908), [English translation at D. S. Lemons and A. Gythiel, Am. J. Phys. 65, 1079 (1997)].
  * Kubo (1957) R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
  * Lau and Lubensky (2007) A. W. C. Lau and T. C. Lubensky, Phys. Rev. E 76, 011123 (2007).
  * Fixman (1981) M. Fixman, Macromolecules 14, 1710 (1981), https://doi.org/10.1021/ma50007a019 .
  * Liu and Dünweg (2003) B. Liu and B. Dünweg, The Journal of Chemical Physics 118, 8061 (2003), https://doi.org/10.1063/1.1564047 .
  * Rotne and Prager (1969) J. Rotne and S. Prager, The Journal of Chemical Physics 50, 4831 (1969), https://doi.org/10.1063/1.1670977 .
  * Dyer and Ball (2017) O. T. Dyer and R. C. Ball, The Journal of Chemical Physics 146, 124111 (2017), https://doi.org/10.1063/1.4978808 .
  * Kirkwood (1954) J. G. Kirkwood, J. Pol. Sci. XII, 1 (1954).
  * Dyer (2019) O. T. Dyer,  _Wavelet Monte Carlo dynamics_ , Ph.D. thesis, University of Warwick, UK (2019).
  * (10) We have rounded up the statistical error on the Gaussian limiting value, as it is based on rather scattered error bars in Fig. 3.
  * Fixman (1986) M. Fixman, The Journal of Chemical Physics 84, 4080 (1986), https://doi.org/10.1063/1.450070 .
  * Uhlenbeck and Ornstein (1930) G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. 36, 823 (1930).

##  III Supplementary Material

##  IV Fixman’s Law in a Memory Medium

Figure 5: Fractional decrease in diffusivity, found by numerically integrating
the data in Fig. 1 for the autocrrelation function $A_{\textrm{F}}(t)$,
plotted against the upper time limit of integration, for various chain lengths
$N$. Master curves shown in gray were obtained by vertical linear shift of the
data in the long time regime, and these are shown in gray behind (and
extending) the curves for each $N$. The shifted master curves then give an
estimate of the eventual long time plateau of the relative decrease in
diffusivity for each curve.

Our general equation of motion is

$\frac{dx}{dt}=\int G(t-t^{\prime})f(x(t^{\prime}))dt^{\prime}+u(t)$ (20)

where for causality $G(t-t^{\prime})=0,\quad t<t^{\prime}$ and we will also
allow ourselves the shorthand notation $f(t)=f(x(t))$. For the multivariate
case, that is where $x$, $f$ and $u$ all become time dependent vectors, the
above displayed equation applies with $G$ being a matrix and we have
$\left\langle
u(t)u^{T}(t^{\prime})\right\rangle=k_{\textrm{B}}T\left(G(t-t^{\prime})+G^{T}(t^{\prime}-t)\right)$.

It proves useful to express the random velocity contributions as

$u(t)=\int K(t-t^{\prime})h(t^{\prime})dt^{\prime}$ (21)

where $\left\langle h(t)h^{T}(t^{\prime})\right\rangle=\delta(t-t^{\prime})I$
is white noise and $K(t-t^{\prime})=0$ for $t<t^{\prime}$ is strictly causal.
Then we can interpret $h(t)$ as the innovation, that is what is new in the
noise at time $t$, with the important consequence that $h(t)$ is uncorrelated
with everything from earlier times, including

$\left\langle h(t)f^{T}(t^{\prime})\right\rangle=0,\;t>t^{\prime}.$ (22)

The defining property of $K$ is that the autocorrelation of $u(t)$
reconstructs correctly, so we require
$k_{\textrm{B}}T\left(G(t-t^{\prime})+G^{T}(t^{\prime}-t)\right)=\int
K(t-t")K^{T}(t^{\prime}-t")dt".$

These equations have a natural time convolution structure, starting with

$dx/dt=G*f+K*h.$ (23)

We will use $R$ to denote negation of time arguments, equivalent in effect to
a transpose of the times as indices. Then we can write

$\left\langle
uu^{T}\right\rangle=k_{\textrm{B}}T\left(G+G^{TR}\right)=K*K^{TR}.$ (24)

The above can be interpreted as an L-U matrix decomposition with respect to
indices corresponding to times, and we will also have need of a partner U-L
decomposition so that

$k_{\textrm{B}}T\left(G+G^{TR}\right)=K*K^{TR}=L^{TR}*L.$ (25)

where $L$ is strictly causal.

###  IV.1 Time Reversal Transformation in terms of Innovation

If we simply substitute all time arguments in terms of $t^{R}=-t$ in our
equation of motion we obtain

$-dx^{R}/dt^{R}=G^{R}*f^{R}+K^{R}*H.$ (26)

Here $x^{R}(t^{R})=x(t)=x(-t^{R})$ and similarly for $f^{R}$, $G^{R}$ and
$K^{R}$ as in the memoryless case, and $H(t_{R})=h(-t^{R})$. However
$G^{R}\equiv G^{R}(t^{R})=G(t)=G(-t^{R})$ and similarly for $K^{R}$ have the
effect that they are propagating backwards in reversed time $t^{R}$, so the
form of the equation of motion is clearly not preserved under this simple time
negation.

The time reversed equation we want is

$dx^{R}/dt^{R}=G^{T}*f^{R}+L^{T}*h^{R},$ (27)

where $G^{T}$ and $L^{T}$ are now propagating forwards in their implied
arguments of reverse time. The transpose on $G$ here is simply the natural
time reverse of the Fluctuation-Dissipation Theorem, and we have chosen our
notation for $L$ to match. This requires

$h^{R}=-(L^{T})^{-1}*\left((G^{T}+G^{R})*f^{R}+K^{R}*H\right),$ (28)

where convolution by $(L^{T})^{-1}$ inverts convolution by $L^{T}$.

The time reversed partner to Eq. (22) is $\left\langle
h^{R}(t^{R})f^{TR}({t^{\prime}}^{R})\right\rangle=0$ for
$t^{R}>{t^{\prime}}^{R}$. Substituting $h_{R}$ from Eq.(28) this leads to

$\begin{split}(L^{T})^{-1}*&K^{R}*\langle Hf^{TR}({t^{\prime}}^{R})\rangle\\\
&=-(L^{T})^{-1}*(G^{T}+G^{R})*\langle
f^{R}f^{TR}({t^{\prime}}^{R})\rangle,\end{split}$ (29)

all for $t^{R}>{t^{\prime}}^{R}$. This can then be expressed in terms of
unreversed times $t<t^{\prime}$ as

$(L^{TR})^{-1}*K*\langle hf^{T}({t^{\prime}})\rangle=-L*\langle
ff^{T}({t^{\prime}})\rangle/(k_{\textrm{B}}T)$ (30)

where on the RHS we have made use of
$(L^{TR})^{-1}*(G+G^{TR})=(L^{TR})^{-1}*K*K^{R}/(k_{\textrm{B}}T)=L/(k_{\textrm{B}}T)$
all from Eq. (25). We can also use (25) to substitute
$(L^{TR})^{-1}*K=L*(K^{TR})^{-1}$ on the LHS of Eq. (30) to obtain

$L*(K^{TR})^{-1}*\langle hf^{T}({t^{\prime}})\rangle=-L*\langle
ff^{T}({t^{\prime}})\rangle/(k_{\textrm{B}}T)$ (31)

for $t<t^{\prime}$.

###  IV.2 The cross correlation of $h$ and $f$

We now seek to bring equations (22) and (31) together as a single expression
for $\left\langle hf^{T}\right\rangle$. Given that (31) is valid for
$t<t^{\prime}$ it follows that it remains valid if we convolve both sides
forwards to a new time $t$ which is still less than $t^{\prime}$, and using
$L^{-1}*$ in particular then leads to the top line in

$(K^{TR})^{-1}*\left\langle
hf^{T}\right\rangle(t,t^{\prime})=-\frac{1}{k_{\textrm{B}}T}\begin{cases}\left\langle
ff^{T}\right\rangle&t<t^{\prime}\\\ 0&t>t^{\prime}\end{cases}.$ (32)

The bottom line of the above follows simply by convolving Eq. (22) by
$(K^{TR})^{-1}$ which being anticausal preserves its time inequality. Finally
it is convenient to denote the braced factor of the RHS above as $\left\langle
ff^{T}\right\rangle_{-}$ and then write the whole relation as

$\left\langle hf^{T}\right\rangle=-K^{TR}*\left\langle
ff^{T}\right\rangle_{-}/(k_{\textrm{B}}T)$ (33)

which can be convolved by $K$ to yield also

$\left\langle uf^{T}\right\rangle=-(G+G^{TR})*\left\langle
ff^{T}\right\rangle_{-}.$ (34)

###  IV.3 Application to correlation function

Now we come to the heart of the matter and consider

$\begin{split}\left\langle\dot{x}\dot{x}^{T}\right\rangle=&\left\langle
uu^{T}\right\rangle+G*\left\langle ff^{T}\right\rangle*G^{TR}\\\
&+\left\langle uf^{T}\right\rangle*G^{TR}+G*\left\langle
fu^{T}\right\rangle\end{split}$ (35)

where the first term direct from the innovation autocorrelation corresponds to
the “free” value of the LHS in the absence of any internal forces $f$. In the
third term we can substituted $\left\langle uf^{T}\right\rangle$ using the
previous subsection IV.2 to give

$\left\langle uf^{T}\right\rangle*G^{TR}=-(G+G^{TR})*\left\langle
ff^{T}\right\rangle_{-}*G^{TR}$ (36)

and the fourth term is this value both transposed and time reversed for which
we need to note $\left(\left\langle
ff^{T}\right\rangle_{-}\right)^{TR}=\left\langle ff^{T}\right\rangle_{+}$
where $\left\langle ff^{T}\right\rangle=\left\langle
ff^{T}\right\rangle_{+}+\left\langle ff^{T}\right\rangle_{-}$, leading net of
some cancellations to

$\begin{split}A_{\textrm{F}}(t-t^{\prime})&=\left\langle\dot{x}\dot{x}^{T}\right\rangle_{\textrm{free}}-\left\langle\dot{x}\dot{x}^{T}\right\rangle\\\
&=G^{TR}*\left\langle ff^{T}\right\rangle_{-}*G^{TR}+G*\left\langle
ff^{T}\right\rangle_{+}*G.\end{split}$ (37)

The instantaneous case introduced by Fixman now serves as a check. For
$G(t-t^{\prime})=M\delta(t-t^{\prime})$, we have $G=G^{R}$ and the above
result simplifies down to Fixman’s Law,

$\left\langle\dot{x}\dot{x}^{T}\right\rangle_{\textrm{free}}-\left\langle\dot{x}\dot{x}^{T}\right\rangle=\left\langle
vv^{T}\right\rangle$ (38)

where $v(t)=M\,f(x(t))$ is the direct contribution to the motion due to the
conservative forces.

##  V Worked Example of Fixman’s Law: Ornstein-Uhlenbeck process in memory
medium

Our Langevin equation (14) becomes linear and directly solvable by Fourier
methods, if we specialise to the case where the conservative forces are linear
in the displacement coordinates,

$f(x(t))=-s.x(t),$ (39)

where $s$ is a real symmetric matrix with non-negative eigenvalues. This
solvable Ornstein-Uhlenbeck process Uhlenbeck and Ornstein (1930) then
provides a direct check on the full memory medium version of Fixman’s law
(19).

The Fourier Transform of the Langevin equation is now

$(i\omega+G(\omega)s)x(\omega)=u(\omega)$ (40)

where $\langle
u(\omega)u^{+}(\omega)\rangle=2k_{\textrm{B}}T\left(G(\omega)+G^{+}(\omega)\right)$
and $G^{+}={G^{*}}^{T}$ is the Hermitian conjugate, in which complex
conjugation takes over the role of reversing time. It then follows directly
that

$\begin{split}\langle
x&(\omega)x^{+}(\omega)\rangle=k_{\textrm{B}}T(i\omega+Gs)^{-1}(G+G^{+})((i\omega+Gs)^{T})^{-1}\\\
&=k_{\textrm{B}}T\left((i\omega+Gs)^{-1}s^{-1}+\textrm{Hermitian
conjugate}\right).\end{split}$ (41)

We can immediately deduce the force correlations and noting that $G$ and hence
$i\omega+G$ is causal we identify

$\langle f(\omega)f(\omega)^{+}\rangle_{+}=k_{\textrm{B}}Ts(i\omega+Gs)^{-1}$
(42)

with $\langle f(\omega)f^{+}(\omega)\rangle_{-}$ given by the Hermitian
conjugate and the subscripts $\pm$ relate back to the time domain versions.

The analogously calculation for
$\langle\dot{x}(\omega)\dot{x}^{+}(\omega)\rangle_{+}$ is slightly more
involved. We start from

$\begin{split}\langle\dot{x}(\omega)\dot{x}(\omega)^{+}\rangle=&k_{\textrm{B}}T\left(G-i\omega
s^{-1}-Gs(i\omega+Gs)^{-1}G\right)\\\ +&\textrm{ Hermitian conjugate
}.\end{split}$ (43)

Then after cancelling terms in $\mp i\omega s^{-1}$ we can read off that

$\langle\dot{x}(\omega)\dot{x}(\omega)^{+}\rangle_{+}=\langle
u(\omega)u^{+}(\omega)\rangle_{+}-G\langle f(\omega)f(\omega)^{+}\rangle_{+}G$
(44)

in full agreement with Fixman’s Law extended to a memory medium (19).