Dynamic Looping of a Free-Draining Polymer

TL;DR

This paper revisits and generalizes the Wilemski-Fixman theory for polymer looping times, introducing a perturbation approach that captures known results and reveals a new term, supported by numerical simulations.

Contribution

A systematic perturbation method is developed to extend the WF theory, accommodating different sink types and revealing a previously unknown term in looping time expressions.

Findings

Reproduces known analytical and asymptotic results for polymer looping.

Identifies a new term in the analytical expression for looping time.

Numerical simulations confirm the theoretical predictions.

Abstract

We revisit the celebrated Wilemski-Fixman (WF) treatment for the looping time of a free-draining polymer. The WF theory introduces a sink term into the Fokker-Planck equation for the $3(N+1)$-dimensional Ornstein-Uhlenbeck process of the polymer dynamics, which accounts for the appropriate boundary condition due to the formation of a loop. The assumption for WF theory is considerably relaxed. A perturbation method approach is developed that justifies and generalizes the previous results using either a Delta sink or a Heaviside sink. For both types of sinks, we show that under the condition of a small dimensionless $\epsilon$, the ratio of capture radius to the Kuhn length, we are able to systematically produce all known analytical and asymptotic results obtained by other methods. This includes most notably the transition regime between the $N^2$ scaling of Doi, and $N\sqrt{N}/\epsilon$ scaling of Szabo, Schulten, and Schulten. The mathematical issue at play is the non-uniform convergence of $\epsilon\to 0$ and $N\to\infty$, the latter being an inherent part of the theory of a Gaussian polymer. Our analysis yields a novel term in the analytical expression for the looping time with small $\epsilon$, which is previously unknown. Monte Carlo numerical simulations corroborate the analytical findings. The systematic method developed here can be applied to other systems modeled by multi-dimensional Smoluchowski equations.

Tables

Table 1. Table 1: Comparison of theoretical results and simulations for selected values of N 𝑁 N and the capture radius ϵ italic-ϵ \epsilon . H. n. (Heaviside numerical) and DS. n. (Delta sink numerical) are obtained from the numerical integration in WF theory using Eq.   26 and Eq.   27 . H. a. (Heaviside analytical) and DS. a. (Delta sink analytical) are analytical results from Eq.   32 and Eq.   33 . SSS is the analytical result from Eq.   13 .

$N$, $ϵ$	Sample points	Simulation	H.n.	H.a.	DS. n.	DS. a.	SSS
50, 0.1	8000	563$\pm$12	677	632	564	530	426
50, 0.5	8000	184$\pm$4	233	200	205	186	85
50, 1.0	8000	113$\pm$2	163	143	141	134	43
75, 0.1	4000	1117$\pm$34	1314	1230	1107	1048	783
75, 0.5	4000	396$\pm$12	498	437	447	414	157
75, 1.0	4000	261$\pm$8	368	330	326	317	783
100, 0.1	2000	1790$\pm$76	2114	1984	1796	1710	1206
100, 0.5	2000	681$\pm$29	858	762	778	733	241
100, 1.0	2000	466$\pm$20	655	595	590	581	121

Equations

$$\frac{\mathrm{d}{\mathbf{R}}_i}{\mathrm{d}t}=-D_0{\nabla }_{{\mathbf{R}}_i}H(\mathbf{R})+{\mathbf{f}}_{i}i=0,1,\dots ,N,⟨{\mathbf{f}}_n(t){\mathbf{f}}_m(t^{\prime })⟩=2D_0\mathbf{I}{\delta }_{nm}{\delta }_{t{t}^{\prime }}$$

$${\theta }_p=\frac{p\pi }{2(N+1)},{\gamma }_p=12{\sin }^2({\theta }_p),$$

$$\displaystyle{\bf X}_{p}=\frac{1}{N+1}\sum_{n=0}^{N}{\bf R}_{n}\cos\Big{(}(2n+1)\theta_{p}\Big{)},\ \ {\bf R}_{n}={\bf X}_{0}+2\sum_{p=1}^{N}{\bf X}_{p}\cos\Big{(}(2n+1)\theta_{p}\Big{)},$$

$$\frac{\mathrm{d}{\mathbf{X}}_p}{\mathrm{d}t}=-{\gamma }_p{\mathbf{X}}_p+{\mathbf{F}}_p,⟨{\mathbf{F}}_p(t){\mathbf{F}}_q(t^{\prime })⟩=\frac{1}{N+1}\mathbf{I}{\delta }_{pq}{\delta }_{t{t}^{\prime }}\text{for}p+q>0$$

$$\displaystyle P_{e}({\bf X})=\prod_{p\ \text{odd}}\Big{(}\frac{(N+1)\gamma_{p}}{\pi}\Big{)}^{3/2}\exp(-(N+1)\gamma_{p}{\bf X}_{p}^{2})$$

$${\mathbf{R}}_{ee}={\mathbf{R}}_0-{\mathbf{R}}_N=\sum _{p\text{odd}}{c}_p{\mathbf{X}}_p,c_p=4\cos ({\theta }_p),$$

$$\displaystyle P_{{\bf r}}({\bf R}_{ee})=\Big{(}\frac{3}{2\pi N}\Big{)}^{3/2}\exp\Big{(}-\frac{3{\bf R}_{ee}^{2}}{2N}\Big{)}$$

$$\frac{\partial P(\mathbf{X},t)}{\partial t}$$

$$\varphi (t)=\frac{⟨{\mathbf{R}}_{ee}(t){\mathbf{R}}_{ee}(0)⟩}{⟨{\mathbf{R}}_{ee}^2⟩}=\frac{2}{N(N+1)}\sum _{p\text{odd}}\frac{1}{{\tan }^2({\theta }_p)}\exp (-{\gamma }_{p}t)$$

$$\displaystyle T({\bf r},t|{\bf r}_{0},0)=\Big{(}\frac{3}{2\pi N(1-\phi^{2}(t))}\Big{)}^{3/2}\exp\Big{(}-\frac{3}{2N}\frac{({\bf r}-\phi(t){\bf r}_{0})^{2}}{1-\phi^{2}(t)}\Big{)}$$

$$\frac{\mathrm{d}r}{\mathrm{d}t}=-D\frac{\mathrm{d}U(r)}{\mathrm{d}r}+\xi (t),⟨\xi (t)\xi (t^{\prime })⟩=2D\delta (t-t^{\prime }),$$

$$\tau (r)={\int }_{ϵ}^r\frac{1}{D\exp (-U(x))}\mathrm{d}x{\int }_x^{\infty }\exp (-U(y))\mathrm{d}y$$

$$\tau=\int_{\epsilon}^{\infty}\tau(r)P_{e}(r){\rm d}r=\int_{\epsilon}^{\infty}\frac{1}{DP_{e}(x)}{\rm d}x\Big{(}\int_{x}^{\infty}P_{e}(y){\rm d}y\Big{)}^{2}$$

$${\tau }_{\text{SSS}}\approx \frac{1}{8\pi ϵP_e(0)}=\frac{\sqrt{\pi }}{6\sqrt{6}}\frac{N\sqrt{N}}{ϵ}$$

$$\displaystyle L_{F}P({\bf X},t)=\nabla\cdot\Big{(}P_{e}({\bf X})\nabla\Big{(}\frac{P({\bf X},t)}{P_{e}({\bf X})}\Big{)}\Big{)}$$

$$\tau ={\int }_0^{\infty }\Sigma (t)\mathrm{d}t=\sum _{i=0}^{\infty }\frac{C_i}{{\lambda }_i^{ϵ}}\approx \frac{1}{{\lambda }_0^{ϵ}}.$$

$${\lambda }_0^{ϵ}=0+8\pi ϵJ+O({ϵ}^2)$$

