\section{Polymer in a double well
\label{sec:doublewell}}
We consider now the case in which the external potential $V(x)$ in Eq.~\eqref{eq:motion} is a double-well. This 
simple model allows us to study some properties of thermally activated barrier 
crossing as, for instance, the dependence of the jump rate $r$ on the physical 
parameters.

The general problem of activated dynamics has been extensively studied since the seminal works
by Kramers, and the reaction-rate theory provides many analytic methods to
compute jump times in different contexts (see Ref.~\cite{hanggi90} and reference therein).
Important results have been derived also for polymeric chains~\cite{poly_su_Barr,poly_dwell, sebastian06}.
We want to stress that our aim here is not to improve such results: we are interested in
analyzing the ability of a data-driven approach to reconstruct an activated dynamics.
In particular, we will focus on the role of 
the chain stiffness on the activated dynamics and, more importantly, 
on its relevance for the description in terms of the RCs.




The external potential reads, in this case,
\begin{equation}
V(x)=\frac{B^2}{4}(x^2-A^2)^2
	\label{eq:biwell}
\end{equation}
where $A$ and $B$ are suitable constants. The typical dynamics of the center of mass, $Q(t)$,
is therefore characterized by jumps over the barrier separating 
the two minima of the potential (two-state model).  For the interaction potential we choose the form 
$U(r)= K(r - \sigma)^2/2$. 
\\
\iffalse
The description of such transition dynamics strongly depends 
on the bond rigidity $K$. We can expect that in the limit of high
$K$ (semi-rigid rod), the dynamics is fairly reproduced by model (M1). 
As soon as the internal interactions become comparable to the one induced by 
the external potential, a single-variable description fails
and one is forced to consider at least the two-dimensional effective 
dynamics (M2).
\fi

In the following, we will always consider the limit $N\sigma \simeq A$, 
i.e. the case in which the equilibrium length of the polymer is 
comparable to the half distance between the well minima: it is reasonable to 
expect that in these conditions the value of the bond rigidity affects the 
qualitative behavior of the chain in a relevant way.

Our aim is to show that for high values of $K$ the model described by Eq.~\eqref{eq:M1}
suffices to reproduce the quantitative macroscopic behaviour of the system; whereas,
as soon as $K$ becomes comparable to $B^2A^2/\sigma$, the evolution of $Q$ is no more
Markovian and any attempt to describe it through model~\eqref{eq:M1} is doomed to fail.
However, if the phase-space is expanded by including a suitable additional RC, it is
still possible that the evolution of the new RCs vector turns out to be Markovian, so that
a dynamical description based on Eq.~\eqref{eq:M2} can be accurate enough.

The validity of such scenario can be tested by using the data-driven 
approach mentioned in the Introduction and detailed in Appendix~\ref{sec:appendixLE}. 
We first perform numerical simulations of the whole system by using a Stochastic 
Runge-Kutta algorithm~\cite{honeycutt92} and measuring the relevant
RCs of the system at every time step. As a first attempt, from long time-series of
such data we build an effective stochastic equation for $Q$ only,
in the form of Eq.~\eqref{eq:M1}; then we apply the extrapolation procedure to
the dynamics of the two-dimensional vector $(Q,L)$, obtaining an M2-like model.
The ``goodness'' of M1 and M2 is tested
by measuring the Kramers' transition times of the reconstructed 
models and comparing the corresponding jump rates to the original ones. 




\subsection{1-variable model}

Before applying the mentioned extrapolation method to infer numerically
the functional form of the terms appearing in Eq.~\eqref{eq:M1}, let us derive 
analytically an effective equation for $Q$ for the 
high-stiffness limit, $K\gg B^2A^2/\sigma$. In this case we can assume that
the position of two consecutive beads is fixed and equal to $\sigma$. Due to the 
simple form of the external potential $V(x)$, the drift term in 
Eq.~\eqref{eq:Qdot1} can be exactly computed in this case:
\begin{equation}
\begin{aligned}
 -\frac{1}{N} \sum_{i=1}^N V'(Q+u_i)&= -\frac{B^2}{N} \sum_{i=1}^N \left[ (Q+u_i)^3 -A^2 (Q+u_i)\right]\\
 &= -\frac{B^2}{N} \sum_{i=1}^N \left[ 3Qu_i^2+Q^3-A^2Q\right]\,.
