Non-Markovian closure kinetics of flexible polymers with hydrodynamic interactions

TL;DR

This paper develops a non-Markovian theoretical model for polymer closure kinetics with hydrodynamic interactions, accurately matching simulations and revealing significant memory effects that influence cyclization times.

Contribution

It introduces a non-Markovian analysis based on preaveraged Zimm dynamics, improving over Markovian approximations and deriving new scaling laws for cyclization times.

Findings

Markovian models overestimate cyclization times by up to a factor of 2.

Scaling law for large N: T ~ N^{3/2} log(N/b^2).

Non-Markovian approach aligns closely with numerical simulations.

Abstract

This paper presents a theoretical analysis of the closure kinetics of a polymer with hydrodynamic interactions. This analysis, which takes into account the non-Markovian dynamics of the end-to-end vector and relies on the preaveraging of the mobility tensor (Zimm dynamics), is shown to reproduce very accurately the results of numerical simulations of the complete non linear dynamics. It is found that Markovian treatments, based on a Wilemski-Fixman approximation, significantly overestimate cyclization times (up to a factor 2), showing the importance of memory effects in the dynamics. In addition, this analysis provides scaling laws of the mean first cyclization time (MFCT) with the polymer size $N$ and capture radius $b$, which are identical in both Markovian and non-Markovian approaches. In particular, it is found that the scaling of the MFCT for large $N$ is given by $T\sim N^{3/2}\ln (N/b^2)$, which differs from the case of the Rouse dynamics where $T\sim N^{2}$. The extension to the case of the reaction kinetics of a monomer of a Zimm polymer with an external target in a confined volume is also presented.

Equations

$${\mathbf{F}}_j=-k(2{\mathbf{x}}_j-{\mathbf{x}}_{j+1}-{\mathbf{x}}_{j-1})=-\sum _{k=1}^N{M}_{jk}{\mathbf{x}}_k,$$

$${\dot{\mathbf{x}}}_i=\sum _{j=1}^N{\mathbb{D}}_{ij}\cdot {\mathbf{F}}_j+{\mathbf{\zeta }}_i(t),$$

$${\mathbf{D}}_{ij}=$$

$$\frac{k_{B}T}{8\pi \eta r_{ij}}\left[\left(1+\frac{2a^2}{3r_{ij}^2}\right)\mathbb{I}+\frac{r_{ij}\otimes r_{ij}}{r_{ij}^2}\left(1-\frac{2a^2}{r_{ij}^2}\right)\right],$$

$${\mathbf{D}}_{ij}$$

$${\mathbf{D}}_{ii}$$

$${\bar{\mathbf{D}}}_{ij}\simeq \frac{k_{B}T}{6\pi \eta }\left[\frac{{\delta }_{ij}}{a}+\frac{(1-{\delta }_{ij})}{l_0}\sqrt{\frac{2}{\pi |i-j|}}\right]\mathbb{I}.$$

