\section{{\bf{Abstract:}}}
\begin{abstract}
Methods that combine collective variable (CV) based enhanced sampling and global tempering approaches are \changes{used in speeding-up} the conformational sampling \changes{and free energy calculation} of large and soft systems with a plethora of energy minima.
In this paper, a new method of this kind is proposed in which the well-sliced metadynamics approach (WSMTD) is united with Replica Exchange with Solute Tempering (REST2) method.
WSMTD employs a divide-and-conquer strategy wherein high-dimensional slices of a free energy surface are independently sampled and combined.
The method enables one to accomplish a controlled exploration of the CV-space with a restraining bias as in umbrella sampling (US), and enhance-sampling of one or more orthogonal CVs using a metadynamics (MTD) like bias.
\changes{The new hybrid method proposed here enables boosting the sampling of more slow degrees of
freedom in WSMTD simulations, without the need to specify associated CVs, through a replica exchange scheme within the framework of REST2}.
The high-dimensional slices of the probability distributions of CVs computed from the united WSMTD and REST2 simulations are subsequently combined using the weighted histogram analysis method (WHAM) to obtain the free energy surface.
We show that the new method proposed here is accurate\changes{,} improve\changes{s} the conformational sampling\changes{, and achieves quick} convergence in free energy estimates.
We demonstrate this by computing the conformational free energy landscapes of solvated alanine tripeptide and Trp-cage mini protein in explicit water.
\end{abstract}
\maketitle

\section{Introduction}
Molecular dynamics (MD) simulations are widely used to study conformational sampling of large biological systems, compute free energetics and identify the mechanism of biochemical processes.\cite{Karplus_McCommon:Nat:2002,Schulten:2009,Dshaw:Science:2010,Shaw:ARB:2012}
Modelling large biological systems pose several challenges, primarily due to large number of conformational states and significant energetic and entropic barriers separating different conformational states.\cite{Schulten:2009}
Advanced enhanced sampling techniques are quintessential for achieving proper sampling of conformational states and obtaining reliable free energy estimates.\cite{Chipot:07,Vanden:2009:jcc,Tuckerman:Book,vanGunsteren:JCC:Rev:2010,Ciccotti:12,Giovanni:Entropy:2014,Pratyush:2016,McCammon:2016,Pietrucci:2017,Peters:Book,shalini:review:2019,vashisth:review:2019}
It is a common practice to use either CV-based or global tempering-based enhanced sampling techniques in combination with MD simulations to improve the sampling.
These approaches are applied to a variety of problems including protein  folding/unfolding, protein-drug binding/unbinding, transport of molecules through membranes, protein aggregation etc.\cite{Adam:AR:2007,Daniel:AR:2011,Wingbermuhle:2020,Acharya:2021}

Biased sampling,
generalized-ensemble or global tempering,
and the combination of \changes{the} both form the three major classes of enhanced sampling methods.
Generalized-ensemble methods achieve a random walk in configurational space by accelerating all the degrees of the freedom either by increasing the temperature of the system or by scaling-down the potential energy.
Several methods in this category use the replica exchange molecular dynamics (REMD)\cite{remd:1999} algorithm.
\changes{Methods such as} parallel tempering\cite{remd:1999},  replica exchange solute tempering (REST),\cite{rest:2005} and the modified version of REST, called REST2\cite{rest2:2011}, belong to this family.
Recently modified versions of REST method have been proposed for improving the sampling.\cite{Sugita:JCP:2018,REHT:2021}
Some other distinct \changes{global tempering methods are}
accelerated molecular dynamics(aMD),\cite{amd:2004}
chemical flooding\cite{chemical_flooding:1995}, integrated tempering sampling.\cite{ITS:2004}
and accelerated weight histogram method (AWH).\cite{AWH:2014,AWH-1:2021}

In other family of methods,
enhanced sampling of certain geometric variables
 of the system, called collective variables (CVs), are carried out by adding a bias potential along CVs or by enhancing the temperature of CVs.
Umbrella sampling (US),\cite{us:1974,us:1977} metadynamics (MTD),\cite{mtd:2002} adaptive biasing force method\cite{Darve:2001}, logarithmic mean-force dynamics (LogMFD),\cite{LogMFD:2012} driven-Adiabatic Free Energy Dynamics/Temperature Accelerated Molecular Dynamics (d-AFED/TAMD)\cite{tamd:1,tamd:2}, Variational Enhanced Sampling\cite{Variational:1}
and
other variants\cite{wt-mtd:2008,be:mtd:1,pbmtd:2015,shalini:2016} are some examples of biased sampling methods.
Such methods require {\em a priori} selection of a set of CVs that describes the process of interest.
The accuracy of the properties computed from the explored conformational ensemble may depend on the quality of the chosen CVs.
Identification of optimal CVs is a challenging task.\cite{Frank_Noe:2017,Henkelman:2017,VPande:2018,Gerhard:2018}
Inclusion of large number of CVs for biased sampling is often required for an efficient exploration of the CV-space and for quick convergence in free energy estimates.
The efficiency of most of the aforementioned
biased-sampling techniques diminishes
with increasing dimensionality of the CV-space.
Sampling of high dimensional free energy landscapes requires advanced techniques such as
d-AFED/TAMD,\cite{tamd:1,tamd:2}
bias-exchange MTD,\cite{be:mtd:1}
parallel-bias MTD,\cite{pbmtd:2015}
unified-free energy dynamics,\cite{ufed:2012}
and
temperature accelerated sliced sampling\cite{tass:2017}.

To take the best out of both the generalized ensemble and the CV based biased-sampling methods, hybrid sampling algorithms are proposed.
Replica exchange umbrella sampling (REUS),\cite{reus:2013}
replica exchange with CV tempering,\cite{rect:2015}
combination of parallel tempering with MTD,\cite{ptmtd:2006,reus-rest}
replica state exchange MTD,\cite{rse-mtd:2015} and multi-scale sampling using temperature accelerated and REMD\cite{mustar-MD} are examples of such methods.
It has been observed that for modelling protein folding/unfolding and protein-drug binding/unbinding, these methods are advantageous.\cite{ptmtd:2006,rect:2015,reus-rest,reus:2013,Applications:1,Application:2, Application:3,Application:4,Application:REDS:2021}

In this paper, we introduce a new method called globally accelerated sliced sampling (GASS) by integrating the WSMTD and REST2 approach.
WSMTD\cite{shalini:2016} was introduced by our group for increasing the efficiency
of MTD method by driving the bias along a specific direction and by controlling the span of the explored CV-space.
\changes{This method is suited to study systems with broad, flat and unbound free energy landscapes.}
This method has been applied to study various problems\cite{Application:wsmtd:2017,Application:wsmtd:2018_1,Application:2018:2,Application:wsmtd:2020:1,Application:wsmtd:2020:2,Application:wsmtd:2021} and has also been extended to deal with high-dimensional CV-space.\cite{tass:2017}
By combining REST2 with WSMTD, we hope to boost the sampling of hidden transverse coordinates while \changes{exploring} the relevant CV space.
The controlled biased sampling feature of WSMTD could thus be extended to the conformational sampling of large biomolecular systems in solution and is expected to be beneficial for studying problems like protein folding.

\changes{Here, we first} introduce the theory behind the GASS method.
Subsequently, we present the results of two applications using the GASS method: computation of \changes{conformational} free energy landscapes for (a) alanine tripeptide, and (b) Trp-cage mini protein in water.
\section{Methods}

\changes{WSMTD is a CV based enhance sampling method
which can help to achieve a controlled exploration of the CV space by combining US and MTD.
While US is ideal for a controlled or directional sampling along one CV, MTD is advantageous in sampling orthogonal CVs in a self guided manner.
In WSMTD, our interest is in computing an $n$-dimensional free energy surface $F(\mathbf s$) as a function of CVs, $\mathbf s(\mathbf R)\equiv \{s_1,\cdots,s_n\}$, where $\mathbf R$ is the set of atomic coordinates.
}
WSMTD\cite{shalini:2016} uses the Lagrangian,
\begin{eqnarray}
      \mathcal L_h^{\rm ws}(\mathbf R ,\dot {\mathbf R})
   = \mathcal L^{0}(\mathbf R ,\dot {\mathbf R}) -
      W^{\rm b}_{h}(s_1)
      - V^{\rm b}({\overline{\mathbf s}},t),
      \label{e:wsmtd:Lagrangian}
\end{eqnarray}
$h=1, ..., M $,
where, $\mathcal L^0$ is the unbiased Lagrangian, \changes{$\dot{\mathbf R}$
is the set of velocities, and ${\overline{\mathbf s}} \equiv \left \{ s_2,\cdots,s_n  \right \}$.}

\changes{In} WSMTD\changes{, we apply a} restraining bias
\begin{eqnarray}
    W^{\rm b}_{h}(s_1) = \frac{1}{2}k_{h}(s_1 - z_h)^2
    \label{e:us:bias}
\end{eqnarray}
placed at $s_1 = z_h$\changes{. T}he parameter $k$ determines the curvature of the biasing potential at $h^{\rm th}$ window.
\changes{Here, $M$ umbrella windows placed from $z_1$ to $z_M$ along  $s_1$ define the span of sampling of that CV.}
This helps in achieving a controlled sampling of the coordinate $s_1$, similar to the conventional US technique.
\changes{In WSMTD, we also apply} a well-tempered MTD bias potential
\begin{equation}
V^{\rm b}( {\overline{\mathbf s}},t)= \sum_{\tau < t}w_\tau \exp\left[- \frac{ \left ( {\overline{\mathbf s}} - {\overline{\mathbf s}}(\tau) \right )^2}{2(\delta s)^2}  \right] \enspace ,
\end{equation}
in order to sample the orthogonal coordinates ${\overline{\mathbf s}}$.
Here, $\delta s$ is the width parameter defining the Gaussian potential and
\begin{equation}
    w_\tau = w_0 \exp \left[- \frac{V^{\rm b}({\overline{\mathbf s}},t)}{k_{\rm B} \Delta T}  \right] \enspace ,