\end{aligned}
\end{equation}
where we have used the fact that, due to the rigidity of the polymer, $\sum_i u_i^3=\sum_i u_i=0$.

Now we substitute the explicit expression for the relative positions of the polymer beads, $u_i=(2i-N-1)\sigma/2$, and after straightforward algebra we get
\begin{equation}
\begin{aligned}
 -\frac{1}{N}\sum_{i=1}^N V'(x_i)&=-B^2 Q\left(Q^2-A^2+\frac{\sigma^2}{4}(N+1)(N-1)\right)\\
 &=-B^2 Q\left(Q^2-A_{\text{eff}}^2\right)
\end{aligned}
\end{equation}
with 
\begin{equation}
 A_{\text{eff}}=A\sqrt{1-\frac{L^2(N+1)}{4A^2(N-1)}}\,,
\end{equation}
where we have used the definition of the polymer length for the rigid case, $L=(N-1)\sigma$.
Let us notice that the above drift corresponds to an effective potential
\begin{equation}
 V_{\text{eff}}(Q)=\frac{B^2}{4}(Q^2-A_{\text{eff}}^2)^2\,,
 \label{eq:poteff}
\end{equation}
i.e. a ``rescaled'' version of the original external potential~\eqref{eq:biwell}.
\iffalse
Let us also remark that, in the limit
$N\gg1$, Eq.~\eqref{eq:closed} (see~\ref{sec:appendixHCA}) is exactly recovered.
\fi
\\
We can now use the theory of escape times~\cite{gardiner} to estimate the jump rate $r$
for the effective potential~\eqref{eq:poteff}. The average waiting time between two consecutive jumps can be computed through the formula~\cite{gardiner}
\begin{equation}
\label{eq:tau_rigid}
\begin{aligned}
 \tau=&\frac{1}{D} \int_{-A}^{A}dy\, e^{V(y)/D} \int_{-\infty}^{y}dz\, e^{-V(z)/D}\\
 =&\frac{1}{D} \int_{-A}^{A}dy\, e^{\frac{B^2}{4D}\left[y^4-2A_{\text{eff}}^2y^2\right]}\int_{-\infty}^{y}dz\, e^{-\frac{B^2}{4D}\left[z^4-2A_{\text{eff}}^2z^2 \right]}
\end{aligned}
\end{equation}
where $D=T/N$ is the diffusivity associated to $Q$. The jump rate $r$ is then found as
\begin{equation}
 r=\frac{1}{\tau}\,.
\end{equation}
The above equations, which are only valid in the rigid-rod limit,
will be a useful touchstone to evaluate the level of accuracy of
the model inferred numerically. 

\begin{figure}
\includegraphics[width=.99\linewidth]{FIG2.pdf}
\caption{\label{fig:potential}
(A): drift term $F(Q)$ of model~\eqref{eq:M1} as reconstructed from data 
(points) and fitted with a 9-th degree odd polynomial (solid lines), for three 
different values of $K$. (B): Effective potential obtained by integration of 
$F(Q)$. Parameters for the 
simulations of the complete system: $A=10$, $B=0.1$, $T=30$, $\sigma=1$, $N=10$, using 
a time-step $dt=10^{-5}$. Simulations on model~\eqref{eq:M1} have been run 
with a time-step $dt=10^{-4}$.} \end{figure}

\begin{figure}
\includegraphics[width=.99\linewidth]{FIG3.pdf}
\caption{\label{fig:markov}
 Checking the validity of model~\eqref{eq:M1}. (A) Autocorrelation $\langle \zeta(t)\zeta(t+\tau) \rangle$ of the ``noise'' term~\eqref{eq:zeta}
vs. $\tau$. The inset is a zoom. It is apparent the fast decay with $\tau$ that can be 
considered effectively a delta-correlation, since the time-step of the
integration algorithm generating the dynamics is $dt=10^{-5}$. Here, $K=14$.
(B): Jump rates as functions of $K$ in the original system (red circles) and in simulations 
of the reconstructed M1 (blue diamonds). The rigid-rod 
approximation (dashed line) is reported as reference. Parameters as in Fig.~\ref{fig:potential}.} 
\end{figure}

