\section{Introduction}
\label{sec:intro}

Sound Event Localization and Detection \textbf{(SELD)} refers to the combined task of Sound Event Detection \textbf{(SED)} and Sound Event Localization \textbf{(SEL)}, whose aim is the recognition of sound sources and their spatial location. In particular, SED requires to identify, instantaneously, the onset and offset of sound events and their correct classification, labeling the event according to the sound class that they belong to. SEL is defined as the estimation of the sound event direction in space with respect to a microphone when an event is active, referred to as Direction-of-Arrival (DOA) estimation.
Formerly, SED and SEL have been explored as two standalone tasks. Only in recent times they started to be considered jointly.
In fact, until 2018, the proposed systems considering SELD as a single task were scarce, with only one method based on deep neural network \cite{hirvonen2015classification}.
In 2018, Advanne et al.~introduced SELDnet \cite{adavanne2018sound}, a convolutional recurrent neural network (CRNN) which simultaneously recognizes, localizes and tracks sound event sources, being the first method to address the localization and recognition of more than two concurrent overlapping sound events.
The system was proposed as baseline for the Detection and Classification of Acoustic Scenes and Events (DCASE) Challenge 2019 Task 3 \cite{dcase2019task} and for the same task of DCASE Challenge 2020 \cite{dcase2020task}. This year's baseline system includes some modifications inspired by the highest ranked architectures of last year's challenge submissions, among which we can highlight \cite{kapka2019sound}, \cite{cao2019polyphonic} and \cite{zhang2019data}.
The network receives as input features log-mel spectral coefficients, together with generalized cross-correlation (GCC) for the microphone array (MIC) format, and the acoustic intensity vector in the First Order Ambisonic (FOA) domain. Regarding the network architecture, the model is initially trained with a SED loss only, and then continued with a joint SELD loss. The localization part of the joint loss is masked with the ground truth activations of each class, so that if an event is not active, it does not contribute to the training of the network.

The methodology proposed in this paper is based on the SELDNet proposed by Adavanne et ~al.~\cite{adavanne2018sound}, including some of the adopted additions in the baseline algorithm of the DCASE 2020 Task 3, such as the use of log-mel spectral coefficients and acoustic intensity vector for the FOA format. However, the system proposed in this paper differs from the baseline system presented for the DCASE 2020 Challenge in several points such as (i) data augmentation, (ii) network architecture and (iii) training loss functions. With respect to (i), -90$^{\circ}$, +90$^{\circ}$ and +180$^{\circ}$ channel rotations of the azimuth angle $\phi$ and its position reflection on the xy and xz planes are used as data augmentation technique, implementing the \textit{16 patterns} spatial augmentation proposed by Mazzon et al. ~\cite{mazzon2019sound}.
Regarding (ii), the network has been increased, adding 2 convolutional layers. Furthermore, the receptive field has been expanded using rectangular filters (instead of squared ones) in order to make the network able to recognize spectral features relevant for the task. With regard to (iii), we used the same loss functions proposed on the baseline system of last year challenge \cite{adavanne2018sound}. Instead of masked mean square error (used in the baseline system of this year), we used binary cross-entropy loss for SED prediction task and mean square error (MSE) loss for DOA estimation.

Results on the development dataset are evaluated considering the metrics proposed by Mesaros et.~al.~\cite{Mesaros_2019_WASPAA}, which take into account the joint nature of localization and detection. These are the metrics used as evaluation criteria for the challenge.

The paper is structured as follows: Section \ref{sec:method} presents the methodology and the architecture of the proposed system. Section \ref{sec:experiments} describes the experiment setup. Section \ref{sec:res} reports the development results compared with the baseline method. Conclusions and future work are presented in Section \ref{sec:conc}.