\end{equation}
\changes{Where} $w_0$ is the Gaussian height parameter, $k_{\rm B}$ is the Boltzmann constant, and $\Delta T$ is a tempering parameter.

The free energy surface, $F(\mathbf s)$, is reconstructed from the probability distribution, $P(\mathbf s)$, of the CVs as,
\begin{equation}
    F(\mathbf s)=-\frac{1}{\beta}\ln P(\mathbf s)
\end{equation}
where $\beta = (k_{\rm B} T)^{-1}$. \changes{In order to obtain $P(\mathbf s)$, a} time independent probability distribution, $P_h^{\rm u}(\mathbf s)$, is obtained by reweighting the MTD time-dependent bias potential as,
\begin{widetext}
\begin{eqnarray}
    P^{\rm u}_{h}(\mathbf s^\prime)=
    \frac{\int_{t_{\rm min}}^{t_{\rm max}}d\tau \, \exp \left[
    \beta \{V^{\rm b}({\overline{\mathbf s}},\tau)-c(\tau)\} \right]
    \prod_{\alpha=1}^n \delta(s_\alpha(\tau)-s_\alpha^\prime)
    }
    { \int_{t_{\rm min}}^{t_{\rm max}}d\tau \, \exp \left[
    \beta \{V^{\rm b}({\overline{\mathbf s}},\tau)-c(\tau)\} \right]
    }  \enspace  ,
    \label{e:prob:unbias:mtd}
\end{eqnarray}
\end{widetext}
\changes{for each umbrella window $h$}.
\changes{In the above,} $c(t)$ is evaluated as\cite{Tiwary-Parinello:2015}
\begin{equation}
    c(t)=\frac{1}{\beta}\ln \left [ \frac{\int d{\overline{\mathbf s}} \exp[-\beta \gamma V^{\rm b}({\overline{\mathbf s}},t)] }{\int d {\overline{\mathbf s}} \exp[-\beta (\gamma-1)  V^{\rm b}({\overline{\mathbf s}},t)]}\right ]
\end{equation}
and $\gamma = {(T+\Delta T)}/{\Delta T}$ .

Now we employ the weighted histogram analysis method (WHAM)\cite{wham:1,wham:2} to combine the distributions $P^{\rm u}_h(\mathbf s)$, $h=1,\cdots,M$ and reweight the restraining bias by self-consistent calculations using,
\begin{equation}
P(\mathbf s) = \frac{\sum^{M}_{h=1}n_h P^{\rm u}_h(\mathbf s)}{\sum_{h=1}^{M}n_h \exp [\beta f_h] \exp[-\beta W_h^{\rm b}(s_1) ]}
\end{equation}
and
\begin{equation}
    \exp [-\beta f_h]=\int d \mathbf s \,  \exp \left [-\beta W_h^{\rm b}(s_1)] \, P(\mathbf s) \right ] \enspace .
\end{equation}
\changes{In the above,} $n_h$ is the number of configurations sampled in the $h^{\rm th}$ \changes{umbrella} window.

In the REST2 method,\cite{rest:2005,rest2:2011} the potential energies of a selected set of atoms in each replica are scaled-down by some parameter $\lambda$ to enhance their sampling, and conformations sampled in these scaled-replica are exchanged with the non-scaled replica to improve the conformational sampling of the latter.
The atoms which are part of the scaled-potential terms (in a molecular mechanics force-field) are referred to be in ``hot'' region and the rest \changes{are said to be} in ``cold'' region.
Pair-wise potential form of the molecular mechanics force-field
makes it easier to selectively scale potential contributions for a set of atoms and
differently scale torsional, electrostatics, and Lennard-Jones terms.
Typically, for a solvated protein system, the potential energy of a replica $m$ is computed from the modified
contributions of the protein-protein (pp), protein-water (pw) and water-water (ww) interactions, as
\begin{equation*}
    U_m^{\rm REST2}({\bf{R}})= \lambda_m \, U_{\rm pp}({\bf{R}})+\sqrt{\lambda_m} \, U_{\rm pw}({\bf{R}})+U_{\rm ww}({\bf{R}}),

\end{equation*}
$ m=0,\cdots,N_{\rm r}-1$,
where
$\lambda_m={\beta_m}/{\beta_0}$,
$N_{\rm r}$ is the number of replicas, and $\beta_m > \beta_0$.
 The $\beta_m$ parameter is \changes{set} different for different replicas.
For $m=0$, i.e., for the unscaled replica, we take $\beta_0 = (k_{\rm B} T_0)^{-1}$, where $T_0$ is the
physical temperature of the system.
Exchange between adjacent replicas is attempted by swapping their atomic coordinates based on the Metropolis exchange criterion
\begin{eqnarray} p(i\rightarrow j)=\min(1,e^{-\Delta_{i,j}}) \enspace, \label{e:detailed:balance}
\end{eqnarray}
\changes{with}
\begin{widetext}
\begin{equation}
\Delta_{i,j}=(\beta_i - \beta_j)\left [ (U_{\rm pp}({{\mathbf R}_j})- U_{\rm pp}({{\mathbf R}_i})) + \frac{\sqrt{\beta_0}}{\sqrt{\beta_i}+\sqrt{\beta_j}}(U_{\rm pw}({\mathbf{R}_j})- U_{\rm pw}({\mathbf{R}_i}))      \right ]
\enspace .
\label{e:rest2:delta}
\end{equation}
\end{widetext}
The absence of solvent-solvent interaction in the above expression boosts the acceptance  between neighbour replicas \changes{compared to conventional parallel tempering simulations}.

\changes{WSMTD is usually used for cases where  the number
of CVs is small.
To improve the efficiency of WSMTD for large number of CVs, TASS method was proposed.\cite{tass:2017}
Like WSMTD, TASS is also a CV-based enhanced sampling method.
However, it is nontrivial to find suitable CVs for enhancing the slow global motions in large soft matter systems like solvated proteins.\cite{Bolhuis:BPJ:10}
As discussed earlier, global tempering methods like REST2 can achieve this by potential energy scaling in replica exchange.
Thus to allow enhanced sampling of global motions of large molecular systems in WSMTD, we introduce the GASS approach by combining WSMTD and REST2.
}

\changes{Like in WSMTD, our aim in GASS simulations is to compute the free energy surface $F(s_1,\cdots,s_n)$ where $n$ is typically not more than three.
We
prefer to achieve a controlled sampling along
$s_1$ together with a self-guided biased sampling along
other $n-1$ CVs.
For a quick and accurate estimation of $F(\mathbf s)$, we also intent to boost the enhance sampling of slow
global conformational changes without defining
 additional CVs.}
In GASS approach, we use the WSMTD Lagrangian (\eref{e:wsmtd:Lagrangian}) for the windows $h=1,\cdots,M$.
For every umbrella $h$, we perform REST2 simulation considering $N_{\rm r}$ replicas, and construct the  probability
density $P^{\rm u}_h(\mathbf s)$ from the REST2 trajectory for the $m=0$ replica using \eref{e:prob:unbias:mtd}.

In order to accommodate the bias-potentials acting on the replica for a given window $h$,
a modified equation for computing exchange probability is used for GASS.
To satisfy the condition of detailed balance  we use
\begin{equation}
    \Delta_{i,j}= \Delta^{(1)}_{i,j} + \Delta^{(2)}_{i,j}
\end{equation}
in \eref{e:detailed:balance}, where $\Delta^{(1)}_{i,j}$ is given by \eref{e:rest2:delta} and
\begin{eqnarray*}
\Delta^{(2)}_{i,j} = \beta_0 \left\{ \left [V^{\rm b}_{i}({\overline{\mathbf s}}_{j}, t) - V_{i}^{\rm b}({\overline{\mathbf s}}_{i},t)]-[V_{j}^{\rm b}({\overline{\mathbf s}}_{j},t) - V_{j}^{\rm b}({\overline{\mathbf s}}_{i},t)  \right ] \right\} \enspace .
\end{eqnarray*}
 Here $i$ and $j$
are the replica indices; ${\overline{\mathbf s}}_i$ specifies the set of CVs  $(s_2,\cdots,s_n)$ corresponding to the coordinates $\mathbf R_i$ of the replica $i$.

Since all the replicas for a given window $h$ have the same umbrella bias (\eref{e:us:bias}),
it is not contributing to the calculation of $\Delta_{i,j}$.
\changes{See Supporting Information for a derivation of the above expression for exchange probability.
}
The GASS method has been implemented using the GROMACS/PLUMED interface.\cite{GROMACS,plumed2:2014}
All simulations in this paper were performed using  GROMACS-2018.6\cite{GROMACS} patched with PLUMED-2.2.6\cite{plumed,plumed2:2014} and HREX.\cite{Bussi_rest2:2014}