We apply now the extrapolation procedure mentioned in the Introduction and 
detailed in Appendix~\ref{sec:appendixLE} to infer the right functional forms for 
the terms of the Langevin equation~\eqref{eq:M1}, assuming that the dynamics of $Q$
is Markovian and a 1-variable description in the form of model~\eqref{eq:M1} holds.
We find that $D$ (not shown here) is always almost constant and equal to $T/N$,
as it would be expected if the process was Markovian, while the drift shows a 
more complex shape (Fig.~\ref{fig:potential}A); we fit the data by a 9-th degree, 
odd polynomial, then we integrate the resulting function in order to get an 
effective potential, which is reported in Fig.~\ref{fig:potential}B for several 
values of $K$. In the large-$K$ limit, as expected, we recover a quartic 
effective potential: terms of higher order become relevant when the bond stiffness is 
low, and their effect is to flatten the potential barrier between the two wells.

As mentioned above, the validity of Eq.~\eqref{eq:M1} relies on the assumption that
the evolution of $Q$ is Markovian, which has to be checked.
First, one can define and measure the following quantity:
\begin{equation}
\label{eq:zeta}
 \zeta(t)=\dot{Q}(t)-F(Q(t))\,,
\end{equation}
which represents the ``noise'' of Eq.~\eqref{eq:M1}, if the dynamics of $Q$ is Markovian. We can compare the autocorrelation time of $\zeta(t)$ and verify that it is much shorter than any characteristic time-scale of the dynamics of $Q$. In our case, $\zeta(t)$ always decorrelates on the scale of the time-step of the integration algorithm, $dt$ (see Fig.~\ref{fig:markov}A).

This first check assures that there is a clear time-scale separation between the dynamics of the center of mass and its ``noise''. However, this does not imply that the original dynamics of $Q$ has to be Markovian: in order to check that, we also have to verify the consistency with the original dynamics. If the one-variable description is able to catch the relevant features of the whole system, we can conclude, \textit{a posteriori}, that the evolution of $Q$ was Markovian also in the complete dynamics; if not, a different description has to be taken into account.

Fig.~\ref{fig:markov}B shows, for several values of 
the rigidity, the jump rates measured in the original dynamics and those 
observed in the reconstructed model, using a standard stochastic integration 
algorithm (the one discussed in~\cite{mannella-palleschi89}, up to order 
$dt^{3/2}$). In the high-$K$ limit the simple rigid-rod 
approximation~\eqref{eq:tau_rigid} holds, there is no dependence on $K$ and the 
agreement between the jump rates of M1 and of system~\eqref{eq:motion} is excellent. As the polymer 
becomes softer, even if a significant improvement on Eq.~\eqref{eq:tau_rigid} 
can be observed, the relative error between M1 and the true dynamics exceeds 
30\%: this is a clear hint that a 1-variable description, even if inferred 
directly from data, cannot reproduce all the relevant features of the dynamics.
This is due to the fact that our implicit assumption on the Markovianity of the
process is wrong.
 


\subsection{2-variables model}


The failure of model~\eqref{eq:M1} for small values of $K$,
revealed by the discrepancies between the reconstructed and the 
original jump rate, suggests to go beyond a 
single variable description in order to achieve satisfactory results.
As discussed in Section~\ref{sec:system}, a reasonable attempt to recover
a Markovian dynamics is to consider the elongation of the polymer, $L$,
as a second $RC$ for our model, and we postulate the validity of an 
evolution equation of the form~\eqref{eq:M2}.
\iffalse
The same also happens under the simplifying
assumption of homogeneous deformations, which makes the problem exactly 
solvable, see Appendix~\ref{sec:appendixHCA}.
\fi