In order to promote reproducibility, the code is available under an open-source license on GitHub\footnote{https://github.com/RonFrancesca/dcase2020-fp}.
\section{METHOD}
\label{sec:method}
The method proposed for the DCASE Challenge Task 3 is based on a enhancement of the system proposed by Advanne et al. \cite{adavanne2018sound}.

\subsection{Feature Extraction}
This method uses Ambisonic format data. Log mel magnitude spectrogram together with acoustic intensity vector are used as input features for the network. Both are represented in the log mel space to better concentrate the input information of the network, as proposed by Cao et. al in \cite{cao2019polyphonic} and also implemented in the baseline system of the challenge.

FOA uses four channels to encode spatial information of a sound field, typically denoted as W, X, Y and Z. The channel W corresponds to an omnidirectional microphone recording the sound pressure. The channels X, Y and Z correspond to directional \textit{figure-of-eight} microphones oriented along the components of the homonym cartesian axes, and measure the acoustic velocity of each directional component. The acoustic intensity vector expresses the energy carried by sound waves per area unit in a direction perpendicular to that area, providing DOA information.
The intensity vector is computed as in \cite{cao2019polyphonic}.

\subsection{The network}
\label{sec:net}
Figure \ref{fig:network} shows the overall architecture of our system with the relative parameter values used in our implementation.

\tikzstyle{rect} = [rectangle, draw, text centered, minimum width = 9.7em, minimum height = 2em, align=center]
\tikzstyle{rect_large} = [rectangle, draw, text centered, minimum width = 20em, minimum height = 2em, align=center]

\begin{figure}[h]
\begin{center}
\begin{tikzpicture}[node distance = 3.5em, auto, align=center]
\centering

    \node (osc)     [rect_large, fill=lightgray] {Feature Extraction};
    \node (conv1)   [rect_large, fill=blue!10, below of = osc] {Conv 2D (64, 1xM) \\
    ReLu + MaxPool2D (5x2) + Dropout (0.2)};
    \node (conv2)   [rect_large,  fill=blue!10, below of = conv1] {Conv 2D (64, 1xM) \\
    ReLu + MaxPool2D (1x2) + Dropout (0.2)};
    \node (conv3)   [rect_large,  fill=blue!10, below of = conv2] {Conv 2D (64, 1xM) \\
    ReLu + MaxPool2D (1x2) + Dropout (0.2)};
    \node (conv4)   [rect_large,  fill=blue!10, below of = conv3] {Conv 2D (64, 1xM) \\
    ReLu + MaxPool2D (1x2) + Dropout (0.2)};
    \node (conv5)   [rect_large,  fill=blue!10, below of = conv4] {Conv 2D (64, 1xM) \\
    ReLu + MaxPool2D (1x2) + Dropout (0.2)};
    \node (GRU1)    [rect_large, fill=yellow!20, below of = conv5] {GRU (128, tahn, bidirectional) \\ Dropout (0.2)};
    \node (GRU2)    [rect_large, fill=yellow!20, below of = GRU1] {GRU (128, tahn, bidirectional) \\ Dropout (0.2)};
    \node (FC1)     [rect, xshift=-5em, fill = red!80!blue!20, below of = GRU2] {FC (128, linear)};
    \node (FC2)     [rect, xshift= 5em, fill=red!80!blue!20, below of = GRU2] {FC (128, linear)};
    \node (fc-sed)  [rect, fill = red!80!blue!20, below of = FC1] {FC (14, sigmoid)};
    \node (fc-doa)  [rect, fill=red!80!blue!20, below of = FC2] {FC (42, tahn)};

    \draw [->] (osc)  ->  (conv1);
    \draw [->] (conv1)  ->  (conv2);
    \draw [->] (conv2)  ->  (conv3);
    \draw [->] (conv3)  ->  (conv4) ;
    \draw [->] (conv4)  ->  (conv5) ;
    \draw [->] (conv5)  ->  (GRU1) ;
    \draw [->] (GRU1)  ->  (GRU2) ;
    \draw [->] (GRU2)  ->  ([xshift=-2em] FC1);
    \draw [->] (GRU2)  ->  ([xshift= 2em] FC2);
    \draw [->] (FC1)  ->  (fc-sed);
    \draw [->] (FC2)  ->  (fc-doa);

\end{tikzpicture}
\end{center}

\caption{The proposed network architecture.}
\label{fig:network}
\end{figure}

The proposed system is based on the SELDNet introduced by  Adavanne et al.~in \cite{adavanne2018sound}, with some modifications. In a similar manner, it features a CRNN network using Gated Recurrent Units (GRU) as recurrent layers. This is followed by two parallel branches of Fully Connected (FC) layers, one for SED and one for DOA estimation, sharing weights along time dimension. The first FC layer of both branches uses linear activation, while the last FC layer of each branch uses a different activation function according to the task. The last FC layer in the SED branch contains 14 nodes using sigmoid activation (one node for each sound event class to be detected), while the last FC layer in the DOA branch consists of 42 nodes using $\tanh$ activation (each of the sound event classes is represented by 3 nodes relative to the 3-dimensional sound event location). We use binary cross-entropy as loss function for the SED branch and mean square error (MSE) loss for DOA estimation branch, keeping the two branches separated.

Regarding the differences respect to the baseline in this implementation, firstly, we added 2 CNN blocks in order to help the network learn more features, increasing the number of CNN blocks from 3 to 5. Each CRNN block consists of a convolutional layer with rectified linear unit (ReLU) activation, batch normalization to normalize the activation output, and MaxPooling along frequency axis to reduce the dimensionality. Although adding layers to a neural network usually helps it to learn more features, it has the disadvantage of leading to possible overfitting, especially when the training dataset size is small as in this case. To prevent overfitting, we use Dropout in each convolutional block, after reducing the dimensionality.

Secondly, we used rectangular filters instead of squared ones, mainly inspired by Pons et al.~\cite{pons2016experimenting}. In the mentioned paper, the authors studied how filter shapes can help to proper model CNN motivated by musical aspects, achieving positive results in music classifications. Filters used in convolutional layers are principally inspired by image processing literature, typically being small and squared (3x3 or 5x5). Considering that, in the audio domain, filter dimensions corresponds to time and frequency axes, wider filters may be capable of learning longer temporal dependencies in the audio, while higher filters may be capable of learning more frequency features.
We propose the same concept, applying it to sound event detection. We used rectangular filters of shape ${1xM}$, being 1 the time dimension and the M the frequency dimension. We hypothesize that setting the time dimension to 1 and increasing the frequency dimension, taking into consideration nearly all the mel-bands, would lead the network to better model the frequency dimensions, helping the system to learn the presence or absence of an event, and consequently improving the metrics for location-aware detection.
Such horizontal filters increase the receptive field only on frequency dimension. Anyway, the temporal information is taken care by the recurrent GRU layers.

The last addition is the use of data augmentation as described in Section 2.3, expanding the number of DOA represented in it and consequently raising the scores related to classification-dependent localization metrics.

\subsection{Data augmentation}
\label{sec:data}
With the aim of improving the results and to reduce system overfitting, the training dataset size has been increased using data augmentation based on channel rotations and reflection on the \textit{xy} plane in the FOA domain.
Among others, this system implements the \textit{16 patterns} technique proposed for the first time by Mazzon et al.~in \cite{mazzon2019sound}, with some small changes. This approach allows to increase DOA combination and correctly compute the corresponding ground truth DOA
labels of the augmented data. Moreover, a relevant advantage of this method is the possibility to be applied regardless of the number of overlapping sound sources \cite{mazzon2019first}, which makes it an easy and straightforward data augmentation method.
We augmented the data following the transformations suggested in \cite{mazzon2019first}, considering only \textit{channel swapping} and \textit{channel sign inversion}. The suggested data manipulations correspond to rotations of 0, -90$^{\circ}$, +90$^{\circ}$, and +180$^{\circ}$ of the azimuth angle $\phi$ and its reflection with respect to the xz plane, leading to 8 rotations around the z axis, and a reflection with respect to the xy plane (considering the opposite elevation angle), for a total of 15 new patterns plus the original one.
Figure \ref{fig:intensity} shows an example of channel rotation on intensity vector, after applying a reflection with respect to xy plane. The figure shows the reflection of channel Z. The reader is referred to \cite{mazzon2019first} for further details.
In \cite{mazzon2019sound}, Mazzon et al.~compute the augmented dataset in time domain and extract the features offline for each of the augmented signals. All the possible transformations are computed offline, and the data generator randomly chooses one of them at each iteration.
In our system we implemented only 15 patterns, without considering the original one as data augmentation pattern. We also implemented the data augmentation offline, but, for memory reasons, instead of computing all the transformations for each audio file, we randomly select one out of the 15 patterns to augment the data during the feature extraction process. The pattern selected for a particular audio file is also used for augmenting the corresponding label. All the new generated files, together with the original ones, are used to train the network.

\begin{figure}[t]
\centering
\includegraphics[width=0.5\textwidth]{intensityvector.png}
\caption{Example of augmented intensity vector. Pattern applied: reflection with respect to xy plane. The image shows the Z-component of the intensity vector, before and after applying a reflection on the Z-channel. The horizontal axis represents time (expressed in frames) and the vertical axis mel-band coefficients.}
\label{fig:intensity}
\end{figure}

\subsection{Hyper-parameters}
All the dataset audio files are sampled at 24kHz. For the STFT, we used a 960 point Hanning window with a 50\% hop size. Considering that the temporal resolution of the label is 100 ms, we interpolated the sub-frames of 20ms as suggested in the baseline system. The number of mel-band filter is set to 64.

The training development set has been increased from 400 to 1200 files using the aforementioned augmentation technique.
With regard to the optimization technique, we used Adam method \cite{kingma2014adam}.

A sound event is considered to be active, and its respective DOA estimation considered, if the SED output exceeds a threshold of 0.5.

\section{Experiments}
\label{sec:experiments}
The dataset used for the evaluation of the system is the one provided for the DCASE 2020 Challenge Task 3: TAU-NIGENS Spatial Sound Events 2020 \cite{politis2020dataset}.
The dataset consists of sound samples belonging to 14 different sound classes, which have been convolved with real spatial room impulse responses (RIRs). The dataset is presented in two spatial recording formats, MIC and FOA, and it is divided between development and evaluation set. The results presented in this paper are based on the FOA development set of the dataset, which consists of 6 predefined splits. We followed the same specification given in the task description, using split 1 for testing, split 2 for evaluation and split 3-6 for training. Further information regarding the dataset can be found at \cite{politis2020dataset}.

The network predictions have been evaluated considering the joint nature of localization and detection, as proposed in \cite{Mesaros_2019_WASPAA}. In particular, $ER\textsubscript{20$^{\circ}$}$ and $F\textsubscript{20$^{\circ}$}$ are related to the SED task and they are location-dependent. A prediction is considered true positive only if it is under a distance threshold of {20$^{\circ}$} from the reference. $LE\textsubscript{CD}$ (localization error) and $LR\textsubscript{CD}$ (localization recall) are related to DOA estimation, being classification-dependent, which means they are going to be calculated only between sounds with the same label in each frame.
All metrics are computed in one-second non-overlapping frames. For more information about the evaluation metrics refer to \cite{Mesaros_2019_WASPAA}.

Several architecture configurations, filter dimensions and data augmentation techniques have been explored before reaching the network architecture described in Section \ref{sec:net}. It is possible to consider the full experiment as the combination of two minor sub-experiments. The first one has concentrated on finding the rectangular filter shape which gives the best performance. The second has been focused on considering different data augmentation methods to select the most appropriate technique.

During the first experiment, we explored different rectangular shapes in order to understand the filters size which gives the best performance and helps the network to properly learn frequency features and correctly detect events. To do so, the results of the baseline system, which uses rectangular filters of shape 3x3, have been compared with the results of the proposed system using rectangular filters of shape $TxM$ of size 1x46, 1x48, 2x48, 3x48, 1x50, 1x52, 1x54, 1x56, and completely dense filters of dimension 1x64 (being 64 the size of mel-band filter). We supposed that the use of rectangular filters would lead the network to better model the frequency features.
We also tested filters of shape $2x48$ and $3x48$ to investigate if increasing the time dimension would facilitate the network to learn temporal information.

After selecting the filter which best performs, different data augmentation techniques have been considered to increase the SELD score and to reduce overfitting. In this second step, the research has concentrated on three data augmentation techniques: time stretching, pitch shifting and channel rotations. Each technique has been independently explored.
Regarding time stretching, we considered a stretch factor  between 0.9 and 1.1. Regarding pitch shifting, we considered the sound in a range between -2 and 2 semitones, considering only integer values. For each augmentation technique, two values within the chosen ranges have been randomly selected for each audio file, increasing the dataset from 400 to 1200 files. More details regarding channel rotation are given in Section \ref{sec:data}.

Across all the experiments, the batch size has been set to 128, and the systems have been trained for 50 epochs at most. An early stopping strategy has been implemented, stopping the training if the validation loss does not improve for 20 epochs.

\section{Results and discussion}
\label{sec:res}

\begin{table}[!t]
 \begin{tabularx}{\columnwidth}{lccccc}
  \toprule
  Filter shape & ER\textsubscript{20$^{\circ}$} & F\textsubscript{20$^{\circ}$} & LE\textsubscript{CD} & LR\textsubscript{CD} & SELD \\
  \midrule
  Baseline (3x3) & 0.72  & 37.4\%    & 22.8$^{\circ}$    & 60.7\%    & 0.47  \\
  1x46          & 0.69  & 40.1\%    & 21.7$^{\circ}$    & 62.3\%    & 0.45   \\
  \textbf{1x48} & \textbf{0.69}  & \textbf{40.3\%}  & \textbf{20.9$^{\circ}$}    & \textbf{62.4\%}    & \textbf{0.45} \\
  1x50         & 0.70  & 40.4\%    & 21.1$^{\circ}$    & 61.1\%    & 0.45 \\
  1x52          & 0.70  & 39.9\%    & 21.6$^{\circ}$    & 61.4\%    & 0.45 \\
  1x54          & 0.72  & 37.1\%    & 21.7$^{\circ}$    & 59.7\%    & 0.47 \\
  1x56          & 0.72  & 38.2\%    & 21.2$^{\circ}$    & 59.4\%    & 0.47 \\
  1x64          & 0.72  & 38.1\%    & 23.1$^{\circ}$    & 61.4\%    & 0.46 \\
  \textbf{2x48}          & \textbf{0.70}  & \textbf{40.6\%}    & \textbf{20.8$^{\circ}$}   & \textbf{61.1\%}    & \textbf{0.45} \\
  3x48          & 0.75  & 36.0\%    & 23.1$^{\circ}$    & 58.6\%    & 0.48 \\
 \bottomrule
 \end{tabularx}
 \caption{Evaluation results on development set for first stage of the experiment using different filters shapes.}
 \label{tab:results_filter}
\end{table}

Tables \ref{tab:results_filter} and \ref{tab:results_data} report the evaluation results of the first and second stage of the study, respectively. In particular, Table \ref{tab:results_filter} shows the evaluation results for the development dataset on the testing split using different rectangular filters shapes. Table \ref{tab:results_data} details the results on the same dataset, using 1x48 filters with different data augmentation methods.
In both tables, the results are compared with the baseline system.

Table \ref{tab:results_filter} reports the results related to the first sub-experiment.
In order to select the filters shapes which give the best results, several sizes has been tested. As it is possible to observe, between all the filters of size $1xM$, the best results are reached with filter shapes of 1x48 and 1x50. Among those, almost all the tested rectangular shapes outperfom the baseline results. The only filter which does not improve the metrics is the 1x54, which is otherwise comparable to it. Between the two filter shapes which better perform, we selected the 1x48 to continue the experiment. In particular, the selection of the 1x48 filters instead of 1x50 is based on the fact that, while the difference on sound event detection metrics is negligible, localization-aware classification scores for the first shape are better than the second one.
The table also report the results of filters's size 2x48 and 3x48 (on the same table for space limitations), which have been tested (after selecting 1x48 filter shape) to explore their relation with temporal information. As it possible to observe, increasing the receptive field also in time dimension does not help improve the results. Comparing 2x48 and 1x48, the difference in terms of results is neglectable, while increasing the time dimension also increase the training time of the network and its size. For this reasons, we selected the 1x48 as the best performant and as the final filter size to follow to study with.
This shows that using rectangular filters increasing the receptive field only in frequency dimension help the network to better model frequency features while GRU layers take care of the temporal information. Moreover, it also highlights the importance of understanding the training datasets in order to properly design the model architecture for a determined task.

The second sub-experiment is based on the exploration of several data augmentation techniques with the aim to increase the score and prevent overfitting. We focused the study on three methods: time stretching, pitch shifting and channel rotation. Results are shown in Table \ref{tab:results_data}. As it possible to observe, channel rotation is the only method that outperforms the baseline results, substantially increasing location-aware detection scores. This might be explained by the fact that channel rotations maintains the physical relations between channels, making the manipulations of sound more realistic. Moreover, the noticeable improvement of SED-related scores ($ER$ and $F-score$) could be explained by the joint nature of localization and detection of the evaluation metrics. In fact, the results suggest that channel rotation augmentation technique increases the DOA information, helping the network to better localize sound sources. It is possible that, without the expansion of the dataset size, the same sound source would have been considered as active, but the system would have been wrong in localizing it, considering the prediction wrong. By augmenting the data and by giving the network more DOA examples for training, the algorithm's results improve in a substantial way.

\begin{table}[!t]
 \begin{tabularx}{\columnwidth}{l@{\hskip 0.2in}c@{\hskip 0.2in}c@{\hskip 0.2in}c@{\hskip 0.2in}c@{\hskip 0.2in}c}
  \toprule
  Method & ER\textsubscript{20$^{\circ}$} & F\textsubscript{20$^{\circ}$} & LE\textsubscript{CD} & LR\textsubscript{CD} & SELD  \\
  \midrule
  Baseline          & 0.72  & 37.4\%    & 22.8$^{\circ}$    & 60.7\%    & 0.47 \\
  TS  & 0.86  & 22.8\%    & 27.9$^{\circ}$    & 51.0\%    & 0.57 \\
  PS    & 0.86  & 22.3\%    & 28.7$^{\circ}$    & 51.0\%    & 0.57 \\
  \textbf{CR}  & \textbf{0.59}  & \textbf{50.6\%}    & \textbf{17.6$^{\circ}$}    & \textbf{66.2\%}    & \textbf{0.38} \\
  \bottomrule
 \end{tabularx}
 \caption{Evaluation results on development for the second stage of the experiment using different data augmentation techniques. TS stands for time stretching, PS stands for pitch shifting and CR stands for channel rotations. }
 \label{tab:results_data}
\end{table}

\section{Conclusions and Future Work}
\label{sec:conc}
This paper describes a system submitted for the DCASE Challenge 2020 Task 3. The method is based on SELDNet presented by Adavanne et. al~ \cite{adavanne2018sound}, with some differences. Data augmentation based on Ambisonic rotations, network architecture and training loss functions are the main changes.
The main contribution of the proposed system is the usage of rectangular filters, highlighting the importance of modeling the network architecture according to the dataset used for training. Data augmentation also helped to increase the evaluation score, especially ER\textsubscript{20$^{\circ}$} and F\textsubscript{20$^{\circ}$} related to the SED task.
The proposed system considerably outperforms the state-of-the-art method presented as baseline, significantly increasing the location-dependent metrics related to SED task.

Future work will include further investigation of different data augmentation techniques based on preserving physical relations between channels, such as random rotation matrices to rotate the acoustic scene of random angles, and explore rectangular filter shapes performances on different dataset formats.
\newpage

\bibliographystyle{IEEEtran}


\section{Introduction}
For more than a century mathematical methods have been developed to
describe the stochastic evolution of complex dynamical systems, with
current applications in physics, chemistry, biology,
engineering and finance~\cite{1996-Nature-Ghashghaie-Dodge,1997-PRE-Naert-Peinke,1999-PRL-Luck-Friedrich,2005-NJP-Hummer,2006-JCP-Ma-Dinner,2009-NJP-Bahraminasab-Friedrich,2010-PNAS-Best-Hummer,2013-JCTC-Comer-Gonzalez,2016-JCIM-Lee-Gumbart,2016-JCTC-Gaalswyk-Rowley}.
Of particular interest in toxicology and pharmacology is the prediction
of kinetic quantities such as transition rates, permeation coefficients,
and mean first-passage times of solutes, including drugs and small
peptides, which provide fundamental understanding of numerous
biochemical transport processes~\cite{2016-JCIM-Lee-Gumbart,2016-JCTC-Gaalswyk-Rowley,2018-JPCB-Badaoui-Rosta}.
As originally described by Hendrik A. Kramers in his seminal work~\cite{1940-Physica-Kramers}, such quantities can be calculated by
employing models of diffusive motion. In the simplest diffusive model,
one can ignore inertia and memory effects and regard the diffusive motion
as the random walk of a particle under a position-dependent potential~\cite{1990-RMP-Hanggi-Borkovec, 2018-PRE-Sicard}.
The dynamics of the reaction coordinate, $X(t)$, is
determined by two functions: the potential of mean force, $W(X)$,
and the diffusivity, $D(X)$, along this coordinate.
In general, $W(X)$ and $D(X)$ are both position-dependent
and are likely to vary substantially in heterogeneous systems.\\

A variety of numerical methods for calculating $W(X)$ based on enhanced
sampling techniques are now well-established within the computational
community. They employ biaising potentials deposited along a few selected
degrees of freedom, named reaction coordinates (RC), to overcome the limits
of Boltzmann sampling in the presence of rare events in which energy
barriers and transition states may be poorly sampled or not sampled at all~\cite{2013-JCP-Sicard-Senet,2015-JCP-Sicard-Manghi,2016-Springer-Tiwary-vandeWalle,2018-Langmuir-Sicard-Striolo,2018-PRE-Sicard,2019-ACSNano-Sicard-Striolo,2018-APX-Camilloni-Pietrucci,2020-PRE-Sicard-Manghi}.
The calculation of $D(X)$, on the other hand, has received less attention.
The Einstein-Smoluchowski relation, which relates the diffusivity to the
mean square deviation of the position of the solute in the one-dimensional
long-time limit, can be used to calculate the diffusion coefficient of a
solute in a homogeneous solution by analysis of a molecular dynamics (MD)
trajectory~\cite{1956-Book-Einstein,2005-NJP-Hummer}. This relationship,
though, offers a very poor approximation of the true diffusivity in
inhomogeneous systems such as a bilayer~\cite{2016-JCTC-Gaalswyk-Rowley}.
In these systems, the variation of the solute diffusivity is large because
the frictional environment varies dramatically as the solute moves from
bulk water through interface, and into the membrane interior, potentially
encountering free energy barriers with heights greater than $k_BT$.
For similar reasons, estimates based on a Green-Kubo relation of the
velocity are expected to be equally poor~\cite{1966-RPP-Kubo,2006-BC-Mamonov-Coalson}.

To circumvent this limitation, Marrink and Berendsen calculated the
diffusivity profile for the permeation of water using the force
autocorrelation function~\cite{1994-JPC-Marrink-Berendsen,1996-JPC-Marrink-Berendsen}.
This method requires the solute to be constrained to a point
along the RC through the modification of the equation of motion of the
MD integration, which makes it relatively difficult to apply. As a result,
it is far more convenient to perform simulations where the solute is
simply restrained to remain near a given position with a biasing potential.
Among the latter, two strategies have been commonly employed for
calculating $D(X)$ using biased MD simulations. The first is based on
the generalized Langevin equation for a harmonic oscillator~\cite{2016-JCTC-Gaalswyk-Rowley,2016-JCIM-Lee-Gumbart}.
The second employs Bayesian inferences on the likelihood of the observed
dynamics of the solute~\cite{2013-JCTC-Comer-Gonzalez,2016-JCIM-Lee-Gumbart}.

The generalized Langevin equation provides straightforward methods to
calculate position-dependent diffusion coefficients from restrained MD
simulations. The solute is restrained by a harmonic potential so that it
oscillates at a point along the RC. The solute is then described as a
harmonic oscillator undergoing Langevin dynamics where the remainder of
the system serves as the frictional bath for the solute. Implicitly,
describing the system as a harmonic oscillator requires the restraining
potential to be dominant over the underlying free energy surface,
\textit{i.e.} the latter is effectively a perturbation on the former.
Values of the spring constant sufficiently large to justify this
assumption also tend to be too large for umbrella sampling (US) simulations,
meaning that it may not be possible to use the same simulation to calculate
$W(X)$ and $D(X)$.
Once a time series of the $X$ position of the solute is collected, the
diffusion coefficient for that point can be calculated from the position
or velocity autocorrelation functions (PACF and VACF, respectively).
These methods were first introduced by Berne and co-workers~\cite{1988-JPC-Berne-Straub}
for the calculation of reaction rates and elaborated by Woolf and Roux
to calculate position-dependent diffusion coefficients~\cite{1994-JACS-Woolf-Roux}.
Hummer proposed a simpler method to calculate diffusion coefficients from harmonically
restrained simulations which avoid the need for multiple numerical
Laplace transforms of the VACF~\cite{2005-NJP-Hummer}. However methods
based on PACF are sensitive to the convergence of the correlation
functions, particularly for heterogeneous bilayer environments
where the PACF does not decay to values near zero at long time scale
due to the lack of ergodicity in sampling~\cite{2016-JCIM-Lee-Gumbart}.
Furthermore, these algorithms are only applicable to restrained MD simulations,
and unbiased simulations cannot be analysed with PACF or VACF.

A fundamentally different approach for the determination of
position-dependent diffusivities employs Bayesian inferences
to reconcile thermodynamics and kinetics. In this method originally
developed by Hummer~\cite{2005-NJP-Hummer} and elaborated thereafter
by Türkcan et al.~\cite{2012-BJ-Turkcan-Masson} and Comer
et al.~\cite{2013-JCTC-Comer-Gonzalez}, no assumptions are made regarding
the form of the free energy landscape on which diffusion occurs.
This approach estimates diffusion coefficient and free energies self-consistently.
The Bayesian scheme uses as parameters the values of the RC, $X$, along the trajectory,
together with the force, $F(X,t)$, which is the sum of the time-dependent bias
and the intrinsic system force, $-\nabla W(X)$. Under the stringent assumption of
a diffusive regime, the motion is propagated using a discretized Brownian integrator.
A Likelihood function can be constructed that gives the probability
of observing exactly the motions along $X$ seen in the simulation runs.
Using Bayes' formula and Metropolis Monte Carlo simulations,
the likelihood of observations given the parametrized diffusive model
is turned into a posterior density of the unknown parameters $F(X,t)$
and $D(X)$ of the diffusive model, given the simulation observations.
Additional prior knowledge and free parameters can be added to improve
the result accuracy. However, they can make it difficult for the Bayesian
scheme to find unique solutions, and greater sampling may be needed than
for simpler models~\cite{2005-NJP-Hummer,2013-JCTC-Comer-Gonzalez}.
Additionally, a crucial component of this scheme is the Brownian
integrator time step, which must be chosen to be much larger than any
correlation time of the motion~\cite{2016-JCIM-Lee-Gumbart}.\\

In the present study, we propose a general, automatic method for estimating
multi-dimensional position-dependent drift and diffusion coefficients
equally valid in biased or unbiased MD simulations. We combine
the dynamic histogram analysis method (DHAM)~\cite{2015-JCTC-Rosta-Hummer},
which uses a global Markov state model (MSM) based on discretized RCs
and a transition probability matrix, with a Kramers-Moyal (KM)
expansion to calculate position-dependent drift and diffusion coefficients
of stochastic processes, from the time series they generate, using their
probabilistic description~\cite{1966-book-Risken,2009-PRE-Anteneodo-Riera,2009-PRE-Petelczyc-Baranowski,2009-NJP-Bahraminasab-Friedrich,2016-JORS-Rinn-Peinke,2019-JOSS-Rydin-Meirinhos}.
DHAM has been successfully employed to compute stationary quantities
and long-time kinetics of molecular systems with significantly high
energy barriers, for which transition states can be poorly sampled or
not sampled at all with unbiased MD simulations~\cite{2018-JPCB-Badaoui-Rosta}.
In this context, MSMs are extremely popular because they can be used
to compute stationary quantities and long-time kinetics from ensembles
of short simulations~\cite{2015-JCTC-Rosta-Hummer,2017-JCTC-Stelzl-Hummer,2017-PNAS-Olsson-Noe,2018-JACS-Husic-Pande,2018-JPCB-Biswas-Stock}.

Here we use MSMs to model MD trajectories generated from biased
simulations. We apply the DHAM method to unbias their transition
probability matrices and use the KM expansion to assess the original
drift and diffusion coefficients.
We apply the formalism to one- and two-dimensional analytic potentials
and data from explicit solvent MD simulations, including
the water-mediated conformations of alanine dipeptide and the permeation of
the Domperidone drug molecule across a lipid membrane.
Our approach neither requires prior assumptions regarding the form of
the free energy landscape nor additional numerical integration scheme.
It can present significant improvement compared to standard methods
when long time scale fluctuations yield a lack of ergodicity in sampling,
as encountered in lipid bilayer systems studied in drug discovery.

\section{Methods}
\subsection{Definition of the Diffusion coefficient}
A wide range of dynamical systems can be described with a stochastic
differential equation, the non-linear Langevin equation~\cite{2009-NJP-Bahraminasab-Friedrich,2009-PRE-Anteneodo-Riera,2014-PRE-Farago-Gronbech,2016-JORS-Rinn-Peinke,2016-SR-Bodrova-Metzler,2016-PRE-Regev-Farago}.
Considering a one-dimensional stochastic trajectory $X(t)$ in time $t$,
the time derivative of the system's
trajectory $dX/dt$ can be expressed as the sum of two complementary contributions:
one being purely deterministic and another one being stochastic, governed by
a stochastic force.
For a stationary stochastic process, the deterministic term is defined
by a function $D^{(1)}(X)$ and the stochastic contribution is weighted
by another function, $D^{(2)}(X)$, which do not explicitly depend on
time, yielding the evolution of the equation of $X$,
\begin{equation}
\label{eqLgv}
\frac{dX}{dt} = D^{(1)}(X) + \sqrt{D^{(2)}(X)} ~\Gamma(t) \,,
\end{equation}
where $\Gamma(t)$ is a zero-average Gaussian white noise, \textit{i.e.}
$\langle \Gamma(t) \rangle = 0$ and $\langle \Gamma(t) \Gamma(t') \rangle
= 2\delta(t-t')$,  with $\delta$ the Dirac function.
According to the Ito's prescription, this is equivalent to the
Fokker-Planck equation~\cite{1966-book-Risken}
\begin{equation}
\frac{\partial P(X,t)}{\partial t}=\Big[-\frac{\partial }{\partial x} D^{(1)}(X) + \frac{\partial^2 }{\partial x^2} D^{(2)}(X) \Big]P(x,t)\,,
\label{eqFP}
\end{equation}
with stationary solution
\begin{equation}
P(X) \propto \frac{1}{D^{(2)}(X)} \exp \Bigg(\int_X dx \frac{D^{(1)}(x)}{D^{(2)}(x)} \Bigg)\,.
\label{solFP}
\end{equation}
Different integration schemes (other than the Ito's prescription) can
be envisaged to integrate Langevin equations with a position-dependent
diffusion coefficient, as in Eq.~\ref{eqLgv}, \textit{e.g.} the Stratonovich
convention. Such schemes lead to different drift coefficients (often referred to
as 'anomalous' drifts) in the Fokker-Planck equation, resulting in different
stationary distributions~\cite{2007-PRE-Lau-Lubensky,2007-PRE-Bussi-Parrinello,2014-PRE-Farago-Gronbech,2014-JSP-Farago-GronbechJensen,2016-PRE-Regev-Farago}. \\

To infer the underlying dynamics of a stationary Markovian stochastic process, we can
evaluate the drift and diffusion coefficients by retrieving the Kramers-Moyal
coefficients from the time series $X(t)$~\cite{1997-PRL-Friedrich-Peinke,1998-PLA-Siegert-Peinke,2011-PR-Friedrich-Tabar}.
The KM coefficients arise from the Taylor expansion of the master equation
describing the Markov process, as
\begin{equation}
\label{eqKM}
D^{(n)}(X) = \lim_{\tau \to 0} \frac{1}{n!~\tau} c^{(n)}(X,\tau)  \,,
\end{equation}
which can be seen as the derivative with respect to $\tau$ of the
$n$-th moment of the stationary conditional probability density function
$\mathit{p}\big(X',t+\tau \mid X,t \big) = \mathit{p}\big(X',\tau \mid X \big)$~\cite{2011-PRE-Honisch-Friedrich}
\begin{equation}
\label{eqConditionalMomentContinuous}
c^{(n)}(X,\tau) = \int dX' [X'-X]^n ~\mathit{p}\big(X',\tau \mid X \big) \,.
\end{equation}
Mathematically the drift and diffusion coefficients are defined as
the first two KM moments, \textit{i.e.} $n=1$ and $n=2$
in Eq.~\ref{eqKM}, respectively.
Under the ergodic hypothesis, the average over the microstates
defined in Eq.~\ref{eqConditionalMomentContinuous} can be equivalently
replaced by the average over time of the trajectory $X(t)$, defined as
\begin{equation}
\label{eqConditionalMomentTimeAverage}
c^{(n)}(X,\tau) = \Big\langle [X(t+\tau)-X(t)]^n \Big\rangle_{X(t)=X} \,.
\end{equation}
Additionally, the deviation from the stochastic Langevin description given
in Eq.~\ref{eqLgv}, \textit{i.e.}  the deviation of the driving noise
$\Gamma(t)$ from a Gaussian distribution, can be tested with the Pawula
theorem~\cite{1966-book-Risken}. To do so, one can compute the fourth-order
coefficient in the KM expansion,
$D^{(4)}(X) = \lim_{\tau \to 0} \frac{1}{4! \tau} c^{(4)}(X,\tau)$
and compare it to the diffusion coefficient, expecting
$D^{(4)}(X) \ll \big(D^{(2)}(X)\big)^2$ at all $X$.
In this work, we use MSMs to model stationary time series
resulting from trajectories generated with a position-dependent drift $D^{(1)}(X)$ and diffusion coefficient $D^{(2)}(X)$.

\subsection{Kramers-Moyal and Diffusion coefficients from MSM}
Typically, in molecular simulations we assume that the state space
of a system evolving stochastically in time can be discretised, and this discretised
data leads to a Markovian process. This usually involves the assumption that we have
a low-dimensional RC, X, which describes the thermodynamical and dynamical
evolution of the system.
In constructing MSMs, we discretise the RC, X, in $N_{\rm bin}$
bins, $\{x_i,\ldots, x_{N_{\rm bin}}\}$, defining the set of states
(also called microstates), which the system can occupy.
Each MD trajectory can then be analysed as a series of microstate assignments
rather than as a series of conformations. The dynamics of the system is regarded
as a memoryless process such that the next state of the system depends only
on its present state.
The number of transitions between each pair of state $i$
and state $j$ in the lagtime $\tau$ can then be counted
and stored as a transition count matrix, $\mathbf{C}_{ji}^{(\tau)}$,
which can be normalised to provide a numerical estimate of
the transition probability matrix $M_{ji}^{(\tau)}$, at lagtime $\tau$~\cite{2010-Methods-Pande-Bowman,2011-JCP-Prinz-Schutte,2014-Springer-Prinz-Noe,2015-JCTC-Rosta-Hummer}.
Subsequently, to test the Markov condition, a spectral decomposition can
be performed to write the Markov matrix $M_{ji}^{(\tau)}$ in terms of
its eigenvalues and eigenvectors which provide information about the dynamics
of the system,~\cite{2018-JCP-Kells-Rosta}
\begin{equation}
\label{eqMMspectral}
M_{ji}^{(\tau)} = \sum_n \psi_n^R(j) \psi_n^L(i) \lambda_n \,,
\end{equation}
where $\psi_n^R$ and $\psi_n^L$ are the right and left eigenvectors, respectively,
corresponding to eigenvalues $\lambda_n$. The latter are ordered such that
$1 = \vert\lambda_1\vert > \vert\lambda_2\vert \geq \dots \geq \vert\lambda_N\vert$.
The second largest eigenvalue (in magnitude) of the MSMs describes the slowest relaxation
process in the system. In practice, the slowest relaxation time,
$\tau^{\textrm{relax}} = -\tau/\ln \lambda_2$, determined from MSMs will have
a functional dependence on the lagtime at which the model is constructed,
as illustrated in Figs.~\ref{fig1}(c), (d), and (e). This poses the question of
what lagtime MSMs should be constructed at.
\begin{figure*}[t]
\includegraphics[width=0.9 \textwidth, angle=-0]{fig1.eps}
 \caption{
 (a) Analytical potentials $V_{\textrm{low}}$ and (b) $V_{\textrm{high}}$
 considered in the 1D Langevin simulations.
 The lagtime dependence of the slowest relaxation time, $\tau^{\textrm{relax}}$,
 in the overdamped Langevin dynamics is shown for (a) the unbiased simulations
 in $V_{\textrm{low}}$, (b) the biased simulations ($K=500$ kcal/mol)
 in $V_{\textrm{low}}$, and (c) the biased simulations ($K=500$ kcal/mol)
 in $V_{\textrm{high}}$ for different numbers of bins, $N_{\textrm{bin}}$.
 The limiting relaxation times,
 $\mu^{\textrm{relax}}$, obtained  from Eq.~\ref{eqRelax} are shown with
 dashed lines for comparison.
 We measured $\mu^{\textrm{relax}} = 5.8\times 10^3$, $5.1\times 10^3$,
 and $1.7\times 10^{10}$ in panels (c), (d), and (e), respectively.
 %
 Uncertainties, as represented in shaded area in panels (c), (d), and (e),
 were estimated from 10 independent runs.
}
\label{fig1}
\end{figure*}
Typically, we can assess the non-Markovian effects by carrying out the
Chapman-Kolmogorov test, verifying that $n$ times the lagtime of
a Markov process is equivalent to raising the transition matrix
to power $n$~\cite{Bowman2014AnIT}. The relaxation time of a Markovian system
must then be invariant under lagtime changes. The smallest lagtime $\tau$ at which
this condition is sufficiently fulfilled is called the Markov timescale,
$\tau_M$~\cite{2009-NJP-Bahraminasab-Friedrich,2011-PRE-Honisch-Friedrich}. \\

Given the one-dimensional time series of microstates $i\in \{1,\ldots, N_{\rm bin}\}$,
one can rewrite the $n$-th conditional moment in its discrete form
\begin{equation}
\label{eqConditionalMoment}
c^{(n)}(x_i,\tau) = \sum_{j=1}^{N_{\rm bin}} (x_j-x_i)^n ~M_{ji}^{(\tau)} \,,
\end{equation}
where $x_i$ is the (discretized) value of the RC in the center of bin $i$.
In many cases, for a given $X$, $c^{(n)}(X,\tau)$ depends linearly on $\tau$
in the range of $\tau$ for which the diffusive regime is satisfied. Consequently the drift
and diffusion coefficients are estimated solely by the quotient between the corresponding
conditional moment and the lagtime $\tau$ in this range.
Integrating Eq.~\ref{eqLgv} within the Taylor-Ito framework yields the stochastic Euler equation~\cite{2014-PRE-Farago-Gronbech,2016-PRE-Regev-Farago}
\begin{equation}
\label{eqEuler}
X(t+\tau) = X(t) + D^{(1)}(X) \tau + \sqrt{D^{(2)}(X) \tau} ~\eta(t) \,,
\end{equation}
where $\eta(t)$ is a Gaussian white noise with the same average and
correlation as $\Gamma(t)$.
Inserting Eq.~\ref{eqEuler} into Eq.~\ref{eqConditionalMomentTimeAverage}
 yields the relations between the conditional moments and the drift and
 diffusion coefficients~\cite{2016-JORS-Rinn-Peinke}
\begin{eqnarray}
c^{(1)}(X,\tau) & \approx & D^{(1)}(X) \tau \,, \label{eqSlope:1} \\
c^{(2)}(X,\tau) & \approx & 2D^{(2)}(X) \tau + \big(D^{(1)}(X) \tau\big)^2 \,. \label{eqSlope:2}
\end{eqnarray}
The latter can be subsequently assessed from the slope of a weighted
polynomial regression of Eqs.~\ref{eqSlope:1} and \ref{eqSlope:2}.\\

Similarly to the one-dimensional case discussed above, the two-dimensional case comprehends
two stochastic variables, $X(t)$ and $Y(t)$, governed by the stochastic differential equation~\cite{2009-NJP-Bahraminasab-Friedrich,2016-JORS-Rinn-Peinke}
\begin{eqnarray}
\label{eqLgv2D}
\frac{d}{dt}  \begin{bmatrix} X \\ Y \end{bmatrix} &=&
\begin{bmatrix} D_1^{(1)}(X,Y) \\ D_2^{(1)}(X,Y) \end{bmatrix} \\
&+& \begin{bmatrix}  g_{11}(X,Y) & g_{12}(X,Y) \\ g_{21}(X,Y) & g_{22}(X,Y) \end{bmatrix}
\begin{bmatrix} \Gamma_1(t) \\ \Gamma_2(t) \end{bmatrix}\,, \nonumber
\end{eqnarray}
where the drift coefficient $\textbf{D}^{(1)} = (D_1^{(1)}, D_2^{(1)})^T$ is
a two-dimensional vector and the diffusion coefficient is a $2 \times 2$ matrix given
by $\textbf{D}^{(2)} = \textbf{g}\textbf{g}^T$, with $\star ^T$ the transpose algebraic
operator. Similar to the one-dimensional case
the integration of Eq.~\ref{eqLgv2D} follows from a simple Euler scheme leading to
\begin{eqnarray}
\label{eqLgv2D-2}
\begin{bmatrix} X(t+\tau) \\ Y(t+\tau) \end{bmatrix} &=& \begin{bmatrix} X(t) \\ Y(t)\end{bmatrix} +
\tau \begin{bmatrix} D_1^{(1)}(X,Y) \\ D_2^{(1)}(X,Y) \end{bmatrix}  \\
&+& \sqrt{\tau} \begin{bmatrix}  g_{11}(X,Y) & g_{12}(X,Y) \\ g_{21}(X,Y) & g_{22}(X,Y)\nonumber \end{bmatrix}
\begin{bmatrix} \eta_1(t) \\ \eta_2(t) \end{bmatrix}\,.
\end{eqnarray}
The trajectory data can then be analyzed with 2D MSMs~\cite{2015-JCTC-Rosta-Hummer}
to construct the two-dimensional Markov transition probability matrix, $\mathbf{M}^{(\tau)}$,
at lagtime $\tau$ and the respective conditional moments yielding the expression
for the drift and diffusion coefficients
\begin{eqnarray}
c^{(1)}(x_i,y_i,\tau) &=&  \sum_{j=1}^{N_{\rm bin}} (x_j-x_i)~M_{ji}^{(\tau)}\,,
\label{eqDriftDiffusion2D:1}\\
c^{(2)}(x_i,y_i,\tau) &=&  \sum_{j=1}^{N_{\rm bin}} (x_j-x_i) (y_j-y_i) ~M_{ji}^{(\tau)}\,, \label{eqDriftDiffusion2D:2}
\end{eqnarray}
which now depend on the initial bin $i$ through the values of the two
discretised RCs, $x_i$ and $y_i$.
Extracting the drift and diffusion coefficients from MSMs and KM expansion
depends on two parameters.
The first parameter is the number of bins, $N_{\textrm{bin}}$, dividing
the RC, at which $D^{(1)}$ and $D^{(2)}$ are estimated. This integer should not
be too large that each bin does no longer include sufficient statistics
of the transition counts and also not too small that no dependence of the drift
and diffusion on the state variable can be observed~\cite{2016-JORS-Rinn-Peinke}.
The second parameter is the range of lagtimes, $\mathbf{L}_{\tau}$, used to
build the Markov transition matrix $\mathbf{M}^{(\tau)}$
and calculate the KM coefficients
for different $\tau$ values in Eq.~\ref{eqKM}. The conditional moments
in Eq.~\ref{eqConditionalMoment} are computed for each bin and for each lagtime.
For each bin, a linear fit is computed for all lagtimes in $\mathbf{L}_{\tau}$.
In practice, we have to carefully consider the range of lagtimes,
as too short lagtimes lead to non-Markovian effects, where the approximations
for MSMs can break down. To address the choice of lagtimes, we can also consider
a more extensive phase space in higher dimensions, with more finely discretized MSMs,
where non-Markovian effects are reduced.\\

The construction of MSMs from biased simulation data has not been traditionally
possible. Biased simulations modify the potential energy function of
the system of interest such that the system is, for example,
harmonically restrained to a given region of the energy landscape. This is
advantageous as it allows sampling of regions which might otherwise not be adequately
visited during the simulation time. However, the kinetic behavior observed is no longer
representative of the true system, and as such, this needs to be accounted for
when constructing the MSMs.
Several recent numerical methods have been proposed in the literature
to estimate unbiased MSMs from biased (e.g. umbrella sampling or replica exchange)
MD simulations~\cite{2015-JCTC-Rosta-Hummer,2016-PNAS-Wu-Noe,2017-JCTC-Stelzl-Hummer,2017-JCP-Donati-Keller,2020-COSB-Kieninger-Keller}.
We used here the dynamic histogram analysis
method (DHAM)~\cite{2015-JCTC-Rosta-Hummer}, which uses a maximum likelihood estimate of the MSM transition
probabilities $M_{ji}^{(\tau)}$ given the observed transition counts during
each biased trajectory. The unbiased estimate $M_{ji}^{(\tau)}$ can then be
inserted into Eq.~\ref{eqConditionalMoment} within the KM framework, which yields
the estimation of the drift and diffusion coefficients.

\section{Results}
In the following we give several illustrative examples of the approach
presented in the Methods section. We used analytical models both in one and
two dimensions within the Brownian overdamped or full inertial Langevin equations.
We also used explicit solvent simulations, including the water-mediated conformations
of alanine dipeptide and the permeation of the Domperidone drug molecule
across a lipid membrane.
\begin{figure}[t]
\includegraphics[width=0.99 \columnwidth, angle=-0]{fig2.eps}
 \caption{
 Evolution of $c^{(2)}(X,\tau)$ in the Brownian unbiased dynamics
 and $V_{\textrm{low}}$ for (a) three different position along the RC and a
 fixed $N_{\textrm{bin}}=400$ and (b) different values of
 $N_{\textrm{bin}}=400$, $800$, and $1200$ and a fixed position $X=0.30$.
 The value of the Markov timescale, $\tau_M \approx 0.04$, obtained from
 the analysis of the relaxation time $\tau^{\textrm{relax}}$ in Fig.~\ref{fig1}
 and $N_{\textrm{bin}}=400$ is indicated.
 The timescale $\tau_L$ indicating the departure from the linear regime of
 $c^{(2)}(X,\tau)$ is shown for comparaison. We measured
 $\tau_L \approx 10^{-2}$, $10^{-3}$, and $5 \times 10^{-4}$ for
 $N_{\textrm{bin}}=400$, $800$, and $1200$, respectively.
 }
\label{fig2}
\end{figure}
\subsection{1D Brownian overdamped Langevin equation}
\begin{figure*}[t]
\includegraphics[width=0.85 \textwidth, angle=-0]{fig3.eps}
 \caption{
 \textbf{Top panels.} Brownian dynamics simulations in $V_{\textrm{low}}$
 and quadratic diffusion profile. The drift (a) and diffusion (c) profiles
 obtained in unbiased simulations are shown along with the diffusion
 profile obtained in US simulations with biasing spring constant
 $K=500$ kcal/mol (e).
 %
 \textbf{Bottom panels.} Brownian dynamics simulations in $V_{\textrm{low}}$
 and steplike diffusion profile. The drift (b) and diffusion (d) profiles
 obtained in unbiased simulations are shown along with the diffusion
 profile obtained in US simulations with biasing spring constant
 $K=500$ kcal/mol (f).
 %
 The diffusion coefficients, $D^{(2)}(X)$, are rescaled to their respective
 minimum values defined in the main text, $D_0 = k_BT/\gamma_0^P$
 or $D_0 = k_BT/\gamma_0^Z$, for clarity.
 %
 The profiles obtained within the DHAM framework (red) agree within standard error
 with the exact profiles obtained from the generalization of Einstein's relation (black).
 %
 Uncertainties, as represented in shaded area, were estimated from 10 independent runs.
}
\label{fig3}
\end{figure*}
As a first example we integrated the 1D Brownian overdamped Langevin equation
\begin{equation}
\frac{dx(t)}{dt} = \frac{F(x(t))}{\gamma(x(t))} + \sqrt{\frac{ k_B T}{\gamma(x(t))}}\eta(t) \,,
\label{eqBrownian}
\end{equation}
with $k_B$ the Boltzmann constant and $T=300~K$ the temperature of the system.
In Eq.~\ref{eqBrownian}, the mass of the system is conveniently set to unity and
$F(x)=-\nabla V(x)$ is the deterministic force derived from the 1D potential
energy $V(x) = V_{\textrm{ref}}(x) + V_{\textrm{bias}}(x)$.
In the following, $V_{\textrm{ref}} = \sum_{n=0}^{6} \alpha(n) x^n$
is defined as a polynomial of degree $6$.
We considered two different choices of the coefficients $\alpha$,
as detailed in the supporting information (SI), leading to
the high ($V_{\textrm{high}}$) and low ($V_{\textrm{low}}$) barrier
potentials plotted in Fig.~\ref{fig1}(a) and (b),
respectively.
The biased potential is defined as
$V_{\textrm{bias}} = \frac{1}{2} K \big(x-x^{(k)}\big)^2$
with $K$ the biasing spring constant and $x^{(k)}$ the center of the harmonic
bias in simulation $k$.
The parameter $\gamma(x)$ represents the position-dependent friction
coefficient with $D(x) = k_BT/\gamma(x)$ the natural generalization
of Einstein's relation defining the position-dependent diffusion coefficient~\cite{2014-JSP-Farago-GronbechJensen}.
We used the It\^o convention to obtain the first order
integrator of the overdamped Langevin equation~\cite{2014-PRE-Farago-Gronbech,2016-PRE-Regev-Farago}.
The total number of timesteps $N_{\textrm{step}}$ in the unbiased or biased
simulation runs was chosen such that the simulation time $\Delta t \times N_{\textrm{step}} =
5 \times 10^4$ was kept constant.
The simulation timestep was set to $\Delta t = 10^{-5}$, which is at least an order
of magnitude smaller than the slow characteristic time scale for the diffusion
in the system, $\tau_{\textrm{diff}}=1/\max_x\big(\gamma(x)\big)$.\\

\textbf{Unbiased simulation.} We studied the evolution of the system
under Eq.~\ref{eqBrownian} with the low barrier analytical potential
$(V_{\textrm{low}})$ shown in Fig.~\ref{fig1}(a) and $V_{\textrm{bias}}=0$.
We considered either a quadratic or step-like position-dependent diffusion
profile with $\gamma(x)$ a parabolic (P) or a $Z$-shaped
membership (Z) function~\cite{1998-IJAR-Medasani-Krishnapuram},
\begin{eqnarray}
    \gamma^{P}(x) &=& \gamma_0^P \Big( 1-\frac{1}{3}(x-x_P)^2 \Big) \,, \\
    \gamma^{Z}(x) &=& \gamma_0^Z \Big( 2+\textrm{zmf}(x,a,b) \Big) \,,
\end{eqnarray}
with $\gamma_0^P = 3000$, $x_P=0.8$, $\gamma_0^Z=1700$, $a=0.5$, $b=1.1$,
and $\textrm{zmf}(x,a,b)$ the sigmoidal membership function
(see details in the SI).
We first assessed the Markov timescale $\tau_M$ associated with the trajectories
generated via Eq.~\ref{eqBrownian} for different binning $N_{\textrm{bin}}=100$,
$200$, and $400$, via the Chapman-Kolmogorov test. Then, we inferred the relaxation
time of the system, in the limit of very long lagtimes, using the fitting procedure
introduced in previous work~\cite{2018-JCP-Kells-Rosta}, which describes
the relaxation time, $\tau^{\textrm{relax}}$, as
\begin{equation}
\label{eqRelax}
\tau^{\textrm{relax}} = \frac{\tau \times \mu^{\textrm{relax}}}{\tau + \epsilon \mu^{\textrm{relax}}} \,.
\end{equation}
In Eq.~\ref{eqRelax}, $\tau$ is the lagtime, and $\mu^{\textrm{relax}}$
and $\epsilon$ are two free parameters which represent the true (limiting)
relaxation timescale and the initial rate of change of the effective relaxation
time, respectively.
As shown in Fig.~\ref{fig1}(c), increasing moderately the number of bins from $100$
to $400$ increases the rate of convergence towards the true relaxation time,
$\mu^{\textrm{relax}}=5.8\times 10^3$,
obtained from Eq.~\ref{eqRelax} with $N_{\textrm{bin}}=400$.
We obtained similar values for the true relaxation times measured
with $N_{\textrm{bin}}=200$ and $100$ (data not shown).
This yields a smaller Markov timescale, $\tau_M$, which, in turn,
gives a more accurate measure of the drift and diffusion coefficients,
in line with the variational principle satisfied by the MSMs~\cite{2019-JCP-Kells-Rosta}.
Increasing further $N_{\textrm{bin}}$, however, each bin would eventually no longer
include the requisite sufficient statistics for the MSM analysis.
Considering $N_{\textrm{bin}}=400$ in Fig.~\ref{fig1}(c), the relaxation time
can be seen to level off in the region of lagtimes greater than $0.2$.
In the analysis that follows, we chose to define $\tau_M \approx 0.4$,
as it is sufficiently large to be insensitive to the precise choice
of the lagtime.
Subsequently, the first and second conditional moments were evaluated
using Eq.~\ref{eqKM} over the range of lagtime $\tau \geq \tau_M$.

In Fig.~\ref{fig2}(a) and (b), the representative evolution
of $c^{(2)}(X,\tau)$ is shown for different positions along the RC
and different number of bins, respectively.
At relatively long lagtime, for $\tau \geq \tau_M$, $c^{(2)}$
follows the linear trend expected from Eq.~\ref{eqSlope:2}.
Close analysis of the evolution of $c^{(2)}$ shows that the linear trend
is satisfied above a lagtime threshold, $\tau_L$, significantly lower
than $\tau_M$.
At sufficiently short lagtime, on the other hand, for $\tau < \tau_L$,
$c^{(2)}$ deviates from the linear trend, which stems from non-Markovian
effects coming from the system discretization.
As shown in Fig.~\ref{fig2}(b), increasing the number of bins from
$N_{\textrm{bin}}=400$ to $1200$ decreases the value of the lagtime threshold,
$\tau_L$, above which the linear trend obtained in Eq.~\ref{eqSlope:2} is satisfied.

Drift and diffusion coefficients were subsequently assessed from the slope
of a weighted polynomial regression following Eqs.~\ref{eqSlope:1}
and \ref{eqSlope:2}.
We also calculated the fourth-order coefficient and evaluated the ratio
$D^{(4)}(X)/\big(D^{(2)}(X)\big)^2 < 5\times 10^{-3}$,
indicating the validity of the condition of the Pawula theorem.
As shown in Fig.~\ref{fig3}(a)-(d),
we observed excellent agreement between the theoretical profiles and
the numerical results for both position-dependent drift and diffusion
coefficients, with higher variability around $X \approx 0.5$
and $\approx 1.2$, where the derivative of the low barrier potential
$V_{\textrm{low}}$, shown in Fig.~\ref{fig1}(a), is maximal.\\

\textbf{Biased simulation.} We extented the previous analysis to the biased evolution
of the system in the low $(V_{\textrm{low}})$ and high $(V_{\textrm{high}})$
barrier potentials shown in Figs.~\ref{fig1}(a) and (b), respectively,
within the US framework.
We ran standard US simulations with the biasing potential
$V_{\textrm{bias}}$ defined above in each umbrella window.
We used $50$ uniformly distributed umbrella windows in the range $[0.25, 1.35]$
along the RC. This number was sufficiently high to obtain accurate sampling,
given the different biasing spring constants considered in this work.
\begin{figure}[t]
\includegraphics[width=0.99 \columnwidth, angle=-0]{fig4.eps}
 \caption{Drift (a) and diffusion (b) profiles obtained in biased
 overdamped Langevin simulations for the high energy barrier potential
 ($V_{\textrm{high}}$). A biasing spring constant $K=500$ kcal/mol and a quadratic
 diffusion profile are used. The diffusion coefficient $D^{(2)}(X)$ is rescaled to
 $D_0 = k_BT/\gamma_0^P$ for clarity.
 %
 The profiles obtained within the DHAM framework (red) agree within standard error
 with the exact profiles obtained from the generalization of Einstein's relation (black).
 %
 Uncertainties, as represented in shaded area, were estimated from 10 independent runs.
}
\label{fig4}
\end{figure}

We first studied the effect of the bias on the reconstruction of the diffusion
coefficient when the system evolves in $V_{\textrm{low}}$.
We set the biasing spring constant to $K=500$ kcal/mol, which was sufficienlty
strong to allow good sampling in both low ($V_{\textrm{low}}$)
and high ($V_{\textrm{high}}$) barrier potentials.
The lagtime dependence of $\tau^{\textrm{relax}}$ is shown in Fig.~\ref{fig1}(d).
We observed a similar behaviour
to the one obtained in the unbiased simulations (Fig.~\ref{fig1}(c)) with the decline of
the non-Markovian effects on the dynamics when the  number of bins,
$N_{\textrm{bin}}$ increases from $100$ to $400$. Most noticeably,
the use of the biasing spring constant decreases significantly
the variability of the relaxation timescale at longer lagtimes.
As shown in Figs.~\ref{fig3}(e) and (f), excellent agreement is observed between the theoretical
profiles and the numerical results.
We also confirmed that $D^{(4)}(X)/\big(D^{(2)}(X)\big)^2 < 5\times 10^{-3}$.
As expected, the use of the biasing spring constant yields better sampling
across the RC.

We complemented this analysis with the reconstruction of the diffusion
coefficient of the system using $V_{\textrm{high}}$ and a biaising
spring constant $K=500$ kcal/mol.
The relaxation timescale is shown in Fig.~\ref{fig1}(e).
We observed a similar behavior to the  unbiased simulations
with the decline of the non-Markovian effects on the dynamics when the
number of bins, $N_{\textrm{bin}}$ increases from $100$ to $400$.
The  relaxation time of the system calculated for $N_{\textrm{bin}}=400$
converges towards a limiting value, $\mu^{\textrm{relax}}=1.7\times 10^{10}$,
significantly higher than the one measured for $V_{\textrm{low}}$,
due to the longer equilibration time needed to cross the high energy barrier.
As shown in Figs.~\ref{fig4}(a) and (b), we observed good agreement
between the theoretical profiles and the numerical results for
the drift and diffusion coefficients.
We also verified that $D^{(4)}(X)/\big(D^{(2)}(X)\big)^2 < 5\times 10^{-3}$.
We noticed, however, higher variability
around $X \approx 0.5$ and $\approx 1.2$, where the derivative of
the high barrier potential $V_{\textrm{high}}$ (Fig.~\ref{fig1}(b)) is maximal.
As shown in the SI, the increase of the biasing spring constant
from $500$ kcal/mol to $2000$ kcal/mol yields
better sampling around the transition states with lower variability on
the reconstruction of the diffusion profile, as already observed
in the low barrier case.

\subsection{1D full inertial Langevin equation}
\begin{figure}[t]
\includegraphics[width=0.99 \columnwidth, angle=-0]{fig5a.eps}
\includegraphics[width=0.99 \columnwidth, angle=-0]{fig5b.eps}
 \caption{
Relaxation times $\tau^{\textrm{relax}}$ as a function of the lagtime, $\tau$,
and the number of bins, $N_{\textrm{bin}}$, obtained in biased 1D full Langevin
simulations for (a) low barrier ($V_{\textrm{low}}$) and (b) high barrier
($V_{\textrm{high}}$) energy  potentials. A biasing spring
constant $K=500$ kcal/mol is used. Drift (c) and diffusion (d) obtained
for US simulations with biasing spring constant $K=500$ kcal/mol and quadratic
diffusion profile. The diffusion coefficient $D^{(2)}(X)$
is rescaled to $D_0 = k_BT/\gamma_0^P$ for clarity.
The profiles obtained within the DHAM framework (red) in panels (c) and (d)
agree within standard error with the exact profiles obtained from the generalization
of Einstein's relation (black).
Uncertainties, as represented in shaded area, were estimated from 10 independent
runs.
}
\label{fig5}
\end{figure}
To take into account the role played by inertia in the reconstruction
of the drift and diffusion coefficients, we extended Eq.~\ref{eqBrownian}
to the full inertial Langevin equation
\begin{equation}
\frac{d^2x(t)}{dt^2} = F(x(t)) - \gamma(x(t))\frac{dx(t)}{dt} + \sqrt{k_B T \gamma(x(t))}\eta(t) \,,
\label{eqFull}
\end{equation}
where the mass of the system is conveniently set to unity.
From this relation, we used the Vanden-Eijnden and Ciccotti algorithm~\cite{2006-CPL-VandenEijnden-Ciccotti} that generalises the Velocity Verlet
integrator to Langevin dynamics, along with the It\^o convention. This scheme takes
the inertial term into account and is accurate to order $\Delta t^2$.

We used the same simulation timestep, $\Delta t = 10^{-5}$,
and biasing spring constants, $K=500$ kcal/mol,  as in the Brownian overdamped
Langevin dynamics.
The lagtime dependence of the relaxation time $\tau^{\textrm{relax}}$,
as shown in Fig.~\ref{fig5}(a) and (b) for $V_{\textrm{low}}$ and $V_{\textrm{high}}$,
respectively, is similar to the one measured in
the Brownian overdamped Langevin simulation, with the decline of the
non-Markovian effects when the  number of bins, $N_{\textrm{bin}}$ increases
from $100$ to $400$, and the convergence towards the limiting relaxation times
$\mu^{\textrm{relax}}=4.4\times 10^3$ and $1.3\times 10^{10}$ for
$V_{\textrm{low}}$ and $V_{\textrm{high}}$, respectively, calculated
for $N_{\textrm{bin}}=400$.

We limited the rest of the analysis to the biased evolution of the system
in $V_{\textrm{high}}$.
The associated drift and diffusion profiles are shown in Fig.~\ref{fig5}(c)
and (d) (see details in the SI). We observed good agreement between theoretical
and numerical profiles, with $D^{(4)}(X)/\big(D^{(2)}(X)\big)^2 < 4\times 10^{-3}$.
As already noted in the overdamped Langevin dynamics, we observed higher variability
around $X \approx 0.5$ and $\approx 1.2$, where the derivative of the potential energy
$V_{\textrm{high}}$ is maximal. Increasing the biasing spring constant
from $K=500$ kcal/mol to $2000$ kcal/mol improved the accuracy of the measure
of the diffusion coefficient (Fig.~S2).

\subsection{2D Brownian overdamped Langevin equation}
We extended the 1D analysis above to the estimation
of the 2D diffusion tensor of the analytical system evolving under
the 2D Brownian overdamped equation~\cite{2000-book-Freund-Poschel}
\begin{equation}
\gamma(t)\frac{d\mathbf{X}}{dt}=\mathbf{F}(\mathbf{X}) + \sqrt{k_BT}~\mathbf{b}(t)~\eta(t)
\label{eq2DLgv1}
\end{equation}
where $\mathbf{X}\equiv \big(X,Y\big)$ is a 2D vector and
$\textbf{F}(\mathbf{X}) = -\nabla V(\mathbf{X})$ is the deterministic force
derived from the 2D potential energy $V(\mathbf{X})$.
In Eq.~\ref{eq2DLgv1}, the force $\textbf{b}(t)$ is related to the friction
tensor $\mathbf{\gamma}(t)$ as follows
\begin{equation}
\mathbf{\gamma}(t) = \textbf{b}(t)\textbf{b}^T(t)\,,
\end{equation}
and $\eta(t)=\big(\eta_1(t), \eta_2(t)\big)$ is a 2D vector of independent,
identical Gaussian white noise sources, \textit{i.e.} $\langle \eta_i(t) \rangle = 0$
and $\langle \eta_i(t)\eta_j(t') \rangle = 2\delta_{i,j}\delta(t-t')$.
From this relation, we used the It\^o convention to obtain the first order
integrator of Eq.~\ref{eq2DLgv1}.
For economy of computational resources, we focused our analysis
on a diagonal friction tensor
\begin{equation}
\gamma(t) = \begin{bmatrix}  \gamma_{11}(\mathbf{X}) & 0 \\ 0 & \gamma_{22}(\mathbf{X}) \end{bmatrix}
\end{equation}
with constant friction $\gamma_{11}(\mathbf{X}) = 300$ and
$\gamma_{2}(\mathbf{X}) = 30$ and the 2D analytical potential
$V(\mathbf{X}) = -3 X^2 + X^4 -3X Y + Y^4$ shown in Fig.~\ref{fig8}(a).
We modeled standard US simulations, evenly positioned along the RC,
with a 1D bias potential
$V_{\textrm{bias}}(\mathbf{X}) = V^{(k)}(X) = \frac{1}{2} K (X-X^{(k)})^2$
in each umbrella
window $k$, with $K$ the biasing spring constant and $X^{(k)}$ the center of
the harmonic bias in window $k$. We used $50$ uniformly distributed umbrella windows
in the range $[0.25, 1.35]$ along the $x$-axis.
We then constructed a discretized 2D grid
to determine the MSMs along the $x$- and $y$-axes.
\begin{figure}[t]
\includegraphics[width=0.99 \columnwidth, angle=-0]{fig8.eps}
 \caption{ (a) Free-energy surface associated with the 2D analytical
 model potential projected along the $x$- and $y$-axes.
 The 2D projections of the deviation of the numerical diagonal elements
 of the diffusion tensor, $\Delta D_{11}$ and $\Delta D_{22}$, are shown in
 panels (b) and (c), respectively.
 }
\label{fig8}
\end{figure}

The diagonal elements of the diffusion tensor were obtained from the analysis
of the KM coefficients defined
in Eqs.~\ref{eqDriftDiffusion2D:1} and \ref{eqDriftDiffusion2D:2}
(see details in the SI).
We also calculated the fourth-order coefficients and verify that the ratio
$D_{11}^{(4)}(X)/\big(D_{11}^{(2)}(X)\big)^2 < 5\times 10^{-3}$ and
$D_{22}^{(4)}(X)/\big(D_{22}^{(2)}(X)\big)^2 < 6\times 10^{-4}$,
indicating the validity of the condition of the Pawula theorem.
To assess the accuracy of the reconstructed diffusion coefficient,
we measured the deviation between the observed (ob) and  theoretical (th)
values
\begin{equation}
    \Delta D_{ii}(\mathbf{X}) = \frac{\big\vert D_{ii}^{(\textrm{ob})}(\mathbf{X}) - D_{ii}^{(\textrm{th})}(\mathbf{X})\big\vert}{D_{ii}^{(\textrm{th})}(\mathbf{X})} \,,
\end{equation}
with $D_{ii}^{(\textrm{th})}(\mathbf{X}) = k_BT/\gamma_{ii}(\mathbf{X})$.
In Fig.~\ref{fig8}(b) and (c), it is shown the numerical measures
for $\Delta D_{11}(\mathbf{X})$ and $\Delta D_{22}(\mathbf{X})$, respectively.
The observations are in good agreement with the theoretical values
$D_{11}^{\textrm{th}}(\mathbf{X}) = 2\times 10^{-3}$ and $D_{22}^{\textrm{th}}(\mathbf{X}) = 2\times 10^{-2}$.
Most noticeably, the lower value of the friction parameter $\gamma_{22}(\mathbf{X})$
allows better sampling, which yields a more accurate measure for $D_{22}(\mathbf{X})$.
Finally, to quantify the accuracy of the reconstruction, we measured the effective
diffusion coefficient along each reaction coordinate
\begin{equation}
D_{ii}^{\textrm{eff}} = \frac{1}{N_{\textrm{bin}}}\int D_{ii}(\mathbf{X})~d\mathbf{X} \,.
\label{Deff2D}
\end{equation}
We measured $D_{11}^{\textrm{eff}}=(1.9 \pm 0.4)\times 10^{-3}$
and $D_{22}^{\textrm{eff}}=(2.0 \pm 0.7)\times 10^{-2}$,
in good agreement with $D_{11}^{\textrm{th}}$ and $D_{22}^{\textrm{th}}$,
respectively.
\begin{figure}[b]
\includegraphics[width=0.99 \columnwidth, angle=-0]{fig6.eps}
 \caption{(a) Ball-and-stick representation of the conformers
 $\alpha$ and $\beta$ of Ala2 along with (b) the free energy
 profile reconstructed along the dihedral angle $\Psi$ $(N1-C\alpha-C2-N2)$ using the DHAM
 in unbiased (red) and biased (blue) simulations. The diffusion profiles obtained in
 (c) unbiased and (d) biased simulations agree within standard errors
 with the one obtained within the PACF method of Hummer (grey)~\cite{2005-NJP-Hummer}.
 Uncertainties, as represented in shaded areas, were estimated from 3 independent
 runs.
}
\label{fig6}
\end{figure}

\subsection{Water-mediated conformations of Alanine Dipeptide in 1D}
The conformational transition between the different conformers of the
solvated alanine dipeptide (Ala2) has been extensively used as a case study
for several theoretical and computational investigations~\cite{2005-JCP-Ren-Weinan, 2010-JPCB-Vymetal-Vondrasek,2011-JCTC-Yonezawa-Nakamura,
2013-PRL-Tiwary-Parrinello,2017-JMM-Cuny-Mineva,2018-PRE-Sicard}.
We studied the transition between the metastable states $\alpha$ and $\beta$ of
Ala2 shown in Fig.~\ref{fig6}(a), which can be differenciated by the values of the
backbone dihedral angle $\Psi$ and are separated by an activation free energy barrier
of $\approx 2$ kcal/mol at the temperature $T=300~K$.
We used a Langevin thermostat to enforce the temperature~\cite{2007-JCP-Bussi-Parrinello},
a time step of $0.2$ fs, AMBER03 force field~\cite{2005-JCC-Case-Woods} with TIP3P water
model~\cite{1983-JCP-Jorgensen-Klein}, and GROMACS 5.1 molecular dynamics
code~\cite{2001-JMM-Lindahl-VanDerSpoel}. A single alanine dipeptide molecule was kept
in a solvated periodic cubic box of edge $\approx 3$ nm. The LINear Constraint Solver (LINCS)
algorithm~\cite{1997-JCC-Hess-Fraaije} handled bond constraints while the particle-mesh Ewald
scheme~\cite{1993-JCP-Darden-Pedersen} was used to treat long-range
electrostatic interactions. The non-bonded van der Waals cutoff radius was $0.8$ nm.

The relatively low energy barriers allow the system to be sampled with both unbiased
and biased simulations, which gives a mean to assess the accuracy of the results
obtained with enhanced sampling methods.
We determined the free energy profile and diffusion coefficient by using either
free simulations ($500$ ns) or US biased simulations and the DHAM approach.
The starting structures for the US simulations were obtained
by pulling the system along the dihedral angle $\Psi$ within the range $[-\pi,\pi]$. During the
simulations, a snapshot was saved every $0.1$ rad generating $60$ windows. Each US window was
subsequently run for $1$ ns to allow equilibration, followed by additional $5$ ns of the production
run using an US force constant of $5$ kcal~$\textrm{mol}^{-1}$~$\textrm{rad}^{-2}$.
In Fig.~\ref{fig6}(b), (c), and (d) are shown the reconstructed FEPs
and diffusion coefficients obtained within the DHAM framework,
either with unbiased (red) or biased (blue) simulations.
Both approaches gave similar results.
In particular, we compared the diffusion profiles obtained with DHAM
with the one obtained within the PACF framework (grey),
\begin{equation}
D(\Psi_k=\langle \Psi \rangle_k) = \frac{\textrm{var}(\Psi)^2}{\int_0^\infty C_{\Psi}(t)~dt}~.
\label{Diffeq}
\end{equation}
In Eq.~\ref{Diffeq}, $\langle \Psi \rangle_k$ is the average of the RC in the US window $k$,
$\textrm{var}(\Psi)=\langle \Psi^2 \rangle - \langle \Psi \rangle^2$ is its variance, and
$C_{\Psi}(t)=\langle \delta \Psi(0) \delta \Psi(t) \rangle$ the PACF calculated directly
from the time series.
Following the original work of Hummer, we increased the strength of the biasing
potential from $\approx 5$ kcal~$\textrm{mol}^{-1}$~$\textrm{rad}^{-2}$
to $\approx 24$ kcal~$\textrm{mol}^{-1}$~$\textrm{rad}^{-2}$ to satisfy
the underlying assumption that the harmonic restraint be sufficiently large to render
the underlying free energy surface a small perturbation on the harmonic
potential~\cite{2005-NJP-Hummer}.
As shown in Fig.~\ref{fig6}(c) and (d), we observed good agreement between the results
derived with DHAM either in the unbiased or biased simulations and
the PACF method. To compare the results with those in the literature,
we measured the effective diffusion coefficient
\begin{equation}
D^{(2)}_{\textrm{eff}} = \frac{1}{N_{\textrm{bin}}}\int_{-\pi}^{+\pi} D^{(2)}(\Psi)~d\Psi\,.
\label{D2eff}
\end{equation}
with $N_{\textrm{bin}}$ the number of bins used in the MSM.
Eq.~\ref{D2eff} yields the effective Diffusion coefficient
$D^{(2)}_{\textrm{eff}} \approx 0.16~\textrm{rad}^2/\textrm{ps}$,
in  agreement with the result obtained by Hummer et al.~\cite{2003-JCP-Hummer-Kevrekidis,2005-NJP-Hummer}
($\approx 0.15~\textrm{rad}^2/\textrm{ps}$) and Ma et al.~\cite{2006-JCP-Ma-Dinner}
($\approx 0.34~\textrm{rad}^2/\textrm{ps}$).

\subsection{Membrane permeation in 1D}
In this section we applied our method to study the permeation of
Domperidone, a dopamine receptor antagonist, across a POPC lipid membrane~\cite{2014-JCR-Eyer-Kramer,2017-JACS-Dickson-Duca,2018-JPCB-Badaoui-Rosta},
as illustrated in Fig.~\ref{fig7}(a).
The simulation data were obtained from previous studies~\cite{2017-JACS-Dickson-Duca,2018-JPCB-Badaoui-Rosta}
and further details of the simulation methods and parameters can be found there.
This system shows long time scale fluctuations inherent to the presence
of the bilayer structure, which is representative of the transport of solute
across heterogeneous media~\cite{2017-JACS-Dickson-Duca}.
As described in the original work of Dickson et al.~\cite{2017-JACS-Dickson-Duca},
the starting structures for the US simulations were obtained by placing the ligand
at the center of a POPC bilayer surrounded by water molecules ($72$ POPC and $60$
waters per lipid). To generate the US windows, the drug was pulled out from the center
of the system to outside the membrane, for a total of $40~\textrm{\AA}$.
Configurations were saved every $1~\textrm{\AA}$, from the center $z=0~\textrm{\AA}$
to $z=40~\textrm{\AA}$ generating $40$ windows.
Each US window was run for $20$ ns to allow equilibration, followed by additional $80$ ns
of the production run using an US force constant of
$2.5~\textrm{kcal}~\textrm{mol}^{-1}\textrm{\AA}^{-2}$.
\begin{figure}[b]
\includegraphics[width=0.99 \columnwidth, angle=-0]{fig7.eps}
 \caption{
 (a) Schematic representation of Domperidone (in brown at the center
 of the image) at the center of a POPC bilayer surrounded by water. The local area
 centered around $Z=10 \AA$, considered in the 2D calculation of the diffusion tensor,
 is shown (solid black rectangle).
 (b) Free energy and (c) diffusion profiles (in semi-log representation) associated
 with the crossing of the drug molecule through the bilayer and reconstructed
 along the distance from the center of the bilayer, $Z$.
 %
Uncertainties, as represented in shaded area, were determined by dividing the data into $5$ equal sections.}
\label{fig7}
\end{figure}

\begin{figure*}[t]
\includegraphics[width=0.99 \textwidth, angle=-0]{fig10.eps}
 \caption{2D projection of the diagonal (a-b) and off-diagonal (c) elements
 of the diffusion tensor, $D_{11}(Z,\Delta z)$, $D_{22}(Z,\Delta z)$,
 and $D_{12}(Z,\Delta z)$, associated with the crossing of Domperidone across
 the POPC membrane bilayer.
 %
 We restrained the exploration of the free energy landscape to the area
 centered around the minimal energy well ($Z \approx 10 \textrm{\AA}$
 in Figs.~\ref{fig7}(a) and (b)).
 }
\label{fig10}
\end{figure*}

The FEP and the diffusion profiles associated with the transport of Domperidone
across the POPC bilayer, as shown in Fig.~\ref{fig7}(b) and (c), respectively,
were obtained using the DHAM.
As the drug reaches the water/bilayer interface, the value of the
diffusion coefficient significantly increases before remaining approximately
constant until the molecule reaches the adjoining area of the two lipid layers,
which is representative of the hydrophobic nature of the compound.
These results can be interpreted within the inhomogeneous solubility-diffusion
model commonly used to provide realistic description of membrane
permeation~\cite{2016-BBA-Shinoda}. In this model, the permeability, Perm,
is written as
\begin{equation}
\frac{1}{\textrm{Perm}}=\int_{-d/2}^{+d/2}\frac{1}{K(Z)D(Z)}=\int_{-d/2}^{+d/2}\frac{\exp\big(\Delta G(Z)/k_B T\big)}{D(Z)}~,
\label{eqPerm}
\end{equation}
where $K(Z)$, $D(Z)$, and $d$ are the position-dependent partition coefficient, the solute diffusion
coefficient, and the membrane thickness respectively. $\Delta G(Z)$ is the free energy difference,
which is related to the partition coefficient $K(Z)=\exp\big(-\Delta G(Z)/k_B T\big)$.
Considering the results shown in  Fig.~\ref{fig7}, we obtained the $\log$ Perm
value associated with the transport of the Domperidone molecule across
the POPC bilayer, $\log \textrm{Perm} = -2.71 \pm 0.15$,  which is in agreement
with experimental results ($\log \textrm{Perm}^{(\textrm{exp})}=-2.6$)~\cite{2014-JCR-Eyer-Kramer,2017-JACS-Dickson-Duca,2018-JPCB-Badaoui-Rosta}.

\subsection{Membrane permeation in 2D}
We extended the 1D study above to the 2D analysis of the permeation of Domperidone
across the POPC lipid membrane. To do so, we  considered the
rotational movement of the drug during its passage across the membrane as
additional degree of freedom.
We constructed a discretized
two-dimensional grid to determine the MSMs along two reaction coordinates
(see details in the SI).
As our first reaction coordinate, we used the distance from the drug center
of mass to bilayer center, $Z$, as already considered in 1D.
For our second coordinate, we used the projection of the molecular
orientation vector onto the axis orthogonal to the bilayer axis, $\Delta Z$,
as a measure of the orientation of the drug with respect to the membrane,
as described in Badaoui et al.~\cite{2018-JPCB-Badaoui-Rosta}.
We restrained the exploration
of the free energy landscape to the area centered around the minimal
energy well ($Z \approx 10 \textrm{\AA}$ in Figs.~\ref{fig7}(a) and (b)).
We measured the diagonal and off-diagonal KM coefficients
defined in Eqs.~\ref{eqDriftDiffusion2D:1} and \ref{eqDriftDiffusion2D:2}
(see details in the SI). Subsequently, the associated elements of the diffusion
tensor were obtained with weighted polynomial regression.

In Fig.~\ref{fig10}(a), (b), and (c), we calculated the 2D projection of the diagonal
and off-diagonal elements of the diffusion tensor, $D_{11}(Z,\Delta z)$,
$D_{22}(Z,\Delta z)$, and $D_{12}(Z,\Delta z)$, respectively.
We observed good agreement in the order of magnitude between the measure
of $D_{11}(Z,\Delta z)$ (Fig.~\ref{fig10}(a)) and the one of $D^{(2)}(Z)$
(Fig.~\ref{fig7}(c)) in the minimal energy well.
Additionally, we measured that the orientation of the drug with respect to the membrane
diffuses an order of magnitude faster than the drug COM through the bilayer.
Finally, we estimated the off-diagonal element of the diffusion tensor,
which is essential for understanding the importance of dynamic coupling
between the coordinates. As shown in Fig.~\ref{fig10}(c), the order of magnitude
of $D_{12}(Z,\Delta z)$ is comparable to the one measured for $D_{11}(Z,\Delta z)$
and $D_{22}(Z,\Delta z)$.

\section{Conclusion}
In the present study, we presented a general, automatic method
for estimating multi-dimensional position-dependent diffusion
coefficients equally valid in biased or unbiased MD simulations.
We combined Markov State Model analysis and Kramers-Moyal  expansion
to link the underlying  stochastic process obtained within the
non-linear Langevin framework to its probabilistic description.
Our approach neither requires prior assumptions regarding the form
of the free energy landscape nor additional numerical integration scheme.
We applied our numerical approach to one- and two-dimensional
analytic potentials and data from explicit solvent molecular dynamics
simulations, including the water-mediated conformations of the alanine
dipeptide.\\

Importantly, we demonstrated the efficiency of our algorithm
in studying the transport of solute across three-dimensional
heterogeneous porous media, which is known to show long time
scale fluctuations potentially breaking ergodicity in sampling~\cite{2016-JCIM-Lee-Gumbart}.
Specifically, our algorithm provided accurate assessment of the diffusion
coefficient associated with the crossing of Domperidone
across a lipid POPC membrane, which was not previously accessible with
standard methods~\cite{2017-JACS-Dickson-Duca}.\\

The method would also allow the measure of the off-diagonal elements
of the diffusion tensor, which is essential for assessing the importance
of dynamic coupling between the reaction coordinates. As discussed by
Ma and coworkers~\cite{2006-JCP-Ma-Dinner}, the presence of significant
off-diagonal element in the diffusion tensor could indeed impact
the determination of the reactive eigenvector at play in the Kramers' theory
of reaction kinetics and then the measure of the rate associated with
the kinetic evolution of the system.
This is particularly important in  toxicology  and pharmacology, where
the prediction of kinetic quantities of solutes, including drugs and small peptides
provides fundamental   understanding of numerous biochemical transport processes
involved in the design of new drugs~\cite{2018-NRC-Amaro-Mulholland}.

\section{Acknowlegments}
We acknowledge use of the research computing facility at
King’s College London, Rosalind (https://rosalind.kcl.ac.uk).
F.S. and E.R. thanks Gerhard Hummer for fruitful discussions and Magd Badaoui
for his help in setting up the membrane bilayer system.
We acknowledge the support of the UK Engineering and Physical Sciences
Research Council (EPSRC), under grant number EP/R013012/1, and ERC project 757850 BioNet.