$$⟨[r_i(t)-r_i(0){]}^2⟩\sim \{\begin{array}{rlrlr}& D_{G}t& \text{for}& t\ll {\tau }_0& \\ & A{t}^{2/3}& \text{for}& {\tau }_0\ll t\ll {\tau }_0{N}^{3/2}& \\ & D_{G}t/\sqrt{N}& \text{for}& {\tau }_0{N}^{3/2}\ll t& \end{array}$$

$$P(\{\mathbf{x}$$

$${\int }_0^{t}d{t}^{\prime }\int d\{{\mathbf{x}}^{\prime }\}f(\{{\mathbf{x}}^{\prime }\},t^{\prime })P(\{\mathbf{x}\},t-t^{\prime }|\{{\mathbf{x}}^{\prime }\}),$$

$$T$$

$${\int }_0^{\infty }dt\left[\int d\Omega P(\{x\},t|{\pi }_{\Omega })-P(\{x\},t|\text{ini},0)\right].$$

$$T_{\text{WF}}={\int }_0^{\infty }\frac{dt}{[1-\varphi (t{)}^2{]}^{3/2}}\left\{e^{-b^2\varphi (t{)}^2/[2\psi (t)]}-\frac{Z(b,\psi (t))}{Z(b,L^2)}\right\},$$

$$\varphi (\tau )=\frac{⟨x_{\mathrm{ee}}(t+\tau )x_{\mathrm{ee}}(t)⟩}{⟨x_{\mathrm{ee}}(t)x_{\mathrm{ee}}(t)⟩},$$

$$L^2=⟨x_{\mathrm{ee}}(t{)}^2⟩=(N-1)l_0^2,$$

$$\psi (t)=\mathrm{Var}(x_{\mathrm{ee}}(t)|x_{\mathrm{ee}}(0)=x_{\mathrm{ee}}^0)=L^2[1-\varphi (t{)}^2],$$

$$\displaystyle\int_{0}^{\infty}\frac{dt}{\psi^{3/2}}\Bigg{\{}$$

$$\displaystyle+\left(\frac{\beta_{i}}{\psi}-\frac{\alpha_{i}}{L^{2}}\right)\frac{b^{3}e^{-b^{2}/2\psi}}{3Z(b,L^{2})}\Bigg{\}}=0.$$

$${\beta }_i(t)=\mathrm{Cov}(x_i(t),x_{\mathrm{ee}}(t)|x_{\mathrm{ee}}(0)=x_{\mathrm{ee}}^0),$$

$${\alpha }_i={\mathrm{Cov}}_{\mathrm{stat}}(x_i,x_{\mathrm{ee}}),$$

$${\mu }_i^{\pi }(t)=⟨{\mathbf{x}}_i(t+t^{\ast })\cdot \frac{{\mathbf{r}}_{\mathrm{ee}}(t^{\ast })}{|{\mathbf{r}}_{\mathrm{ee}}(t^{\ast })|}⟩,$$

$$R_{\pi }(t)=⟨{\mathbf{r}}_{\mathrm{ee}}(t+t^{\ast })\cdot \frac{{\mathbf{r}}_{\mathrm{ee}}(t^{\ast })}{|{\mathbf{r}}_{\mathrm{ee}}(t^{\ast })|}⟩={\mu }_N^{\pi }(t)-{\mu }_1^{\pi }(t),$$

$${\partial }_t{\mu }_i^{\pi }(t)=-\sum _{j,k=1}^N{D}_{ij}{M}_{jk}{\mu }_k^{\pi }(t),{\mu }_i^{\pi }(0)=m_i^{\pi }.$$

$$\frac{T}{L^3}={\int }_0^{\infty }\frac{dt}{{\psi }^{3/2}}\left[e^{-R_{\pi }^2/2\psi }-\frac{Z(b,\psi )}{Z(b,L^2)}\right].$$

$$T$$

$$\varphi (t)\simeq \sum _{p\text{odd}}\frac{8}{p^2{\pi }^2}{e}^{-p^{3/2}t/{\tau }_1},$$

$${\tau }_1=\frac{3\eta l_0^3}{\sqrt{\pi }{k}_{B}T}{N}^{3/2}.$$

$$\psi (t)=⟨[x_{\mathrm{ee}}(t)-x_{\mathrm{ee}}(0){]}^2⟩\simeq \frac{8\Gamma (1/3)}{3^{2/3}{\pi }^{5/3}}{\left(\frac{k_{B}Tt}{\eta }\right)}^{2/3}$$

$$T={\tau }_1{\int }_0^{\infty }d\tau \left[\frac{e^{-{\tilde{b}}^2\Phi (\tau {)}^2/2\Psi (\tau )}}{\Psi (\tau {)}^{3/2}}-1\right]\equiv {\tau }_{1}f(\tilde{b})$$

$$\kappa =8\Gamma (1/3)/{\pi }^2.$$

$$f(\tilde{b})\simeq {\int }_0^{\frac{ϵ}{{\tilde{b}}^3}}du\frac{e^{-\frac{1}{2\kappa u^{2/3}}}}{{\kappa }^{3/2}u}+{\int }_{ϵ}^{\infty }d\tau \left(\frac{1}{[1-{\Phi }^2{]}^{3/2}}-1\right)$$

$$f(\tilde{b})$$

$$+{\int }_0^{\infty }d\tau \left(\frac{1}{[1-\Phi (\tau {)}^2{]}^{3/2}}-1-\frac{\theta (1-\tau )}{{\kappa }^{3/2}\tau }\right)$$

$$f(\tilde{b})\simeq \frac{3}{{\kappa }^{3/2}}(-\ln \tilde{b}+\mathrm{0.721...})$$

Full text

I Introduction

Reactions involving polymers are ubiquitous in nature. Among them, reactions of closure of linear chains are of particular interest since they are involved in a number of chemical and biological processes. Examples cover gene regulation by the formation of RNA-hairpinsBonnet et al. (1998); Wang and Nau (2004) or DNA-loopsWallace et al. (2001); Allemand et al. (2006), the folding of polypeptidesLapidus et al. (2000); Möglich et al. (2006), as well as the appearance of cycles in synthetic polymers Gooden et al. (1998); Zheng et al. (2011); Burgath et al. (2000). In the diffusion-controlled regime, the kinetics of contact formation strongly depends on the complex dynamics of the reactive monomers. As a result of the collective dynamics of all the monomers in the chain, the motion of a single monomer is often subdiffusiveDoi and Edwards (1988) and presents non-Markovian featuresPanja (2010) (i.e. memory effects), which lead to nontrivial reaction kineticsDe Gennes (1982) even for the simplest models of polymers.

On the theoretical level, various approaches have been proposed to quantify the kinetics of polymer closureWilemski and Fixman (1974a); Wilemski and Fixman (1974b); Szabo et al. (1980); Sokolov (2003); Likthman and Marques (2006); Guerin et al. (2012); Friedman and O’Shaughnessy (1989); Chakrabarti (2012); Amitai et al. (2012); Toan et al. (2008). An important step has been provided by Wilemski and FixmanWilemski and Fixman (1974a); Wilemski and Fixman (1974b), who made a local equilibrium assumption, thereby replacing the non-Markovian problem by an effective Markovian approach. Other theories include the so-called SSS analytical approachSzabo et al. (1980) , which neglects important aspects of the polymer dynamics, and the perturbative renormalization group theoryFriedman and O’Shaughnessy (1989), which provides results at leading order in the parameter $\epsilon=4-d$, with $d$ the spatial dimension. Recent works have improved these approaches in the case of the Rouse chain (i.e. flexible without hydrodynamic interactions) by introducing several methods: (i) a refined way to take into account the memory of the initial configurationsSokolov (2003), (ii) an exact formal iterative resolution scheme in one dimension Likthman and Marques (2006), (iii) a strong localized perturbation analysis Amitai et al. (2012), and (iv) an approach based on the calculation of the distribution of chain configurations at the instant of cyclizationGuerin et al. (2012, 2013a, 2013b); Bénichou et al. (2015), which has a strong influence on the contact kinetics.

The Rouse model, however, provides an incorrect description of polymer chain dynamics in dilute solutions, where hydrodynamic interactions are long-ranged and deeply modify the motion of the polymerZimm (1955); Doi and Edwards (1988) as well as the closure kineticsFriedman and O’Shaughnessy (1989); Chakrabarti (2012). Since hydrodynamic interactions are non-linear, most of available analytical treatments involve approximations. A standard approach consists in using a pre-averaged form introduced by ZimmZimm (1955), which is known to be accurate Doi and Edwards (1988). In the context of closure kinetics, Ortiz-Repiso et al.Ortiz-Repiso et al. (1998) found by simulations that the Wilemski-Fixman treatment of the cyclization kinetics of Zimm chains systematically overestimates the reaction times, in a way that cannot be attributed to the simplified treatment of these interactions. This suggests that non-Markovian effects, which have not been described before in the presence of hydrodynamic interactions, play an important role in the cyclization kinetics.

This paper is devoted to the theoretical description of the cyclization kinetics of chains with hydrodynamic interactions in the diffusion controlled regime. This problem has already been discussed in the Wilemski-Fixman approachChakrabarti (2012); Ortiz-Repiso et al. (1998) and the perturbative renormalization group theoryFriedman and O’Shaughnessy (1989). However, a non-Markovian description of the closure kinetics has not been considered so far. The first goal consists in describing these effects by adapting a theory proposed recently Guerin et al. (2012, 2013a, 2014), which was so far restricted to the case of flexible and semi-flexible free draining chains. One could expect that the treatment of pre-averaged interactions is not valid for the problem of chain closure, since interactions averaged over equilibrium configurations could differ from the actual (non-preaveraged) interactions for closed chains, which could then play an important role for the contact kinetics. We will show however that this effect is relatively small: indeed, calculating the average shape of the polymer at the instant of cyclization for pre-averaged interactions yields a reaction time that is in quantitative agreement with the results of Brownian dynamic simulations performed without pre-averaging. The second goal of the paper is to derive new scaling laws for the mean cyclization time in the limit of small capture radius and long chains. In particular, we prove that for long chains, the mean cyclization time $T$ scales as $T\sim N^{3/2}\log(N/b^{2})$, with $N$ the number of monomers and $b$ the capture radius, for both the Wilemski-Fixman approach and the non-Markovian theory, with however a different prefactor. In addition, we predict scaling laws for intermolecular reactions, by extending the formalism to the case of the search of an external target of size $r$ in a confined volume $V$. We show that in this case, the mean reaction time scales as $T\sim V\log(r/\sqrt{N})$ for long chains.

The outline of the paper is as follows. In section II, we recall briefly some important results of the Zimm model. In section III, we show how to adapt the recent non-Markovian theory Guerin et al. (2012, 2013a, 2014) of cyclization kinetics to the Zimm model. We compare the numerical solution of the model to both Brownian dynamics simulations and to the results of a Wilemski-Fixman approach. Section IV is devoted to the derivation of scaling laws for the mean cyclization time with the chain length and the capture radius, for both the Markovian and the non-Markovian theories. Finally, we present in section V the direct application of our calculations to the case of the reaction of one monomer with an external target in confinement and derive explicit asymptotic formulas for the reaction time for this problem.

II Polymer dynamics with hydrodynamic interactions

We consider the dynamics in a three-dimensional ($3d$) space of a polymer represented by $N$ beads at positions $(\mathbf{x}_{1},...,\mathbf{x}_{N})$ (in this paper quantities in bold represent vectors or tensors in the $3d$ space). The beads are linked by springs of stiffness $k$, so that the force $\mathbf{F}_{j}$ exerted on a bead $j$ by the neighboring beads is

[TABLE]

where we used the convention $\mathbf{x}_{0}=\mathbf{x}_{1}$ and $\mathbf{x}_{N+1}=\mathbf{x}_{N}$ and $M$ is the (Laplacian) tridiagonal connectivity matrix. This force $\mathbf{F}_{j}$ is balanced by the force due to the fluid solvent, and from Newton’s third law a force $\mathbf{F}_{j}$ is exerted on the fluid at the position of bead $j$. This force generates a velocity field in the fluid that influences the motion of all other beads. We introduce a non-isotropic mobility tensor $\mathbf{D}_{ij}$ to describe these hydrodynamic interactions, such that the average velocity of the fluid at the position of the bead $i$ is $\mathbf{D}_{ij}\mathbf{F}_{j}$. The dynamics in the presence of hydrodynamic interactions then follows the Langevin equation in the overdamped limit

[TABLE]

where $\bm{\zeta}_{i}(t)$ is a stochastic Gaussian white noise term whose amplitude follows from the fluctuation-dissipation relation, $\langle\bm{\zeta}_{i}(t)\otimes\bm{\zeta}_{j}(t^{\prime})\rangle=2k_{B}T\mathbb{D}_{ij}\delta(t-t^{\prime})$, with $k_{B}T$ the thermal energy. We also introduce $l_{0}=\sqrt{k_{B}T/k}$ the typical bond length and $\tau_{0}=6\pi\eta a/k$ the single bond characteristic relaxation time, where $\eta$ is the fluid viscosity and $a$ the monomer radius.

Different choices of mobility matrix $\mathbf{D}_{ij}$ exist. The simplest one is the Oseen tensor, describing the fluid velocity induced by a punctual force in the fluid, and therefore valid for quasi-punctual monomers. However, for very small distances between beads, the Oseen tensor may not be positive definite, making the dynamics (2) unphysical due to the appearance of negative relaxation times Ermak and McCammon (1978). Instead, for numerical purposes, we use the so-called Rotne-Prager tensor, which does not have this problem. We denote by $\mathbb{r}_{ij}=\mathbf{x}_{j}-\mathbf{x}_{i}$ the vector from the particle $i$ to the particle $j$ and $r_{ij}=|\mathbb{r}_{ij}|$. The elements of the Rotne-prager tensor are then given as followsGeyer and Winter (2009): when $i\neq j$, and when the beads $i,j$ do not overlap ($r_{ij}>2a$) ,

[TABLE]

while for $i\neq j$ and $r_{ij}<2a$,

[TABLE]

and finally, when $i=j$

[TABLE]

where $\mathbb{I}$ is the $3\times 3$ identity tensor. Eqs. (1,2,3,4,5) define the polymer dynamics that will be analyzed in this paper. For numerical purposes, we will make use of the choice $a=0.25\,l_{0}$ Osaki (1972); Osaki et al. (1971).

The non-linear dependence of the mobility matrix $\mathbf{D}_{ij}$ on the positions $\mathbf{x}_{j}$ makes the Langevin equation (2) very difficult to solve. In particular it does not admit Gaussian solutions. In his pioneering workZimm (1955), Zimm overcame this difficulty by replacing the mobility matrix $\mathbf{D}_{ij}$ by its average value over the equilibrium distribution, thereby making the equation linear. It has been shown that this approximation catches the main physics of hydrodynamic interactions Zimm (1955); Doi and Edwards (1988). In the model introduced above, it would be therefore relevant to use a pre-averaged version of the Rotne-Prage tensor rather than the Oseen tensor. However, it is shown in Appendix A that these tensors differ at most by a few percents for the parameters that we use. Hence, it is sufficient to use for analytical calculations the pre-averaged form of the Oseen tensor, given by Zimm (1955); Doi and Edwards (1988):

[TABLE]

Note that the pre-averaged Oseen tensor is isotropic.

The effects of hydrodynamic interactions on polymer dynamics have already been studied in details. In particular, it is well-known that any monomer of a long polymer chain performs a subdiffusive motion at intermediate time scales, which results from the collective dynamics of all the monomers. More precisely, for $N\gg 1$ we have three different regimes for the mean-square displacement of a monomer Teraoka (2002) :

[TABLE]

where $D_{G}=k_{B}T/\eta l_{0}$, $A=(k_{B}T/\eta)^{2/3}$ (see Eq. (24)) and all the numerical constants have been omitted. We can already note that in the subdiffusion regime, the dimension of the random walk of a monomer (defined such that $\langle[r(t)-r(0)]^{2}\rangle\sim t^{2/d_{w}}$) is given by $d_{w}=3$. The process in $3d$ is therefore marginally recurrent, whereas it is transient for both short and long time scales (see ref.Condamin et al. (2007)). This will play a key role in the determination of the scaling laws that we present below.

The goal of this paper is to characterize the Mean First Cyclization Time (MFCT, denoted $T$), defined as the average first time for which the distance $|\mathbf{x}_{N}-\mathbf{x}_{1}|$ between the two chain ends becomes smaller than a capture radius $b$. Initially the polymer is in an equilibrium configuration, with the constraint that $|\mathbf{x}_{N}-\mathbf{x}_{1}|>b$. The effect of the complex monomer dynamics resulting from the hydrodynamic interactions will be analyzed with both numerical simulation methods (which we describe in App. D) and analytical means (see section III).

III Theories of cyclization kinetics

We present here a theoretical approach that enables the determination of the mean cyclization time $T$, which can be seen as the mean first passage time of a non-Markovian problem. The approach consists in adapting a recent theory, developed so far for free-draining chainsGuerin et al. (2013b, a, 2014), in which the MFCT is expressed in terms of the conformational distribution of chains at the very instant of cyclization, whose moments are computed by solving a set of self-consistent equations.

The main steps of this approach can be summarized as follows. We introduce the joint probability density $f(\{\mathbf{x}\},t)$ that the contact is made for the first time at time $t$ and that, at this first passage event, the macromolecule has a configuration described by the set of positions $\{\mathbf{x}\}=(\mathbf{x}_{1},\mathbf{x}_{2}...)$. We partition the trajectories that lead to a configuration $\{\mathbf{x}\}$ (in which the contact condition is satisfied) into two steps, the first step consisting in reaching the target for the first time at $t^{\prime}$, and the second step consisting in reaching the final configuration $\{\mathbf{x}\}$ in a time $t-t^{\prime}$. The mathematical formulation of this decomposition of events is

[TABLE]

where $d\{\mathbf{x}\}\equiv d\mathbf{x}_{1}d\mathbf{x}_{2}...d\mathbf{x}_{N}$, $P(\{\mathbf{x}\},t|\{\mathbf{x^{\prime}}\})$ is the probability of $\{\mathbf{x}\}$ at $t$ starting from $\{\mathbf{x^{\prime}}\}$ at $t=0$ while $P(\{\mathbf{x}\},t|\mathrm{ini},0)$ is the probability of $\{\mathbf{x}\}$ at $t$ starting from the initial conditions at $t=0$ (in which the chain is at equilibrium, with the condition that the reactive monomers are not necessarily in contact). Next, taking the Laplace transform of (8) and expanding for small values of the Laplace variable, we obtainGuerin et al. (2013a)

[TABLE]

In this equation, $\Omega$ represents the angular directions parametrized by $\theta$ and $\varphi$ in spherical coordinates, with the normalization $\int d\Omega=1$; $\pi_{\Omega}(\{\mathbf{x}\})$ represents the probability distribution of configurations at the very instant of cyclization, given that the angular direction of the end-to-end vector at first contact is $\Omega$, and $P(\{\mathbb{x}\},t|\pi_{\Omega})$ is the probability of $\{\mathbf{x}\}$ at $t$ starting from the distribution $\pi_{\Omega}$ at initial time. Equation (9) is exact and does not depend on the particular hypotheses of chain dynamics. It is derived in details in Ref.Guerin et al. (2013a). However, this equation cannot be solved explicitly to the best of our knowledge and approximations have to be introduced.

The first simplifying step is to approximate the dynamics by a Gaussian dynamics in order to be able to evaluate the propagators appearing in (9); hence we consider only the Zimm dynamics with the pre-averaged mobility tensor. Next, the simplest approach is to make a Markovian approximation, which consists in neglecting any memory effect by assuming that the distribution $\pi_{\Omega}$ is the equilibrium distribution conditional to $\mathbf{x}_{N}-\mathbf{x}_{1}=b\hat{\mathbb{u}}_{r}(\Omega)$ (with $\hat{\mathbb{u}}_{r}(\Omega)$ the unit vector pointing in the direction $\Omega$). This corresponds to the so-called Wilemski-Fixman approximationPastor et al. (1996); Wilemski and Fixman (1974a); Guerin et al. (2013a). Introducing this approximation into Eq. (8), integrating over all configurations and taking the long time limit lead to the estimate $T_{\mathrm{WF}}$ of the MFCT Guerin et al. (2013a)

[TABLE]

where $\phi(t),L^{2},\psi(t)$ characterize the dynamics of the end-to-end vector $\mathbf{r}_{\mathrm{ee}}=\mathbf{x}_{N}-\mathbf{x}_{1}$, which is assumed to be Gaussian, and $Z(b,x)=\int_{b}^{\infty}dR_{0}e^{-R_{0}^{2}/2x}$. The function $\phi$ is the normalized temporal auto-correlation function of any coordinate of the end-to-end vector $\mathbf{r}_{\mathrm{ee}}$,

[TABLE]

where the spatial coordinates of $\mathbf{r}_{\mathrm{ee}}$ are $(x_{\mathrm{ee}},y_{\mathrm{ee}},z_{\mathrm{ee}})$, $L^{2}$ is the equilibrium mean-square length

[TABLE]

and $\psi(t)$ is the Mean-Square-Displacement of any coordinate of $\mathbf{r}_{\mathrm{ee}}$ when the initial value of $\mathbf{r}_{\mathrm{ee}}$ is fixed, and is easily shown to be related to $\phi$ by

[TABLE]

where we denote $\mathrm{Var}(y|B)$ the variance of the variable $a$ given the event $B$ is realized.

Going beyond the Wilemski-Fixman approximation requires a more precise description of the distribution $\pi_{\Omega}$. Here, our key hypothesis is to assume that $\pi_{\Omega}$ is a multivariate Gaussian distribution. For symmetry reasons, the average of $\mathbb{x}_{i}$ over $\pi_{\Omega}$ is along $\hat{\mathbb{u}}_{r}(\Omega)$. We denote this quantity $\mathbb{E}_{\pi_{\Omega}}(\mathbf{x}_{i})=m_{i}^{\pi}\hat{\mathbb{u}}_{r}(\Omega)$. We make the additional assumption that the covariance matrix of $\pi_{\Omega}$ is the same covariance matrix that characterizes the equilibrium distribution of chains with the constraint $\mathbf{r}_{\mathrm{ee}}=b\hat{\mathbb{u}}_{r}(\Omega)$. Such an approximation is not necessary but largely simplifies the calculations. Previous studies on free-draining chainsGuerin et al. (2013b, 2012) have revealed that releasing this "stationary covariance approximation" only slightly improves the estimate of the cyclization time at an important calculation cost. A set of self-consistent equations that defines the $m_{i}^{\pi}$’s is found by multiplying (9) by $\mathbf{x}_{i}\,\delta(\mathbf{x}_{N}-\mathbf{x}_{1}-\mathbf{r}^{*})$ (for a fixed $\mathbf{r}^{*}$ satisfying $|\mathbf{r}^{*}|<b$) and integrating over all configurations. Adapting existing calculations for free-draining chainsGuerin et al. (2013a), the resulting equation is (see the complete derivation in Appendix C)

[TABLE]

Here, we have denoted by $\mu_{i}^{\pi}(t)$ the average of $\mathbf{x}_{i}$ at time $t$ after the instant of cyclization in the direction $\hat{\mathbf{u}}(\Omega)$, $R_{\pi}=\mu_{N}^{\pi}-\mu_{1}^{\pi}$, $\beta_{i}(t)$ is the covariance between $x_{i}(t)$ and $x_{\mathrm{ee}}(t)$ when $x_{\mathrm{ee}}(0)$ is fixed, while $\alpha_{i}$ is the covariance between $x_{i}$ and $x_{\mathrm{ee}}$ at equilibrium. In other words, $\alpha_{i}$ and $\beta_{i}$ characterize the dynamics and the equilibrium of the motion starting from a constrained equilibrium state,

[TABLE]

(note that $\beta_{i}$ does not depend on $x_{\mathrm{ee}}^{0}$) and $\mu_{i}^{\pi},R_{\pi}$ characterize the motion of the chain in the future of the first contact,

[TABLE]

where $t^{*}$ is the first cyclization time. The time evolution of $\mu_{i}^{\pi}$ follows from the Langevin equation (2):

[TABLE]

Here, the unknowns $m_{i}^{\pi}$ are contained in the $\mu_{i}^{\pi}$ as initial conditions of the dynamical system (19), while $\alpha_{i},\beta_{i},\phi,\psi$ characterize the dynamics of the chain and are analytically calculated in Appendix B. We stress that equation (14) is fully general for a $3d$ isotropic Gaussian non-Markovian process, and does not depend on the particular structure of the mobility matrix $D$ or the connectivity matrix $M$. Then, the expression of the mean reaction time is obtained by multiplying (9) by $\delta(\mathbf{x}_{N}-\mathbf{x}_{1})$ and integrating over all configurations, leading to

[TABLE]

To conclude, Eq. (19) provides a system of $N$ equations for the unknowns $\mu_{i}^{\pi}(t)$. Note that this system is of rank $N-2$, because the polymer center of mass can be set arbitrarily, and the constraint $m_{N}^{\pi}-m_{1}^{\pi}=b$ must hold. Solving this system allows us to compute the function $R_{\pi}(t)$, and then to calculate the mean first-passage time $T$ from (20).

We have solved numerically these equations for different values of the number of monomers $N$ and capture radius $b$. The results are shown in Fig. 1. We clearly see that the non-Markovian theory accurately predicts the mean-cyclization time, whereas the Wilemski-Fixman approach quickly fails when the capture radius is large enough. More precisely, the theory is in quantitative agreement with the simulations performed with preaveraging. Remarkably, the theory is also in quantitative agreement with the simulations performed with the exact (non pre-averaged) interactions, in which the full non-linear stochastic dynamics is taken into account when the capture radius is not small.

The average positions $m_{i}^{\pi}$ of the monomers at the instant of cyclization are shown in Fig. 2, which shows that the Zimm chain is significantly more elongated in the direction of the end-to-end vector at the reaction than in an equilibrium looped configuration. This allows the chain to perform cyclization more rapidly than if it had to equilibrate. The monomers neighboring the reactive ones are on average outside the reactive region. Note that the small shift between the theoretical prediction and the simulation indicates that the theory is not exact. The shape of the curve of the $m_{i}^{\pi}$ is however significantly closer to the simulations for the non-Markovian theory than in the case of the Wilemski-Fixman approximation.

For small capture radius, the cyclization times obtained numerically without pre-averaging are slightly larger than those obtained numerically with pre-averaging, as well as those predicted by the non-Markovian theory. This can be understood from the fact that, for very small capture radius, the monomer motion at short time scales plays a key role. For small $t$, the MSD of the end-to-end distance reads $\psi(t)\simeq 4(D_{11}+D_{1N})\,t$, and the motion is diffusive. With the pre-averaging, $D_{1N}\sim 1/\sqrt{N}$ is negligible compared to $D_{11}$. However, for chains that are close to form a loop, the actual mobility tensor reads $D_{1N}\sim D_{11}$ since the distance between the end-beads is small. Hence, the effective diffusion constant at small time scales is not correctly estimated with the pre-averaging procedure, and we can expect discrepancies for the cyclization times due to the pre-averaging for small capture radius, as observed in Fig. 1.

IV Asymptotic behavior of the cyclization time

IV.1 Scalings of the MFCT in the Wilemski-Fixman theory

Let us now derive the different scaling behaviors of the MFCT with $N$ and the capture radius $b$, first in the Wilemski-Fixman approximation. The MFCT are then computed with Eq. (10). We first consider the limit of small capture radius, $b\rightarrow 0$ at fixed chain length $N$. If $b=0$, the integral (10) is divergent due to the linear behavior of the MSD function $\psi$ at short times. Introducing the effective short time diffusion coefficient $D_{s}$ (such that $\psi\simeq_{t\rightarrow 0}2D_{s}t$), and considering that the integral (10) is governed by the short time regime, we get

[TABLE]

As discussed above, the actual value of $D_{s}$ could be quantitatively underestimated by the pre-averaging procedure, which computes it for equilibrium rather than looped configurations. Nevertheless, $D_{s}$ remains of the order of $k_{B}T/(6\pi\eta l_{0})$. The scaling (21) is similar to that appearing for free-draining flexible (Rouse) chains Amitai et al. (2012); Pastor et al. (1996); Guerin et al. (2012, 2013a). It is known that such scaling does not depend on the structure of the chain since it is the same for semi-flexible chains Guerin et al. (2014) and for hyperbranched structures; here it is clear that it also appears in the presence of hydrodynamic interactions. Eq. (21) means that in the small capture radius regime, the MFCT is (up to a prefactor) the time needed for a diffusive particle to find a target of size $b$ in a confining volume $L^{3}$ Condamin et al. (2007); Bénichou and Voituriez (2014), and does not result from the collective dynamics of the monomers.

We now consider the scaling of the MFCT in the limit of long chains. In this limit, we use the commonly used dominant diagonal approximationDoi and Edwards (1988); Chakrabarti (2012), in which the orthogonal matrix $Q$ that diagonalizes the Laplacian matrix $M$ is assumed to diagonalize also the product $DM$. One can check in Fig. 3 that the MSD function $\psi$ calculated with this approximation is very close to its exact value even for moderately large $N$. In this approximation, the correlation function in the continuous limit can be shownDoi and Edwards (1988); Chakrabarti (2012) to be

[TABLE]

where $\tau_{1}$ is the slowest relaxation time scale of the chain,

[TABLE]

In the limit $t/\tau_{1}\ll 1$, the sum (22) can be replaced by a continuous integral, leading to the identification of the short time subdiffusive behavior

[TABLE]

which shows no dependence on the bond size $l_{0}$ or the chain length $N$. Here $\Gamma(\cdot)$ represents the Gamma function. Defining $\tau=t/\tau_{1}$, and introducing the rescaled functions $\phi(t)=\Phi(t/\tau_{1})=\Phi(\tau)$ and $\psi(t)=l_{0}^{2}N\Psi(\tau)$, with $\Psi$ and $\Phi$ independent on $N$ in the continuous limit, a simple change of variable in Eq. (10) leads to

[TABLE]

with $\tilde{b}=b/(l_{0}\,\sqrt{N})$ and $f$ a dimensionless function. Equation (25) finally gives the behavior of the mean cyclization time for continuous Zimm chains.

We focus now on the asymptotics of $f$ for small rescaled capture radius $\tilde{b}$, which is highly dependent of the short time behavior of $\Psi\simeq\kappa\tau^{2/3}$, where $\kappa$ is deduced from Eqs. (23,24) and reads

[TABLE]

Because of this subdiffusive behavior, replacing $\tilde{b}$ by zero in (25) leads to a divergent integral, meaning that $f$ diverges with $\tilde{b}$ for small capture radii. The reason for this divergence is that, when $\tau\sim b^{3}$, the term $\tilde{b}^{2}/\Psi$ becomes of order $1$ and cannot be replaced by [math]. This suggests to introduce an intermediate "time" scale $\varepsilon$ such that $\tilde{b}^{3}\ll\varepsilon\ll 1$, and to split the integral (25) into two contributions, leading to

[TABLE]

in which we have used the short time expression (24) for $\tau<\varepsilon$, and simply set $\tilde{b}=0$ in the contribution coming from the large times. In the joint limit $\varepsilon\rightarrow 0$ and $\varepsilon/\tilde{b}^{3}\rightarrow\infty$, the integral (27) can be recast under the form

[TABLE]

with $\theta(\cdot)$ the Heaviside step function. All integrals appearing in this equation are convergent, and their numerical evaluation leads to the scaling form

[TABLE]

In other words, we have identified the asymptotic scaling for the mean cyclization time of long Zimm chains

[TABLE]

We note that the weak logarithmic dependence on the size of the target is due to the fact that the motion of a monomer is a marginally compact process in this regime (the dimension of the walk satisfies $d_{w}=d=3$) Condamin et al. (2007); Bénichou and Voituriez (2014). The scaling $T\sim N^{3/2}\log N$ is in particular consistent with the renormalization group approach of Friedman and O’Shaughnessy (1989); Chakrabarti (2012). The effect of hydrodynamic interactions is clearly visible, since this scaling of $T$ with $N$ is very different from the case of flexible chains without hydrodynamic interactions, where $T$ for long chains scales as $N^{2}$ and is independent of the capture radius. Finally, let us note that a non-Markovian analysis reproduces the scaling of the MFCT of Eq. (30), with however a different numerical prefactor.

V Case of an external target

The formalism presented above can be adapted to the case of an external target, which models an inter-molecular reaction between one monomer at the end of the chain and a fixed target of radius $b$ in a confined volume $V$. This case has already been studied for a Rouse chain Guerin et al. (2013b), and we here briefly describe how the case of a Zimm chain can be treated. Assuming that the target is at the origin, the quantity of interest is $\mathbb{R}=\mathbb{R}_{1}$ if we assume that the reactive monomer is the first one. We present below the scaling of the mean reaction-time $T$ with the target, for a chain confined in a large volume $V$ in Wilemski-Fixman approximation. Eq.(9) remains valid, and by integrating it over all configurations such that $\mathbb{x}_{1}=\mathbf{0}$, one gets :

[TABLE]

The propagators appearing above are a priori the propagators in a confined volume. In the limit of large volume, we argue that $P_{\text{stat}}(\mathbb{x}_{1}=0)$ is simply equal to $1/V$. The other propagators are well defined in infinite volume. The Wilemski-Fixman approximation then consists in assuming that the reactive conformations given by $\pi$ are simply the equilibrium configurations such that $|\mathbb{x}_{1}|=b$. Hence, one has :

[TABLE]

where $R_{0}$ is the initial position of the reactive monomer and $\psi(t)=\langle[x_{1}(t)-x_{1}(0)]^{2}\rangle$ the MSD of the reactive monomer, whose asymptotic behavior is

[TABLE]

where $D_{\mathrm{cm}}$ is the large time center-of-mass diffusion coefficient. For large $N$, the following expressions are found Teraoka (2002) :

[TABLE]

The expression of the mean reaction time $T$ in the Wilemski-Fixman approximation can then be deduced from Eq. (32) and yields in the limit of chains starting far from the target ($R_{0}\rightarrow\infty$)

[TABLE]

where $\nu$ is a numerical coefficient that reads

[TABLE]

As in the case of cyclization dynamics, the scaling of the mean reaction time with $N$ is different from the case of a Rouse chain. In particular, for intermediate capture radii $r$, the very weak (logarithmic) dependence on $N$ leads to reaction times that can be significantly shorter than for a Rouse polymer, for whichGuerin et al. (2013b) $T\propto V\sqrt{N}$. The non-Markovian analysis of the reaction time with an external target closely follows the steps developed above for the determination of the MFCT and is not described here. As in the case of the cyclization dynamics, the scaling laws with $r$ and $N$ would be the same as predicted by the Wilemski-Fixman treatment (up to numerical prefactors), given in Eq. (36).

VI Conclusion

In conclusion, we have presented in this paper a theoretical analysis of the closure kinetics of a polymer with hydrodynamic interactions. We have provided both a Markovian (Wilemski-Fixman) and a non-Markovian analytical approach, based on a recent method introduced in Guerin et al. (2012, 2013a, 2013b); Bénichou et al. (2015). Although our theory relies on the preaveraging of the mobility tensor (Zimm dynamics), it is shown to reproduce very accurately the results of numerical simulations of the complete non linear dynamics. It is found that the Markovian treatment significantly overestimates cyclization times (up to a factor 2), showing the importance of memory effects in the dynamics. Such non-Markovian effects can be understood by analyzing the distribution of the polymer conformations at the instant of reaction, which is found to significantly depart from the equilibrium distribution.

In addition, we derived asymptotic expressions of the mean cyclization time with the polymer size $N$ and capture radius $b$, which are identical in both Markovian and non-Markovian approaches, but with different prefactors. We computed the precise values of the prefactors in the case of the Wilemski-Fixman approach. In particular, it is found that the scaling of the MFCT for large $N$ is given by $T\sim N^{3/2}\ln(N/b^{2})$, whereas for the Rouse chain one has $T\sim N^{2}$ (see Ref. Bénichou et al. (2015) for review). Hydrodynamic interactions therefore change both the dependence on $N$ and $b$. This difference comes from the fact that subdiffusive exponent that characterize the monomer dynamics at intermediate length scales are different in both models.

The present work demonstrates that the physics of cyclization kinetics in realistic models of polymers can be described by taking into account non-Markovian effects, which turn out to be much more important than the errors due to the approximate treatment of hydrodynamic interactions.

Acknowledgements.

N.L., T.G., O.B., and R.V. acknowledge the support of the Campus France (project No. 28252XE) and of the European Research Council starting Grant No. FPTOpt-277998. M.D. and A.B. acknowledge the support of the DAAD through the PROCOPE program (project No. 55853833) and of the DFG through Grant No. Bl 142/11-1 and through IRTG Soft Matter Science (GRK 1642/1).

Appendix A Pre-averaged Rotne-Prager tensor

We give in this section some details about the pre-averaging of the Rotne-Prage tensor, in order to quantify the differences obtained with the pre-averaged Oseen tensor. Let us first remind that the hydrodynamic interactions do not modify the equilibrium state of the polymer, which is that of a flexible Gaussian chain, hence the equilibrium probability density of the vector $\mathbf{r}_{ij}$ is

[TABLE]

where $\sigma_{ij}=l_{0}^{2}|i-j|$ is the variance of each spatial coordinate of $\mathbf{r}_{ij}$ at equilibrium. The pre-averaged Rotne-Prager tensor is defined as

[TABLE]

We note that the average over all orientations of the non-isotropic tensor appearing in the expression of $\mathbf{D}_{ij}$ is

[TABLE]

Using this relation, the integral over rotational degrees of freedom can be performed in Eq. (39); using Eqs. (3),(4) we obtain an integral over radial components only:

[TABLE]

with $H$ the Heaviside step function. The result of this integral in the limit of small $a$ is

[TABLE]

Note that for $a=0$ one recovers the pre-averaged Oseen tensor (6). The corrections are of order $a^{2}/(3l_{0}^{2}|i-j|)$; with our choice of parameter $a=l_{0}/4$ these corrections are always less than $\simeq 2\%$, which justifies to chose to take the averaged Oseen tensor instead of the averaged Rotne-Prager tensor in the analysis.

Appendix B Explicit expressions of the dynamic quantities $\alpha_{i},\beta_{i},\phi,\psi$

We consider the Fokker-Planck equation describing the chain dynamics with the pre-averaged mobility tensor. Since this tensor is isotropic, we consider the dynamics of a single spatial component, say $x_{i}$, of the vector positions $\mathbf{x}_{i}$. To the Langevin equation (2) we can associate the Fokker-Planck equation

[TABLE]

By replacing the mobility tensor by its pre-averaged version, Gaussian solutions do exist. The evolution of the mean vector $\mu=(\mu_{1},...,\mu_{N})$ and the covariance matrix $\Gamma_{ij}$ of $x_{i},x_{j}$ with time can be found by multiplying (43) by $x_{i}$ or $x_{m}x_{n}$, integrate over all $x_{k}$ and use the divergence theorem, resulting inVan Kampen (1992)

[TABLE]

where $K=DM$ is the product of the pre-averaged mobility matrix $D$ by the connectivity matrix $M$. To solve these equations, it is convenient to introduce modes, which diagonalize $K$ (which ones coincide with Rouse modes if $D_{ij}\propto\delta_{ij}$). We define an invertible matrix $P$ such that $PKP^{-1}=\text{Diag}(\nu_{1},...,\nu_{N})$ is diagonal, with $0=\nu_{1}<\nu_{2}<\nu_{3}<...$. The vanishing eigenvalue is associated to the translational motion of the chain, while all other eigenvalues describe the internal conformational degrees of freedom. We consider the amplitudes of the Zimm modes, defined by

[TABLE]

and the associated mean vector $\mu_{i}=\langle a_{i}\rangle$ and covariance matrix $\kappa_{ij}=\mathrm{Cov}(a_{i},a_{j})$ read

[TABLE]

Note that $P$ is not an orthogonal matrix. The evolution of $u$ and $\kappa$ reads

[TABLE]

The solutions are :

[TABLE]

At large times, the covariance matrix reaches the stationary value (for $i,j\geq 2$):

[TABLE]

We introduce the set of coefficients $\tilde{c}$,

[TABLE]

where we remind that $h=(-1,0,....,0,1)^{t}$. The vector $\tilde{c}$ is such that

[TABLE]

Note that $\tilde{c}_{1}=0$, meaning that the motion of the drift center is not involved in the evolution of the internal variable $x_{\mathrm{ee}}$. The equilibrium end-to-end distance is :

[TABLE]

Similarly, we obtain for $\alpha_{i}$

[TABLE]

Consider the correlation matrix $C_{ij}(t)=\langle a_{i}(t)a_{j}(0)\rangle$, starting from an equilibrium configuration. This matrix is a solution of

[TABLE]

The solution is straightforward: $C_{ij}=\kappa_{ij}^{s}e^{-\nu_{i}t}$. Given that $\langle x_{\mathrm{ee}}(t)x_{\mathrm{ee}}(0)\rangle=\sum_{ij}\tilde{c}_{i}\tilde{c}_{j}C_{ij}$, we obtain for the function $\phi$ [defined in Eq. (11)]

[TABLE]

Next, we denote $\kappa^{*}$ the covariance matrix of $a_{i},a_{j}$ at equilibrium with the constraint of fixed $x_{\mathrm{ee}}$, which from (67) reads

[TABLE]

Consider now

[TABLE]

which is related to $\beta_{i}$ by

[TABLE]

Taking the covariance matrix (59) as an initial condition for the dynamics (49), we obtain

[TABLE]

Using (58) and (61), we obtain for $\beta_{i}$

[TABLE]

Finally, the function $\psi=\beta_{N}-\beta_{1}$ reads

[TABLE]

where we have again used (58). Hence, the expressions of all the dynamical quantities $\alpha_{i},\beta_{i},\phi,\psi$ are given explicitly in this section.

Appendix C Derivation of the self-consistent equations (14)

We present here a derivation of the set of equations (14) for the moments $m_{i}^{\pi}$, which adapts the method used for Rouse polymersGuerin et al. (2013a). First, we multiply (9) by $x_{iz}\delta(\mathbf{x}_{N}-\mathbf{x}_{1}-R_{f}\hat{\mathbf{e}}_{z})$ (where $\hat{\mathbf{e}}_{z}$ is the unit vector in the vertical direction, and $R_{f}$ is fixed and satisfies $0<R_{f}<b$) and integrate over all configurations, to get :

[TABLE]

with $p$ representing the probability distribution function of the end-to-end vector, $p(R_{f}\hat{\mathbf{e}}_{z},t|\pi_{\Omega},0)$ is the probability that $\mathbf{r}_{\mathrm{ee}}=R_{f}\hat{\mathbf{e}}_{z}$ at $t$ with an initial distribution $\pi_{\Omega}$, $\mathbb{E}(x_{iz},t|R_{f}\hat{\mathbf{e}}_{z},t;\pi_{\Omega},0)$ is the conditional average of the $z$ coordinate of $\mathbf{x}_{i}$ at $t$ given that $\mathbf{r}_{\mathrm{ee}}=R_{f}\hat{\mathbf{e}}_{z}$ at the same time $t$ and starting from $\pi_{\Omega}$ initially; other notations are similar. The mention “$P_{\text{ini},\Omega}$" means the equilibrium distribution with $\mathbf{r}_{\mathrm{ee}}=R_{0}\mathbb{u}_{r}(\Omega)$ initially, we perform the average over initial end-to-end distances at the end of the calculation.

Let us remind the following formula for the conditional mean of a Gaussian variable $X$ given that a second Gaussian variable $Y$ takes the value $Y_{0}$:

[TABLE]

A similar general formula for the conditional covariances of Gaussian variables is

[TABLE]

Now, we consider a fixed angular direction $\Omega$ and define $\theta$ the angle with the vertical direction. We write

[TABLE]

with $x_{ir}$ the component of $\mathbf{x}_{i}$ in the direction $\hat{\mathbf{u}}(\Omega)$. We note that conditioning the end-to-end vector to have the value $R_{f}\hat{\mathbf{e}}_{z}$ imposes that its component in the direction $\hat{\mathbf{u}}(\Omega)$ takes the value $R_{f}\cos\theta$. Then, applying (66) and using the definitions of $R_{\pi},\beta_{i},\psi$, we get:

[TABLE]

Applying the same reasoning in the direction normal to $\mathbf{u}(\Omega)$, we get

[TABLE]

Inserting these expressions into (68) leads to

[TABLE]

Let us pose now

[TABLE]

where initial equilibrium conditions (apart from the constraint for $x_{\mathrm{ee}}(0)$) are understood. A reasoning similar to that leading to Eq. (71) gives

[TABLE]

Similarly, from (66) and (16), we obtain

[TABLE]

Finally, the propagators for the end-to-end distance read :

[TABLE]

All the terms appearing in (65) have been evaluated, the self-consistent equation becomes

[TABLE]

This equation should be verified for any $R_{f}$ between [math] and $b$. We choose to write it in the limit $R_{f}\rightarrow 0$ (more precisely, we develop both expressions at first order in $R_{f}$). Noting that $\hat{\mathbf{u}}_{r}\cdot\hat{\mathbf{e}}_{z}=\cos\theta$, the integration over $\theta$ can be performed, leading to

[TABLE]

We can eliminate $T$ by replacing its expression :

[TABLE]

Let us now find the dependence of $R$ and $\mu_{i}^{\text{stat},R_{0}}$ with $R_{0}$. From (72) and (67), we deduce

[TABLE]

Furthermore, applying the formula for the conditional covariances in (15), we get

[TABLE]

From which we obtain

[TABLE]

and

[TABLE]

The last step consists in adapting the reasoning to the case of initial end-to-end distance that is distributed $p_{r}(R_{0})=R_{0}^{2}e^{-R_{0}^{2}/2L^{2}}/\int_{b}^{\infty}r^{2}e^{-r^{2}/2L^{2}}dr$. This can be achieved by keeping an average over $R_{0}$ at all steps of the derivation. Averaging (80) with respect to $R_{0}$ [by using (84,85)] leads to the final form of the self-consistent equation (14) for the moments given in the main text.

Appendix D Numerical methods

We performed the numerical integration of the stochastic equation (2) by using either the full Rotne-Prager tensor or its pre-averaged form. At each time step of size $\Delta t$, the positions of the monomers evolve according to the algorithm of Ermak and McCammonErmak and McCammon (1978):

[TABLE]

with the same notations as before. $\bm{\xi}_{i}(\Delta t)$ is a random Gaussian noise with zero mean and covariance $\left\langle\bm{\xi}_{i}(\Delta t)\bm{\xi}_{j}(\Delta t)\right\rangle=2\mathbf{D}_{ij}(t)\Delta t$. The generation of the $3N$ random numbers $\bm{\xi}_{i}$ requires to find a Cholesky decomposition of the mobility tensor, which can be done in $N^{3}$ operations. In the case of simulations using the pre-averaged mobility tensor (6), this decomposition needs to be performed only once, whereas simulations using the full Rotne-Prager tensor require to perform a Cholesky decomposition at each time step, resulting in considerably longer computational times which prevented us to explore the same range of parameters. Initial configurations are generated from the equilibrium Gibbs-Boltzmann Gaussian distribution $P_{\mathrm{eq}}(\{\mathbf{x}\})\propto e^{-k\sum_{i=1}^{N-1}(\mathbf{x}_{i+1}-\mathbf{x}_{i})^{2}/2k_{B}T}$, and are rejected if the condition $|\mathbf{x}_{N}-\mathbf{x}_{1}|>b$ is not satisfied. Once an equilibrium configuration is generated, the positions evolve through (86) until reaching a configuration $|\mathbf{x}_{N}-\mathbf{x}_{1}|\leq b$, the cyclization time $t_{i}$ for this run is recorded. The MFCT is finally found by ensemble averaging $t_{i}$ over many runs. The time step is chosen as suggested by Pastor et al. Pastor et al. (1996). Noting $R$ the end-to-end distance :

[TABLE]

if $R^{2}<b^{2}+3l_{0}$, and

[TABLE]

if not.

We choosed $\Delta_{\text{high}}=\sqrt{N}10^{-5}$ and $\Delta_{\text{low}}=10^{-5}$, and we controlled that convergence was reached. The results of simulations are represented on Fig. 1.