The requirement of a variable accounting for the elongation of the polymer
can be easily understood by looking at Fig.~\ref{fig:scatter} reporting
three scatter plots of the original dynamics in the $(Q,L)$ plane,
for different bond rigidity, where $L = x_N - x_1$.
When $K$ is high, and the system is well approximated by the rigid-rod model,
the region of the phase space explored by the dynamics is a thin strip
around the equilibrium value $L\simeq (N-1) \sigma$. 
As soon as the rigidity condition
is relaxed, and the system is allowed to vary its length in a significant way,
a two-lobe distribution takes place: $L$ tends to be smaller than the rest 
length of the chain when the polymer occupies one of the two minima of the double-well
potential, while it significantly increases during the transition across the barrier.
This particular shape of the scatter plot indicates that 
the typical pathways in the space $(Q,L)$ include a non-negligible 
deformation in $L$, which can be straightforwardly interpreted as follows: when the rigidity $K$ 
is low, the transition across the barriers of the polymer
occurs with a concomitant stretching of the bonds, presumably those that instantaneously lay on top the barrier. 
As a consequence, any Markovian effective description involving only the center of mass 
is completely insufficient to fairly approximate the dynamics of the system. 
\begin{figure}[h!]
\includegraphics[width=0.5\linewidth]{FIG4.png}
\caption{\label{fig:scatter}
Scatter plot of the end-to-end distance $L$ vs. $Q$ (total integration time: $10^6$). Green dots: $K=4$; blue dots: $K=20$; red dots: $K=600$. Other parameters as in Fig.~\ref{fig:potential}.}
\end{figure}
Following again the strategy discussed in Appendix~\ref{sec:appendixLE}, 
we provide numerical values for $F_Q$, $D_Q$, $F_L$ and $D_L$ in the $(Q,L)$ 
space,
which have to be fitted using suitable functional forms. Due to the symmetries
of the system, $F_Q(Q,L)$ has to be odd with 
respect to the variable $Q$, while $F_L(Q,L)$ should be even. 
Fig.~\ref{fig:fit_3d} shows the results obtained by fitting the following 
polynomial form:
\begin{equation}
\begin{aligned}
F_Q(Q,L)& = Q\left[c_{10}^{(Q)}   + c_{12}^{(Q)}L^2   + c_{13}^{(Q)}L^3\right] \\
	& + Q^3\left[c_{30}^{(Q)} + c_{32}^{(Q)}L^2 + c_{33}^{(Q)}L^3\right]\\
F_L(Q,L)& = c_{00}^{(L)}    + c_{01}^{(L)}L      + c_{03}^{(L)}L^3
	  + c_{21}^{(L)}Q^2L\,.
\end{aligned}
\label{eq:FQ_FL}
\end{equation}

\begin{figure}[h!]
\includegraphics[width=.8\linewidth]{FIG5.pdf}
\caption{\label{fig:fit_3d}
	Reconstructed drift terms of $Q$ and $L$ in the model \eqref{eq:M2}, 
	case $K=4$. Points are extrapolated from data (see Appendix~\ref{sec:appendixLE}); the 
	surface is obtained by fitting the polynomial~\eqref{eq:FQ_FL}. Other parameters are as in Fig.~\ref{fig:potential}.}
\end{figure}

The agreement between the actual data and the proposed functional form is good enough
to hope that the guessed model catches the most relevant features of the dynamics.
The diffusivity terms $D_Q$ and $D_L$ are again fitted by constant functions.
Once model~\eqref{eq:M2} is determined, 
we can check the reliability of its stochastic evolution
by a direct comparison with the original dynamics.

\begin{figure}[h!]
\includegraphics[width=.5\linewidth]{FIG6.pdf}
\caption{\label{fig:pdfs}
Pdfs in the 2-variable phase space of the original system, first column, 
and in the reconstructed 2-variables model~\eqref{eq:M2}, second column. 
Stiffness: $K=4$ (top), $K=20$ (center), $K=600$ (bottom).  Other parameters
are the same as Fig.~\ref{fig:potential}.}
\end{figure}