\begin{figure*}[t]
\centering
    \includegraphics[width=0.5\textwidth]{Ala_tri_pep_GASS_3.eps}
    %
    \caption{(a) Ball and stick representation of alanine tripeptide (Color codes: H - white, C - black, O - red, and N - blue).
    %
     Free energy surface  $F(\phi_1, \phi_2)$
     computed using
     GASS (b), WSMTD (c), and REST2 (d) simulations of alanine tripeptide in water are given.
     %
    Convergence of free energy barriers for $\mathbf P \rightarrow \mathbf Q$, $\mathbf Q \rightarrow \mathbf P$, $\mathbf R \rightarrow \mathbf Q$, and $\mathbf Q \rightarrow \mathbf R$ as a function of simulation time in
    GASS (e), WSMTD (f), and REST2 (g) simulations are shown.
    %
    Here the angles are in radians and free energies are in kcal~mol$^{-1}$.
    }
    \label{fig:fes_alatripep}
\end{figure*}

\section{Results and Discussion}

The accuracy and efficiency of the GASS method was first investigated for
the exploration of the  conformational free energy landscape of alanine tripeptide in explicit water.
The performance of the method is compared with REST2 and WSMTD.
Thereafter, it was tested for computing the free energy landscape of unfolding of Trp-cage protein in water.
\subsection{Alanine Tripeptide in Water}

\begin{table*}[htbp]
\centering
\begin{tabular}{||c| c | c | c | c | c | c | c ||}
 \hline
 Method & Simulation Time (ns) & \multicolumn{4}{|c|}{$\Delta F^\ddagger$} & \multicolumn{2}{|c||}{$\Delta F$} \\ \cline{3-8}
        &                      &  $\mathbf P \rightarrow \mathbf Q$    & $\mathbf Q \rightarrow \mathbf P$ & $\mathbf R \rightarrow \mathbf Q$ & $\mathbf Q \rightarrow \mathbf R$ & $\mathbf P - \mathbf Q$ & $\mathbf R - \mathbf Q$
     \\
 \hline\hline
 GASS & 10 & 4.4 & 6.0 & 4.4 & 6.0 &   1.6 &   1.6  \\
 WSMTD & 10 & 4.4& 6.0 & 4.4 & 6.0 &   1.6  & 1.6  \\
 REST2 & 1000 & 4.6 & 6.1 & 5.1& 6.6 & 1.5  & 1.5 \\[1ex]
 \hline
\end{tabular}
\caption{Free energy barriers ($\Delta F^\ddagger$)  and free energy differences ($\Delta F$), in kcal~mol$^{-1}$, computed from the free energy landscape $F(\phi_1,\phi_2)$ of alanine tripeptide in water using various methods.}
\label{table:tab_alatripep}
\end{table*}

Alanine tripeptide in water was modeled using {\tt AMBER14SB}\cite{ff14SB} and TIP3P\cite{TIP3P} force-fields.
We took the amino acid sequence ACE-ALA-ALA-NME where the terminal residues ACE and NME were acetyl and N-methyl amide, respectively, and the total number of water molecules in the system is 785.
All the bonds in the system were constrained and the equations of motion were integrated using a time step of 1~fs.
Long range electrostatics was treated using the particle-mesh Ewald technique, as available in GROMACS.\cite{PME}
Temperature of the system was controlled using the stochastic velocity rescaling thermostat.\cite{v-scaling:thermostat}
Density was converged in 100~ps of $NPT$ simulation and the
equilibrated cell volume \changes{achieved was} $34 \times 36 \times 30 $~{\AA}$^{3}$.
GASS simulations were then carried out in the $NVT$ ensemble using the equilibrated density.
In REST2, all the peptide atoms were taken in the hot region and 5 replicas were considered with the values of $\lambda_m$ ranging from 1.0 to 0.3 following a geometric distribution.
Exchanges between the neighboring replicas were attempted every 1000 MD steps.
 The two Ramachandran angles $\phi_1$ and $\phi_2$ were chosen as CVs for biasing; see \fref{fig:fes_alatripep}a.
Umbrella bias was applied along $\phi_1$ from $-\pi$ to $+\pi$ at an interval of 0.2 radians with $\kappa_h=1.2 \times 10^{2}$ kcal~mol$^{-1}$ rad$^{-2}$ and \changes{a well-tempered MTD} bias was applied along $\phi_2$.
The MTD bias was updated every 500~fs, and the bias parameters $w_0=0.6$~kcal~mol$^{-1}$,  $\delta s=0.05$~radians, and $\Delta T = 900~$K were taken.
\changes{It is noted that the choice of types of biases}
applied on \changes{$\phi_1$ and $\phi_2$}
\changes{were arbitrary}.

For demonstrating the performance of GASS, we carried out independent WSMTD, and REST2 simulations with identical setups. %

 Figure ~\ref{fig:fes_alatripep} shows the free energy surface $F(\phi_1,\phi_2)$ computed using GASS, WSMTD, and REST2 and the convergence of free energy barriers for these  surfaces.
The free energy barriers $\mathbf Q \rightarrow \mathbf P$ and $\mathbf Q \rightarrow \mathbf R$
are expected to be nearly the same due to the symmetry, and the same holds true for the corresponding reverse barriers.
From \fref{fig:fes_alatripep}e, it is clear that very accurate predictions of free energy barriers are possible even after 2~ns of GASS simulation.
As anticipated, the $\mathbf Q \rightarrow \mathbf P$ and  $\mathbf Q \rightarrow \mathbf R$ barriers are found to be equal, and the same is true for $\mathbf P \rightarrow \mathbf Q$ and $\mathbf R \rightarrow \mathbf Q$ barriers.
A smooth free energy surface was obtained after 10~ns of GASS simulation, and the system  swept through the entire two-dimensional landscape.

The free energy barriers computed from \changes{independent} WSMTD are agreeing very well with the GASS results.
Interestingly, the convergence of free energy barriers \changes{in their runs} was as quick as that \changes{using} GASS.
However, the higher energy region\changes{s near } $\phi_1=2.4$, $\phi_2=2.4$ (radians) \changes{were} not well explored in WSMTD (\fref{fig:fes_alatripep}c).
The \changes{results of independent } REST2 simulations  (\fref{fig:fes_alatripep}d,g) show that conformations far from the free energy minima are not well explored even after \changes{1~$\mu$s} and the estimates of the free energy barriers \changes{($\Delta F^\ddagger$)} are not well converged.
On the other hand, free energy differences \changes{($\Delta F$)} computed using GASS, WSMTD, and REST2 are
agreeing very well; see \tref{table:tab_alatripep}.

These results give us the confidence that GASS method is able to \changes{provide accurate free energy estimates and could} achieve quick conformational sampling.
The method is evidently performing better than both WSMTD and REST2 methods.

\subsection{Conformational Landscape of Trp-Cage in Water}

\begin{figure*}[htbp]
\centering
    \includegraphics[width=0.5\textwidth]{Figure2_updated.eps}
    \caption{ Conformational free energy landscape of Trp-cage protein computed along the RMSD and Rg CVs using (a) GASS (b) WSMTD (c) REST2. (d) Convergence of free energy difference
    ($\Delta F = F ( {\mathbf{U}} - {\mathbf{F}})$ )
    between the unfolded ({\bf U}) and the folded state ({\bf F}) and convergence of free energy barrier for going from {\bf F} to the intermediate {\bf I}  ($\Delta F^\ddagger ({\mathbf{F}\rightarrow\mathbf{I}})$) are plotted as function of simulation time.

    (e)  Representative conformations of {\bf F}, {\bf I} and {\bf U} states are shown. The NMR structure (green) is overlapped with {\bf F} for a comparison.  Number of other conformations were also observed in {\bf I} and {\bf U} states; See SI Section 4.}

    \label{fig:fes_trpcage}
\end{figure*}

As next, we investigated the unfolding/folding landscape of Trp-cage in explicit water using the GASS method.
This is an ideal problem to study using GASS as a controlled sampling along a ``folding/unfolding coordinate'' using the restraining bias could drive the conformations from folded to unfolded states, or vice-versa.
At the same time, accelerated sampling of several orthogonal \changes{coordinates} are essential to boost the exploration of the conformational states \changes{for such systems} which can be achieved by \changes{using} MTD bias and REST2.

Trp-\changes{c}age is a 20 amino acid mini-protein (NLYIQ WLKDG GPSSG RPPPS)  designed by Neidigh et al.\cite{trp:2002:original} from 39 amino acid extendin-4 peptide.
It is a fast-folding protein and is \changes{considered as} an ideal model system for testing \changes{computational} methods developed for protein-folding problems.
It contains a short $\alpha$-helix (residues 2-9), a $3_{10}$-helix (residues 11-14), and a C-terminal polyproline-II helix.

The initial Trp-cage structure for our simulations was built based on the folded NMR structure PDB ID 1L2Y.\cite{trp:2002:original}
The protein was first solvated in a periodic $80 \times 80 \times 80$~\AA$^3$
TIP3P water box \changes{with} 1~g~cm$^{-3}$ density.
\changes{We used} AMBER ff14SB\cite{ff14SB} force-field \changes{for the protein.} \changes{L}ong-range electrostatic interactions were evaluated using the Particle Mesh Ewald method.\cite{PME}
All the bonds in the system were constrained using the LINCS algorithm.\cite{LINCS}
A time step of 2~fs was used to integrate the equation of motion.
Temperature of the system was controlled by stochastic velocity rescaling thermostat.\cite{v-scaling:thermostat}
A 1~ns long $NPT$ run was initially carried out using Parrinello-Rahman barostat\cite{Parrinello:Rehman} to obtain converged density and the equilibrated cell volume \changes{was found to be} $83 \times 83 \times 83$~{\AA}$^3$.
Using the equilibrated \changes{system}, we performed  10~ns equilibration in $NVT$ ensemble to generate initial structure\changes{s} for the GASS simulation.
 %

 All the protein atoms were chosen in the hot region for the REST2 \changes{calculations.}
 %
 We chose 20 replicas with $\lambda_m$ values ranging from 1.0-0.3 to obtain a good exchange probability.
 %
 The exchange of coordinates was attempted \changes{after} every 1000~MD steps.
 %