$$\frac{\partial Q(\mathbf{X},t)}{\partial t}=L_{F}Q(\mathbf{X},t)-\kappa S(\mathbf{X})Q(\mathbf{X},t),Q(\mathbf{X},0)=P_e(\mathbf{X})$$

$$Q(\mathbf{X},t)=P_e(\mathbf{X})-\kappa {\int }_0^t\mathrm{d}{t}^{\prime }\int \mathrm{d}{\mathbf{X}}^{\prime }G(\mathbf{X},t|{\mathbf{X}}^{\prime },t^{\prime })S({\mathbf{X}}^{\prime })Q({\mathbf{X}}^{\prime },t^{\prime }).$$

$$\int Q(\mathbf{X},t)S(\mathbf{X})\mathrm{d}\mathbf{X}=$$

$$\times \int \mathrm{d}{\mathbf{X}}^{\prime }G(\mathbf{X},t|{\mathbf{X}}^{\prime },t^{\prime })S({\mathbf{X}}^{\prime })Q({\mathbf{X}}^{\prime },t^{\prime })$$

$$C(t)=⟨S(\mathbf{r},0),S(\mathbf{r},t)⟩=\int \mathrm{d}\mathbf{r}\int \mathrm{d}{\mathbf{r}}^{\prime }S(\mathbf{r})T(\mathbf{r},t|{\mathbf{r}}^{\prime },0)S({\mathbf{r}}^{\prime })P_{\mathbf{r}}({\mathbf{r}}^{\prime })$$

$$g(t,\kappa )P_0=P_0-\kappa {\int }_0^{t}g(t^{\prime },\kappa )C(t-t^{\prime })\mathrm{d}{t}^{\prime },P_0=\int P_{\mathbf{r}}(\mathbf{r})S(\mathbf{r})\mathrm{d}\mathbf{r}.$$

$$P_0={\int }_0^{t}q(t^{\prime })C(t-t^{\prime })\mathrm{d}{t}^{\prime }.$$

$$\tau$$

$$\displaystyle=\hat{I}(0)=\int_{0}^{\infty}{\rm d}t\Big{(}\frac{C(t)}{C(\infty)}-1\Big{)}=\lim_{s\rightarrow 0}\frac{\hat{C}(s)}{C(\infty)}-\frac{1}{s}.$$

$$\tau =\frac{1-C(0)}{C(0)}{\int }_0^{\infty }\frac{C(t)-C(\infty )}{C(0)-C(\infty )}\mathrm{d}t.$$

$$C(t)=$$

$$\displaystyle\times\sinh\Big{(}\frac{3\phi rr^{\prime}}{N(1-\phi^{2})}\Big{)}/\frac{3\phi rr^{\prime}}{N(1-\phi^{2})}$$

$$\displaystyle I_{H}(t)=\frac{C(t)}{C(\infty)}-1\approx\Big{(}\frac{1+\frac{4}{5}x_{0}}{1-\phi^{2}(t)+\frac{4}{5}x_{0}}\Big{)}^{3/2}-1$$

$$\displaystyle C(t)=\frac{12x_{0}/N\pi}{\phi\sqrt{1-\phi^{2}}}\exp\Big{(}\frac{-2x_{0}}{1-\phi^{2}}\Big{)}\sinh\Big{(}\frac{2x_{0}\phi}{1-\phi^{2}}\Big{)},C(\infty)=\frac{24x_{0}^{2}}{N\pi}\exp(-2x_{0})$$

$$\displaystyle I_{DS}(t)=\frac{C(t)}{C(\infty)}-1=\frac{\exp\Big{(}2x_{0}\phi/(1+\phi)\Big{)}-\exp\Big{(}-2x_{0}\phi/(1-\phi)\Big{)}}{4x_{0}\phi\sqrt{1-\phi^{2}}}-1.$$