A first, important benchmark is given by the ability of the model to reproduce 
the static properties of the system, namely the joint probability distribution 
in the $(Q,L)$ space. This test is reported in Fig.~\ref{fig:pdfs} for different 
values of $K$, showing a reasonable qualitative agreement even in the non-trivial
case of low bond stiffness: in particular, the stretching occurring when the polymer
crosses the barrier is clearly reproduced.
\begin{figure}[h!]
\includegraphics[width=0.5\linewidth]{FIG7.pdf}
\caption{\label{fig:rates}
Relative errors of the jump rates $r$ in the two reconstructed models,
varying $K$. For other parameters, see caption of Fig.~\ref{fig:potential}.}
\end{figure}
The improvement of our effective description when also $L$ is taken into 
account can be fully appreciated by looking at dynamical observables
as the jump rate $r$. 
Fig.~\ref{fig:rates} displays the relative errors
between the values of $r$ obtained in the reconstructed models M1, M2 and
in the original dynamics. As already discussed, 
the 1-variable model fails when the polymer is soft, while the accuracy
of M2 does not seem to be affected in this limit.
Let us notice, on the other hand, that for $K\gg A^2B^2/\sigma$ the reconstructed model~\eqref{eq:M2} is less reliable than the 1-variable version: this is probably
a consequence of the larger number of parameters involved, which leads to a lower
degree of precision on their determination with the discussed method.

\subsection{Remarks on the structure of the effective equations}
The procedure for the reconstruction of a model in the form~\eqref{eq:M2} that 
we used in the previous subsection is based on a ``dynamical'' analysis, which 
builds the coefficients of the guessed stochastic model by looking at the 
time evolution of suitable observables of the original system (see 
Appendix~\ref{sec:appendixLE}). 
One could wonder whether such procedure is really 
needed in order to get a realistic description of the studied process; for 
example, in analogy with many statistical mechanical problems, one may expect
that the following recipe works:
\begin{enumerate}
\item measure the stationary p.d.f. $\rho(Q,L)$ from long time series of data;
\item deduce an effective 2-variable ``potential'', also known as ``potential
of mean force'' in chemical and biophysical contexts:
\begin{equation}
\label{eq:pot_stat}
W_S(Q,L)= -\log[\rho(Q,L)]\,;
\end{equation}
\item define $F_Q\equiv-c\partial_Q W_S$ and $F_L\equiv-c\partial_L 
W_S$, where the constant $c$ has to be determined.
\end{enumerate}

Let us notice that in many practical situations the multidimensional free-energy landscapes obtained by simulations or experiments are assumed to be generated by a system with a gradient (or gradient-like) structure and, using this hypothesis, Langer's formula~\cite{langer69,hanggi90} is applied to derive the transition rates over saddles.

The above procedure is a completely ``static'' analysis,
because it involves only quantities measured in equilibrium conditions.
Leaving apart the problem of finding the right multiplicative constant $c$
and the noise terms, this approach has a 
major issue: it is not sure, \textit{a priori}, whether the dynamics of 
the chosen reaction coordinates can be described by a potential, 
even if the complete system actually can; in other words,
there is in general no reason to expect that the model is a gradient system in the reduced phase space.
\iffalse
This gradient structure, however, is not guarateed \textit{a priori} and it needs to be carefully checked: the case discussed in this section is a clear example in which this assumption fails and one cannot rely on the mentioned approach.
\fi

The knowledge of $W_S$ alone is not sufficient to obtain the drift; additional informations on the dynamics need to be taken into account. A possibility is the theoretical approach discussed in Ref.~\cite{grabert06}, which also involves the transport coefficients. Our numerical method, instead, exploits the dynamical information by computing suitable conditioned moments of the RCs, as discussed in Appendix~\ref{sec:appendixLE}.