We chose to apply restraining bias along the backbone root mean square deviation (RMSD) CV  in order to drive the unfolding of the protein from the folded starting structure;
See SI~Section~1.1 for the definition of the coordinate.
 %
As the reference structure for computing RMSD, we took the backbone structure of all the heavy atoms \changes{for the residues 1 to 15} in the PDB 1L2Y.
 %
This CV will serve as the ``unfolding coordinate'' which will guide the system from folded state to unfolded state in a controlled manner.
 %
A total of 45 umbrella windows were placed every {0.20}~{\AA}, starting from 0.20~{\AA}.
 %
MTD bias was applied along the radius of gyration (Rg) in all the windows; See SI~Section~1.2 for the definition of Rg.
This bias potential could enhance the conformational sampling by triggering the breakage of $\alpha$-helices, thereby promoting the system to visit unfolded states.
Since GASS approach enables us to have different transverse coordinates for different windows, we
 chose the $\alpha$-helicity (SI~Section~1.3) CV together with Rg for applying MTD bias for the windows placed between RMSD values 6.0~{\AA} and 8.0~{\AA}.
For the MTD bias, we choose the Gaussian parameters $w_0=0.6$~kcal~mol$^{-1}$, and $\delta s = 0.02$.
We took $\Delta T= 2700$~K and $5700$~K when we used one-dimensional and two-dimensional MTD biases, respectively.
We scaled-up the width of the Gaussian by a factor 10 in the direction of the $\alpha-$helicity coordinate because of the
large amplitude fluctuation of this CV compared to Rg.
The starting umbrella window at RMSD=$0.20$~{\AA}  pertained to a folded structure, resembl\changes{ing} well with the NMR structure (PDB~ID: 1L2Y).
Starting structure of all the other windows with higher value of RMSD was taken as the equilibrated structure of the \changes{preceding} window.

We reconstructed the free energy surface in the RMSD-Rg space after 120~ns of the simualation (\fref{fig:fes_trpcage}a).
Free energy surfaces constructed after 50~ns, 75~ns, 100~ns of GASS simulation are presented in SI~Section~2.
The \fref{fig:fes_trpcage}a clearly indicates folded $(\mathbf F)$, semi-unfolded/intermediate $(\mathbf I)$ and unfolded $(\mathbf U)$ states of Trp-cage.
Representative structures of various conformational minima found on the free surface are obtained by clustering the trajectory of the corresponding GASS window (\fref{fig:fes_trpcage}f).
They indicate that the conformational space, including the unfolded states, is very well explored.
The free energy barrier for $\mathbf F \rightarrow \mathbf I$, and free energy difference between $\mathbf F$ and $\mathbf U$ are computed as a function of simulation time (\fref{fig:fes_trpcage}d).
This plot indicates that the estimate of free energy difference was converged by 100~ns and free energy barrier was converged by 60~ns.
After 90~ns, new unfolded states were explored, resulting in a
small drift in the convergence curve of $\Delta F$.

\changes{
Detailed studies on the mechanism of folding have been done experimentally\cite{WountersenJPCB:2013,trp:2002:original,Qiu:2002,Ahmed:2005,Neuweiler:2005,Iavarone:2005,Bunagan:2006,Streicher:2007,Iavarone:2007,Mok:2007,Hudaky:2008,Barua:2008,Culik:2011,Rovo:2011,Ronge:2012,Anna:2012,Rovo:2013,ADAMS:2006,Tucker:JPCL:2020}
and computationally.\cite{WountersenJPCB:2013,Vijaypande:2002,Adrian:2002,CHOWDHURY:2003,Zhou:PNAS:2003,Bolhuis:PNAS:2006,Bolhuis:2008,HU:2008,Lindorff:2011,Shao:2012,WU:2011,Marino:2012,Lai:2013,XU:2008,Bandyopadhyay:PCCP:2016,Ferguson:2018,Kim:2016,Hatch:2014,Kannan:2014,Alexander:2005,Zhou:Proteins:2003,Levy:JPCB:2013,Patriksson:2007,Pitera:2003,trp:2002,be:mtd:1,be:mtd:2,trp:2011,trp:2012,trp:2013,trp:2014,trp:2015,trp:2015:amd,trp:2017,trp:2018,REHT:2021}
The protein folding landscape of Trp-cage shows intermediate states, different to the classical  picture with only two states.
The presence of metastable intermediate states have been seen observed computationally~\cite{Zhou:PNAS:2003,Lindorff:2011} and experimentally.\cite{Sauer:PNAS:2005}
Two major folding pathways have been identified for this protein.\cite{Bolhuis:PNAS:2006, Debenedetti:JCP:2015,Laio:PLOS:2009,Levy:JPCB:2013}
The global minimum on the GASS free energy surface
(\fref{fig:fes_trpcage}a) is located at RMSD$=1.1$~{\AA},
and Rg$=6.8$~{\AA}, and it corresponds to the folded state ($\mathbf{F}$).
The same was also observed in independent 120~ns WSMTD (\fref{fig:fes_trpcage}b) and 3~$\mu$s REST2 (\fref{fig:fes_trpcage}c) simulations.
The intermediate ($\mathbf{I}$) and the unfolded ($\mathbf{U}$) states are located on the (RMSD, Rg) free energy surface at ($5.7$~{\AA}, $7.5$~{\AA}), and ($6.2$~{\AA}, $11.0$~{\AA}), respectively.
The intermediate state {\bf I} and the unfolded state {\bf U} comprised of large number of protein conformations (\fref{fig:fes_trpcage}e).
The representative structures presented in \fref{fig:fes_trpcage}e were obtained by backbone RMSD based clustering of the biased GASS trajectories.
We also observed the intermediates {\bf SB-I}, {\bf LOOP-I}, and {\bf HLX-I} reported by Kim~et~al.\cite{Debenedetti:JCP:2015}; see Supporting Information.
However, it is noted that proper characterization of the intermediate states and reaction pathways require analysis of the unbiased GASS trajectories, which was not carried out here.
The free energy barrier $\mathbf{F}\rightarrow\mathbf{I}$ is $2.3$ kcal~mol$^{-1}$ using GASS, while it is  $2.4$ kcal~mol$^{-1}$ and $3.2$ kcal~mol$^{-1}$ using WSMTD and REST2 simulations, respectively.
From the GASS-computed free energy surface, we calculated the barrier for going from {\bf F} to {\bf U} state as $3.0$~kcal~mol$^{-1}$.
The {\bf I} and the {\bf U} states are only $0.0$ and $0.7$~kcal~mol$^{-1}$ higher than the {\bf F} state, respectively.

The free energy data and the conformational landscape reported in our study are in excellent  qualitative and quantitative agreement with various reports in the literature.\cite{trp:2015:amd,Gracia:2010}
Barrier for ${\bf F} \rightarrow {\bf I}$
 also agrees with that reported in  Refs.\cite{REHT:2021,Lindorff:2011}
The free energy difference between the unfolded and the folded state was experimentally measured as 0.77 kcal~mol$^{-1}$ at 298K by Streicher and Makhatadze,\cite{Makhatadze:BioChem:2007} which is in remarkably good agreement with our estimate ($0.7$~kcal~mol$^{-1}$).

Although the folded (${\mathbf{F}}$) and intermediate (${\mathbf{I}}$) states were reasonably sampled in WSMTD, the unfolded (${\mathbf{U}}$) states were not (\fref{fig:fes_trpcage}b).
The computed free energy difference between the {$\mathbf{U}$} and ${\mathbf{F}}$ states is  not converged in 120~ns of WSMTD (\fref{fig:fes_trpcage}d).
Sampling of the unfolded states was insufficient even after 3~$\mu$s of REST2; See \fref{fig:fes_trpcage}c.
The extent of exploration of the RMSD-Rg space is much less in REST2 compared to that in GASS.
Further, the free energy surface computed from REST2 trajectory is noisy.
Thus we think that REST2 free energy estimates are likely to be not converged.
These results underline the importance of the GASS approach.
}

\section{Conclusions}
It has been shown that the combination of the global tempering REST2 method with the CV-based biased sampling technique WSMTD, as done in the proposed GASS method, is an efficient way to study conformational sampling and compute free energies of large soft matter systems in solution.
A directed conformational sampling achieved by a restrained-bias along a CV, in tandem with an exhaustive sampling of orthogonal coordinates by MTD bias and REST2, makes the GASS method different to other sampling methods.
This is a much needed feature for studying various biochemical processes, such as protein folding.
Test calculations performed on solvated alanine tripeptide show that the GASS method provides accurate prediction of conformational free energy landscapes.
The method has been applied to study the unfolding/folding free energy landscape of Trp-cage protein in water, wherein a controlled unfolding of the protein was accomplished by applying restraining bias along the RMSD coordinate.
Orthogonal coordinates were enhanced sampled by REST2, in concert with the explicit MTD bias on Rg and $\alpha$-helicity CVs.
A quick convergence in free energy estimates was observed and a good conformational sampling of the unfolded states was noted.
The free energy landscape projected on the RMSD-Rg space has three distinct free energy minima corresponding to fully folded, partially unfolded (intermediate), and unfolded states.
On the landscape, the free energy barrier to go from folded to the intermediate state is $2.3$~kcal~mol$^{-1}$ and the folded to the unfolded state is $3.0$~kcal~mol$^{-1}$.
The intermediate and the unfolded states are only $0.0$ and $0.7$~kcal~mol$^{-1}$ higher than the folded state, respectively.
We hope that the new method proposed here will be very useful to study mechanism and free energetics of complex biochemical processes such as protein folding, drug-binding, and  diffusion of molecules through membrane.
\changes{The reweighting codes and the input files for performing GASS simulations reported here are  available online.\cite{github:gass} }

