Wet and dry internal friction can be measured with the Jarzynski equality

TL;DR

This paper proposes a new protocol using the Jarzynski equality to distinguish and measure wet and dry internal friction in polymers through repeated stretching experiments, validated by analytical and simulation models.

Contribution

It introduces a simple, experimentally feasible method to quantify wet and dry internal friction coefficients in polymers using work measurements and the Jarzynski equality.

Findings

The protocol accurately estimates internal friction coefficients.

Analytical calculations confirm the method's validity.

Simulations demonstrate applicability to complex polymer models.

Abstract

The existence of two types of internal friction wet and dry is revisited, and a simple protocol is proposed for distinguishing between the two types and extracting the appropriate internal friction coefficient. The scheme requires repeatedly stretching a polymer molecule, and measuring the average work dissipated in the process by applying the Jarzynski equality. The internal friction coefficient is then estimated from the average dissipated work in the extrapolated limit of zero solvent viscosity. The validity of the protocol is established through analytical calculations on a one-dimensional free-draining Hookean spring-dashpot model for a polymer, and Brownian dynamics simulations of: (a) a single-mode nonlinear spring-dashpot model for a polymer, and (b) a finitely extensible bead-spring chain with cohesive intra-chain interactions, both of which incorporate fluctuating hydrodynamic interactions. Well-established single-molecule manipulation techniques, such as optical tweezer-based pulling, can be used to implement the suggested protocol experimentally.

Tables

Table 1. Table 1: Typically observed lower and upper bounds on optical tweezer parameters.

Parameter	Lower bound	Upper bound
Trap stiffness (pN/nm)	$0.0002$ [ref. 51]	$0.9$ [ref. 35]
Pulling velocity (nm/s)	$10$ [ref. 35]	$13560$ [ref. 56]
Stretching distance (nm)	$10$ [ref. 35]	$8000$ [ref. 13]

Table 2. Table 2: Internal friction coefficients estimated using the protocol described in Fig.  8 , for various values of the molecular, and control parameters. The error associated with the protocol is calculated as, % error = 100 × | ( K BD − K ) / K | \%\,\textrm{error}=100\times\left|\left(K_{\textrm{BD}}-K\right)/K\right| . An ensemble size of N = 1 × 10 4 𝑁 1 superscript 10 4 N=1\times 10^{4} and a time-step width of Δ ​ t ∗ = 1 × 10 − 3 Δ superscript 𝑡 1 superscript 10 3 \Delta t^{*}=1\times 10^{-3} was used to obtain the results. For the K = 100.0 × 10 − 9 𝐾 100.0 superscript 10 9 K=100.0\times 10^{-9} kg/s case (Marko-Siggia), a smaller time-step, Δ ​ t ∗ = 1 × 10 − 4 Δ superscript 𝑡 1 superscript 10 4 \Delta t^{*}=1\times 10^{-4} , was used. For the K = 4.0 × 10 − 9 𝐾 4.0 superscript 10 9 K=4.0\times 10^{-9} kg/s case (FENE), a larger ensemble, N = 1 × 10 5 𝑁 1 superscript 10 5 N=1\times 10^{5} , was used.

	Spring force law:	Marko-Siggia
	Input,	$h^{\ast }=0.0$			$h^{\ast }>0.0$
Parameters	$K[\times {10}^9$ kg/s]	$K_{\text{BD}}[\times {10}^9$ kg/s]	$\%$ error	$h^{\ast }$	$K_{\text{BD}}[\times {10}^9$ kg/s]	$\%$ error
$b=200,l_H=150\mathrm{nm},\mathrm{\hspace{0.17em}2}{l}_H\to 3l_H$	$1.0$	$1.04\pm 0.01$	$4.42$	$0.3$	$0.97\pm 0.02$	$3.11$
$b=400,l_H=350\mathrm{nm},\mathrm{\hspace{0.17em}4}{l}_H\to 6l_H$	$1.0$	$0.99\pm 0.03$	$1.36$	$0.16$	$1.04\pm 0.02$	$4.29$
$b=200,l_H=150\mathrm{nm},\mathrm{\hspace{0.17em}2}{l}_H\to 3l_H$	$10.0$	$10.04\pm 0.05$	$0.47$	$0.3$	$9.87\pm 0.09$	$1.24$
$b=400,l_H=350\mathrm{nm},\mathrm{\hspace{0.17em}4}{l}_H\to 6l_H$	$10.0$	$9.922\pm 0.008$	$0.78$	$0.16$	$9.95\pm 0.09$	$0.52$
$b=800,l_H=500\mathrm{nm},\mathrm{\hspace{0.17em}5}{l}_H\to 7l_H$	$3.0$	$2.95\pm 0.02$	$1.74$	$-$	$-$	$-$
$b=800,l_H=500\mathrm{nm},\mathrm{\hspace{0.17em}5}{l}_H\to 7l_H$	$6.0$	$6.00\pm 0.07$	$0.11$	$-$	$-$	$-$
$b=800,l_H=500\mathrm{nm},\mathrm{\hspace{0.17em}5}{l}_H\to 7l_H$	$100.0$	$97.4\pm 1.4$	$2.56$	$-$	$-$	$-$
	Spring force law:	FENE
	Input,	$h^{\ast }=0.0$			$h^{\ast }>0.0$
Parameters	$K[\times {10}^9$ kg/s]	$K_{\text{BD}}[\times {10}^9$ kg/s]	$\%$ error	$h^{\ast }$	$K_{\text{BD}}[\times {10}^9$ kg/s]	$\%$ error
$b=400,l_H=350\mathrm{nm},\mathrm{\hspace{0.17em}2}{l}_H\to 3l_H$	$2.0$	$-$	$-$	$0.23$	$1.993\pm 0.009$	$0.35$
$b=600,l_H=250\mathrm{nm},\mathrm{\hspace{0.17em}3}{l}_H\to 4l_H$	$3.0$	$-$	$-$	$0.09$	$2.985\pm 0.009$	$0.50$
$b=800,l_H=500\mathrm{nm},\mathrm{\hspace{0.17em}3}{l}_H\to 5l_H$	$4.0$	$3.992\pm 0.007$	$0.18$	$-$	$-$	$-$
$b=800,l_H=500\mathrm{nm},\mathrm{\hspace{0.17em}3}{l}_H\to 5l_H$	$6.0$	$6.02\pm 0.02$	$0.23$	$-$	$-$	$-$

Table 3. Table 3: Parameter values for the two representative cases for which the free energy differences are evaluated using Jarzynski’s equality and numerical integration.

	Parameter sets
	1	2
$b$	$50$	$80$
$c_1$	$20$	$15$
$c_2$	$1$	$15$
$𝝌_{}^{\ast }{}_1{}^{}$	$(0,0,0)$	$(0,0,0)$
${𝝌}_2^{\text{(i)}\ast }$	$(1,0,0)$	$(4,0,0)$
${𝝌}_2^{\text{(f)}\ast }$	$(3,0,0)$	$(5,0,0)$

Table 4. Table 4: A comparison of the free-energy differences calculated using numerical integration [Eq. ( 104 )], and BD simulations using Jarzynski’s equality [Eq. ( 1 )], over N = 1 × 10 5 𝑁 1 superscript 10 5 N=1\times 10^{5} trajectories. Simulation data reported for freely-draining dumbbells with no internal friction [ h ∗ = 0.0 , φ = 0.0 ] delimited-[] formulae-sequence superscript ℎ 0.0 𝜑 0.0 [h^{*}=0.0,\,\varphi=0.0] . The error is quantified as, % error = 100 × | ( Δ A ∗ − Δ A num ∗ ) / Δ A num ∗ | \%\,\text{error}=100\times\left|\left(\Delta A^{*}-\Delta A^{*}_{\text{num}}\right)/\Delta A^{*}_{\text{num}}\right| .

	Parameter set 1 :	$\Delta A_{\text{num}}^{\ast }=2.11504$
$v^{\ast }$	$\Delta A^{\ast }$	$\%$ error
$0.001$	$2.1151\pm 0.0002$	$0.0006$	$0.0017\pm 0.0003$
$0.005$	$2.1149\pm 0.0004$	$0.008$	$0.0084\pm 0.0006$
$0.01$	$2.1147\pm 0.0006$	$0.02$	$0.0169\pm 0.008$
$0.02$	$2.1144\pm 0.0008$	$0.03$	$0.033\pm 0.001$
$0.05$	$2.116\pm 0.001$	$0.04$	$0.081\pm 0.002$
$0.1$	$2.114\pm 0.002$	$0.06$	$0.155\pm 0.003$
$0.2$	$2.116\pm 0.003$	$0.05$	$0.281\pm 0.004$
$0.5$	$2.116\pm 0.004$	$0.04$	$0.508\pm 0.005$
$1.0$	$2.115\pm 0.005$	$0.01$	$0.673\pm 0.006$
	Parameter set 2 :	$\Delta A_{\text{num}}^{\ast }=5.55479$
$v^{\ast }$	$\Delta A^{\ast }$	$\%$ error
$0.001$	$5.5551\pm 0.0002$	$0.005$	$0.0031\pm 0.0004$
$0.005$	$5.5559\pm 0.0006$	$0.02$	$0.0156\pm 0.0008$
$0.01$	$5.5546\pm 0.0008$	$0.003$	$0.031\pm 0.001$
$0.02$	$5.553\pm 0.001$	$0.02$	$0.063\pm 0.002$
$0.05$	$5.554\pm 0.002$	$0.01$	$0.155\pm 0.003$
$0.1$	$5.551\pm 0.003$	$0.07$	$0.306\pm 0.004$
$0.2$	$5.549\pm 0.005$	$0.08$	$0.593\pm 0.006$
$0.5$	$5.54\pm 0.01$	$0.26$	$1.38\pm 0.01$
$1.0$	$5.50\pm 0.03$	$0.95$	$2.42\pm 0.03$

Equations

$$\begin{array}{rl}⟨\exp \left[-W/k_{B}T\right]⟩& =\exp \left[-\Delta A/k_{B}T\right];\\ ⟨W_{\text{dis}}⟩& =⟨W⟩-\Delta A\end{array}$$

$$\mathcal{H}$$

$$\frac{d{r}_2^{\ast }}{d{t}^{\ast }}=-\frac{E{r}_2^{\ast }}{\theta }+\frac{c_2{\chi }^{\ast }(t^{\ast })}{4\theta }+\frac{1}{\theta }\xi (t^{\ast })$$

$$r_2^{\ast }(t^{\ast })=r_2^{\ast }(0)G(t^{\ast })+\frac{1}{\theta }{\int }_0^{t^{\ast }}d{t}_1^{\ast }$$