\begin{figure}[h!]
\includegraphics[width=.5\linewidth]{FIG8.pdf}
\caption{\label{fig:stat}
Static analysis. (A): effective potential $W_S(Q,L)$ from the empirical pdf in 
the $(Q,L)$ space. Panels below compare the drift computed with this static approach,
$\nabla W_s$, and the one reconstructed with our dynamical approach, $(F^{(2)}_Q, F^{(2)}_L)$,
to the measured coefficients $F_Q, F_L$. (B): $|\partial_Q W_s/F_Q|$. (C): 
$|F^{(2)}_Q/F_Q|$. (D): $|\partial_L W_s/F_L|$. (E): $|F^{(2)}_L/F_L|$.} 
\end{figure}
In order to show that in our case the simplifying assumption of a gradient structure 
for the $(Q,L)$ dynamics is wrong, let us determine $W_S$ using Eq.~\eqref{eq:pot_stat}. 
In Fig.~\ref{fig:stat}A we fit such ``potential'' with a 2-variable polynomial
(4-th order in both $Q$ and $L$), getting a nice superposition.

If the 2-variables system were gradient, $F_Q(Q,L)$ should be equal to $c\partial_Q W_S$,
for some value of $c$; therefore the ratio $\partial_Q W_S/F_Q$
should be constant all over the phase-space. Fig.~\ref{fig:stat}B shows instead that, in our case,
such function strongly depends on $Q$ and $L$. 
Just for comparison, Fig.~\ref{fig:stat}C displays the ratio 
$$
\epsilon =  F^{(2)}_Q/F_Q
$$ 
where $F^{(2)}_Q$ represents here the expression from the fitting 
\eqref{eq:FQ_FL}: not surprisingly, this function
is constant and equal to one almost everywhere, meaning only that we have made
sensible choices for the polynomial functions to use in the fit.
In Fig.~\ref{fig:stat}D,~\ref{fig:stat}E  the same comparison is done for the drift of the 
end-to-end distance $L$.

This simple check clearly shows that the simplified 2-variables system is not gradient,
and therefore a static analysis
is not sufficient to infer reasonable effective equations for its evolution.



\section{Conclusions}
\label{sec:conclusions}
We studied the dynamics of a 1-dim polymer whose monomers are subjected to a
double-well external potential. Under certain conditions, the phenomenology is characterized
by thermally activated barrier crossing (classical Kramers' problem). 
We addressed the issue of describing this complex high-dimensional 
dynamics in terms of few suitable observables, the reaction 
coordinates (RCs). 
In our system the center of mass, $Q$, and the end-to-end distance, $L$, are the 
most natural candidates.

These RCs evolve according to effective Langevin equations that are 
generally difficult to be derived via a systematic procedure.
The proper reconstruction of the effective stochastic equations for $Q$ and 
$L$ is achieved via a data-driven numerical method that extrapolates 
the drift and diffusion terms from a long trajectory of the original system.
Let us stress that this method allows us to find nonlinear terms for the
reconstructed Langevin equation; this is an important difference, e.g., from
the standard Mori-Zwanzig approach, where the complexity of the problem
is shifted into the shape of the memory kernels.~\cite{zwanzig61}
The reliability of the reconstructed dynamics is tested on the way it 
fairly reproduces some essential properties of the original dynamics, 
such as a correct estimate of the jump rate over the barrier.

From our study it emerges that the description level in terms of RCs strongly 
depends on the bond stiffness $K$. More precisely: if the bonds are 
rigid enough, we are allowed to consider the evolution of $Q$ only, given by 
model \eqref{eq:M1}, to fully characterize the jump dynamics.
However, when $K$ decreases, the internal motion of the 
chain cannot be neglected and also the dynamics of $L$ has to be considered, 
so that a satisfactory description can only be obtained in the plane $(Q,L)$ 
through the model \eqref{eq:M2}. On a more mathematical perspective, lowering $K$ can be regarded as a loss of Markovianity of the one-variable description. A second coordinate is necessary to recover the Markov property.
\iffalse
The emergence of such qualitative scenario can be also shown analytically,
under the simplifying assumption of HCA.
\fi

Our work shows how subtle is the procedure of reducing the dynamics of a
many-dimensional system to a low-dimensional model, even in the simplest cases 
where the physical intuition leaves little ambiguity to the choice of the RCs.
In fact, the choice of an RC can be correct in certain regimes but not 
sufficient in others. Specifically in our case, we can only know 
a posteriori that it is the stiffness of the polymer to determine 
whether one or two RCs are needed.

 
 