\begin{acknowledgments}
Authors gratefully acknowledge the discussions with Prof. Ricardo L. Mancera (Curtin University). ABK thanks IIT Kanpur and Curtin University for the PhD fellowship and travel support. Authors thank IIT Kanpur for providing computing resources at the HPC2013 cluster.
\end{acknowledgments}

\section*{Supporting Information}
Additional supporting information may be found in the online version of this article.


\section{Introduction}

Adversarial examples were firstly introduced by Szegedy et al. who demonstrated that by applying a certain hardly perceptible perturbation to an image, one neural network is driven to misclassify the image.
The reason is that deep neural networks learn input-output mappings that are fairly discontinuous to a significant extent. The attack can be performed by maximizing the network’s prediction error.
In addition, the specific nature of these perturbations is transferable. The same perturbation can cause the same input to be misclassified by a different network, even when it was trained on a different dataset \cite{szegedy2013intriguing}.

Most methods used to generate adversarial perturbations rely either on gradient-based information or on confidence scores such as class probabilities. For instance, DeepFool \cite{moosavi2016deepfool} relies on an iterative linearization of the classifier to generate minimal perturbations that are sufficient to change classification labels.
However, black-box attacks that are solely based on the final model decision have been recently shown feasible \cite{brendel2017decision}. One strategy consists
in training a local model to substitute for the target network,
using inputs synthetically generated by an adversary and
labeled by the target network \cite{papernot2017practical}. The local substitute  is used to craft adversarial examples, that found to be misclassifed by the targeted deep network.

Adversarial machine learning imposes a threat in real world applications concerning people's safety. For example, an adversarial attack that fools traffic sign classification in autonomous cars may cause devastating consequences.
It is risky to apply neural networks in security-critical areas today without adequate defense techniques and adversarial training.
That's why defense and detection of adversarial attacks is currently a research hotspot.
Several defense mechanisms have been proposed to this date such as regularization and data augmentation. A new
variant of distillation was proposed in  \cite{papernot2016distillation} to provide for defense training.
The basic idea is to smooth the network representative function by transferring knowledge encoded into soft probability vectors.  The soft labels are obtained through training with the same network architecture and "Temperature" parameter.

More similar to our work is Feature Squeezing \cite{xu2017feature}.
Feature squeezing coalesces samples that correspond
to many different feature vectors in the original space into a single sample. This can be achieved for instance by reducing the
color bit depth of each pixel or spatial smoothing. This technique aims at reducing the search space available to an adversary. By comparing the prediction of two DNN models: the original one and the input-squeezed one,
adversarial examples are shown to be detected with high accuracy and few false positives.
A similar idea appeared in \cite{zhang2018detecting} based on saliency maps and is connected to the concept of explainability in machine learning.
At visual level, saliency is the regions of the image
that were particularly influential to the final classification.
The saliency maps of adversarial samples appear to be different than those derived from benign samples. The authors train a binary classifier benign/adversarial based on saliency maps.

Multi Layer - Leave One Out (ML-LOO) \cite{yang2020ml} proposes another idea using feature (pixel) attribution rather than saliency. Adversarially crafted samples manifest differences in feature attribution histograms when compared to those of original samples. These differences can be captured in form of Inter Quartile Range (IQR) values. A binary classifier can be further trained to flag suspicious input samples.
The approach can work in black box mode where only class probabilities are used, or in white box (Multi-Layer) mode where feature attribution with respect to internal units' values are used as well. The results of ML-LOO are promising based on an extensive set of experiments, but we noticed a significant computational cost.

In this paper, we revisit ML-LOO and propose to support more coarse-grained features by using classical image segmentation techniques. We evaluate the segmentation techniques and compare to the original work in \cite{yang2020ml}.

The rest of the paper is organized as follows: in Section \ref{related} related work is surveyed. We revisit ML-LOO and present our implementation in Section \ref{revisit}. We present our "Segmentation Fault" approach in Section \ref{segfault}. Finally Section \ref{concl} concludes the paper and outlines future work.

\section{Related Work}
\label{related}

Neural networks, deep convolutional ones in particular, are currently the best machine learning algorithms at object classification. However, these models have been shown to be vulnerable to many adversarial deceptions.
By slightly changing the pixel values of an image, one can lead a convolutional neural network to make wrong predictions \cite{szegedy2013intriguing}. Attacks can be divided into model-specific (white box) and model-agnostic (black box). Model-specific assumes access to the gradients of the learning model \cite{kurakin2016adversarial}. Model-agnostic techniques transact with classifiers as black boxes and do not require internal knowledge of the model or the training data \cite{papernot2017practical}.
Some adversarial techniques lead to stealthy, unperceived to the human observer, samples. For example, changing only one pixel in the input image is shown to be enough in \cite{su2019one}. However, the feature attribution of this pixel will be significant. Recently, more interest is put on physical samples such as adversarial patches and adversarial 3D printed artifacts \cite{athalye2018synthesizing}.
We showed how adversarial patches may help shoplifting smart stores in \cite{nassar2019shoplifting}.
A recent review of adversarial attacks and defenses in the physical world is in \cite{renadversarial}. A nice summary of the history of adversarial machine learning is in \cite{biggio2018wild}. ML-LOO \cite{yang2020ml}, a recently proposed defense technique, shows applicability over a wide set of attacks and datasets. We carry on the same direction of ML-LOO and propose image segmentation to boost the performance of this approach.

\section{ML-LOO: Background \& Implementation}
\label{revisit}
The idea of Leave One Out (LOO) \cite{yang2020ml} is simply to take an $m \times n$ pixels image and sample $m \times n$ images out of it by blacking out pixels one by one. Sometimes, $m \times n \times 3$ are sampled by blacking R, G and B channels of each pixel. The sample images are forward-propagated through the neural network classifier in question. The output values of of a predefined set of network nodes are recorded. Several operational modes are available:
\begin{description}
    \item[1D] Only the probability of the predicted class for the original image is recorded:
    $f(x_i)_c$ where $x$ is the original image, $x_i$ is the image with  blacked pixel $i$, $f$ is the neural network representative function, and $c= \argmax_{j\in C} f(x)_j$. $C$ is the set of output classes.
    \item[Output layer] Only the output layer values are recorded.
    \item[Multi Layer (ML-LOO)] The values for the output layer and a randomly pre-selected set of hidden nodes from different intermediate layers are recorded.