$$\varphi (t)\approx \{\begin{array}{ll}1-\frac{6t}{N}& \text{Short Timescale}\\ 1-\frac{4}{N}\sqrt{\frac{3t}{\pi }}& \text{Median Timescale}\\ \frac{8}{{\pi }^2}\exp (-{\gamma }_{1}t)& \text{Long Timescale}\end{array}$$

$$D_{eff}(t)=\frac{⟨({\mathbf{R}}_{ee}(t)-{\mathbf{R}}_{ee}(0){)}^2⟩}{6t}=\frac{N(2-2\varphi (t))}{6t}$$

$$\bar{I}_{H}(t)=\begin{cases}\Big{(}\frac{1+\frac{4}{5}x_{0}}{\frac{12t}{N}+\frac{4}{5}x_{0}}\Big{)}^{3/2}&\text{Short Timescale}\\ \Big{(}\frac{1+\frac{4}{5}x_{0}}{\frac{8}{N}\sqrt{\frac{3t}{\pi}}+\frac{4}{5}x_{0}}\Big{)}^{3/2}&\text{Median Timescale}\\ \frac{96}{\pi^{4}(1+\frac{4}{5}x_{0})}\exp(-\frac{6\pi^{2}t}{(N+1)^{2}})&\text{Long Timescale}\end{cases}$$

Full text

{AMS}

82D60, 60G20, 60J70, 35P15, 47A55

1 Introduction

The theory of polymer dynamics is one of the most successful stochastic-process models in chemical science [13, 12]. Polymer dynamics in an aqueous solution is naturally stochastic. A polymer molecule dissolved in the so-called theta solvent, known also as free-draining, can be mathematically represented in terms of a multi-dimensional Ornstein-Uhlenbeck (OU) process. While the general theory of an OU process is well developed (see [27] and the references cited within), explicitly analytical results on the kinetics of the formation of a end-to-end loop are still highly sought after in theoretical chemistry and biochemistry. The problem is essentially a perturbation of a linear operator [18, 37, 7].

In the literature, “polymer looping” has been a classic problem with important biological applications. It is an essential step in various intra-cellular processes, including gene regulation, DNA replication, protein and RNA folding [29, 22, 25]. In some biological systems, two reacting molecules may bind distal target sites along a single polymer chain and a biological function requires the polymer to “spontaneously” form a loop structure that brings the two in contact [20]. This polymer-loop-mediated phenomena is due to thermal fluctuations and the loop formation consists of two steps: the “encounter” of the two ends of the polymer defined by a small capture distance, followed by a bond formation between the two ends when they are within the capture radius. When the time scale of the former is much longer than the latter, the loop formation is called diffusion controlled: The diffusive Brownian motion of the chain is then the rate determining step.

As an in-depth study of a multi-dimensional stochastic process, the problem was first discussed by Jacobson and Stockmayer in 1950 even before the dynamic Rouse model [16, 28]. They formulated the free energy cost of looping in terms of an equilibrium distribution. It is quantified as the ratio of respective equilibrium constants for intermolecular and intramolecular synapsis reactions when studying the polymer melts, which is called the J factor. The J factor is directly proportional to the statistical fraction of polymer conformations within the capture radius $a$, which is the ratio of the rate constants of the cyclization reaction, $k_{\text{coil}\rightarrow\text{loop}}/k_{\text{loop}\rightarrow\text{coil}}$. The J factor is mainly discussed in the more complicated case where the effect of bending and twisting rigidity is taken into account and so the polymer is no longer flexible [30]. The rate $k_{\text{coil}\rightarrow\text{loop}}$ is the inverse of the mean time for the two ends of a polymer to meet within the capture radius $a$ (or simply looping time $\tau$) from some configurations.

The study of the kinetics of loop formation is also a classical problem in polymer dynamics. It has led to intense theoretical and numerical research, even for the simplest case, of a Rouse chain. Though seemingly simple, the problem is actually challenging [23]. The mathematical formulation of the kinetics of this looping process involves a non-simple boundary-value problem of a linear PDE in a high dimensional space.

There were two major theoretical approaches in the early stage of the polymer looping investigation; they led to seemingly contradicting results. First, Wilemski and Fixman (WF) [39, 40] and Doi [11] estimated that the looping time for long polymer chains scales as $\sim N^{2}$, where $N$ is the number of monomers. Surprisingly, according to this result, the looping time does not depend on the capture radius $a$ (more details of this method will appear in later sections). On the other hand, Szabo, Schulten and Schulten (SSS) [33] estimated the looping time scales as $\sim N^{3/2}/a$. They approximated the dynamics of the end-to-end distance of the polymer as a single Brownian particle diffusing in a potential of mean force (this theory will be discussed in Section 3). Recent numerical simulations, however, indicated these time scales are showing up in different time regimes but the prediction from WF is more accurate than that from SSS theory [6, 23]. In fact, more recent theoretical advances show that the looping time may follow a mixed scaling law with $\tau\approx c_{1}N\sqrt{N}/a+c_{2}N^{2}$. Toan et al [36] proposed an effective space-dependent diffusion coefficient in SSS theory, such that both theories could be unified into one and the mixed scaling law was revealed. However, it has been difficult to verify whether this space dependency is really a mathematical approximation of a multi-dimensional OU process, or a different mathematical model all together. Very recently Amitai et al [1] also discovered a similar scaling law, based on the expansion of the eigenvalues of the Fokker-Planck equation in the limit of small capture radius. The first order perturbation in $a$ of the largest eigenvalue matches the SSS result and the next order perturbation roughly scales as $N^{2}$. But their theory works well only when $N$ is relatively small ($N<64$), and the contribution from other eigenvalues is unclear. Guèrin et al [34] and Bènichou et al [3] adapted the renewal equation method on the non-Markovian process of the end-to-end vector and reviewed other previous approaches. The numerical result they obtained is better than previous theories but an analytic scaling law is currently out of reach.

Within the OU process framework and in terms of analytical results, the mean first passage time is a solution to a boundary value problem associated with a high dimensional backward equation. Significant progress has been made in recent years on the asymptotic of passage-time problem associated with three-dimensional diffusion processes with a small exit boundary [7, 4, 8], crossing an unstable limit-cycle barrier [10], and rigorous lower bound on the density of the passage-time of the OU process [35]. The WF theory distinguishes itself from these other works in dealing with a high-dimensional OU process, which is equivalent to a low-dimensional non-Markovian process. With this in mind, the WF theory is currently the only analytically feasible approximation for the non-Markovian process of the end-to-end vector.

In the more than 40 years since the WF theory was first proposed, many efforts have been made to improve this theory numerically and theoretically. A comprehensive theoretical discussion is given in [38]. In particular, WF theory was found to agree well with simulation results for the case of small capture radius. However, no explicit analytical expression was provided and instead numerical integration was employed. Pastor et al [23] observed that this integral can be evaluated in the small time range in the limit of small capture radius $a$ and one finds the same limit as SSS. Moreover, they claimed the looping time must combine both kinds of behavior, $N^{3/2}/a$ and $N^{2}$. This shows the potential to extract the mixing scaling law from WF theory. Sokolov [31] and Likthman et al [21] further improved the WF theory to accommodate different initial conditions and used iterative methods to find exact solutions (for both cases numerical calculations were needed).

This paper is organized as follows. To make the presentation self contained, in Section 2 we review the Rouse model and derive the end-to-end vector as the summation of independent, non-identical OU processes. We express the looping time in terms of a boundary-value problem for a $3(N+1)$-dimensional Fokker-Planck equation. In Section 3, we review the Markovian approximation for the end-to-end vector and show that SSS theory corresponds to the first order perturbation of the largest eigenvalue for the Fokker-Planck equation. In Section 4, we discuss the WF theory rigorously and articulate the assumptions behind it. Specifically, we show that the assumption for the Delta sink is weaker than the Heaviside sink, so that the Delta sink gives better approximation results. The integrand appearing in WF theory in the expression for the looping time is not a genuine survival probability and we give a probabilistic interpretation for the Heaviside sink. The theory presented in this section is much more general and could be applied to other types of boundary-value problems. In Section 5, we extract the postulated mixed scaling laws directly from WF theory in the limit of small dimensionless $\epsilon$, the ratio of capture radius $a$ to the Kuhn length $b$, with the help of perturbation methods. The mixed scaling law we find includes an extra term $N\sqrt{N}$ which will be dominated by $N^{2}$ term as $N\rightarrow+\infty$. In Section 6, we use Monte Carlo simulations to estimate the looping time for various $N$ and the capture radius $\epsilon$. In general, WF theory overestimates the looping time but agrees very well with the simulation in the small $\epsilon$ regime (irrespective of the value of $N$). We also verify the mixed scaling law numerically through regression. Last but not least, the survival probability function is numerically obtained by deconvolution and for small $\epsilon$, it is a single exponential function and agrees remarkably well with numerical simulation.

2 Rouse model

The Rouse model is the simplest flexible polymer model which uses beads and springs to represent the polymer chain in a viscous fluid. It assumes beads have no excluded volume and no hydrodynamic interactions among them. Actually, the notion of a theta solvent is defined as an exact cancellation between these two opposing effects. Two adjacent beads are connected with a harmonic spring with the same spring constant $k$ [12, 19] and in total there are $N+1$ beads and $N$ springs. The potential energy of the entire mechanical system, $H(\vec{{\bf R}})$, is given by $H(\vec{{\bf R}})=\frac{1}{2}k\sum_{n=1}^{N}({\bf R}_{n}-{\bf R}_{n-1})^{2}$, where ${\bf R}_{n}$ is the position of $n$th bead in three-dimensional space. Let $b$ denote the Kuhn segment length or effective bond length, that is defined as the mean square length of the bond. Then, if $k=3/b^{2}$ in $k_{B}T$ units, the Boltzmann distribution of the conformation and the end-to-end vector distribution will be the same as the ones for the Gaussian chain. In addition to the Hookean force from its connecting springs, a bead experiences a frictional force when it moves. Each bead is assumed to have the same friction coefficient $\zeta$. Inertia is negligible since the motion is overdamped. The diffusion coefficient of a monomer is $D_{0}=1/\zeta$ in $k_{B}T$ units as dictated by the Einstein relation. The dynamics are described by Langevin equations Eq. 1 with Gaussian white noise [17],

[TABLE]

Through a normal mode transformation (2), the Langevin equations Eq. 3 for the modes ${\bf X}_{p}$ are decoupled and each mode evolves independently but not identically to the others.

[TABLE]

Here $\gamma_{p}$ are eigenvalues of the Laplacian matrix of the Rouse chain. The zeroth mode of the chain represents the dynamics of the center of mass for the chain which is Brownian motion. The dynamics is also called a reversible Ornstein-Uhlenbeck process [27] and its equilibrium distribution is the Boltzmann distribution $P_{e}({\bf X})$ in Eq. 4 [26].

[TABLE]

Nondimensionalize the system with the Kuhn length $b$ as the characteristic length and $b^{2}/D_{0}$ as the characteristic time, one can define a dimensionless parameter as the ratio of capture radius to the Kuhn length, $\epsilon=a/b$. It specifies the capture radius in the unit of the Kuhn length so we will call $\epsilon$ as capture radius as well. From now on, all quantities and equations are nondimensionalized unless otherwise specified.

The end-to-end vector ${\bf R}_{ee}$ can be expressed as the linear combination of the odd order terms of ${\bf X}_{p}$ with coefficient $c_{p}$ (see Equation 5). All even order modes do not contribute if the chain has a homogenous spring constant. The equilibrium distribution of the end-to-end vector is also a Gaussian with variance $N$. Then the J factor has the simple expression as scaling of $N$, $J\propto(3/2\pi N)^{3/2}$.

[TABLE]

The relaxation time $\tau_{p}$ for the odd modes ${\bf X}_{p}$ is $\tau_{p}=1/\gamma_{p}$. Specifically, $\tau_{1}$ is the largest, $\tau_{1}\approx\frac{(N+1)^{2}}{3\pi^{2}}$. So the relaxation time for the end-to-end vector ${\bf R}_{ee}$ is dominated by $\tau_{1}$.

The corresponding Fokker-Planck equation for the Langevin equation Eq. 3 describes the time evolution of the probability density $P({\bf X},t)$ [14]. This linear operator is denoted as $L_{F}.$ The Green function of the Fokker-Planck equation is the probability propagator, $G({\bf X},t|{\bf X}^{0},0),$ and is a Gaussian function.

[TABLE]

in which $\Gamma$ is a diagonal matrix with elements $\gamma_{p}$. But the dynamics of the end-to-end vector ${\bf R}_{ee}$, being a projection of a $N$-dimensional OU process, is non-Markovian. One could sample from Boltzmann (canonical) distribution polymer configurations with fixed ${\bf R}_{ee}={\bf r}_{0}$ and study the average end-to-end vector dynamics. These dynamics follow a non-Markovian Gaussian process with the conditional probability $T({\bf r},t|{\bf r}_{0},0)$. The probability at time $t$ on location ${\bf r}$ is also w.r.t. the canonical ensemble. It is determined by the time correlation function of the end-to-end vector, $\phi(t)$, which is the summation of the correlation functions of the odd modes. The conditional probability $T({\bf r},t|{\bf r}_{0},0)$ looks like a probability propagator, but it doesn’t satisfy the Kolmogorov-Chapman equation, i.e, $\int T({\bf r},t|{\bf r}^{\prime},t^{\prime})T({\bf r}^{\prime},t^{\prime}|{\bf r}_{0},0){\rm d}{\bf r}^{\prime}\neq T({\bf r},t|{\bf r}_{0},0)$, for $t>t^{\prime}>0$.

[TABLE]

Once the two ends of the chain get close (within the capture distance $\epsilon$), the reaction may start. So this reactive surface is described as a tubular neighborhood $S_{\epsilon}=\{{\bf X}:\|{\bf R}_{ee}\|=\epsilon\}$. In the ideal case, every time the two ends get closer than $\epsilon$, the reaction happens and this corresponds to the absorbing boundary condition $P({\bf X},t)=0$ on the surface $S_{\epsilon}$. Once the reaction term is imposed, the modes are no longer independent.

In the literature, when addressing the looping problem, it is assumed that the initial distribution is the Boltzmann distribution $P_{e}({\bf X})$. But the initial condition doesn’t fulfill the boundary condition. Physically all configurations that are in a looped state react immediately, so the initial condition becomes the normalized Boltzmann distribution outside the reactive surface and zero inside the surface. The difference between these two initial conditions are negligible since the probability inside the surface is very small compared with outside. It is worth noting that the polymer chain sampled from the Boltzmann distribution doesn’t imply all internal modes of the chain are fully relaxed.

If the PDE Eq. 7 with the boundary condition specified could be solved, one could integrate the solution $P({\bf X},t)$ over the configuration space ${\bf X}$ and get the probability that the chain hasn’t formed a loop at time $t$, the survival probability $\Sigma(t)$. Then $1-\Sigma(t)$ is the cumulative probability distribution for the random time of capturing, and the mean looping time $\tau$ is the time integral of survival probability through the time domain with some mild assumptions, which can be fully determined. $\tau$ is a function of the number of beads $N$ and the dimensionless capture radius $\epsilon$. One intuitively expects that $\tau\to\infty$ as $N\to\infty$, or $\epsilon\to 0$.

3 Markovian approximation

A simple approach is to consider that the two ends of the polymer behave as if they were diffusing apart from the rest of the chain, such that no memory effect is expected. We can project all these internal modes of the polymer into a 1D single variable, the end-to-end distance $r=\|{\bf R}_{ee}\|$ (after angular components are averaged). This is the theory first developed by Szabo, Schulten and Schulten [33]. Suppose that the end-to-end distance obeys a reduced Langevin equation of the form with diffusion coefficient D in the unit of $D_{0}$,

[TABLE]

where the effective potential $U(r)$ of the mean field is given by $U(r)=-\ln P_{e}(r)$, where P_{e}(r)=4\pi r^{2}\Big{(}\frac{3}{2\pi N}\Big{)}^{3/2}\exp\Big{(}-\frac{3r^{2}}{2N}\Big{)}. One can calculate the mean first passage time $\tau(r)$ starting from the fixed end-to-end distance $r$ and integrate this space dependent function w.r.t. the equilibrium distribution to get the looping time $\tau$.

The dynamics can be considered those of a single particle diffusing in the effective potential. The looping problem becomes a standard first passage time problem which can be addressed by solving the adjoint operator equation, D\Big{(}\frac{{\rm d}^{2}\tau(r)}{{\rm d}r^{2}}-\frac{{\rm d}U}{{\rm d}r}\frac{{\rm d}}{{\rm d}r}\tau(r)\Big{)}=-1, with boundary condition $\tau(\epsilon)=0$. One can assume a reflective boundary condition when the chain is fully stretched, i.e. $\frac{{\rm d}}{{\rm d}r}\tau(N)=0.$ In fact, the difference with the natural boundary condition here is negligible since the probability of the chain reaching that far is very small for the flexible chain. It is easy to show the solution of the adjoint equation is

[TABLE]

Averaging over the equilibrium distribution of the end-to-end distance and find,

[TABLE]

Equation 12 is known as the SSS theory of polymer looping. It can be simplified through the Laplace method if we assume that the potential barrier at $r=\epsilon$ is very high. In SSS theory, the diffusion coefficient $D$ is chosen as the relative diffusion coefficient of the two end monomers, $D=2.$ (in Section 5, it is shown that this choice of diffusion coefficient is related to a short time approximation for the end-to-end correlation function $\phi(t)$). For the Rouse model, the looping time becomes

[TABLE]

Specifically, $P_{e}(0)$ is exactly proportional to the J factor of the Gaussian chain. Another way to derive the SSS result is to coarse grain the polymer chain to a dumbbell model with the effective spring constant $k_{\text{eff}}=3/N.$

It turns out that the equilibrium assumption about the chain is not an accurate one. Not all internal modes of the polymer chain have relaxed before the looping process starts. So when $\tau_{1}\gg\tau_{\text{SSS}}$, the effective potential is no longer time independent. The motion of all internal modes has be taken into consideration in the problem of looping. If $\epsilon\ll 1/\sqrt{N}$, the looping time $\tau$ is longer than the relaxation time of the chain, so all internal modes are fully relaxed before the two ends meet and SSS theory works well.

Amitai et al [1] calculate the perturbation of the spectrum for the Fokker-Planck equation Eq. 7 and estimate the looping time accordingly. If we assume the natural boundary condition, the system is detailed balanced with the Boltzmann distribution as unique stationary density [27]. Also, in this case, the operator $L_{F}$ in the linear PDE Eq. 7 is hermitian in the Hilbert space $L^{2}({\bf X})$ with the inner product defined as $\langle f,g\rangle=\int_{\Omega}\frac{f({\bf X})g({\bf X})}{P_{e}({\bf X})}{\rm d}{\bf X}$ [17]. The linear operator can be rewritten as

[TABLE]

Therefore it has a discrete real spectrum $\{-\lambda_{i}^{0}\}$ with complete orthogonal eigenfunction basis $\{\phi_{i}^{0}\}$. The spectrum consists of nonnegative integer-weighted sums of eigenvalues of the Laplace matrix, $-\gamma_{p},$ and the eigenfunctions are products of the corresponding Hermite polynomials. The top eigenvalue $\lambda_{0}^{0}=0$ and the corresponding eigenfunction is the stationary density function. The next two are $\lambda_{1}^{0}=\gamma_{1}$ and $\lambda_{2}^{0}=\gamma_{2}$.

But with the absorbing boundary condition, the spectrum $\{-\lambda_{i}^{\epsilon}\}$ is shifted due to removing the small tubular neighborhood $S_{\epsilon}$ and there is no longer a stationary density. If $\epsilon$ is small enough, the spectrum will be still distinct discrete and real. The time dependent solution can still be expressed as an expansion in the eigenfunctions $\{\phi_{i}^{\epsilon}\}$ with coefficients $c_{i}$, $P({\bf X},t)=\sum_{i=0}^{\infty}c_{i}\phi_{i}^{\epsilon}({\bf X})\exp(-\lambda_{i}^{\epsilon}t)$. If the initial condition is the previous stationary density, i.e. the Boltzmann distribution, the eigenfunction $\phi_{0}^{\epsilon}$ is close to the Boltzmann distribution so that most of the energy is contributed by this zeroth mode, i.e, $c_{0}\approx 1$. The survival probability $\Sigma(t)=\int_{\Omega-S_{\epsilon}}P({\bf X},t){\rm d}{\bf X}$, is in fact the summation of infinitely many exponential distributions. If the survival probability decays to 0 sufficiently fast, i.e, $\lim_{t\rightarrow+\infty}t\Sigma(t)=0,$ then the mean first passage time $\tau$ is simply the integral of the survival probability from 0 to $+\infty,$

[TABLE]

Amitai et al [1] discovered that the first passage time roughly follows one single exponential distribution for relatively small $N$ ($N<64$) by numerical simulation, so the zeroth mode in Equation 15 is enough. We will discuss the distribution of the first passage time in Section 6. When $\epsilon$ is small enough, the first order perturbation of the zeroth eigenvalue is proportional to the ratio of the partition function of the closed polymer chain to the whole configuration space. This ratio is again exactly the J factor and the proportionality spatial factor is $8\pi\epsilon$.

[TABLE]

The series expansion of the eigenvalues follows results in [7]. The J factor in the scaling of $N$ appears in the expression was discussed in [15]. If one uses Eq. 16 in Eq. 15, then the first order perturbation result exactly matches the SSS Markovian approximation, i.e, the looping time is exactly Eq. 13.

In summary, the SSS theory describes the kinetics of loop formation as a diffusion process in an effective potential of mean force that is derived from the equilibrium distribution for the end-to-end distance $P_{e}(r)$. It approximates the non-Markovian dynamics of the end-to-end distance $r$ by these simple Markovian dynamics. It assumes that the internal modes have relaxed before the looping process starts (this is called the local equilibrium assumption). The condition for this to hold is $\epsilon\ll 1/\sqrt{N}$. As far as the spectrum of the linear operator is concerned, this Markovian approximation corresponds to the first order perturbation of the zeroth mode. In Section 6, it is shown that the SSS result significantly underestimates the looping time (it can also be proved through a variational principle [24]). Roughly speaking, the Markovian estimate ignores the case that two ends meet each other due to fluctuation, but the center of mass of the chain may be far away from the two ends. Intuitively, this case is more likely to happen for relatively large $\epsilon$ or large $N$. This drawback brings us to a more comprehensive method, Wilemski-Fixman theory.

4 WF theory

One could incorporate a distance-dependent reaction term in the Fokker Planck equation Eq. 7 with a microscopic rate constant $\kappa$ and relate the looping time with the time integral of a normalized sink-sink time correlation function. This is the celebrated Wilemski-Fixman (WF) approximation [39, 40]. Instead of solving a Fokker-Planck equation Eq. 7 with a complicated boundary condition, the modified equation is of convection-diffusion-reaction type in free space (see Eq. 17) and one can express the solution in terms of the Green function $G({\bf X},t|{\bf X}^{\prime},0)$ for Equation 7. This method can also be applied to study the Fokker-Planck equation with other difficult boundary conditions, like diffusion-limited catalytically-activated chemical reactions with a special catalytic subvolume [2], facilitated diffusion and looping with a heterogenous Rouse chain. The microscopic rate constant $\kappa$, which measures the effectiveness of the reaction, is not the same as the coefficient in the partially absorbing boundary condition in Collins and Kimball’s kinetic theory. In fact, Szabo et al [32] shows that the partially absorbing boundary condition is equivalent to the use of a delta sink on the reactive surface in conjunction with a reflecting boundary condition. However, as we will show later, if we let $\kappa\rightarrow+\infty$ in the WF approximation, the absorbing boundary condition can be recovered under some assumptions.

[TABLE]

Two popular choices for the sink function are the Heaviside function and the radial delta function. For a given configuration ${\bf X}$, the Heaviside function is $S_{1}({\bf X})=\begin{cases}1,&\text{if}\ \|{\bf R}_{ee}\|\leq\epsilon\\ 0,&\text{Otherwise}\end{cases}$. Intuitively, it corresponds to the Brownian particle starting to react once inside the reactive surface $S_{\epsilon}$. The delta function is $S_{2}({\bf X})=\frac{1}{4\pi\epsilon^{2}}\delta(\|{\bf R}_{ee}\|-\epsilon)$. It is also called the Smoluchowski reaction sink. The reactive spherical surface is in three-dimensional space and it separates the whole space into two regions, inside and outside of the reactive surface. But this sink doesn’t allow the Brownian particle to pass the surface. Both types of sinks have a three-dimensional representation in terms of the end-to-end vector ${\bf r}$.

The solution of Eq. 17 starting with the equilibrium distribution $P_{e}({\bf X})$ can be expressed with the help of Dyson’s formula as

[TABLE]

The inner space integral gives the probability density at time $t$ in the original Fokker-Planck equation Eq. 7 if the initial distribution at time $t^{\prime}$ is $S({\bf X})Q({\bf X},t^{\prime})$. The meaning of Eq. 18 is as follows: the probability of a chain to have the configuration ${\bf X}$ at time $t$, is the probability to observe the configuration ${\bf X}$ without sink, minus the probability of reaching the configuration ${\bf X}$ at time $t$ starting in the sink at some point $t^{\prime}$ between 0 to $t.$ The advantage of this approach is that the Green function for Eq. 7 is known. The survival probability $\Sigma(t)$ is given by $\Sigma(t)=1-\kappa\int_{0}^{t}{\rm d}t^{\prime}\int{\rm d}{\bf X}^{\prime}S({\bf X}^{\prime})Q({\bf X}^{\prime},t^{\prime})$.

To solve analytically for $Q({\bf X},t)$ is still very difficult. The WF theory takes advantage of the conditional probability $T$ to reduce the dimensionality from $\mathbb{R}^{3N+3}$ to $\mathbb{R}^{3}$. One can multiply $S({\bf X})$ on both sides of Eq. 18 and integrate w.r.t. ${\bf X}$ to get Eq. 19. If one multiplies with a sink form other than $S({\bf X}),$ the result will be an unbalanced sink-sink correlation.

[TABLE]

If $S({\bf X})Q({\bf X},t)$ is the canonical ensemble with the Boltzmann distribution for the end-to-end vector ${\bf r},$ one could rewrite Eq. 19 in terms of $T({\bf r},t|{\bf r}_{0},0).$ Specifically, for the Heaviside sink, the space and time dependencies of $Q({\bf X},t)$ are separated inside the sink, $Q({\bf X},t)=P_{e}({\bf X})g(t,\kappa)$ when $\|{\bf r}\|\leq\epsilon$; for the delta sink, $Q({\bf X},t)$ only requires the same space and time separation on the reactive surface, i.e, $Q({\bf X},t)=P_{e}({\bf X})g(t,\kappa)$ when $\|{\bf r}\|=\epsilon$. The condition for the delta sink is weaker than for the Heaviside sink, so one would expect that the delta sink performs better than the Heaviside sink. Note that both are weaker than the original WF assumption which requires time and space separation in the whole space. In either case, $Q({\bf X},t)$ is not homogenous on the reactive surface. In fact, WF theory allows non-Markovian dynamics for the end-to-end vector of the chain but doesn’t capture the full non-Markovian effect.

With the help of the sink-sink correlation function $C(t)$, Equation 19 can be rewritten into a more compact form Eq. 20 and Eq. 21, if the separation condition is satisfied.

[TABLE]

One observes that $C(\infty)=(P_{0})^{2}$ for both sinks. For the delta sink, $C(t)$ has a singularity at $t=0$; for the Heaviside sink, $C(0)=P_{0}$, the probability of the looped state at equilibrium. In general, there is no closed form expression for $C(t)$.

Assume that on the reactive surface, the time function $g(t,\kappa)$ is asymptotically $q(t)/\kappa+O(1/\kappa^{2})$. If we let the microscopic rate constant $\kappa\rightarrow+\infty$, then $Q({\bf X},t)\rightarrow 0$ on the reactive surface asymptotically and the absorbing boundary condition is recovered. Under this assumption, the survival probability will be \Sigma(t)=1-\Big{(}\int_{0}^{t}P_{0}q(t^{\prime}){\rm d}t^{\prime}\Big{)}. In this limit, the left hand side of Eq. 21 vanishes and we find that

[TABLE]

Equation 22 is valid inside the sink for the Heaviside sink and on the reactive surface for the delta sink respectively, and looping time only requires the solution in these regions. This defines a deconvolution problem with the kernel given by $C(t).$ This is an inverse problem and one way to solve it is to use the Laplace transform. The Laplace transform of $q(t)$ and $\Sigma(t)$ are, $\hat{q}(s)=\frac{P_{0}}{s\hat{C}(s)}$ and $\hat{\Sigma}(s)=\frac{1}{s}-\frac{(P_{0})^{2}}{s^{2}\hat{C}(s)}$ resepctively. The looping time is given by $\tau=\lim_{s\rightarrow 0}\hat{\Sigma}(s).$ Since $C(t)$ is a decreasing function to $C(\infty)$, in the transform variable we have $\hat{C}(s)>C(\infty)/s.$ If one approximates $\hat{\Sigma}(s)$ in the denominator by $sC(\infty),$ then $\tau\leq\lim_{s\rightarrow 0}\frac{\hat{C}(s)}{C(\infty)}-\frac{1}{s}.$ In fact, we can show the inequality is an equality as follows.

Define $I(t)=\frac{C(t)}{C(\infty)}-1.$ The improper integral of $I(t)$ from 0 to $+\infty$ is finite, and it is equivalent with $\hat{I}(0)<\infty$ in Laplace transform. The looping time $\tau$ is given by

[TABLE]

Finally, the integral \int_{0}^{\infty}{\rm d}t\Big{(}\frac{C(t)}{C(\infty)}-1\Big{)} is the famous WF approximation formula which one can evaluate numerically. But the integrand $\frac{C(t)}{C(\infty)}-1$ is clearly not the survival probability $\Sigma(t)$ because this small alteration in the Laplace domain would change the integrand function significantly in the time domain.

We can provide a probabilistic interpretation of the WF formula for the Heaviside sink. $C(0)=P_{0}$ is the fraction of the polymer that has already formed the loop initially. If one only takes this fraction of the polymers to start the process, the polymer will start to un-loop gradually. We assume that once they are un-looped, these polymers will not form the loop again under the time scale of interest. After sufficiently long time, the distribution of end-to-end vector will be the Gaussian distribution with normalization factor $P_{0}$ and the fraction of polymers that is still in the loop state will be $C(\infty)=P_{0}^{2}.$ So C(t) describes the fraction of the polymer that is still in the loop state at time $t$ if the process starts with the looped polymer. The process of the polymers un-looping themselves is intimately related to the looping problem. We can rewrite the integrand as follows

[TABLE]

The expression $\frac{C(t)-C(\infty)}{C(0)-C(\infty)}$ is taken as the approximation of the survival probability for the un-looping process. The integral calculates the expected un-looping time. If one considers this looping process as the dynamical process in a bi-stable system, the ratio $\frac{1-C(0)}{C(0)}$ is the relative stability of the two potential wells, which builds the connection from un-looping time to looping time. A disadvantage with this argument is that it does not work for the delta sink case.

In summary, the approximations made are:

(i) The solution $Q({\bf X},t)$ of Eq. 17 has the space and time separation inside the sink for the Heaviside sink and on the reactive surface for the delta sink. This approximation is the key step for the construction of this conditional probability w.r.t. canonical ensemble and the reduction in dimension. It is a plausible assumption because the canonical ensemble of polymers will behave as quasi-stationary, at least for the small capture radius case. In order to go beyond the WF approximation, one has to relax this condition.

(ii) The time function $g(t,\kappa)$ is analytic with respect to the variable $1/\kappa$, such that the absorbing boundary is recovered. Although the assumption is difficult to verify, it is reasonable physically. As the microscopic rate constant $\kappa$ increases, the reaction is more and more likely to happen once the two ends are within the capture radius. In the limit, the reaction will happen immediately which corresponds to the absorbing condition.

(iii) Instead of solving the deconvolution problem numerically, the WF theory approximates the denominator in the Laplace domain, and it gives a semi-analytical form for the looping time $\tau$ directly, but the analytical form for the survival probability is unknown. We will address this issue in Section 6 and solve numerically for the survival probability $\Sigma(t)$ there.

5 Perturbation theory

Although the WF theory provides a good estimate of the looping kinetics, the numerical integration neither provides a reduced model nor gives the scaling law in the two parameters $N$ and $\epsilon.$ We are going to apply perturbation techniques on the time integral to extract asymptotic estimates of the looping time.

The sink-sink correlation function $C(t)$ can be expressed as a double integral over two radial variables r and $r^{\prime},$ after averaging out the azimuthal and polar angles. It uses the fact that $\int{\rm d}{\bf r}\int{\rm d}{\bf r}_{0}\exp({\bf r}{\bf r}^{\prime})=\int{\rm d}r4\pi r^{2}\int{\rm d}r^{\prime}4\pi r^{\prime 2}\frac{\sinh(rr^{\prime})}{rr^{\prime}}$.

[TABLE]

We introduce the small dimensionless quantity $x_{0}=\frac{3\epsilon^{2}}{2N}\ll 1.$ $x_{0}$ is a small quantity since the nature of the looping assumes the capture radius $\epsilon$ is much smaller than the average end-to-end distance $\sqrt{N}$. Using $x_{0}$ we can find explicit expressions for the integrand $I(t).$

For the Heaviside sink, the double integral is evaluated by expanding in powers of $x_{0}$, C(t)=\frac{16x_{0}^{3}}{9\pi(1-\phi^{2})^{3/2}}\Big{(}1-\frac{6x_{0}}{5(1-\phi^{2})}+\dots\Big{)}. To match the similar form in [11], Pastor et al [23] proposed the closed form C(t)\approx\frac{16x_{0}^{3}}{9\pi}\Big{(}1-\phi^{2}+\frac{4}{5}x_{0}\Big{)}^{-3/2}, which matches the first two order of the expansion in $x_{0}$. Then

[TABLE]

The time integral of $I_{H}(t)$ is roughly approximated in two different time scales. First, in the short time scale, $\phi\approx 1$, the denominator is approximated by $2(1-\phi)+\frac{4}{5}x_{0}.$ The integrand in this time scale is much larger than 1 since $x_{0}$ is sufficiently small. The approximation of $I_{H}(t)$ is, I_{H}(t)\approx\Big{(}\frac{1+\frac{4}{5}x_{0}}{2(1-\phi(t))+\frac{4}{5}x_{0}}\Big{)}^{3/2}. Second, in the long time scale, $\phi\approx 0$, one could estimate the asymptote as the integrand goes to 0 by exploiting that the quantity $\frac{\phi^{2}(t)}{1-\phi^{2}(t)+\frac{4}{5}x_{0}}$ is small. So the approximation of $I_{H}(t)$ under this time scale is, $I_{H}(t)\approx\frac{3}{2}\frac{\phi^{2}(t)}{1+\frac{4}{5}x_{0}}$.

For the delta sink, the double integral is evaluated exactly

[TABLE]

Similarly to the Heaviside sink case, it can be approximated in two different time scales. First, in the short time scale, $\phi\approx 1$, so the approximation of $I_{DS}(t)$ is given by $I_{DS}(t)\approx\frac{\exp(x_{0})-\exp(-2x_{0}/(1-\phi))}{4\sqrt{2}x_{0}\sqrt{1-\phi}}.$ Second, in the long time scale, $\phi\approx 0,$ one can use a Taylor expansion to estimate the limit of the integrand. So the approximation of $I_{DS}(t)$ is $I_{DS}(t)\approx\frac{3}{2}\phi^{2}$.

Since the time correlation function for the end-to-end vector $\phi(t)$ doesn’t have a closed analytical form, the integral is still not feasible analytically. We will approximate $\phi(t)$ in three different time scales: First, when the time scale is within the relaxation time for the largest mode, $t\leq\tau_{N}=\frac{1}{12}$, take the approximation $\exp(-\gamma_{p}t)=1-\gamma_{p}t+O(t^{2})$. Then $\phi(t)$ can be approximated by $\phi(t)\approx 1-\frac{6t}{N}$. Second, when time $t$ between two time scales, $\tau_{N}\ll t\ll\tau_{3}=\frac{N^{2}}{27\pi^{2}}$, the approximation will be based on $\frac{1}{\tan^{2}(\theta_{p})}\approx\frac{1}{\theta_{p}^{2}}$ and $\gamma_{p}\approx p^{2}\gamma_{1}$. Both approximations work well on small $p$ that also contribute the most to the correlation function, \phi(t)\approx 1-\frac{8}{\pi^{2}}\sum_{p\ \text{odd}}\frac{1}{p^{2}}\Big{(}1-\exp(-p^{2}\gamma_{1}t)\Big{)}. One has to further approximate by turning the summation into the integral. Define $x=p\sqrt{\gamma_{1}t}$, the coefficient $N^{-1}\ll\sqrt{\gamma_{1}t}\ll 1$ is a small quantity, $\phi(t)\approx 1-\frac{4}{\pi^{2}}\sqrt{\gamma_{1}t}\int_{0}^{+\infty}\frac{1}{x^{2}}(1-\exp(-x^{2})){\rm d}x=1-\frac{4}{N}\sqrt{\frac{3t}{\pi}}$. Last, when the time $t\gg\tau_{3},$ all other modes are relaxed except the first mode. Then it can be approximated by $\phi(t)\approx\frac{8}{\pi^{2}}\exp(-\gamma_{1}t)$. In summary, the end-to-end vector correlation function $\phi(t)$ has the following analytical approximation

[TABLE]

The effective diffusion coefficient $D_{eff}$ for the end-to-end vector is defined by the end-to-end correlation function $\phi(t)$,

[TABLE]

In the short timescale, $D_{eff}=2$ which is time homogeneous. It also verifies the choice of diffusion coefficient in SSS theory. SSS theory doesn’t capture the behavior from other timescales, so it cannot reproduce the mixed scaling law.

Numerical simulations show that all three asymptotic results perform very well under the appropriate timescale. Specifically, for $N=75,$ the short timescale approximation works better than the median for $t\leq t_{1}=\frac{4}{3\pi}\approx 5.1\tau_{N}$ and the long timescale approximation is better than the median for $t\geq 2.8\tau_{3}$. Notice when $t\approx\tau_{3},$ the correlation function $\phi$ is about 0.76, so one can still treat $\phi\approx 1$ under the median timescale. Then $I_{H}$ and $I_{DS}$ have analytical approximations $\bar{I}_{H}$ and $\bar{I}_{DS}$ as follows,

[TABLE]

[TABLE]

Specifically, for $N=75$ and $\epsilon=0.75$, the analytical approximations perform very well for both sinks. The time range within which they deviate most is on the boundary of the short timescale and the median timescale. In terms of numerical integration, most of deviation is contributed by the median timescale and its boundaries as expected (figure not shown). The advantage of the approximation is that the integral which provides an estimate of the looping time can be evaluated analytically. Set the two timescale break points as $t_{1}=\frac{4}{3\pi}$ and $t_{2}=3\tau_{3}$. The second breakpoint between the median and long timescales depends on the number of beads, $N.$ We roughly choose $3\tau_{3}$ as the second breakpoint but this will not change the behavior of scaling law qualitatively.

For the Heaviside sink, one can compute the mixed scaling law explicitly,

[TABLE]

When the capture radius $\epsilon\ll 1,$ this integral is $\int_{0}^{\infty}\bar{I}_{H}(t){\rm d}t\approx h_{1}\frac{N\sqrt{N}}{\epsilon}+h_{2}N\sqrt{N}+h_{3}N^{2}$ where $h_{1},h_{2}$ and $h_{3}$ are constants. In the limit of $N\rightarrow+\infty$, we can ignore the middle term $h_{2}N\sqrt{N}$ and recover the mixed scaling law hypothesized before.

For delta sink, the integral is much more complicated,

[TABLE]

[TABLE]

Similarly when $\epsilon\ll 1$, one can use Taylor expansions for $\exp(x)$ and $\operatorname{erf}(x)$ functions at $x=0$. The integral is roughly $\int_{0}^{\infty}\bar{I}_{DS}(t){\rm d}t\approx d_{1}\frac{N\sqrt{N}}{\epsilon}+d_{2}N\sqrt{N}+d_{3}N^{2}$ as well, where $d_{1},d_{2}$ and $d_{3}$ are constants. In the limit $N\rightarrow+\infty$, it recovers the mixed scaling law again. Specifically, $d_{1}$ is mostly contributed by the short timescale approximation in $I_{1}$ and is about $\frac{\sqrt{\pi}}{6\sqrt{6}}$. It exactly matches the SSS result as predicted in [23]. It implies that the SSS and WF theories give the same asymptotic result $\frac{N\sqrt{N}}{\epsilon}$ in the limit $\epsilon\rightarrow 0$ for given $N$. In practice, this asymptotic result is realized when $\epsilon$ is extremely small. $d_{2}$ and $d_{3}$ are mostly contributed by the median and long timescale approximations. This is also predicted by Doi in [11] and Doi provided a dynamical explanation as well. So both $N\sqrt{N}$ and $N^{2}$ term are considered as the next order approximation result when $\epsilon$ is still relatively small. One would insightfully rewrite the looping time with the mixed scaling law as \tau=\frac{N\sqrt{N}}{\epsilon}\Big{(}d_{1}+(d_{2}+d_{3}\sqrt{N})\epsilon+O(\epsilon^{2})\Big{)}. Note in the derivation of this scaling law, we assume $N$ as a constant is much larger than 1, i.e, $N\gg 1$. From the scaling law, two seemingly contradicting results from Doi and SSS are the consequence of the non-uniforming convergence of $\epsilon\to 0$ and $N\to+\infty$. Such results are typical in singular perturbation problems.

6 Numerical Simulation

The Monte Carlo simulation algorithm is based on [23] but with more sample points and smaller time step. The simulation results are also within the range of their paper. The numerical solver is the Euler-Maruyama method, and the initial condition is sampled from the equilibrium distribution. The trajectory is terminated once the end-to-end distance is within the capture radius and record the passage time. The time step is chosen adaptively: when the end-to-end distance is within 2 Kuhn lengths, much finer time step is picked to prevent overshooting; but a larger time step is permitted when the end-to-end distance is outside of the range. The time step defined is

[TABLE]

with $\Delta_{l}=10^{-6}$ and $\Delta_{h}=10^{-3}$. This choice of time step is finer than that in [23]. For each parameter set, there are at least $n=2000$ samples. The passage time is considered to follow a single exponential distribution in [23], so the 95% confidence interval for the mean passage time is $\bar{\tau}\pm\frac{1.96}{\sqrt{n}}\bar{\tau}$. We include this confidence interval for comparison. In fact, one can use the Lilliefors test, which is an improvement of the Kolmogorov-Smirnov test, to return a test decision for the null hypothesis that the passage time comes from an exponential distribution. The result is 1 if the test rejects the null hypothesis at 1% significant level, 0 otherwise [9]. We applied the test to our simulation data for various parameters (we used $N=20,25,35,45,50,65,75,85,100,115,125,150$ and $\epsilon=0.1,0.25,0.5,0.75,1,2$ for a total of 72 parameter combinations). The test shows: when $\epsilon=0.1$ for all $N$ and $\epsilon=0.25$ with $N\geq 100,$ the test does not reject the null hypothesis at 1% significant level; for other parameter ranges, the test rejects the null hypothesis and that means the data do not fit an exponential distribution. Our result is different from Amitai’s result for large $N,$ where they claim that the passage time does not follow a single exponential distribution when $N>64.$ In fact, we can visualize the survival probability for the large $N$ case in Figure 4.

Table 1 and Figure 2 show that the WF theory overestimates the mean passage time, as can also be proved by using a variational principle [24]. They also verify our argument that the delta sink should perform better than the Heaviside sink. The analytical results underestimate the numerical integration at various points in the parameter range, by less than 15 percent for Heaviside sink and by less than 10 percent for delta sink. From Table 1 and Figure 2, we see that the delta sink results are in remarkable agreement with the simulation for small capture radius $\epsilon$. This implies that the space and time separation approximation relies mostly on the small capture radius assumption and not on the large number of beads. So, it is reasonable to use the WF theory as the simulation results when $\epsilon\ll 1$. At the same time, both analytical results predict that the WF theory has the mixed scaling law with $N\sqrt{N}$, $\frac{N\sqrt{N}}{\epsilon}$ and $N^{2}$. We use the WF theory for both sinks under the parameter range $0.1\leq\epsilon\leq 0.15$ and $100\leq N\leq 150$ to fit the scaling law. One could use simulation results to fit but the computational cost is enormous. The coefficient $d_{1}$ is estimated as 0.1225, very close to the analytically predicted result $\frac{\sqrt{\pi}}{6\sqrt{6}}\approx 0.1206$.

[TABLE]

With the scaling law at hand, we use the regression equation Eq. 35 and Eq. 36 to test other parameter ranges given by $\epsilon=0.25,0.75$ and $50\leq N\leq 100,150\leq N\leq 200$ (see Figure 3). When $\epsilon=0.25$, which is relatively small, the fit agrees remarkably well with the WF approximation for both sinks and in both ranges of $N$. However, when $\epsilon=0.75$, the fit starts to deviate from the WF theory for both sinks and the difference grows with $N.$ Since we know that the WF theory can overestimate the looping time and the predictions from the regression fit are larger than those of the WF theory, we conclude that the predictions of the regression fit are not accurate for $\epsilon=0.75.$ For $\epsilon=0.25$, if one omits the new term $N\sqrt{N}$, Figure 3 shows that for both sinks the regression lines deviate from the numerical integration points in this parameter range and also the values of the Akaike information criterion (AIC) are much higher. The AIC is a measure of the relative quality of a statistical model for a given set of data and and the lower AIC model is better [5]. This indicates that the $N\sqrt{N}$ term needs to be included in the model.

In the previous section, it was mentioned that WF approximates the denominator in the Laplace transform domain to get a semi-analytic form for the looping time and, as a result, the integrand is not a survival probability $\Sigma(t)$. However, the Volterra integral equation Eq. 22 can be solved numerically by the trapezoidal method. Correspondingly, the survival probability is found numerically and each moment of passage time can be calculated. The most useful moment of passage time, of course, is the mean.

Rewrite Eq. 22 to get

[TABLE]

The kernel $C(t)/C(\infty)$ is the renormalized sink-sink correlation function and $P_{0}q(t)$ is the function to solve for. The survival probability is $\Sigma(t)=1-\int_{0}^{t}P_{0}q(t^{\prime}){\rm d}t^{\prime}$. There is a difficulty because the kernel has a singular point at $t=0$ and the numerical integration is very stiff. With the perturbation result in Eq. 30, we know the order of the singularity is 1/2, so careful handling of the kernel function at short time range $t$ is needed.

We plot in Figure 4 the survival probability $\Sigma(t)$ predicted by the WF theory for two different values of $N$ and two different capture radius values $\epsilon$. As expected, the time integral of the survival probability recovers looping time from the WF theory in Eq. 23. In addition, we compare with the survival probability from the simulation data. This is computed by creating the histogram for the passage time and calculating the cumulative probability function for it. Then, the survival probability is one minus the cumulative probability function of the passage time. If the passage time follows a single exponential distribution, the survival probability will be an exponential function with the exponent given by the inverse of the looping time. For small capture radius $\epsilon=0.1,$ the survival probability function from WF agrees with the one from the simulation extremely well. It is interesting to note that the survival probability function obtained through simulation is an exponential function, even for large $N.$ This was also verified previously. The exponent of the exponential function is the inverse of the WF looping time. However, for large capture radius $\epsilon=1,$ the survival probability function from the WF theory is clearly no longer an exponential function and it deviates from the simulation results which also do not show exponential behavior.

7 Discussion and future work

We have reviewed the existing analytical theories for polymer looping in the Rouse model, including the SSS theory, the eigenvalue perturbation method, and the WF theory. Specifically, we reviewed the linkage between the J factor and the kinetics of looping by showing that J factor as the scaling of $N$ appears in the first order perturbation of the largest eigenvalue of the Fokker-Planck equation. We have rigorously investigated the conditions under which each asymptotic method holds and articulated the parameter range in which each of them is accurate. In particular, the main assumption for WF has been considerably relaxed, i.e, the solution has the space and time separation inside the sink for the Heaviside sink and on the reactive surface for the delta sink. Our results show that the WF theory as a non-Markovian method can reproduce other existing analytical theories and it has the largest parameter range where it remains accurate. In fact, as far as the capture radius is small, i.e, $\epsilon$ is small, the WF theory can predict the looping time extremely well. Moreover, under the same small $\epsilon$ condition, we are able to derive that the looping time follows the mixed scaling laws asymptotically and it is considered as the next order approximation of the Markovian approximation in powers of the capture radius $\epsilon$. Our findings have also been numerically verified.

We have to admit that further relaxing the time and space separation assumption systematically is not feasible analytically. In order to fully capture this non-Markovian effect, a considerably more involved method will be required. However, the WF theory as a general method of dimension reduction deserves more attention. It could be applied to any reversible OU processes with complicated absorbing boundary conditions. The WF theory can also be used to extract semi-analytical or even analytical expressions for the mean first passage time and this may provide guidelines for numerical simulations and physical experiments. Natural extensions beyond the Rouse model include the hydrodynamical interaction or heterogeneous spring coefficients. Both models can be studied through independent but not identical normal modes with different drift and diffusion terms. Then the end-end correlation function $\phi(t)$ will be different, but the rest of the analysis carries through the same.

Acknowledgments

FY thanks particularly Attila Szabo (NIH) and Nathan Baker (PNNL) for helpful discussions and suggestions, and continuous encouragement in the initial phase of this project. The work of FY was partially supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program, Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4), under Award Number DE-SC0009280. The work of PS was partially supported by the Pacific Northwest National Laboratory Laboratory Directed Research and Development (LDRD) Project “Multiscale modeling and uncertainty quantification for complex non-linear systems”.