$$\displaystyle\times\Biggl{(}\dfrac{c_{2}\chi^{*}(t_{1}^{*})}{4}$$

$$W$$

$$={\int }_0^{\tau }{c}_{2}H\left(\chi (t)-r_2\right)\dot{\chi }(t)dt$$

$$=k_{B}T\left[c_2{\int }_0^{{\tau }^{\ast }}\left({\chi }^{\ast }(t^{\ast })-r_2^{\ast }\right){\dot{\chi }}^{\ast }(t^{\ast })d{t}^{\ast }\right]$$

$$W^{\ast }=\frac{c_2}{2}\left[{\chi }^{2\ast }\left({\tau }^{\ast }\right)-{\chi }^{2\ast }\left(0\right)\right]-c_2{\int }_0^{{\tau }^{\ast }}d{t}^{\ast }{\dot{\chi }}^{\ast }(t^{\ast })r_2^{\ast }$$

$$A^{\ast }\left({\chi }^{\ast }\right)=\left[\frac{c_2}{2(c_2+1)}\right]{\chi }^{2\ast },$$

$$⟨W^{\ast }⟩=\Delta A^{\ast }+\frac{c_2^2}{c_2+1}{\int }_0^{{\tau }^{\ast }}d{t}^{\ast }$$

$$\displaystyle\times G\left(t^{*}-t_{1}^{*}\right)\dot{\chi}^{*}(t_{1}^{*})\Bigr{]}$$

$${\sigma }^2=\frac{2c_2^2}{\left(c_2+1\right)}{\int }_0^{{\tau }^{\ast }}d{t}^{\ast }{\int }_0^{t^{\ast }}d{t}_1^{\ast }{\dot{\chi }}^{\ast }(t^{\ast })G(t^{\ast }-t_1^{\ast }){\dot{\chi }}^{\ast }(t_1^{\ast })$$

$$⟨W^{\ast }⟩=\Delta A^{\ast }+\frac{{\sigma }^2}{2}$$

$$⟨W_{\text{dis}}^{\ast }⟩\equiv ⟨W^{\ast }⟩-\Delta A^{\ast }=\frac{{\sigma }^2}{2}$$

$$P^{\ast }(W^{\ast })=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{{\left(W^{\ast }-⟨W^{\ast }⟩\right)}^2}{2{\sigma }^2}\right]$$

$$⟨\exp \left(-W^{\ast }\right)⟩$$

$$=\exp \left(-\Delta A^{\ast }\right)$$

$$\Delta A^{\ast }=\frac{c_2}{2\left(c_2+1\right)}\left[{\left({\chi }^{\text{(f)}\ast }\right)}^2-{\left({\chi }^{\text{(i)}\ast }\right)}^2\right],$$

$$⟨W_{\text{dis}}^{\ast }⟩$$

$$+16\left(\frac{c_2^2}{{\left(c_2+1\right)}^3}\right){\theta }^2{v}^{\ast 2}\left[\exp \left(-\frac{(c_2+1)d^{\ast }}{4v^{\ast }\theta }\right)-1\right]$$

$${⟨W⟩}_{\text{dis}}$$

$$+\frac{c_2^2}{{\left(c_2+1\right)}^3}$$

$${⟨W⟩}_{\text{dis},{\eta }_{\text{s}}\to 0}$$

$$+\frac{c_2^2}{{\left(c_2+1\right)}^3}$$

$$\frac{{⟨W⟩}_{\text{dis},{\eta }_{\text{s}}\to 0}}{d}=Kv,$$

$$\mathcal{H}$$

$$W={\int }_0^{\tau }\left(\frac{\partial \mathcal{H}}{\partial {\mathbf{\chi }}_2}\right)\cdot \mathbf{v}dt$$

$$W^{\ast }=c_2\sum _{j=1}^{N_{\text{t}}}{\left({\chi }_{2x}^{\ast }-r_{2x}^{\ast }\right)}_j{v}_x^{\ast }\Delta t_j^{\ast }$$

$$H=\frac{k_{B}T}{l_{\text{H}}^2}=\frac{4.142\text{pN nm}}{{\left(500\right)}^2{\left(\text{nm}\right)}^2}=1.657\times 10^{-5}\text{pN/nm}$$