\end{description}
The $m \times n$ recorded values for a node are collected. More precisely we record the absolute difference of the LOO values and the original single value $\phi(x)_{i,\text{node}_j} = \vert f(x_i)_{\text{node}_j} -  f(x)_{\text{node}_j} \rvert$. The original value corresponds to the value of the node taken from the forward propagation of the original image. This difference is called \textit{feature attribution}.
We abstract the group of these recorded differences for all pixels $i$ and a single node $j$ in one statistical measure, namely the Interquartile range (IQR).
$$ IQR (\phi(x)_{\text{node}_j}) = Q_{\phi(x)_{\text{node}_j}}(0.75) - Q_{\phi(x)_{\text{node}_j}}(0.25) $$
where: $$ Q_{\phi(x)_{\text{node}_j}}(p) = \min \{\beta:\frac{\#\{i :\phi(x)_{i,\text{node}_j}< \beta\}}{m \times n} \geq p \} $$
The IQR is useful for measuring the dispersion of feature attributions. The authors of \cite{yang2020ml} observed that adversarial samples cause a larger dispersion. The intuition is that in normal images a few "foreground" pixels have larger feature attributions compared to a majority "background" pixels. These pixels with high feature attribution serve as an explanation of the classifier decision. The feature attributions are more dispersed in adversarial images, which can be captured by a wider IQR value.  We are a bit more cautious in generalizing this observation as a common effect for all adversarial techniques.
Still, we postulate that in multi-dimensional IQR, interquartile values computed over neurons across the output layer and many hidden layers in the network, may serve as an interesting  feature vector representation of an image.  The IQR feature vectors for a large and varied labeled dataset of normal and adversarial samples can be used to train a simple machine learning binary classifier.
The binary classifier task is to take an input image, as represented by its IQR feature vector, and decides whether it is benign or malicious. Experiments show that ML-LOO efficiently filters the majority of adversarial samples.

We however notice that ML-LOO is not cheap in terms of computational performance as it requires $m \times n$ forward-propagations per image. We discuss our implementation of this proposed  technique and findings in this section and differ our approach to boost its computational speed to the next section.

We experimented with CIFAR10 dataset \cite{krizhevsky2009learning} consisting of 32x32 tiny images  (e.g.
\vcenteredinclude{pics/cifar10airplaneEx.png}, \vcenteredinclude{pics/cifar10CatEx.png},
\vcenteredinclude{pics/cifar10TruckEx.png}), 10 classes (airplane, automobile, bird, cat, deer, dog, frog, horse, ship, truck), and 60,000 examples
separated into training (50,000) and testing (10,0000) subsets. Our CIFAR10 pretrained classifier is a deep convolutional neural network with 26 layers and 309,290 parameters\footnote{Skeleton code from \url{https://www.kaggle.com/c/cifar-10/discussion/40237}}. The model has 94\% training accuracy and 89\% testing accuracy.  We used Foolbox \cite{rauber2017foolboxnative,rauber2017foolbox} to generate adversarial samples. Four attack algorithms have been evaluated: Fast Gradient Sign Method (FGSM) \cite{goodfellow2014explaining}, PGD \cite{madry2017towards}, DeepFool \cite{moosavi2016deepfool},  and C\&W \cite{carlini2017towards}.
The efficiency of these attacks for different parameters is shown in Table \ref{tab:adv}. For C\& W we left epsilons as None and shorten the number of steps to 500 for lessening the generation time.

\begin{table}[tbp]
\centering
\caption{Neural Networks Model Classification Accuracy for One Random Batch of 128 Adversarial Images \label{tab:adv}}
\begin{tabular}{|l|c|c|c|}
\hline
 Attack $\downarrow \backslash$ $\epsilon \xrightarrow{}$ & $\epsilon=0.02$ & $\epsilon=0.06$ & $\epsilon=0.1$\\
  \hline
    FGSM  & 41.41\%  & 13.28\% & 11.72\% \\
    PGD & 32.03\% & 1.56\% &  0.0\%\\
    DeepFool & 40.62\% & 11.72\%&   10.94\%\\
    \hline
    C\&W  &  \multicolumn{3}{|c|}{17.19\%}\\
    \hline
    No attack &  \multicolumn{3}{|c|}{95.31\%} \\
    \hline
\end{tabular}
\end{table}

We implemented 3 flavors of ML-LOO:
\begin{description}
\item [1D] Only taking the IQR of the predicted-class output node,
\item [10D]Taking the IQR vector for the 10 nodes of the output layer,
\item [4K-D] The IQR vector has 4K values: 200 nodes are randomly chosen from each layer in the set of the last 20 layers of the network (since the last layer is only 10 nodes, the real number is 3810 values).
\end{description}

A scatter plot for a batch of 128 original images and 128 adversarial images obtained by FGSM ($\epsilon=0.1$) is shown in Figure \ref{fig:IQR1D}. The attack was successful in reducing the classifier accuracy to around 10\% for the adversarial samples, making it no better than a random classifier with 10 classes. Scatter plots for the dimensions in IQR 10D are show in Figure \ref{fig:IQR10D}.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{pics/IQR1D.png}
    \caption{IQR-1D: Adversarial images show larger IQR values in average}
    \label{fig:IQR1D}
\end{figure}

\begin{figure*}[t]
    \centering
    \includegraphics[width=\textwidth]{pics/iqr10.png}
    \caption{IQR-10D: Adversarial images show larger IQR values in average for most of the dimensions}
    \label{fig:IQR10D}
\end{figure*}

Furthermore, our results for the three aforementioned IQR models, adversarial techniques and off-the-shelf binary classifiers are presented in Table \ref{tab:MLLOO_Res}. For each attack/experiment, we use 30 batches of benign images and their counterpart 30 batches of adversarial samples ($30*128=3840$ images) generated using any misclassification as criteria and multiple epsilons. Each third of adversarials (10 batches) is using a different epsilon: 0.02, 0.06, or 0.1. For C\&W we keep epsilons free to choose by the algorithm and limit the number of steps to 500. 80\% of the data were used for training the classifier and only 20\% for testing. The values shown correspond to the classifiers performance over the testing data.

\newcommand{\fcell}[1]{%
\FPeval{\grayshade}{(1.5-#1)*1.1}%
\xdef\tempa{\grayshade}%
\cellcolor[gray]{\tempa} #1}

\begin{table*}[t]
    \centering
    \caption{Off-the-shelf classifiers detection AUC for the different IQR representations}
    \begin{tabular}{|c|c|c|c||c|c|c||c|c|c||c|c|c|}
    \cline{2-13}
          \multicolumn{1}{c|}{}& \multicolumn{12}{|c|}{Attacks}  \\
         \cline{2-13}
         \cline{2-13}
          \multicolumn{1}{c|}{}& \multicolumn{3}{|c||}{\textbf{FGSM}} & \multicolumn{3}{|c||}{\textbf{PGD}} & \multicolumn{3}{|c||}{\textbf{DeepFool}} & \multicolumn{3}{|c|}{\textbf{C\&W}} \\
         \hline
         Classifier & \emph{1-D} & \emph{10-D} & \emph{4K-D} & \emph{1-D} & \emph{10-D} & \emph{4K-D} & \emph{1-D} & \emph{10-D} & \emph{4K-D} &\emph{1-D} & \emph{10-D} & \emph{4K-D} \\
         \hline
         \textbf{Random Forest} & \fcell{0.633} & \fcell{0.765} & \fcell{0.854} & \fcell{0.693} & \fcell{0.782} & \fcell{0.685} & \fcell{0.671} & \fcell{0.792} & \fcell{0.649} & \fcell{0.917} &  \fcell{0.929} & \fcell{0.660} \\
         \textbf{SVC} & \fcell{0.738} & \fcell{0.768} & \fcell{0.802} & \fcell{0.720} & \fcell{0.735}& \fcell{0.744} & \fcell{0.745} &\fcell{0.755} & \fcell{0.746}& \fcell{0.926} & \fcell{0.930 }& \fcell{0.819}\\
         \textbf{XGB} & \fcell{0.737}  & \fcell{0.810}& \fcell{0.838} & \fcell{0.740} & \fcell{0.785} & \fcell{0.786} & \fcell{0.750} &\fcell{0.811}  & \fcell{0.836}& \fcell{0.946} & \fcell{0.956} &  \fcell{0.956}\\
         \hline
    \end{tabular}
    \label{tab:MLLOO_Res}
\end{table*}

The best detection accuracies are recorded against the C\&W attack.
Results shows that increasing the IQR representation dimensionality improves the detection accuracy consistently when XGB is used, even though the random forest classifier and sometimes SVC badly handle the case of high dimensional 4K features.
The shown numbers belong to the default parameters of the classifiers. However, we noticed that these numbers can be significantly enhanced by a grid search for the best parameters, especially for the high-dimensional (4K) case. For example, a grid search for the best $C$ and $\gamma$ parameters of the SVC classifier showed an improvement of 3\% in detection accuracy for the FGSM attack and 4K-D features.

\section{Segmentation Fault}
\label{segfault}
It is worth noting that the time taken to compute the feature attribution values for a batch of 128 images is in the order of one minute for IQR-1D and IQR 4k-D. The time is a bit less for IQR-10D since no post-processing is required. This time is independent of whether the batch is benign or adversarial. We used a Tesla P100-PCIE-16GB GPU from Google Colab. At each iteration we black-out one pixel for all the images in the batch and forward propagate the batch at once. Also note the considerable memory requirements for IQR-4K-D which is in the order of \verb|BATCH_SIZE*W*H*4000|. For CIFAR10 we need around 2GB to store the feature attribution vectors for a batch of 128 images, before reducing it to IQR vectors of size 4000 each (\verb|BATCH_SIZE*4000|).
Our approach consists on improving upon these processing and memory requirements by using segmentation. The idea is to simply to black out one segment at a time, reducing the number of blackouts from $w \times h$ to \verb|nb_segments|, hence reducing the processing and memory requirements by a ratio of 10 to 100 in the case of CIFRA10 tiny images, and even more in the case of a larger image dataset. However, it is not clear which segmentation technique would be more successful on preserving a good detection accuracy, hence the motivation of our work.

We consider the following segmentation methods, and used their scikit-image implementation \cite{scikit-image}.

\begin{table}[htbp]
\centering
\caption{Examples of Felzenszwalb’s  Segmentation\label{tab:Felzen}}
\begin{tabular}{c|cccc}
  \hline
\textbf{Scale}   & 1 & 10 & 100 & 1000 \\
  \hline
  \textbf{Label} & \multicolumn{4}{c}{Truck} \\
  \textbf{\#Segments}  & 15 & 15 & 16 & 7 \\
  \textbf{Image} & \vcenteredinclude{pics/truck1.png} & \vcenteredinclude{pics/truck1.png} & \vcenteredinclude{pics/truck1.png}& \vcenteredinclude{pics/truck1.png} \\
  \hline
    \textbf{Label} & \multicolumn{4}{c}{Automobile} \\
  \textbf{\#Segments}  & 18 & 18 & 19 & 12 \\
  \textbf{Image} & \vcenteredinclude{pics/automobile1.png} & \vcenteredinclude{pics/automobile2.png} & \vcenteredinclude{pics/automobile3.png}& \vcenteredinclude{pics/automobile4.png} \\
  \hline
  \textbf{Label} & \multicolumn{4}{c}{Deer} \\
  \textbf{\#Segments}  & 11 & 11 & 11 & 5 \\
  \textbf{Image} & \vcenteredinclude{pics/deer1.png} & \vcenteredinclude{pics/deer2.png} & \vcenteredinclude{pics/deer3.png}& \vcenteredinclude{pics/deer4.png} \\
  \hline
\end{tabular}
\end{table}
\begin{table}[htbp]
\centering
\caption{ Examples of QuickShift  Segmentation\label{tab:quickshift}}
\begin{tabular}{|l|cc|cc|cc|}
\multicolumn{7}{c}{\textbf{\#Segments} } \\
\hline
 labels $\xrightarrow{}$ & \multicolumn{2}{c|}{Airplane} & \multicolumn{2}{c|}{Bird} & \multicolumn{2}{c|}{Ship}\\
 \hline
$ \sigma \downarrow \backslash m \xrightarrow{}$     & 10 & 20  & 10 & 20  & 10 & 20 \\
  \hline
    0  & 61& 12& 224& 24 & 170 & 25\\
    1 & 35& 12& 87& 19 & 56& 11 \\
    2 & 16& 6& 30 & 11 & 22& 6 \\
    \hline
\end{tabular}
\end{table}
\begin{table}[htbp]
\centering
\caption{Examples of SLIC Segmentation\label{tab:slic}}
\begin{tabular}{c|cccc}
  \hline
\textbf{n\_segments}  & 1 & 32 & 64 & 128  \\
  \hline
  \textbf{Label} & \multicolumn{4}{c}{Ship} \\
  \textbf{Actual \#Segments}  & 1 & 18 & 48 & 98 \\
  \textbf{Image} & \vcenteredinclude{pics/ship-slic1} & \vcenteredinclude{pics/ship-slic2} & \vcenteredinclude{pics/ship-slic3}& \vcenteredinclude{pics/ship-slic4} \\
  \hline
    \textbf{Label} & \multicolumn{4}{c}{Automobile} \\
  \textbf{Actual \#Segments}  & 1 & 13 & 33 & 79 \\
  \textbf{Image} & \vcenteredinclude{pics/auto-slic1.png} & \vcenteredinclude{pics/auto-slic2.png} & \vcenteredinclude{pics/auto-slic3.png}& \vcenteredinclude{pics/auto-slic4.png} \\
  \hline
  \textbf{Label} & \multicolumn{4}{c}{Frog} \\
  \textbf{Actual \#Segments}  & 1 & 19 & 38 & 91 \\
  \textbf{Image} & \vcenteredinclude{pics/frog1.png} & \vcenteredinclude{pics/frog2.png} & \vcenteredinclude{pics/frog3.png}& \vcenteredinclude{pics/frog4.png} \\
  \hline
\end{tabular}
\end{table}
\begin{description}
\item [Felzenszwalb’s Segmentation \cite{felzenszwalb2004efficient}] This algorithm has a single parameter called \textit{scale} (default value = 1). Higher scale means
less and larger segments. The actual segment size and number of segments can vary depending on local contrast. Some examples for random images from CIFAR10 are shown in Table \ref{tab:Felzen}.

\item[Quickshift Segmentation \cite{vedaldi2008quick}]
This algorithm computes a hierarchical segmentation over five dimensions (i.e. 3D color, width and height). It has two main parameters: \textit{sigma} (default value = 0) controls the scale of the local density approximation, and \textit{max\_dist} (default value = 10) selects a level in the hierarchical segmentation that is produced. Higher max\_dist and sigma means less segments. Some examples for random images from CIFAR10 and different values of max\_dist ($m$) and sigma ($\sigma$) are shown in Table \ref{tab:quickshift}.

\item[SLIC \cite{achanta2012slic}]
Simple Linear Iterative Clustering (SLIC) performs K-means clustering in the 5d space of color information (3 dimensions) and image location (x,y). The main parameter is n\_segments which corresponds to the number of centroids for K-means. Some examples for random images from CIFAR10 and different values of n\_segments and  are shown in Table \ref{tab:slic}.
\end{description}

\begin{table*}[t]
    \centering
    \caption{Off-the-shelf XGBoost classifier detection AUC for different IQR representations and different segmentation techniques}
    \begin{tabular}{|c|l|c|c|c||c|c|c||c|c|c||c|c|c|}
    \cline{3-14}
          \multicolumn{2}{c|}{}& \multicolumn{12}{|c|}{Attacks -- epsilons $=[0.02, 0.06, 0.1]$}  \\

         \cline{3-14}
         \cline{3-14}
          \multicolumn{2}{c|}{}& \multicolumn{3}{|c||}{\textbf{FGSM}} & \multicolumn{3}{|c||}{\textbf{PGD}} & \multicolumn{3}{|c||}{\textbf{DeepFool}} & \multicolumn{3}{|c|}{\textbf{C\&W}} \\
         \hline
         \multicolumn{2}{|c|}{Segmentation $\downarrow$}  & \emph{1-D} & \emph{10-D} & \emph{4K-D} & \emph{1-D} & \emph{10-D} & \emph{4K-D} & \emph{1-D} & \emph{10-D} & \emph{4K-D} &\emph{1-D} & \emph{10-D} & \emph{4K-D} \\
         \hline
         \multirow{4}{*}{\textbf{Felzenszwalb}} & scale = 1 &  \fcell{0.707} & \fcell{0.782} & \fcell{0.851} & \fcell{0.688} & \fcell{0.749} & \fcell{0.804} & \fcell{0.746} & \fcell{0.790}  & \fcell{0.822} & \fcell{0.877}&\fcell{0.936}  &\fcell{0.930} \\
        & scale = 10 & \fcell{0.705} & \fcell{0.797} & \fcell{0.846} & \fcell{0.692}  & \fcell{0.753} & \fcell{0.766} & \fcell{0.762}  & \fcell{0.816} & \fcell{0.838} &   \fcell{0.887} & \fcell{0.927}&\fcell{0.916}  \\
          & scale = 100 &\fcell{0.695}  & \fcell{0.814} &\fcell{0.834}  &  \fcell{0.688} & \fcell{0.754} & \fcell{0.795} & \fcell{0.762}  & \fcell{0.811}  & \fcell{0.807} & \fcell{0.887}& \fcell{0.926} & \fcell{0.925} \\
          & scale = 1000 & \fcell{0.632} & \fcell{0.782} & \fcell{0.830} & \fcell{0.593} & \fcell{0.665} & \fcell{0.763} & \fcell{0.672} &\fcell{0.735}  & \fcell{0.775} &\fcell{0.797} & \fcell{0.825}& \fcell{0.825}\\
         \hline
         \hline
         \multirow{6}{*}{\textbf{QuickShift}} & $\sigma=0, m=10$  & \fcell{0.726} & \fcell{0.817} & \fcell{0.853} & \fcell{0.744} & \fcell{0.781} & \fcell{0.771} & \fcell{0.776} & \fcell{0.839} & \fcell{0.835} &  \fcell{0.944}& \fcell{0.946} & \fcell{0.940}\\
          & $\sigma=0, m=20$  & \fcell{0.694} & \fcell{0.797} &  \fcell{0.842} & \fcell{0.703} & \fcell{0.782} & \fcell{0.782} & \fcell{0.761} & \fcell{0.808} &  \fcell{0.827} &  \fcell{0.923}& \fcell{0.937} & \fcell{0.932} \\
          & $\sigma=1, m=10$ & \fcell{0.738} & \fcell{0.798} & \fcell{0.829} & \fcell{0.711} & \fcell{0.776} & \fcell{0.765} & \fcell{0.770} & \fcell{0.810} & \fcell{0.801} & \fcell{0.897}& \fcell{0.923}&  \fcell{0.923}\\
          & $\sigma=1, m=20$ & \fcell{0.710} & \fcell{0.816} & \fcell{0.849} & \fcell{0.637} & \fcell{0.759} & \fcell{0.767} & \fcell{0.753} & \fcell{0.805} & \fcell{0.792} &\fcell{0.877} & \fcell{0.941}& \fcell{0.925}\\
         & $\sigma=2, m=10$ &  \fcell{0.720} & \fcell{0.808} & \fcell{0.849} & \fcell{0.664} & \fcell{0.736}   & \fcell{0.757} &\fcell{0.770} & \fcell{0.806}  & \fcell{0.834} & \fcell{0.895} & \fcell{0.922} & \fcell{0.925} \\
         & $\sigma=2, m=20$ & \fcell{0.692} & \fcell{0.799}& \fcell{0.839} & \fcell{0.615} & \fcell{0.694} & \fcell{0.744}  & \fcell{0.706} & \fcell{0.774} & \fcell{0.790} & \fcell{0.833}& \fcell{0.897}& \fcell{0.879} \\
         \hline
         \hline
          \multirow{4}{*}{\textbf{SLIC}} & n\_segments = 32 &  \fcell{0.727} & \fcell{0.814} & \fcell{0.864} &\fcell{0.682}  & \fcell{0.755} & \fcell{0.755} &  \fcell{0.774} &\fcell{0.805}   & \fcell{0.805} &\fcell{0.895} &\fcell{0.934} & \fcell{0.923}\\
        & n\_segments = 64 & \fcell{0.730} & \fcell{0.809} & \fcell{0.865}  & \fcell{0.723} & \fcell{0.758} & \fcell{0.794} & \fcell{0.782} & \fcell{0.826} & \fcell{0.820} &  \fcell{0.929} & \fcell{0.947}& \fcell{0.941} \\
          & n\_segments = 128 &  \fcell{0.761}&  \fcell{0.812}& \fcell{0.855} & \fcell{0.733} & \fcell{0.773} & \fcell{0.777} & \fcell{0.773} & \fcell{0.826} & \fcell{0.831} & \fcell{0.945}& \fcell{0.959}& \fcell{0.946}\\
         \hline
         \hline
    \multicolumn{2}{|c|}{\textbf{ML-LOO}} & \fcell{0.728} & \fcell{0.829}& \fcell{0.871}& \fcell{0.748}& \fcell{0.793}& \fcell{0.806}& \fcell{0.793}&\fcell{0.829} & \fcell{0.763}& \fcell{0.918}& \fcell{0.967} & \fcell{0.984} \\
    \hline
    \end{tabular}
    \label{tab:Seg_Res}
\end{table*}

The results of experimentation with different IQR representations, different segmentation techniques and different adversarial techniques are shown in Table \ref{tab:Seg_Res}. For each attack/experiment, we use 30 batches of benign images and 30 batches of adversarial samples  ($30*128=3840$ images) generated using any misclassification as criteria and multiple epsilons (0.02, 0.06, and 0.1). 80\% of the data were used for training the classifier and only 20\% for testing. The XGB classifier is used off-the-shelf for these experiments. The values shown correspond to the XGB performance over the testing data.
It can be observed that generally 4K-D improves over 10-D and 10-D improves over 1-D. In general, the accuracy values found are very similar to what can be achieved with ML-LOO, even though slightly lower. The more segments are generated, the more accuracy is obtained, yet SLIC segmentation method looks superior in the sense that it can obtain better accuracy given that equal average number of segments is generated.

The time and memory comparison of ML-LOO versus our approach is shown in Table \ref{tab:comparison}. The time value corresponds to the time for extracting the feature attributions for one batch of 128 images. The size memory value is the size of the object for storing all the feature attributions for one batch of 128 images before being reduced to IQR values. The order of magnitude of the values are represented in grayscale. For time, darker cells correspond to longer processing times. For memory, darker cells correspond to larger memory space required. Results show that for most parameters, our segmentation approach requires much less memory and runs much faster than ML-LOO. We notice that 10-D is consistently the fastest to execute since it does not require any post-processing beyond forward-propagation through the classifier. 1-D has an additional step for fetching the value corresponding to the predicted class for the original image. 4K-D has to loop through the layers to fetch the values of the selected nodes. In terms of memory, 10-D requires 10 times more memory than 1-D and 4K-D requires 381 times more memory than 10-D.

Figure \ref{scatter1} shows an example of the performance of different segmentation techniques and their AUC scores. The scatter plot corresponds to the experiments with the C\&W attack and the 10-D representation. Each point represents a segmentation technique and particular parameters.
Results show that SLIC segmentation is generally better in reducing processing time while scoring AUC scores that are practically the same as those of ML-LOO.
\begin{figure}[t]
\centering
\begin{tikzpicture}
\begin{axis}[
    legend pos=north west,
    ylabel=Time(s),
    xlabel=AUC,
	enlargelimits=0.05
]
\addplot[
    scatter/classes={Felzenszwalb={blue}, QuickShift={green}, SLIC={black}, MLLOO={red}},
    only marks,
    scatter,
    mark=halfcircle*,
    mark size=4pt,
    scatter src=explicit symbolic,

    ]
table[meta=Class]{scattered.dat};
\legend{Felzenszwalb, QuickShift, SLIC, ML-LOO};
\end{axis}
\end{tikzpicture}
\caption{Trade-off Time vs. Detection Accuracy \label{scatter1}}
\end{figure}

Our ML-LOO implementation and the source code for all the experiments in this paper are available in Google Colab notebooks format at \url{https://github.com/mnassar/segfault}.

\begin{table*}[t]
\newcommand{\fcellsec}[1]{%
\FPeval{\grayshade}{(100-#1)/100}%
\xdef\tempa{\grayshade}%
\cellcolor[gray]{\tempa} #1s}

\centering
\caption{Comparison of time and memory requirements for representing one batch of 128 images \label{tab:comparison}}
       \begin{tabular}{|c|l|c|c|c||c|c|c|}

         \cline{3-8}
         \cline{3-8}
          \multicolumn{2}{c|}{}& \multicolumn{3}{|c||}{\textbf{Time}} &  \multicolumn{3}{|c|}{\textbf{Memory}}\\
         \hline
         \multicolumn{2}{|c|}{Segmentation $\downarrow$}  & \emph{1-D} & \emph{10-D} & \emph{4K-D}  & \emph{1-D} & \emph{10-D} & \emph{4K-D}\\
         \hline
         \multirow{4}{*}{\textbf{Felzenszwalb}} & scale = 1 & \fcellsec{1.79} & \fcellsec{1.17} & \fcellsec{1.77}  & \cellcolor[gray]{0.9} 9.8KB & \cellcolor[gray]{0.9} 98KB & \cellcolor[gray]{0.8} 37MB \\
        & scale = 10 & \fcellsec{2.06} & \fcellsec{1.27} & \fcellsec{2.01}  & \cellcolor[gray]{0.9} 10.8KB & \cellcolor[gray]{0.9} 108KB & \cellcolor[gray]{0.8} 41MB \\
          & scale = 100 & \fcellsec{2.24}  & \fcellsec{1.63} & \fcellsec{2.21} & \cellcolor[gray]{0.9} 12.4KB & \cellcolor[gray]{0.9} 124KB & \cellcolor[gray]{0.8} 47MB  \\
          & scale = 1000 &  \fcellsec{1.09}& \fcellsec{0.7} & \fcellsec{1.07} & \cellcolor[gray]{0.9} 5.2KB& \cellcolor[gray]{0.9} 52KB & \cellcolor[gray]{0.8} 19.8MB\\
         \hline
         \hline
         \multirow{6}{*}{\textbf{QuickShift}} & $\sigma=0, m=10$  & \fcellsec{44.8} & \fcellsec{28.4} & \fcellsec{48} & \cellcolor[gray]{0.9} 270KB & \cellcolor[gray]{0.8} 2.7MB & \cellcolor[gray]{0.5} 1GB \\
          & $\sigma=0, m=20$  & \fcellsec{13.9} & \fcellsec{8.86} & \fcellsec{13.8} & \cellcolor[gray]{0.9} 69KB & \cellcolor[gray]{0.9} 690KB &\cellcolor[gray]{0.8}  263MB \\
          & $\sigma=1, m=10$ & \fcellsec{30.5} & \fcellsec{18.7} & \fcellsec{30.4} & \cellcolor[gray]{0.9} 173KB & \cellcolor[gray]{0.8} 1.7MB & \cellcolor[gray]{0.8} 659MB   \\
          & $\sigma=1, m=20$ &  \fcellsec{7.77} & \fcellsec{5.41}  & \fcellsec{7.29} & \cellcolor[gray]{0.9} 32KB & \cellcolor[gray]{0.9} 317KB &  \cellcolor[gray]{0.8} 120MB\\
         & $\sigma=2, m=10$ & \fcellsec{21.3} & \fcellsec{13.2} & \fcellsec{20.9} & \cellcolor[gray]{0.9} 117KB  & \cellcolor[gray]{0.8} 1.2MB  & \cellcolor[gray]{0.8}  446MB \\
         & $\sigma=2, m=20$ & \fcellsec{4.3} & \fcellsec{3.46} & \fcellsec{4.3}  & \cellcolor[gray]{0.9} 13KB & \cellcolor[gray]{0.9} 128KB & \cellcolor[gray]{0.8}  49MB\\
         \hline
         \hline
          \multirow{4}{*}{\textbf{SLIC}} & n\_segments = 32 &  \fcellsec{2.17} & \fcellsec{1.42} & \fcellsec{3.24} &  \cellcolor[gray]{0.9} 11KB & \cellcolor[gray]{0.9} 113KB & \cellcolor[gray]{0.8}  43MB \\
        & n\_segments = 64 & \fcellsec{5} & \fcellsec{3.26} & \fcellsec{5.14} & \cellcolor[gray]{0.9} 28KB & \cellcolor[gray]{0.9} 287KB & \cellcolor[gray]{0.8} 109MB \\
          & n\_segments = 128 & \fcellsec{9.75} & \fcellsec{5.96}  & \fcellsec{9.98} &\cellcolor[gray]{0.9} 58KB & \cellcolor[gray]{0.9} 579KB&  \cellcolor[gray]{0.8}  220MB\\
         \hline
         \hline
         \multicolumn{2}{|c|}{\textbf{ML-LOO}} & \cellcolor[gray]{0.4} \texttildelow 1min20s &  \texttildelow \fcellsec{45}  &\cellcolor[gray]{0.4} \texttildelow 1min22s &  \cellcolor[gray]{0.9} \texttildelow 524 KB & \cellcolor[gray]{0.8}  \texttildelow 5.24 MB & \cellcolor[gray]{0.4}  \texttildelow 2GB  \\
         \hline
    \end{tabular}

\end{table*}

\section{Conclusion}
\label{concl}
ML-LOO approach is a simple defense to wipe out a majority of adversarial samples that was extensively tested against  different attack techniques and mixtures of parameters.   Nevertheless, the computational cost and memory requirements cannot be disregarded. In this paper, we proposed a resource optimization approach based on segmentation. Our approach permits a considerable improvement in resource efficiency while maintaining the original main property of ML-LOO: a simple and efficient adversarial detection that quickly wipes out most adversarial samples. In future work, we aim at applying our approach and experiment on other datasets, image datasets with higher resolution in particular. We also want to experiment with more adversarial techniques and study the transferability properties of our approach.
\section*{Acknowledgements}
This work was supported partially by a grant from the university research board of the American  university of Beirut (URB-AUB-2020/2021). The authors would like to thank Joshua Peter Handali, Johannes Schneider and Pavel Laskov from University of Liechtenstein for the early discussion of the paper topic.

\bibliographystyle{ACM-Reference-Format}