\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. 