$$\displaystyle\scriptsize\dfrac{U_{\text{SDK}}(r)}{k_{B}T}=\left\{\begin{array}[]{ll}4\left[\left(\dfrac{\sigma_{\text{{s}}}}{r}\right)^{12}-\left(\dfrac{\sigma_{\text{s}}}{r}\right)^{6}+\dfrac{1}{4}\right]-\epsilon;&r\leq 2^{1/6}\sigma_{\text{s}}\\[15.0pt] \dfrac{\epsilon}{2}\left[\cos\left(m_{1}\left(\dfrac{r}{\sigma_{\text{s}}}\right)^{2}+m_{2}\right)-1\right];&2^{1/6}\sigma_{\text{s}}\leq r\leq r_{\text{c}}\\[15.0pt] 0;&r\geq r_{\text{c}}\\ \end{array}\right.$$

$$\bm{{r}^{*}}_{\mu}\left(t^{*}+\Delta t^{*}\right)=\bm{{r}^{*}}_{\mu}\left(t^{*}\right)+\dfrac{1}{4}\sum_{\nu=1}^{N_{\text{b}}}\Biggl{[}\bm{\Upsilon^{*}}_{\mu\nu}\cdot\bm{{F}}_{\nu}^{\bm{*}(\phi)}-\boldsymbol{\Upsilon}^{\bm{*}}_{\!{}_{1\nu}}\cdot\left[c_{1}\left(\bm{r^{*}}_{1}-\bm{\chi}^{*}_{1}\right)\right]-\boldsymbol{\Upsilon}^{\bm{*}}_{\!{}_{N_{\text{b}}}\nu}\cdot\left[c_{2}\left(\bm{r^{*}}_{N_{\text{b}}}-\bm{\chi}^{*}_{2}\right)\right]\Biggr{]}\Delta t^{*}+\dfrac{1}{\sqrt{2}}\sum_{\nu=1}^{N_{\text{b}}}\bm{\mathscr{B}^{*}}_{\mu\nu}\cdot\Delta\bm{\mathscr{W}}^{\bm{*}}_{\nu}$$

$$⟨W_{\text{dis}}⟩$$

$$\frac{⟨W_{\text{dis}}⟩}{{⟨W_{\text{dis}}⟩}_{\theta }}-1=\frac{{\gamma }_{\text{G}}-{\gamma }_{\text{s}}}{{\gamma }_{\text{s}}}$$

Full text

I Introduction

Conformational transitions in polymer molecules are impeded by solvent molecules, and sometimes additionally by intramolecular interactions. The dissipation caused by the latter are termed as internal friction Kuhn and Kuhn (1945); de Gennes (1979); Manke and Williams (1985); Sagnella et al. (2000); Hagen (2010); Hridya and Mukherjee (2018), and examples of such dissipation include the damping of protein folding Ansari et al. (1992); Qiu and Hagen (2004); Cellmer et al. (2008); Wensley et al. (2010); Soranno et al. (2017), the modulation of stretching transitions in polysaccharides Khatri et al. (2007), and the enhancement of dissipated work in the stretch-relaxation of polymers Murayama et al. (2007); Alexander-Katz et al. (2009); Schulz et al. (2015). While the microscopic origin of internal friction is manyfold Khatri et al. (2007); Murayama et al. (2007); Alexander-Katz et al. (2009); Schulz et al. (2015); Echeverria et al. (2014); de Sancho et al. (2014); Sashi et al. (2016); Jain et al. (2016); Ameseder et al. (2018); Jas et al. (2001); Soranno et al. (2017), it has been broadly classified as being either of the wet or dry kind Soranno et al. (2012); Hagen (2010).

The transition of a protein from an unfolded to its native folded state is commonly interpreted as a diffusive search process over a rugged energy landscape Onuchic et al. (1997), and the internal friction associated with landscape roughness is typically considered to be of the wet type Wensley et al. (2010); de Sancho et al. (2014); Neupane et al. (2017). An analytical solution for the diffusion coefficient in one dimension was derived by Zwanzig Zwanzig (1988) who observed that the effective friction $\gamma_{\text{eff}}$, is related purely multiplicatively to the solvent friction $\gamma_{\text{s}}$, as $\gamma_{\text{eff}}=\gamma_{\text{s}}\exp\left[\left(\delta/k_{B}T\right)^{2}\right]$, where $k_{B}$ is Boltzmann’s constant, $T$ is the absolute temperature, and $\delta^{2}$ is the variance of the heights of the normally distributed undulations. Since $\gamma_{\text{s}}\sim\eta_{\text{s}}$, where $\eta_{\text{s}}$ denotes the solvent viscosity, the effective friction is higher than that expected due to solvent friction alone, at any finite value of $\eta_{\text{s}}$. It is also clear that the internal friction would vanish in the extrapolated limit of zero solvent viscosity ($\eta_{\textrm{s}}\to 0$), which is a characteristic feature of wet internal friction.

On the other hand, experimental measurements of the dependence of the reconfiguration time of small proteins on $\eta_{\text{s}}$ find a finite value in the extrapolated limit of $\eta_{\text{s}}\to 0$, indicating the presence of a solvent-viscosity-independent resistance to folding Qiu and Hagen (2004); Cellmer et al. (2008); Soranno et al. (2012). Hagen Hagen (2010) has proposed a heterogeneous reaction friction model with the unfolded and native states separated by two consecutive barriers, one for each mode of friction, as an explanation for the presence of dry internal friction.

The commonly accepted operational definition for dry internal friction Qiu and Hagen (2004); Avdoshenko et al. (2017); Soranno et al. (2017) as the reconfiguration time in the limit $\eta_{\text{s}}\to 0$ is not a direct quantitative measure of the internal friction coefficient, nor does it apply to the case of wet internal friction. In principle, the effective friction could be calculated from experimental measurements of the diffusion coefficient of biomolecules hopping between native and unfolded states Chung and Eaton (2013); Chung et al. (2015); Neupane et al. (2016, 2018). However, such measurements would not determine if the internal friction were of the wet or the dry type unless they were performed at multiple solvent viscosities, followed by extrapolation to the limit $\eta_{\text{s}}\to 0$. There is clearly a need for a protocol that can directly estimate the internal friction coefficient and distinguish between the two types.

The internal friction coefficient can be obtained directly from dissipated work, as demonstrated by the experiments on condensed DNA by Murayama et al. (2007), and simulations of polypeptide-stretching by Schulz et al. (2015). Indeed in the protocol proposed by Netz and coworkers Alexander-Katz et al. (2009); Schulz et al. (2015), the work required to stretch a macromolecule is separated into two parts: reversible free energy increase due to the extension of the molecule, and irreversible work required to overcome rate-dependent restoring forces arising from solvent and internal friction. Within this framework, they show that the average dissipated work in the limit of $\eta_{\text{s}}\to 0$ scales with the number of hydrogen bonds, which are considered to be responsible for internal friction Schulz et al. (2015).

Here we propose a novel application of the Jarzynski equality (JE) Jarzynski (1997, 2007) and show that by focussing on measuring the dissipation associated with stretching a macromolecule rather than on obtaining the free-energy difference, a quantitative measure of the internal friction can be obtained. The JE has been routinely employed for reconstructing the free energy landscape of biomolecules from experiments Liphardt et al. (2002); Harris et al. (2007); Gupta et al. (2011) and simulations Hummer and Szabo (2010); Hodges et al. (2016), while dissipation has largely been ignored (except for estimating the accuracy of the free-energy difference Ritort et al. (2002); Jarzynski (2006); Yunger Halpern and Jarzynski (2016)). In the proposed protocol, multiple realizations of the pulling experiment are performed and the JE is used to extract both the free-energy difference and the average dissipated work at finite pulling rates,

[TABLE]

where the $\left<...\right>$ in Eq. (1) represents an average with respect to the probability distribution of work values. Prior studies Murayama et al. (2007); Alexander-Katz et al. (2009); Schulz et al. (2015) estimate $\Delta A$ from the work done in the quasi-static limit Callen (1985) and calculate $\left<W_{\text{dis}}\right>$ at finite pulling rates by subtracting $\Delta A$ from the total work done, rather than estimating both components of work simultaneously, as is done here.

In essence, the proposed protocol consists of calculating $\left<W_{\text{dis}}\right>$ at fixed values of both the pulling velocity $v$ and distance $d$ over which the molecule is stretched, but at various values of $\eta_{\text{s}}$. The value in the limit $\eta_{\text{s}}\to 0$, $\left<W_{\text{dis}}\right>_{\eta_{\!{}_{\,\text{s}}}\to\,0}$, is then obtained by extrapolation. By repeating this process for a number of values of $v$ and plotting the ratio $\left<W_{\text{dis}}\right>_{\eta_{\!{}_{\,\text{s}}}\to\,0}/d$ as a function of $v$, the internal friction coefficient can be determined from the slope of the linear region at sufficiently small velocities. Clearly, dry internal friction corresponds to cases where $\left<W_{\text{dis}}\right>_{\eta_{\!{}_{\,\text{s}}}\to\,0}$ is non-zero, while wet friction is indicated when it is zero. In the latter case, the protocol measures the enhancement in friction at any finite value of $\eta_{\text{s}}$.

The validity of the proposed protocol is established for both types of internal friction using coarse-grained polymer models. Additionally, since hydrodynamic interactions (HI) are known to affect the dynamic response of polymers Kailasham et al. (2018); Prakash (2019), the effect of HI on dissipated work is also examined.

For dry internal friction, a spring-dashpot model Kuhn and Kuhn (1945); Manke and Williams (1985); Khatri and McLeish (2007); Samanta and Chakrabarti (2016) is considered where the molecule is represented as massless beads connected by a spring and dashpot in parallel with each other. The spring accounts for entropic elasticity, while dissipative effects due to internal friction are captured by the dashpot Bird et al. (1987). The drag on the beads is responsible for solvent friction. By its very construction, this model describes dry internal friction, as the dashpot contributes to dissipation even in the limit of $\eta_{\text{s}}\to 0$. Within this framework, two examples are considered. In the first case, the work distribution for a free-draining Hookean spring-dashpot model subjected to constant-velocity pulling is analytically calculated. In the second case, pulling simulations on a nonlinear-spring-dashpot model with fluctuating HI are performed using Brownian dynamics (BD). In both these cases, it is demonstrated that the internal friction coefficient estimated from $\left<W_{\text{dis}}\right>$ in the limit $\eta_{\text{s}}\to 0$ is identical to the damping coefficient of the dashpot, which is a model input parameter, thereby establishing the validity of the proposed protocol. It is also shown that HI does not affect $\left<W_{\text{dis}}\right>$.

For wet internal friction, a bead-spring chain with cohesive interactions between the beads is considered. A similar model was used by Netz and coworkers Alexander-Katz et al. (2009) to compare simulated values of internal friction with experimental data on force-induced unraveling of collapsed DNA Murayama et al. (2007). By using Zwanzig’s interpretation Zwanzig (1988) to estimate energy landscape roughness due to cohesive interactions, they implicitly assume wet internal friction. Using our protocol, it is established directly that the internal friction due to cohesive interactions in this coarse-grained polymer model is wet in nature. Further, it is observed that while HI reduces the total resistance to pulling, the enhancement in the friction coefficient remains unaffected.

The rest of the paper is organized as follows. In Sec. II, analytical calculations and Brownian dynamics simulations on a spring-dashpot model are presented. Sec. III covers the application of the protocol to pulling simulations on a single polymer chain with cohesive interactions between the beads. A discussion of the results and concluding remarks are provided in Sec. IV. Appendix A contains the derivation of the governing equation of motion for a one-dimensional Hookean spring-dashpot tethered at one end and pulled at the other. In Appendix B, additional details pertaining to the derivation of the moments of the work distribution in the Hookean spring-dashpot model is presented. In Appendix C, the derivation of the Fokker-Planck and stochastic differential equations for a spring-dashpot with a nonlinear force-extension relation and fluctuating hydrodynamic interactions is presented, along with the details of the solver algorithm. In Appendix D, the validation of the code for the dumbbell case is presented, by comparing the free energy difference obtained from the Jarzynski equality and that obtained from direct numerical integration. In Appendix E, the probability of the work distribution determined numerically for a Rouse chain tethered at one end and subjected to constant velocity pulling is compared against the analytical results derived by Dhar Dhar (2005), thereby validating the code for the single chain case.

II Dry internal friction

II.1 One-dimensional free draining Hookean spring-dashpot

The simple analytically tractable dumbbell model is shown in Fig. 1, with one bead, at $r_{1}$, fixed at the origin ($r_{1}=0$), and the other bead, at $r_{2}$, connected to a bead at $\chi$ which is indicative of the cantilever of an atomic force microscope (AFM), or the location of the optical trap. The dumbbell is suspended in an incompressible, Newtonian solvent of viscosity $\eta_{\text{s}}$. The bead radius is taken to be $a$, and its associated friction co-efficient given by $\zeta=6\pi\eta_{\text{s}}a$. The bead at $\chi$ is manipulated using a pre-determined protocol, given by $\chi=\chi(t)$. All the springs considered in the present model are Hookean: the spring in parallel with the dashpot has a spring constant of $H$, whereas the spring connecting the spring-dashpot setup to the driven bead has a spring constant of $c_{2}H$, where $c_{2}$ is an arbitrary positive constant. The damping coefficient of the dashpot is denoted by $K$. It is evident that the only degree-of-freedom in the system is $r_{2}$, which is allowed to execute stochastic motion. The Hamiltonian of the system is then written as

[TABLE]

For ease of algebra, it is convenient to work with non-dimensional variables until the need for dimensional variables arises. Using $l_{\text{H}}\equiv\sqrt{k_{B}T/H}$ and $\lambda_{\text{H}}\equiv\zeta/4H$ as the length- and time-scales, respectively, and denoting dimensionless quantities using an asterisk as superscript, the governing stochastic differential equation can be derived to be

[TABLE]

where $E=\left[\left(c_{2}+1\right)/4\right]$, $\theta=\left[1+\left(K/\zeta\right)\right]$, and the noise term, $\xi(t^{*})$, obeys $\left<\xi(t^{*})\right>=0$ and $\left<\xi(t^{*})\xi(t_{1}^{*})\right>=\left({\theta}/{2}\right)\delta(t^{*}-t^{*}_{1})$. The steps for the derivation of Eq. (3) have been detailed in Appendix A. The solution to Eq. (3) can be written as

[TABLE]

where $G(t^{*})=e^{-Et^{*}/\theta}$.

The work done during one realization of the pulling performed in the interval $\left[0,\tau\right]$ is Hodges et al. (2016); Chaki and Chakrabarti (2018)

[TABLE]

and the dimensionless work, $W^{*}=W/k_{B}T$, is then

[TABLE]

Upon substituting the expression for $r^{*}_{2}$ from Eq. (II.1) into Eq. (6), one obtains the complete expression for $W^{*}$. It is clear that the distribution of $W^{*}$ ought also be Gaussian, since $W^{*}$ is linear in $r^{*}_{2}(0)$ and $\xi(t^{*})$, both of which are Gaussian variables. It therefore suffices to evaluate the mean and variance of $W^{*}$ in order to completely determine the distribution. There is only one bead that is allowed to move freely in this problem, and as a result, the dimensionless free-energy, $A^{*}=-\ln Z^{*}$, can be obtained once the dimensionless partition function, $Z^{*}$ is known. Details of all the intermediate steps in the calculation of $\left<W^{*}\right>$, $\left<(W^{*}-\left<W^{*}\right>)^{2}\right>$ and $A^{*}$ are given in Appendix B. For ease of exposition, only the salient results are reproduced below.

The dimensionless free-energy can be shown to be given by the expression,

[TABLE]

while the expression for the average work is,

[TABLE]

where $\Delta A^{*}\equiv A^{*}\left[\chi^{*}(\tau^{*})\right]-A^{*}\left[\chi^{*}(0)\right]$. The variance of the work distribution, $\sigma^{2}\equiv\left<(W^{*}-\left<W^{*}\right>)^{2}\right>$, can be shown to be

[TABLE]

From Eqs. (8) and (9), it is readily seen that

[TABLE]

and the average dissipated work is given by

[TABLE]

It follows that the probability distribution of work is given by

[TABLE]

The quantities $\left<W^{*}\right>,\,\Delta A^{*},\,\text{and}\left<W^{*}_{\text{dis}}\right>$ have been calculated analytically without explicit recourse to the Jarzynski equality. This is a consequence of the governing equation being linear in the position variable and the noise term, resulting in a Gaussian distribution of the work trajectories. The Jarzynski’s equality is satisfied trivially for such systems, since

[TABLE]

which has also been reported previously Jarzynski (1997); Speck and Seifert (2004).

The development so far is applicable to any arbitrary pulling protocol, $\chi^{*}(t^{*})$. The particular value of the average dissipation is dependent on the protocol used to transition the system between its initial and final states. We focus attention on the special case of a constant-velocity pulling protocol commonly encountered in single-molecule force spectroscopy Harris et al. (2007); Gupta et al. (2011). Within this framework, $\chi^{*}(t^{*})=\chi^{\text{(i)}*}+\left(d^{*}t^{*}/\tau^{*}\right)$, implying that the last bead is moved across a distance $d^{*}$ over a time $\tau^{*}$, with a dimensionless pulling velocity given by $\dot{\chi}^{*}(t^{*})=v^{*}=d^{*}/\tau^{*}$. $\chi^{\text{(i)}*}$ and $\chi^{\text{(f)}*}$ represent the position of the last bead at $t^{*}=0$ and $t^{*}=\tau^{*}$ respectively. Under this protocol, the free energy difference, $\Delta A^{*}$, is given by

[TABLE]

and the average dissipated work can be evaluated to be

[TABLE]

It is now appropriate to obtain the expression for the average dissipated work in the limit of zero solvent friction. Since the solvent friction is absorbed into the definition of the timescale, it is first necessary to convert all quantities in Eq. (II.1) to their dimensional form, before taking the limit. Using the non-dimensionalization scheme discussed previously, we obtain

[TABLE]

In the extrapolated limit of zero solvent friction, $\zeta\to 0$ by definition of the bead-friction coefficient. Upon taking this limit in Eq. (II.1),

[TABLE]

In the limit of high pulling trap stiffness ($c_{2}\gg 1$), the second term on the RHS of Eq. (II.1) vanishes, and the parenthesized prefactor in the first term asymptotically tends to unity, leading to

[TABLE]

Clearly, the proposed protocol for determining the internal friction coefficient based on the Jarzynski equality recovers the damping coefficient $K$, establishing its validity in the case of the simple analytical model considered here.

II.2 Non-linear spring-dashpot with hydrodynamic interactions

II.2.1 Model description

In the more general case, a dumbbell model with fluctuating internal friction and hydrodynamic interactions is considered, as shown in Fig. 2 (a). The beads, each of radius $a$, are joined by a spring, with maximum stretchability $Q_{0}$ and a Hookean spring constant $H$, in parallel with a dashpot of damping coefficient $K$. The Marko-Siggia force expression Marko and Siggia (1995), widely employed to model the force-extension relationship in synthetic polymer molecules Black et al. (2017), as well as biopolymers Latinwo et al. (2014); Raman et al. (2014); Sunthar and Prakash (2005); Sasmal et al. (2017), is used to describe the entropic elasticity in the dumbbell. The positions of the two beads are $\bm{r}_{1}$ and $\bm{r}_{2}$, the connector vector joining the two beads is denoted by $\bm{Q}\equiv\bm{r}_{2}-\bm{r}_{1}$, and the position of the centre of mass by $\bm{R}\equiv(1/2)\left(\bm{r}_{1}+\bm{r}_{2}\right)$. Note that while the bead co-ordinates are allowed to sample the entirety of the three-dimensional coordinate space, the pulling is restricted to the $x$-axis alone. An alternative protocol in which the pulling direction is also in general three-dimensional space can be implemented, but this would not alter the analysis and arguments presented here. The positions of the two beads can be manipulated using optical traps, modelled here as harmonic potential wells. The trap stiffnesses are denoted by $H_{1}=c_{1}H$, and $H_{2}=c_{2}H$ (in units of the dumbbell spring constant), and the co-ordinates of the minimum of the wells are represented by $\bm{\chi}_{1}$ and $\bm{\chi}_{2}$, respectively. A temperature of $T=300\,K$ is considered in all our simulations, as a matter of convenience. The viscosity of the solvent at this temperature is taken to be $\eta_{\textrm{s},0}=0.001$ kg/m s, which is close to the viscosity of water at room temperature. In this protocol, values of solvent viscosity which are multiples of $\eta_{\textrm{s},0}$ will be considered.

The finite extensibility parameter, $b$, is defined as $b\equiv Q^{2}_{0}/l^{2}_{\textrm{H}}$. The internal friction parameter, $\varphi\,(\coloneqq 2K/\zeta)$, is defined as the ratio of the internal friction coefficient to the bead friction coefficient. The hydrodynamic interaction parameter is given by $h^{*}=a/\sqrt{\pi}l_{\textrm{H}}$, and $h^{*}=0$ corresponds to the free-draining case.

In Fig. 2 (b), the pulling protocol employed in this study is depicted. Without any loss of generality, $\bm{\chi}_{1}$ is chosen as the origin of our frame of reference. In all pulling simulations throughout this work, the first trap is held stationary, and the second trap is moved from its initial position, $\bm{\chi}_{2}^{\text{(i)}}\equiv(\chi^{\textrm{(i)}}_{2x},0,0)$, to its final position, $\bm{\chi}_{2}^{\text{(f)}}\equiv(\chi^{\textrm{(f)}}_{2x},0,0)$. The notation “$\chi^{\textrm{(i)}}_{2x}\to\chi^{\textrm{(f)}}_{2x}$” represents such a pulling event. The stretching distance is denoted by $d\equiv\left[\chi^{\textrm{(f)}}_{2x}-\chi^{\textrm{(i)}}_{2x}\right]$, the time interval for stretching by $\tau$, and the pulling velocity by $\bm{v}\equiv(v_{x},0,0)$, where $v_{x}=d/\tau$.

In the more general case, the Hamiltonian, $\mathcal{H}$, of the system for any value of the trap position, $\bm{\chi}_{2}$, is given by

[TABLE]

where $U_{\text{MS}}(\bm{Q})$ represents the potential energy in the Marko-Siggia spring. The generalized Jarzynski work corresponding to the pulling protocol discussed above is, in this case, given by

[TABLE]

where $\bm{v}=d\boldsymbol{\chi}_{2}/dt$. The average dissipation associated with the stretching process is calculated using Jarzynski’s equality as shown in Eq. (1).

II.2.2 Simulation details

The system of stochastic differential equations governing the time evolution of the non-dimensional connector vector $\bm{Q}^{*}$ and centre of mass coordinates $\bm{R}^{*}$, subjected to the pulling protocol depicted in Fig. 2, can be derived as shown in Appendix C. Brownian dynamics simulations are used here to solve these equations. To begin with, the initial values of $\bm{Q^{*}}$ and $\bm{R^{*}}$ are picked from a Gaussian distribution of zero mean and unit variance. With the first trap held at the origin, and the second at $\bm{\chi}^{\text{(i)}*}_{2}=(\chi^{\text{(i)}*}_{2x},0,0)$, the dumbbell is equilibrated for a duration of fifty-five dimensionless times. Equilibration is ascertained by checking that $\left<Q^{*2}\right>$ has reached a steady value with respect to time. Then, the pulling is commenced (at $t^{*}=0$), by varying the position of the second trap linearly, as $\chi^{*}_{2x}=\chi^{\text{(i)}*}_{2x}+v_{x}^{*}t^{*}$, till $t^{*}=\tau^{*}$. The window $[0,\tau^{*}]$ is uniformly divided into $N_{\text{t}}$ intervals, such that $\Delta t_{j}^{*}\equiv t^{*}_{j}-t^{*}_{j-1}=\tau^{*}/N_{\text{t}}$, where $j=1,2,...,(N_{\text{t}}+1)$.

The dimensionless equivalent of the work done by the mobile trap during one realization of the pulling event is calculated using a simple rectangular quadrature as follows,

[TABLE]

where the subscript $j$ on the first term indicates that it is evaluated at time, $t_{j}^{*}$, and $r_{2x}^{*}$ refers to the $x$-coordinate of the position of the dumbbell bead subjected to pulling. For representative values of the molecular and control parameters, which are discussed in more detail in the next subsection, the average dissipated work is computed using the time-step widths $\Delta t^{*}=\left\{10^{-5},10^{-4},10^{-3}\right\}$. The values of the average work calculated at all the time-steps concur within statistical error bars of the simulation, and the largest of the three time-step widths, i.e; $\Delta t^{*}=1\times 10^{-3}$, is used for all the cases where $c_{2}=1000$. For $c_{2}=100$, $\Delta t^{*}=1\times 10^{-2}$ is found to suffice, whereas $c_{2}=10000$ requires $\Delta t^{*}=1\times 10^{-4}$. The choice of the time-step width is also affected by the values of the internal friction and the finite extensibility parameter as higher values of these parameters necessitate the use of smaller time-steps.

The protocol proposed here involves pulling the molecule over a pre-determined distance at the same dimensional velocity but different solvent viscosities. In this context, it is essential to note that the timescale varies linearly with the solvent viscosity, $\lambda_{\text{H}}\propto\eta_{\text{s}}$. In order to maintain the same dimensional pulling time ($\tau=\tau^{*}\lambda_{\text{H}}$) across simulations with differing solvent viscosity, the dimensionless pulling time ($\tau^{*}$) is scaled by $1/\eta_{\text{s}}$ as the solvent viscosity is increased.

II.2.3 Molecular and control parameters

The parameters used in the present work are broadly classified into molecular and control parameters. Molecular parameters pertain to the polymer that is being stretched, whereas control parameters are set by the experiments or the simulations used in the study of stretching the molecule.

The choice of molecular parameters is based on the $\lambda$-phage DNA (48.5 kbp) used in Murayama et al.’s Murayama et al. (2007) work, which has a contour length, $L_{\text{c}}$, of $16.5\,\mu\text{m}$, and Kuhn segment length, $b_{\text{K}}$, of approximately 88 nm. In order to model this molecule as a dumbbell, the model parameters, $b$, and $l_{\text{H}}$, are chosen such that the contour length and the radius of gyration of the model and the DNA molecule are the same. Following the procedure for parameter selection described in detail in ref. 54, we obtain $b=811.25$ and $l_{\text{H}}=580.95\,\text{nm}$. We round down both these values, and use $b=800$ and $l_{\text{H}}=500\,\text{nm}$ as the parameters to model $\lambda$-phage DNA.

The Hookean spring constant of the model, $H$, is then found using

[TABLE]

The choice of the bead radius, $a$, is motivated by Alexander-Katz et al.’s Alexander-Katz et al. (2009) work, where it is suggested that the monomeric radius may be taken as the persistence length of the molecule. We choose $a=30\,\text{nm}$, which is identical to the choice made by Alexander-Katz et al. Alexander-Katz et al. (2009) for comparing the results of BD simulations against experiments on DNA.

Other values of $b$ and $l_{\text{H}}$, of the same order-of-magnitude as obtained for the $\lambda$-phage DNA case, have been used in this study and in addition to $a=30\,\text{nm}$, bead radii of $80\,\text{nm}$ and $100\,\text{nm}$ have also been used.

Fig. 3 provides a snapshot of the $x$-$y$ projection of the positions of an ensemble of beads of the dumbbell, obtained after an equilibration of fifty-five dimensionless times at the initial trap positions $\bm{\chi}^{*}_{1}=\left(0,0,0\right),\,\bm{\chi}^{\text{(i)}*}_{2}=\left(5,0,0\right)$ and final trap states $\bm{\chi}^{*}_{1}=\left(0,0,0\right),\,\bm{\chi}^{\text{(f)}*}_{2}=\left(7,0,0\right)$, as a function of the optical trap stiffness. It is clearly seen that the strength of the trap ($c_{1}$ or $c_{2}$) determines its ability to confine the bead near the position of its minimum. Since it is intended to hold the position of the first bead fixed at the origin, a trap strength of $c_{1}=1000$ is used throughout our simulations.

In Figs. 4, the effect of the mobile trap stiffness ($c_{2}$) on the dissipation is shown for a fixed value of the stationary trap stiffness, for freely-draining dumbbells. In Fig. 4 (a), the average dissipation as the dumbbell is pulled at a constant velocity and a fixed pulled distance, is plotted as a function of the mobile trap stiffness. It is seen that the dissipation reaches a plateau at $c_{2}\approx 100$. This behavior is in agreement with the dissipation calculated for the analytical model [Eq. (II.1)], which reaches $\sim 98\%$ of its asymptotic value at $c_{2}\approx 100$. In all our simulations, we set $c_{1}=c_{2}=1000$ (unless specified otherwise), in order to operate in a regime where the dissipation is independent of the trap stiffness.

In Fig. 4 (b), the average dissipation is plotted as a function of the pulling velocity, over the same fixed distance, for three different values of the mobile trap stiffness. At lower values of the trap stiffness, the dissipation grows linearly before scaling sub-linearly with the velocity. The onset of the non-linear regime is pushed to higher velocities as the trap stiffness is increased. In the asymptotic limit of high trap stiffness, the non-linear regime vanishes, and a linear scaling of the dissipation with the pulling velocity is observed for over two orders of magnitude in the velocity. There is a good agreement between the trends predicted by the simple analytical model using a Hookean spring-dashpot discussed in Sec. II.1 and BD simulations on a model with a nonlinear force law. Remarkably, the effect of the non-linear force law is more perceptible at lower trap stiffnesses, vanishing as the mobile trap stiffness is increased, until a quantitative agreement is obtained between the simulation and the analytical results.

The Marko-Siggia force expression is linear at low values of the extension of the molecule, and diverges as the fractional extension approaches unity. For initial extensions in the linear regime of the force-extension profile, a trap stiffness of $c_{2}=1000$ is sufficient to make the bead track the position of the trap, as shown in Fig. 3. Higher trap stiffnesses are found to be required for operating in the non-linear regime of the force-extension curve.

The Jarzynski equality is strictly exact only in the limit of an infinite number of work trajectories, $N\to\infty$. In applications of the JE, the number of trajectories required to accurately recover the free-energy difference increases with average dissipated work in the process, as discussed in refs. 38; 39; 40.

In Figs. 5, the effect of internal friction and pulling velocity on the probability distribution of the work trajectories is plotted for freely-draining dumbbells. The vertical green lines in the figures indicate the free-energy difference, $\Delta A^{*}$, obtained by taking an error-weighted mean of the values of the free energy difference obtained at dimensionless pulling velocities $v^{*}\leq 0.02$. Interestingly the work-distributions are normally distributed, as seen by the good agreement between the histogram data and the Gaussian fit. That it is so, despite the equations of motion for the system being non-linear, is an observation previously made by Speck and Seifert Speck and Seifert (2004). The Gaussianity of the distribution is attributed to the slow rate of the driving protocol ($1/\tau$) with respect to the molecular relaxation rate ($1/\lambda_{\text{H}}$). For all the data-sets plotted in Figs. 5, the driving rate is at least five times slower than the molecular relaxation rate.

From Fig. 5 (a), it is seen that increasing the pulling velocity at a fixed value of the internal friction parameter increases the average dissipated work, and the width of the distribution. An identical trend is observed in Fig. 5 (b), where an increase in the internal friction parameter at a fixed pulling velocity causes the work distribution to shift rightwards, and results in an increased dissipation. Thus, the dissipation in our model is directly correlated with the pulling velocity, and the internal friction in the system. Under such conditions of high dissipation, the estimates for $\Delta A$ are dominated by rare realizations that occur near the tail of the work distribution, necessitating the use of a larger number of trajectories to obtain an accurate estimate of the free energy difference.

In Fig. 6, we illustrate the above point using an alternative representation. Without prior knowledge of the ensemble size required for the simulations, an initial guess of $N=1\times 10^{4}$ was chosen. The green horizontal line represents the free-energy difference obtained by taking an error-weighted mean of the values of the free energy difference obtained at pulling velocities $v^{*}\leq 0.02$ using Jarzynski’s equality. An ensemble size of $N=1\times 10^{4}$ is sufficient for an accurate estimation of the free energy difference at lower velocities (dissipation), but is found to become inadequate at $v^{*}\geq 0.04$. Upon increasing the ensemble size for the higher velocity cases empirically, the accuracy of the estimated free energy difference improves. The ensemble size is plotted as a function of the average dissipated work in the inset of Fig. 6, which shows that the choice of the pulling velocity and the ensemble size are mutually related.

In Fig. 7, the average dissipated work (scaled by the thermal energy $k_{B}T$) calculated for a variety of molecular and control parameters is plotted against the magnitude of the dimensionless pulling velocity, $v^{*}$. It is seen that the average dissipated work varies linearly over the entire range of the pulling velocity, $v^{*}=0.001-0.1$. The velocity range in dimensional units would depend on the molecular parameters.

The experimental feasibility of the proposed protocol can be discussed in the context of the molecular parameters used for the dataset represented by filled circles in Fig. 7. For these set of parameters, the pulling velocities explored in Fig. 7 vary from $v=29.3\,\textrm{nm/s}\,(v^{*}=0.001)$ to $v=2.93\,\mu\textrm{m/s}\,(v^{*}=0.1)$. The molecule is stretched over a distance of $1\mu\textrm{m}$. The stiffness of this molecule is $H=1.657\times 10^{-5}$ pN/nm. In order to operate in a regime where the dissipated work is independent of the trap strength, as discussed earlier, the stiffness of the trap must be at least a hundred times that of the molecule, which implies $H_{\textrm{trap,min}}=1.657\times 10^{-3}$ pN/nm.

In Table 1, based on a survey of the literature, the range of trap stiffnesses, pulling velocities, and stretching distances typically accessible by optical tweezers is given. Additionally, the position and force resolution limits of optical traps are inversely correlated: stiffer traps improve the spatial resolution but also introduce large fluctuations in the measured force. A rough estimate of these resolution limits may be obtained using the equipartition theorem, as explained in refs. 57; 58. Most commercial optical tweezer setups are equipped with filtering mechanisms that aid in improving the precision in the measurements, by reducing the resolution limits Gupta et al. (2011). A detailed discussion of the resolution offered by optical tweezers can be found in refs. 57 and 58.

From Table 1, it is clear that the values of $v$, $d$, and $H_{\textrm{trap,min}}$ for the representative case lie well within the range of values explored experimentally.

Ritort et al. Ritort et al. (2002) have established from computer simulations of mechanical unfolding that the number of trajectories required to obtain estimates for free energy difference within an error of $\mathcal{O}\left(k_{B}T\right)$ increases exponentially with the average dissipation associated with the unfolding process. They predict that for dissipation less than $4k_{B}T$, around 100 trajectories would suffice, and for a dissipation of $5k_{B}T$, about 1000 trajectories would be required. These predictions agree well with the average dissipation and ensemble sizes encountered in optical-tweezer-based pulling experiments. For example, Liphardt et al. Liphardt et al. (2002) stretch RNA hairpins using optical tweezers, and estimate $\left<W_{\text{dis}}\right>=2\,\text{-}\,3\,k_{B}T$ with $N=47$. Similarly, for pulling experiments on DNA hairpins performed by Gupta et al. Gupta et al. (2011), $\left<W_{\text{dis}}\right>=1.1\pm 0.7k_{B}T$ for $N=99$, and $\left<W_{\text{dis}}\right>=4.9\pm 0.3k_{B}T$ for $N=1293$.

For simulations on the single-mode spring dashpot, the statistical error in the free-energy difference is maintained to be $\sim\mathcal{O}\left(0.01k_{B}T\right)$, in order to obtain a sufficiently accurate estimate of the average dissipated work that enables the internal friction coefficient to be extracted reliably. By restricting the regime of operation to the boxed region in Fig. 7, with $v^{*}\leq 0.02$ and $\left<W_{\textrm{dis}}\right>\sim k_{B}T$, it is found that $N=1\times 10^{4}$ trajectories are sufficient to obtain the free energy difference within the desired error limits. It is possible to operate at higher values of dissipation, outside the boxed regime, provided that the ensemble size is suitably increased, as shown in the inset of Fig. 6.

All the results prior to Fig. 7 have been presented for the case of freely-draining dumbbells. In Fig. 7, it is observed that the inclusion of fluctuating HI does not affect the dissipated work values in a single-mode spring dashpot. Speck and coworkers Speck (2017) have shown in the context of colloidal suspensions that the inclusion of hydrodynamic interactions does not alter the dissipation along a single trajectory. The effect of fluctuating hydrodynamic interactions on the dissipation is markedly different in the case of wet internal friction, as will be discussed in greater detail in Sec. III.4.

II.2.4 Dry friction in equals dry friction out

The methodology to extract the internal friction coefficient is illustrated using a molecule with parameters: $\left\{b=800,l_{\mathrm{H}}=500\,\mathrm{nm},K=3.0\times 10^{-9}\,\text{kg/s},h^{*}=0.0\right\}$ as an example. An ensemble of such molecules is pulled from an initial trap position of ${\chi}_{2x}^{\text{(i)}}=5\,l_{\text{H}}$ to a final trap position of ${\chi}_{2x}^{\text{(f)}}=7\,l_{\text{H}}$ at different pulling velocities. At each value of the pulling velocity, the average dissipated work is calculated at several values of the solvent viscosity in the range, $\eta_{\textrm{s}}=\eta_{\textrm{s},0}$ to $\eta_{\textrm{s}}=10\,\eta_{\textrm{s},0}$. In an experimental setting with water as the solvent, suitable viscogens, such as glucose or sucrose, may be added to the solvent in order to realize an approximately four-fold increase in its viscosity Jas et al. (2001); Qiu and Hagen (2004). In experiments that study the kinetics of intrachain contact formation in polypeptides Bieri et al. (1999) suspended in a solvent mixture of ethanol and glycerol, the solvent viscosity was varied over two orders of magnitude by adjusting the proportion of glycerol in the mixture.

As shown in Fig. 8 (a), for each value of the pulling velocity used, the average dissipated work in the hypothetical limit of zero solvent viscosity, $\left<W_{\text{dis}}\right>_{\eta_{\!{}_{\,\text{s}}}\to\,0}$, is obtained from a linear fit to the average dissipated work at finite solvent viscosities. Since the extrapolated value is finite, it is clear signature of the presence of dry internal friction. In Fig. 8 (b), the extrapolated values of the average dissipated work in the limit of zero solvent viscosity (divided by the stretching distance $d$), is plotted against the pulling velocity. The slope of the graph $(K_{\textrm{BD}})$ represents the internal friction coefficient extracted from simulations.

Table 2 shows a comparison between the value of the internal friction coefficient used as an input parameter in the Brownian dynamics simulations, and the corresponding value extracted from the dissipated work using the protocol proposed here, for various molecular and control parameters. Our protocol recovers the input internal friction coefficient to within $5\%$ accuracy, and is insensitive to the choice of the spring-force law, as shown for the Marko-Siggia and the FENE force laws in Table 2. Further, values of $K_{\textrm{BD}}$, for models with and without HI, lie close to each other, indicating that HI does not affect the dissipated work due to dry internal friction.

The validity of the proposed protocol for the extraction of the dry internal friction coefficient has thus been established. In the next section, the application of the protocol to a model with wet internal friction is discussed, in the context of a coarse-grained polymer model with cohesive intra-chain interactions.

III Wet internal friction

In this section, the case of wet internal friction is investigated in the context of the force-induced unraveling of a coiled globule which has previously been studied experimentally by Murayama et al. (2007) and with simulations by Alexander-Katz et al. (2009). The problem is revisited here with the goal of understanding the role of solvent viscosity and HI, both of which have not been considered previously.

III.1 Model description

A bead-spring model with $N_{\text{b}}$ beads connected by FENE springs, each stretchable up to a maximum length of $Q_{0}$, is considered. The excluded volume interactions between beads are modelled using the Soddemann-Dünweg-Kremer (SDK) potential Soddemann et al. (2001), whose functional form is given by

[TABLE]

The minimum of the potential occurs at at $r=2^{1/6}\sigma_{\text{s}}$, and the value of the potential at this location is $U_{\text{SDK}}=-\epsilon\,k_{B}T$. As seen from Eq. (25), the repulsive part of the pair-potential is modeled after the Weeks-Chandler-Andersen (WCA) potential, while the attractive part is constructed using a cosine function which smoothly approaches zero at the cut-off distance, $r_{\text{c}}$. A detailed comparison of the SDK potential with the Lennard-Jones and the WCA potential has been performed recently Santra et al. (2019). The parameter values for $m_{1}$ and $m_{2}$ depend on the choice of $r_{\text{c}}$. They have been refined in ref. 62 in comparison to the values in the paper by Soddemann et al. (2001) in order to be more applicable to BD simulations, as discussed below and in Sec. III.3. Setting $\epsilon=0$ in Eq. (25) leads to purely repulsive inter-bead interactions, and corresponds to the athermal limit of solvent quality. Increasing the value of $\epsilon$ beyond zero results in a decrease in the solvent quality. A special feature of the SDK potential Santra et al. (2019) is that modifying the value of $\epsilon$ allows one to tune the attractive interactions selectively, without affecting the repulsive branch of the pair-potential. This is in stark contrast to the more commonly used Lennard-Jones (LJ) potential for which changing the well-depth affects both the attractive and the repulsive branches. Furthermore, the exact truncation of the SDK potential at the cutoff distance results in an increased computational efficiency Soddemann et al. (2001) in comparison to the LJ potential whose long attractive tail approaches zero only asymptotically, at large inter-bead separations.

The force on the $\mu$th bead due to bonded and non-bonded excluded volume interactions is denoted by $\bm{F}^{\text{(c)}}_{\mu}$, and $\bm{F}^{\text{(E)}}_{\mu}$, respectively, and the notation $\bm{F}_{\mu}^{\phi}\equiv\bm{F}^{\text{(c)}}_{\mu}+\bm{F}^{\text{(E)}}_{\mu}$ is used to denote the total force on a bead due to the FENE and SDK potentials.

As discussed previously, the terminal beads of the chain are subjected to harmonic trap potentials. One of the traps is held fixed at $\boldsymbol{\chi}_{1}=(0,0,0)$, and the other is moved from $\bm{\chi}_{2}^{\text{(i)}}\equiv(\chi^{\textrm{(i)}}_{2x},0,0)$, to $\bm{\chi}_{2}^{\text{(f)}}\equiv(\chi^{\textrm{(f)}}_{2x},0,0)$ at a constant velocity, $\bm{v}\equiv\left(v_{x},0,0\right)$.

The hydrodynamic interaction between any pair of beads $\mu$ and $\nu$ is accounted for by defining the diffusion tensor $\bm{\Upsilon}_{\mu\nu}=\delta_{\mu\nu}\boldsymbol{\delta}+\zeta\boldsymbol{\Omega}_{\mu\nu}$, where $\delta_{\mu\nu}$ is the Kroenecker delta, and the hydrodynamic interaction tensor, $\boldsymbol{\Omega}_{\mu\nu}$, is approximated using the Rotne-Prager-Yamakawa expression, as shown in Eqs. (65) — (67). For notational convenience, we define $\bm{\mathcal{D}}$, a block matrix of size $N_{\text{b}}\times N_{\text{b}}$, whose each element is the $3\times 3$ matrix, $\bm{\Upsilon}_{\mu\nu}$. Additionally, the block matrix $\bm{\mathcal{B}}$ is defined as $\bm{\mathcal{B}}\cdot\bm{\mathcal{B}}^{T}=\bm{\mathcal{D}}$.

The stochastic differential equation describing the time-evolution of the position of the $\mu$th bead is given by

[TABLE]

where $\bm{\mathscr{B}^{*}}_{\mu\nu}$ is the $\left(\mu,\nu\right)-$th element of $\bm{\mathcal{B}^{*}}$, and $\Delta\bm{\mathscr{W}^{*}}$ is a dimensionless Wiener process of zero mean and variance $\Delta t^{*}$.

In Fig. 9 (a), a schematic representation of the pulling is shown, and in Fig. 9 (b), snapshots from BD simulations on a ten-bead chain is presented.

III.2 Simulation details

The stochastic differential equation governing the pulling of a single polymer chain [Eq. (26)] is solved numerically using Brownian dynamics simulations. The initial bead positions are picked from the equilibrium distribution function corresponding to the FENE force law. The chain is then equilibrated at the initial state, with $\chi^{*}_{2x}=\chi^{*\textrm{(i)}}_{2x}$, for fifty Rouse times. Equilibration is ascertained by checking that the mean-squared value of the dimensionless radius of gyration, $\left<R^{2*}_{\text{g}}\right>$, has reached a steady value with respect to time. Pulling is then commenced at $t^{*}=0$, by changing the position of the mobile trap linearly in time, as $\chi^{*}_{2x}=\chi^{*\textrm{(i)}}_{2x}+v^{*}_{x}t^{*}$, till $t^{*}=\tau^{*}$. The work done in one realization of the pulling event is given by Eq. (21), with $r^{*}_{2x}$ replaced with $r^{*}_{N_{\text{b}}x}$, where $r^{*}_{N_{\text{b}}x}$ refers to the $x$-coordinate of the last bead in the chain, and the remaining symbols retain their original meaning as defined in Sec. II.2.2. A timestep width of $\Delta t^{*}=1.0\times 10^{-4}$ is used after ascertaining, for a ten-bead chain with representative parameter values, that the average work obtained for the $\Delta t^{*}=1.0\times 10^{-4}$ and $\Delta t^{*}=1.0\times 10^{-5}$ cases agree within error bars. The calculated work is then used to estimate the free-energy difference and the average dissipation as shown in Eq. (1).

III.3 Parameter space specification

A value of $b=50$ is used throughout this section for the dimensionless FENE parameter. For a free chain (without confining potentials acting on the terminal beads), with this particular value of $b$, the choice of $r_{c}=1.82\sigma$, $m_{1}=1.5306333121$, and $m_{2}=1.213115524$ have been shown to lead to the correct scaling predictions in poor, $\theta$, and good solvent conditions for the radius of gyration with the number of beads in the chain, as described in detail in ref. 62, and hence have been used for all simulations on single chains discussed in this section. The $\theta$-temperature for this system is observed to be at $\epsilon\approx 0.45$, where $\left<R^{2}_{\text{g}}\right>\sim\left(N_{\text{b}}-1\right)$ and the second virial coefficient, $B_{2}=0$ Santra et al. (2019). This value of cohesive strength is denoted as $\epsilon_{\theta}$. Poor solvent scaling, namely, $\left<R^{2}_{\text{g}}\right>\sim\left(N_{\text{b}}-1\right)^{0.67}$, is observed for $\epsilon\geq 0.55$. The bead radius, $a$, is defined on the basis of $\sigma_{\text{s}}$, as $a=0.5\sigma_{\text{s}}$. We set $\sigma_{\text{s}}=l_{\text{H}}$ in all our simulations.

A trap stiffness of $c_{1}=c_{2}=1000$ is used for both the stationary and the mobile traps. The initial and final positions of the mobile traps are chosen as $\chi^{\textrm{(i)}}_{2x}=0.1\,L_{\text{c}}$ and $\chi^{\textrm{(f)}}_{2x}=0.3\,L_{\text{c}}$, respectively, where $L_{\text{c}}\equiv\left(N_{\text{b}}-1\right)Q_{0}$ is the contour length of the chain. A chain size of $N_{\text{b}}=10$ has been used for all the simulation results reported in this section.

III.4 Globule unraveling is wet

Netz and coworkers Alexander-Katz et al. (2009) have measured the internal friction associated with collapsed homopolymers by measuring the work dissipated in the force-induced unfolding of a single polymer chain. As mentioned previously, they estimate the free energy difference, denoted by $W_{\text{eq}}$, as the work done in the quasi-static pulling limit, that is, $W_{\text{eq}}\left(\epsilon\right)\equiv W\left(\epsilon,v\to 0\right)$. The dissipated work at any finite pulling velocity is then calculated as the difference between the average work done at that velocity, and the reversible work. This is exactly the definition of dissipation that we have adopted throughout this paper, as indicated in Eq. (1). However, rather than using the work done in the quasi-static limit to estimate the free-energy difference, we use the Jarzynski equality to evaluate the same quantity.

In Fig. 10, the total dimensionless work done in pulling a ten-bead chain with a representative set of parameters is plotted as a function of the dimensionless pulling velocity, for cases with and without hydrodynamic interactions. Horizontal lines represent the error-weighted mean of the total work at the lowest four values of the pulling velocity, and is therefore a measure of the free-energy difference in the classical sense. Solid lines correspond to the freely-draining case, while dashed lines indicate the case with hydrodynamic interactions. It is seen that the free-energy differences for cases with and without HI concur within $~{}10^{-1}\,k_{B}T$, in agreement with the expectation that the free-energy difference, a static equilibrium property, remains unaffected by hydrodynamic fluctuations. Furthermore, it is also seen from the figure that the work done in the quasi-static limit agrees, within error bars, with the free-energy difference estimated using the Jarzynski equality (square symbols), thus establishing the validity of our approach for the estimation of the free-energy difference and the dissipation.

In Fig. 11, the dimensionless average dissipation is plotted as a function of the dimensionless pulling velocity. It is seen that the data for the phantom chain (no excluded volume interactions), athermal chain (purely repulsive interactions) and the chain under $\theta-$conditions superimpose. This indicates that for uncollapsed chains, the work dissipated during pulling is expended entirely against the solvent friction, and there is no internal friction associated with these cases. Furthermore, increasing the well-depth beyond the $\theta-$point results in an increased dissipation at a fixed value of the pulling velocity, which is due to the additional work needed to unravel the globule as the polymer chain is stretched. The average dissipation scales linearly with the pulling velocity over the entire range of the latter quantity.

In the linear regime, for a collapsed globule ($\epsilon>\epsilon_{\theta}$), the following relationship for the dissipation can be written,

[TABLE]

where $\gamma_{\text{G}}$ is the globule friction coefficient. The same expression is valid in the uncollapsed case, with $\gamma_{\text{G}}$ replaced with $\gamma_{\text{s}}$, the solvent friction coefficient. As $\epsilon\to\epsilon_{\theta}$, the dissipation is entirely due to the solvent, and $\gamma_{\text{G}}\to\gamma_{\text{s}}$.

Fig. 12 elucidates the effect of hydrodynamic interactions on the dissipation more clearly. As shown in Fig. 12 (a), for chains under both $\theta$ and athermal conditions, the dissipation due to the solvent decreases identically upon the inclusion of fluctuating hydrodynamic interactions. In Fig. 12 (b), it is seen that the enhancement in the dissipation due to cohesive interactions, measured as the difference between the total dissipation and the dissipation due to solvent alone, remains practically unaffected by hydrodynamic interactions.

In Fig. 13, the average dissipated work is plotted for two values of the well-depth, for one value of the dimensional pulling velocity, as a function of the solvent viscosity. At any finite value of the solvent viscosity, the work dissipated for the collapsed globule case is greater than that for that in $\theta-$condition. However, in the extrapolated limit $\eta_{\text{s}}\to 0$, the dissipated work goes to zero. Clearly, this suggest that the additional dissipation due to cohesive interactions between the beads corresponds to the case of wet internal friction, which is not clear a priori, and can only be established following the protocol proposed here.

In Figs. 14, the friction coefficients calculated using Eq. (27) have been plotted as a function of the solvent viscosity. From Fig. 14 (a), it is seen that friction coefficient due to the solvent scales linearly with the solvent viscosity, and that hydrodynamic interactions reduce the friction coefficient in comparison to the freely draining case. From Fig. 14 (b), it is seen that the inclusion of cohesive interactions results in an enhancement in the friction coefficient at all finite values of the solvent viscosity. In the extrapolated limit of zero solvent viscosity, however, the friction coefficient is also seen to tend to zero, as is typical of wet friction.

By taking a ratio of the slopes of the dissipation versus pulling velocity for the collapsed and the uncollpased states, for identical pulling distances in the linear regime, one gets

[TABLE]

This quantity represents the rescaled excess contribution to the dissipation due to internal friction, and enables an investigation of the relationship between $\gamma_{\text{G}}$ and $\gamma_{\text{s}}$, and comparison with Zwanzig’s prediction Zwanzig (1988), as was done previously by Alexander-Katz et al. (2009) in the absence of hydrodynamic interactions. For a Brownian particle moving in a corrugated one-dimensional potential of the form $U(x)=\left(\omega/2\right)\left(\epsilon-\epsilon_{\theta}\right)\sin\left(\pi x/a\right)$, the effective friction and the solvent friction are related by Alexander-Katz et al. (2009); Zwanzig (1988)

[TABLE]

where $I_{0}(..)$ is the modified Bessel function of zeroth order, and $\omega$ is a fitting parameter.

In Fig. 15, the rescaled excess dissipation due to internal friction, for models with and without fluctuating hydrodynamic interactions, is plotted as a function of the well-depth relative to the $\theta-$condition. It is seen from the figure that a good qualitative agreement is observed between the simulation results and Eq. (29), but with a fitting parameter that depends on whether hydrodynamic interactions are incorporated in the simulations. This agreement between simulations and theory, also observed by  Alexander-Katz et al. (2009), suggests that Zwanzig’s formulation, albeit based on a one-dimensional energy landscape, satisfactorily captures the scaling of internal friction with the strength of cohesive interactions in force-spectroscopy simulations on single molecules. The subject of diffusion on rugged energy landscapes of dimension higher than one has been treated rigorously in refs. 63 and 64.

IV Conclusions

In summary, we have introduced a simple and novel protocol based on the Jarzynski equality for determining both dry and wet internal friction coefficients of macromolecules that can be implemented experimentally using optical tweezers. Using Brownian dynamics simulations on a spring-dashpot model for a polymer, we establish proof-of-principle by recovering the dry internal friction coefficient which is used as a model input, and show that a bead-spring chain with cohesive interactions is an example of wet friction. It is conceivable that some real polymer chains might possess both wet and dry internal friction, and modeling such molecules would require the use of multi-bead-spring-dashpots with cohesive interactions. We envisage that the scheme proposed here may be applicable to a variety of macromolecules, and would enable a succinct characterization of the dissipative properties of the molecule.

Acknowledgements

We thank Burkhard Dünweg, Subhashish Chaki, Ranjith Padinhateeri and Gareth McKinley for enlightening discussions. The work was supported by the MonARCH and SpaceTime-2 computational facilities of Monash University and IIT Bombay, respectively. CPU time grants from the National Computational Infrastructure (NCI) facility hosted by Australian National University is also acknowledged. R.C. acknowledges IRCC-IIT Bombay for funding (Project No. RD/0518-IRCCAW0-001). We also acknowledge the funding and general support received from the IITB-Monash Research Academy.

Appendix A Governing equation for the one-dimensional spring-dashpot

The Fokker-Planck equation for the configurational distribution function $\Psi\left(r_{2},t\right)$ can be derived by following the procedure commonly used in polymer kinetic theory, i.e., by combining a force balance on the beads with an equation of continuity in probability space Bird et al. (1987). The force balance essentially states that the (i) the internal friction force due to the dashpot, (ii) the restoring force due to the finitely extensible spring, (iii) external forces (like the force due to the optical traps, in the present case), (iv) the random Brownian force due to bombardment by solvent molecules, and (v) the hydrodynamic force responsible for the solvent-mediated propagation of momentum on each bead, must sum up to zero.

The force balance over the free bead may be written as

[TABLE]

Upon simplification,

[TABLE]

Grouping together the terms containing $\llbracket\dot{r}_{2}\rrbracket$, and defining $\theta=\left[1+\left(K/\zeta\right)\right]$, where $\varphi\coloneqq 2K/\zeta$ is the internal friction parameter,

[TABLE]

The equation of continuity for the probability density, $\Psi\left(r_{2},t\right)$, is written as

[TABLE]

Substituting the expression for $\llbracket\dot{r}_{2}\rrbracket$ from Eq. (32) into Eq. (33), one obtains

[TABLE]

which can be rewritten as

[TABLE]

The stochastic differential equation corresponding to Eq. (35) is given by Öttinger (1996)

[TABLE]

where $W_{t}$ represents a Wiener process, and has dimensions of $[\text{time}]^{1/2}$. Eq. (36) can be recast in the Langevin form as

[TABLE]

where $\left<f_{\text{B}}(t)\right>=0$ and $\left<f_{\text{B}}(t)f_{\text{B}}(t^{\prime})\right>=\delta(t-t^{\prime})$, and $f_{\text{B}}(t)$ has dimensions of $[\text{time}]^{-1/2}$. Setting $f^{*}_{\text{B}}(t^{*})=f_{\text{B}}(t)\sqrt{\lambda_{\text{H}}}$, where $\left<f^{*}_{\text{B}}(t^{*})\right>=0$ and $\left<f^{*}_{\text{B}}(t^{*})f^{*}_{\text{B}}(t_{1}^{*})\right>=\delta(t^{*}-t_{1}^{*})$, the variables in Eq. (37) are cast into their dimensionless form as,

[TABLE]

Multiplying Eq. (A) throughout by $\left(\lambda_{\text{H}}/l_{\text{H}}\right)$, and simplifying, the dimensionless governing equation is obtained as

[TABLE]

Eq. (39) is rewritten as Eq. (3), after further simplification.

Appendix B Work distribution for freely draining one-dimensional Hookean spring-dashpot

The intermediate steps in the derivation of the mean and variance of the probability of work distribution is provided in this section.

Firstly, the complete expression for the work done during one realization of pulling is obtained by substituting the expression for $r^{*}_{2}$ from Eq. (II.1) into Eq. (6), as

[TABLE]

where the noise term, $\xi(t^{*})$, obeys $\left<\xi(t^{*})\right>=0$ and

[TABLE]

The dimensionless partition function of the system can be derived to be,

[TABLE]

and the dimensionless free-energy, $A^{*}=-\ln Z^{*}$, is then simply

[TABLE]

after ignoring constant prefactors. This equation is reproduced as Eq. (7) in Sec. II. The probability distribution function for the position of the bead is given by

[TABLE]

and the moments of the distribution are

[TABLE]

The expression for the average work is obtained by taking an ensemble average of Eq. (40),

[TABLE]

Substituting Eq. (44) into the second term on the RHS of Eq. (46),

[TABLE]

The underlined term is simplified as

[TABLE]

Integrating expression on the RHS of Eq. (48) by parts,

[TABLE]

one obtains,

[TABLE]

which can then be written as

[TABLE]

Substituting Eq. (51) into Eq. (47) and simplifying,

[TABLE]

one obtains

[TABLE]

Recognizing that the first term on the RHS of Eq. (53) is the free-energy difference, $\Delta A^{*}\equiv A^{*}\left[\chi^{*}(\tau^{*})\right]-A^{*}\left[\chi^{*}(0)\right]$, with $A$ given by Eq. (43), the expression for the average work can be rewritten as

[TABLE]

This equation is reproduced as Eq. (8) in Sec. II. Using Eqs. (40) and (46), the variance of the work distribution, $\sigma^{2}=\left<\left(W^{*}-\left<W^{*}\right>\right)^{2}\right>$, is written as,

[TABLE]

Eqs. (45) and (41) can be used to simplify the underlined terms in Eq. (55), and we obtain,

[TABLE]

Subsequent integration over $\tilde{t}^{*}_{1}$ in the second integral yields

[TABLE]

The following identity

[TABLE]

can be used used to simplify the second term on the RHS of Eq. (57), resulting in

[TABLE]

This equation is reproduced as Eq. (9) in Sec. II.

Appendix C Governing equation and simulation details for the dumbbell obeying the Marko-siggia force law

The dumbbell model described in Sec. II.2.1 is suspended in a fluid where the velocity field at any location $\bm{r}_{\text{f}}$ is given by $\bm{v}_{\text{f}}(\bm{r},t)\equiv\bm{v}_{0}+\boldsymbol{\kappa}(t)\cdot\bm{r}_{\text{f}}$, where $\bm{v}_{0}$ is a constant vector, and $\boldsymbol{\kappa}\equiv\left(\nabla\bm{v}_{\text{f}}\right)^{\bm{T}}$ is the transpose of the velocity gradient tensor. In the present work, both $\bm{v}_{0}$ and $\boldsymbol{\kappa}$ are set to $\bm{0}$ as the pulling experiments are simulated in a quiescent fluid. However, these terms have been included in the governing equations for the sake of generality. The configurational distribution function, $\Psi\left(\bm{Q},\bm{R},t\right)$, denotes the probability of finding the dumbbell at a position between $\bm{R}$ and $\bm{R}+d\bm{R}$, with an extension that lies between $\bm{Q}$ and $\bm{Q}+d\bm{Q}$, at any time $t$. The total restoring force in the connector vector, $\bm{F}^{\text{(c)}}$, has contributions from both the spring, $\bm{F}^{\text{(s)}}\equiv\partial U_{\text{MS}}/\partial\bm{Q}$, and the dashpot, $\bm{F}^{\text{(d)}}$

[TABLE]

with

[TABLE]

where the potential energy of the spring, $U_{\text{MS}}$, is given by

[TABLE]

and

[TABLE]

where $\llbracket{\dot{\bm{Q}}}\rrbracket$ is the momentum-averaged rate-of-change of the connector vector, $\bm{Q}$.

The force balance on the beads can be solved to obtain the following equations of motion for the position vectors of the beads Kailasham et al. (2018),

[TABLE]

and

[TABLE]

where $\bm{\Omega}$ is the hydrodynamic interaction tensor, defined as

[TABLE]

with $h=\left(3/4\right)a$ and $\bm{M}=2\left(\boldsymbol{\delta}-\zeta\boldsymbol{\Omega}\right)$. The terms $A$ and $B$ depend on the choice of the expression for the HI tensor, as shown in ref. 42. Here we choose the Rotne-Prager-Yamakawa (RPY) expression for the HI tensor Rotne and Prager (1969); Yamakawa (1971), in which the variables $A$ and $B$ are defined as follows

[TABLE]

The quantity $\beta$ that appears in Eqs. (C) and (C) is defined as

[TABLE]

Using $\bm{r}_{1}=\bm{R}-(1/2)\bm{Q}\,;\,\bm{r}_{2}=\bm{R}+(1/2)\bm{Q}$, and the chain-rule for partial differentiation to operate on $\partial\Psi/\partial\bm{r}_{1}$ and $\partial\Psi/\partial\bm{r}_{2}$, leads to,

[TABLE]

and

[TABLE]

By adding and subtracting Eqs. (C) and (C) suitably, we obtain

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

and both $\bm{X}$ and $\bm{Y}$ have dimensions of force.

The equation of continuity in terms of $\bm{R}$ and $\bm{Q}$ is given by Bird et al. (1987),

[TABLE]

Substituting Eqs. (C) and (C) into the above expression leads to the Fokker-Planck equation that governs the configurational distribution function $\Psi\left(\bm{Q},\bm{R},t\right)$,

[TABLE]

We define the following dimensionless quantities,

[TABLE]

In terms of these non-dimensional variables, the Fokker-Planck equation assumes the following form,

[TABLE]

where

[TABLE]

with $\alpha=\left(3/4\right)\sqrt{\pi}h^{*}$ and $h^{*}=a/\left(\sqrt{\pi}l_{\text{H}}\right)$. The quantity $\beta^{*}$ is dimensionless and is defined as,

[TABLE]

Using the following identity,

[TABLE]

where $\bm{L}$ is a tensor and $f$ is a scalar, the second and fourth terms on the right-hand-side of Eq. (C) can rewritten in a way that renders the Fokker-Planck equation amenable to Itô’s interpretation.

Since $\boldsymbol{\hat{\Omega}}\left(\bm{Q^{*}}\right)$ is independent of $\bm{R^{*}}$, and

[TABLE]

the Fokker-Planck equation can be rewritten as follows,

[TABLE]

where

[TABLE]

The Fokker-Planck equation (Eq. (C)) can be written in the following compact form,

[TABLE]

where the definitions of the quantities $\bm{e},\tilde{\bm{e}},\bm{g}$ and $\tilde{\bm{g}}$ are clear by comparison of Eqs. (C) and (82).

It is convenient to define a collective variable, $\bm{C}$, which is a six-element vector containing the components of $\bm{R^{*}}$ and $\bm{Q^{*}}$, such that $\bm{C}\equiv\left[R^{*}_{x},R^{*}_{y},R^{*}_{z},Q^{*}_{x},Q^{*}_{y},Q^{*}_{z}\right]$. Similarly, a six-element vector $\bm{j}$ can be defined, containing the components of $\bm{e}$ and $\bm{g}$, along with the definition of a $2\times 2$ block matrix $\bm{D}$, whose off-diagonal elements are $\bm{0}$, and the diagonal elements are the matrices $\tilde{\bm{e}}$ and $\tilde{\bm{g}}$ (each of size $3\times 3$). With these definitions, the Fokker-Planck equation in Eq. (82) can be written as

[TABLE]

The stochastic differential equation (SDE) corresponding to Eq. (83) can be obtained using the Itô interpretation, as

[TABLE]

where $\bm{w}_{t}$ is a Wiener process and $\tilde{\bm{b}}\cdot\tilde{\bm{b}}^{T}=\bm{D}$. The SDE (Eq. (84)) is solved using a semi-implicit predictor-corrector scheme Öttinger (1996), as discussed in the following subsection.

Solver details

With reference to Eq. (83), $\bm{D}$ is a $6\times 6$ matrix, and its square root, $\tilde{\bm{b}}_{j}$, at any time $t^{*}_{j}$ is found using Cholesky decomposition Press et al. (2007). Although Eq. (84) is written in terms of the collective variable $\bm{C}$, the equation for $\bm{R^{*}}$ is solved purely explicitly, whereas the equation in $\bm{Q^{*}}$ is solved by treating only the spring force term implicitly. For the sake of clarity, the predictor and corrector equations for $\bm{R^{*}}$ and $\bm{Q^{*}}$ are presented separately. It is useful to define another six-element vector, $\Delta\bm{S}_{j}$, as

[TABLE]

where $\bm{w}_{j}$ is a vector of six independent Wiener processes, each of mean zero and variance $\Delta t^{*}_{j}$. The first three elements of $\Delta\bm{S}_{j}$, denoted by $\Delta\bm{S}^{\,(R)}_{j}$, contain the noise contribution to $\bm{R}^{*}$, and the next three elements, denoted by $\Delta\bm{S}^{\,(Q)}_{j}$, contribute to the noise in $\bm{Q^{*}}$. In the following discussion, Eqs. (86)—(93) are in their dimensionless form, but the asterisk has been dropped from these equations for the sake of notational simplicity. Predictor step

[TABLE]

[TABLE]

where

[TABLE]

and the notations $\bm{Q}_{j}$ and $\bm{Q}(t_{j})$ have been used interchangeably to refer to the same quantity.

Corrector step

[TABLE]

[TABLE]

where

[TABLE]

By setting the length of the vector on the RHS of Eq. (90) to be $L$, and the length of $\bm{Q}(t_{j+1})$ to be $\Phi$, the following cubic equation is obtained:

[TABLE]

where

[TABLE]

Eq. (92) has three roots—two complex and one real—and the real root is obtained using the Newton-Raphson scheme Press et al. (2007). Note that the equations are solved in their dimensionless form, and the dimensional quantities are obtained by a suitable multiplication with the scaling factors, as explained in the discussion surrounding Eq. (75).

Appendix D Code validation for the dumbbell case

The total Hamiltonian of the dumbbell and trap system is written as

[TABLE]

The expression for $\mathcal{H}^{*}$ can be rewritten in terms of $\bm{Q^{*}}$ and $\bm{R^{*}}$ as,

[TABLE]

The steady-state configurational distribution function can be written as

[TABLE]

where $\mathcal{Z}$ is the partition function of the system, given by

[TABLE]

Substituiting the definition of $\mathcal{H}^{*}$ from Eq. (95) into Eq. (97) yields the following equation,

[TABLE]

where

[TABLE]

The inner integral in Eq. (98) can be evaluated using the following identity Bird et al. (1987) for Gaussian integrals,

[TABLE]

resulting in

[TABLE]

Upon simplification, one obtains

[TABLE]

where $k=(c_{1}c_{2})/2(c_{1}+c_{2})$, and $\bm{s^{*}}=\bm{\chi^{*}}_{2}-\bm{\chi^{*}}_{1}$. The integral in Eq. (102) can be evaluated by converting to spherical co-ordinates, recognising that $Q^{*}_{x}=Q^{*}\sin\theta\cos\phi$, $Q^{*}_{y}=Q^{*}\sin\theta\sin\phi$, $Q^{*}_{z}=Q^{*}\cos\theta$. Therefore,

[TABLE]

The integral in Eq. (LABEL:eq:norm_spher) does not have an analytically closed-form solution, and is evaluated numerically using MATLAB. The free-energy difference in going from the initial state to the final state is then given by,

[TABLE]

where the subscript ‘num’ indicates that the free energy difference has been calculated numerically.

Fig. 16 shows a comparison between the free energy difference obtained from Brownian dynamics simulations of $N=1\times 10^{5}$ trajectories using Jarzynski’s equality, and that obtained from Eq. (104), for the two parameter sets indicated in Table 3.

Appendix E Code validation for the single chain case

A chain of $N_{\text{b}}$ beads connected by $N_{\text{s}}$ springs is considered, with the stiffness of each spring denoted by $H$. The first bead is held fixed at the origin, and the last bead is subjected to a harmonic trap of stiffness $c_{2}H$. The trap is moved from an initial position of $\chi^{\text{(i)}}$ to a final position of $\chi^{\text{(f)}}$, over a time $\tau$. The distance traveled by the trap is denoted by $d\equiv\chi^{\text{(f)}}-\chi^{\text{(i)}}$, and the pulling velocity, $v$, given by $v=d/\tau$.

By following the procedure proposed by Dhar Dhar (2005), and using the non-dimensionalization scheme introduced in this paper [see Eq. (75), for example], the free energy difference associated with the stretching process is evaluated to be

[TABLE]

and $P^{*}\left(W^{*}\right)$ is found to be a Gaussian of the form,

[TABLE]

whose variance and mean are given by

[TABLE]

[TABLE]

where the notation $\left\{...\right\}_{ij}$ refers to the $\left(ij\right)^{\text{th}}$ matrix element, and $\bm{E}$ is a symmetric, tridiagonal $N_{\text{s}}\times N_{\text{s}}$ matrix of the following form

[TABLE]

The pulling of Hookean chains was simulated using Brownian dynamics, with an ensemble size of $\mathcal{O}(10^{5})$. The first bead was held fixed at the origin by means of a stiff harmonic trap of strength $c_{1}=1000$. In Fig. 17, the probability distribution of work values obtained from BD simulations is compared against the analytical solution for two sample cases. The good agreement between the two suggests the validity of the code.