Uniform description of polymer ejection dynamics from capsid with and without hydrodynamics

TL;DR

This study investigates polymer ejection from a capsid using stochastic rotation dynamics, revealing how hydrodynamics influence ejection speed, polymer expansion, and the universal scaling of waiting times, with implications for understanding ejection times.

Contribution

It introduces a universal scaling function for polymer ejection dynamics and compares hydrodynamic effects using SRD and Langevin simulations, providing new insights into ejection time scaling.

Findings

Hydrodynamics speed up ejection and allow closer to equilibrium expansion.

Waiting times grow with ejected monomers as a sum of two exponents.

Ejection time scales linearly with polymer length for long polymers.

Abstract

We use stochastic rotation dynamics to examine the dynamics of the ejection of an initially strongly confined flexible polymer from a spherical capsid with and without hydrodynamics. The results obtained using SRD are compared to similar Langevin simulations. Inclusion of hydrodynamic modes speeds up the ejection but also allows the part of the polymer outside the capsid to expand closer to equilibrium. This shows as higher values of radius of gyration when hydrodynamics are enabled. By examining the waiting times of individual polymer beads we find that the waiting time $t_w$ grows with the number of ejected monomers $s$ as a sum of two exponents. When $\approx 63 \%$ of the polymer has ejected the ejection enters the regime of slower dynamics. The functional form of $t_w$ vs $s$ is universal for all ejection processes starting from the same initial monomer densities. Inclusion of hydrodynamics only reduces its magnitude. Consequently, we define a universal scaling function $h$ such that the cumulative waiting time $t = N_0 h(s/N_0)$ for large $N_0$. Our unprecedently precise measurements of force indicate that this form for $t_w(s)$ originates from the corresponding force towards the pore decreasing super exponentially at the end of the ejection. Our measured $t_w(s)$ explains the apparent superlinear scaling of the ejection time with the polymer length for short polymers. However, for asymptotically long polymers $t_w(s)$ predicts linear scaling.

Figures

Tables

Table 1. Table 1: The exponents β 𝛽 \beta of the fits to the apparent relation τ ∼ N 0 β similar-to 𝜏 superscript subscript 𝑁 0 𝛽 \tau\sim N_{0}^{\beta} plotted in Fig.  4  (a),(b), and (c) for different initial monomer densities ρ 0 subscript 𝜌 0 \rho_{0} . The errors of the fits are of the order 0.04.

${\rho }_0$	HD	noHD	Langevin
0.50	1.30	1.36	1.37
0.75	1.26	1.33	1.33
1.00	1.25	1.29	1.28
1.25	1.22	1.30	1.29
1.50	1.22	1.29	1.27

Table 2. Table 2: The exponents α 𝛼 \alpha of the fits τ ∼ ρ 0 − α similar-to 𝜏 superscript subscript 𝜌 0 𝛼 \tau\sim\rho_{0}^{-\alpha} plotted in Figs.  4  (d),(e), and (f) for different polymer lengths N 0 subscript 𝑁 0 N_{0} . The errors of the fits are of the order 0.1.

$N_0$	HD	noHD	Langevin
25	0.77	0.89	0.83
35	0.85	0.89	0.82
50	0.94	0.96	0.90
71	0.98	1.00	0.97
100	0.96	1.00	1.03
141	0.99	1.08	1.01
200	0.96	1.05	1.07
283	1.00	0.99	0.98

Equations

$$U_F=-\frac{K}{2}{r}_{\max }^2\ln \left(1-{\left(\frac{r}{r_{\max }}\right)}^2\right),r<r_{\max },$$

$$\displaystyle U_{LJ}=\left\{\begin{array}[]{ccl}4.8\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\ \right]+1.2\epsilon&,&r_{ij}\leq\sqrt[6]{2}\sigma\\ 0&,&r_{ij}>\sqrt[6]{2}\sigma\end{array},\right.$$

$${\mathbf{r}}_i(t+\Delta t)={\mathbf{r}}_i(t)+{\mathbf{v}}_i(t)\Delta t,$$

$${\mathbf{v}}_i(t+\Delta t)={\mathbf{v}}_{\mathrm{cm}}(t)+\Omega \left[{\mathbf{v}}_i(t)-{\mathbf{v}}_{\mathrm{cm}}(t)\right],$$

$${\rho }_0=\frac{N_0}{\frac{4}{3}\pi R_0^3}.$$

$$\frac{\mathrm{d}{\mathbf{p}}_i}{\mathrm{d}t}(t)=-\xi {\mathbf{p}}_i(t)+{\mathbf{\eta }}_i(t)+{\mathbf{f}}_i(t),$$

$$v=\frac{f_{\mathrm{drag}}}{\xi m_b{N}_0},$$

$$\frac{V_{\mathrm{DNA}}}{V_{\mathrm{capsid}}}=\frac{48502\pi (1\mathrm{nm}{)}^{2}0.34\mathrm{nm}}{\frac{4}{3}\pi (27.5\mathrm{nm}{)}^3}\approx \mathrm{0.6.}$$

$$t_w(s)=t(s)-t(s-1),$$

$$t_w(s)=A\left[\exp \left(\frac{2.8}{N_0}s\right)+\exp \left(\frac{10.8}{N_0}\left(s-0.625N_0\right)\right)\right],$$

$$h\left(\frac{s}{N_0}\right)=$$

$$+\frac{1}{10.8}\exp \left[\frac{10.8}{N_0}\left(s-0.625N_0\right)\right]-\text{const.}\}.$$

$${\rho }_e=\frac{N}{V_e}=\frac{N}{\frac{4}{3}\pi (R_0+ϵ{)}^3},$$

$$t_f(s)=\frac{kT}{\xi m}{\int }_0^s\exp \left(-\frac{U(y)}{kT}\right){\int }_y^s\exp \left(\frac{U(z)}{kT}\right)\mathrm{d}z\mathrm{d}y,$$

$$U(s)=-{\int }_0^{s}f(s)\mathrm{d}s=\frac{C}{B}\left[N_0\exp \left(B\frac{N_0-s}{N_0}\right)+Bs\right]$$

Full text

I Introduction

Packaging and ejection of macromolecules in confinements are of high interest due to the potential technological and medical applications, such as drug delivery and gene therapy Glasgow and Tullman-Ercek (2014). The basic understanding of these processes is important also due to the relevant fundamental biological processes, the most prominent of which is the viral packaging in and ejection from bacteriophages Muthukumar (2001); Forrey and Muthukumar (2006); Smith et al. (2001); Grayson et al. (2007); Ali et al. (2006); Ghosal (2012); Cacciuto and Luijten (2006a); Sakaue and Yoshinaga (2009); Linna et al. (2014). The theoretical treatments of polymer ejection are strictly based on fully flexible chains Muthukumar (2001); Cacciuto and Luijten (2006a); Sakaue and Yoshinaga (2009), in spite of the experimental studies being almost solely done on semiflexible double-stranded DNA. This makes sense, since theoretically ejection from confinements is most intriguing when the spring force of the semiflexible polymer does not dominate over the more subtle mechanisms. Besides, the ejection of fully flexible polymers is highly relevant outside the purely theoretical realm due to many important polymers, such as proteins, single-stranded DNA, and RNA, belonging to this class.

In our previous study we showed that the blob-scaling picture used as a basis for analyzing the ejection dynamics from strong confinement is not valid. The blob picture presumes semidilute conditions, which does not hold for in vivo encapsulated polymers. In computer simulations polymers are far too short to justify the blob-scaling assumption. In spite of these shortcomings the blob-scaling has been used to explain the apparent scaling of the ejection time $\tau$ with the length of the polymers $N_{0}$. Indeed, if only the ejection time $\tau$ as a function of $N_{0}$ is measured, $\tau$ seemingly scales superlinearly with $N_{0}$. However, by inspecting the waiting times $t_{w}(s)$ we showed that there, in fact, is no such scaling. $t_{w}(s)$ is defined as the time it takes for a monomer labeled $s$ to translocate after the previous monomer $s-1$ has translocated. From our simulations using a hybrid method consisting of molecular dynamics (MD) and stochastic rotation dynamics (SRD) Malevanets and Kapral (1999, 2004) we obtained $t_{w}$ that grows essentially exponentially with $s$ Piili and Linna (2015).

Due to this exponential growth of $t_{w}(s)$ conflicting the theoretical predictions, some suspicion was cast on the simulation method. Accordingly, to be conclusive a confirmation using a well established method is called for. To this end, we have implemented an identical capsid model in our Langevin Dynamics (LD) algorithm Allen and Tildesley (2006). LD algorithm is a numerical implementation of a stochastic differential equation describing Brownian motion of particles, so it can be regarded as the most fundamental method available for a dynamical simulation of polymer ejection. We present here a thorough comparison of polymer ejection dynamics obtained using LD and SRD.

Our main objective in the present paper is to establish precise forms for $t_{w}(s)$ in the absence and presence of hydrodynamic interactions, that is, to determine how the inclusion of hydrodynamic modes changes the ejection dynamics. To achieve this we use SRD where hydrodynamic modes can easily be switched on or off. Measuring $t_{w}(s)$ is the most precise way of gaining detailed information on ejection dynamics. The form of $t_{w}(s)$ reflects the form of the force $f$ reduced to the pore during the ejection. We determine $f(s)$ with high precision using LD that due to being computationally more effective than SRD allows us to gain much better statistics than what was possible in our previous study using only SRD Piili and Linna (2015). This study reveals the surprisingly strong influence of the local effects in the vicinity of the pore on the overall ejection dynamics. Indications of this were observed already in our earlier study on capsid ejection, where a force applied to aligning a polymer close to the pore was found to give an effective bias to the ejection Linna et al. (2014).

The paper is organized as follows. The computational models are described in Section II. The procedure for matching the models based on LD and SRD is described at the end of this section. Results are presented in Section III. In this section we first compare the results obtained by SRD and LD after which we extract precise forms for $t_{w}(s)$ from each model thus pinpointing the effect of hydrodynamics. We then present our measurements of the radius of gyration of the polymer segment outside the capsid and the force at the pore. By analyzing these measurements we are able to give an accurate account of the prevailing mechanisms during polymer ejection. Finally, in Section IV we summarize our results and present the conclusions based on them.

II The Computational Models

Here, we describe the computational models. The primary computational method is stochastic rotation dynamics (SRD), also called multi-particle collision dynamics, that allows for the inclusion of hydrodynamics Malevanets and Kapral (1999, 2004). We validate our SRD model by closely comparing it to the identically implemented model, see Fig. 1, in our Langevin dynamics (LD) algorithm. LD is based on a thoroughly analyzed and understood stochastic differential equation and hence serves as a perfect reference for our SRD model in the case where hydrodynamics is switched off Allen and Tildesley (2006). LD also has the benefit of being computationally much more efficient than SRD.

II.1 The polymer model

Polymers are modeled as chains of point-like beads of mass $m_{b}$. Adjacent beads are connected via the Finitely Extensible Nonlinear Elastic (FENE) potential

[TABLE]

where $r$ is the distance between adjacent beads and $K$ and $r_{\rm max}$ are potential parameters describing the strength and maximum distance limit of adjacent beads. Each bead interacts with all other beads via the (shifted and truncated) Lennard-Jones potential

[TABLE]

where $\epsilon$ and $\sigma$ are potential parameters and $r_{ij}$ is the distance between beads $i$ and $j$. The potential is truncated at $r=\sqrt[6]{2}\sigma$ in order to model a good solvent. The potential parameters are chosen as $\sigma=1.0$, $\epsilon=1.0$, $K=30/\sigma^{2}$, and $r_{\rm max}=1.5\sigma$ in reduced units.

II.2 The solvent and polymer dynamics

The polymer is immersed in a solvent modeled by stochastic rotation dynamics (SRD) Malevanets and Kapral (1999, 2004). The SRD method was chosen because it allows taking both hydrodynamics and Brownian motion directly into account in a computationally feasible way. A particular benefit of the method is the possibility to switch off hydrodynamics to better understand its effects. This also allows us to verify the polymer escape in SRD against that in Langevin dynamics.

The SRD solvent consists of point-like particles whose dynamics can be divided into the streaming and the collision steps. In the streaming step the solvent particles are moved ballistically,

[TABLE]

where $t$ is the simulation time, $\Delta t$ is the SRD time step, $\mathbf{r}_{i}$ is the position, and $\mathbf{v}_{i}$ is the velocity of solvent particle $i$. If in this step the solvent particle hits the capsid wall, it is bounced back to the direction of incidence and its velocity is reversed. In other words, the capsid wall constitutes a no-slip boundary for the particle. No-slip boundary conditions ensure that the flow velocity in the surface of a wall is zero Lamura et al. (2001).

In the collision step the simulation space is divided into a grid of cubic cells whose edges are of length $1.0$. The interactions between particles are modeled by rotating the random part of particle velocities within each cell by the equation

[TABLE]

where $\mathbf{v}_{\rm cm}(t)$ is the center-of-mass velocity of the particles in the cell and $\Omega$ is a rotation of angle $\theta$ around a randomly chosen axis. The rotation axis is drawn randomly for each cell each time step. The rotation angle is chosen as $\theta=3\pi/4$. It can be used to adjust the viscosity of the solvent. The solvent is kept at constant temperature of $kT=1.0$ by scaling the random part of particle velocities such that the equipartition theorem holds at each time step Frenkel and Smit (2001). The density of the SRD solvent was chosen such that on average there are 5 particles per unit volume. When hydrodynamics is included $\eta=4.67$ is obtained for the solvent viscosity in reduced units with the chosen parameter values Kikuchi et al. (2003).

The polymer is coupled with the solvent in the collision step where the velocities of polymer beads are updated similarly as those of solvent particles, see Eq. (6). The collisions retain the total momentum and energy within each cell. In order to maintain Galilean invariance, the grid is shifted randomly at each time step Ihle and Kroll (2001). Hydrodynamic interactions can be switched off by randomly permuting the solvent particles’ velocities after each collision step.

The polymer performs molecular dynamics (MD). In the velocity Verlet algorithm used for polymer dynamics, the time step is chosen as $\delta t=0.0002$. The SRD time step $\Delta t=0.5$. A relatively small $\delta t$ was used because with larger time steps the numerical errors accumulate inside the capsid when the polymer is tightly packed. The MD and SRD steps are performed in turns such that after $\Delta t/\delta t=2500$ velocity Verlet steps a single SRD step is performed (including polymer in collision step of Eq. (6)). The mass of the polymer beads $m_{b}=16$ and the mass of the SRD particles $m_{s}=4$.

II.3 The simulation geometry and initial polymer conformations

The simulation geometry is depicted in Fig. 1. A polymer ejects from inside a spherical capsid shell through a narrow pore of radius 0.4. The inside of the capsid is referred to as the cis side and the outside of the capsid is referred to as the trans side. The thickness of the capsid shell is 3. The radius $R_{0}$ of the inner shell of the capsid depends on the chosen initial monomer density $\rho_{0}$ and the initial number of polymer beads $N_{0}$ inside the capsid via

[TABLE]

Notice that in some publications volume fraction $\phi_{0}=4/3\pi\rho_{0}$ is used, instead. Also the beads have often a hard sphere potential. In simulations using molecular dynamics, as in the present study, soft sphere potentials must be used. Hence, the values are not directly comparable. In effect, the largest densities used here supersede the densities used in most earlier studies (see  Piili and Linna (2015)).

The capsid geometry is created using constructive solid geometry technique Wyvill and Kunii (1985), which we have implemented for use with the SRD and LD. In the method intersections with the polymer particles’ trajectories and capsid walls are traced and collisions are handled by slip boundary conditions. As stated before, for the solvent particles no-slip boundary conditions are applied, instead. We use a pore of radius 0.8 for the solvent particles which is twice as wide as that for the polymer beads. The larger pore for the solvent allows for a smoother fluid flow in the pore while the narrow pore for the polymer prevents hairpinning.

The initial conformations are created by injecting polymers inside the capsid through the pore with a large enough packing force within the pore. Force is ramped up until the polymer is packed. Different conformations result from using different initial random generator seeds. Creating an initial conformation also includes thermalizing the polymer by scaling the bead velocities so as to have the polymer reside at temperature $kT=1$. Before ejection, a new SRD solvent is initialized for the created polymer conformation and the polymer is allowed to equilibrate for 2000 time units before the ejection is allowed to start. This way we create an ensemble of random initial conformations. Inevitably, conformations created this way may include knots that have been shown to affect the ejection rate Matthews et al. (2009); Marenduzzo et al. (2013). Identification of knots is beyond the scope of the present study.

To estimate the effect of initial conformations on ejection dynamics we performed simulations where polymers having a bending potential ben with persistence length of $\sim 20$ were packed. This resulted in fundamentally different, spooled conformations. After this we removed the bending potential, which made the chains flexible and released them for ejection. The ejection times measured and averaged over multiple ejections starting from random and spool conformations differed only slightly. This indicates that in our simulations the details of initial conformations do not significantly affect ejection dynamics. It is in place to point out that the packing method, for example whether allowing for intermittent relaxation or not, does have some effect on initial conformations and potentially to ejection dynamics. To our knowledge these effects have not been thoroughly investigated. The packing method we use here falls into the category of generally used methods and consequently in this respect our study is directly comparable to previous studies on polymer ejection.

II.4 Matching the LD and SRD models via friction

Our LD algorithm is implemented as derived by Ermak Ermak and Buckholz (1980). LD is a stochastic method where solvent particles are not explicitly simulated but the polymer resides in a Brownian heath bath satisfying the Langevin equation

[TABLE]

where $\xi$, $\mathbf{p}_{i}(t)$, and $\boldsymbol{\eta}_{i}(t)$ are the friction constant, momentum and random force of the bead $i$, respectively. $\mathbf{f}_{i}(t)$ is the sum of all forces exerted on the bead $i$. $\boldsymbol{\eta}_{i}(t)$ is a zero mean delta correlated Gaussian process, with $\left\langle\boldsymbol{\eta}_{i}(t)\cdot\boldsymbol{\eta}_{i}(t^{\prime})\right\rangle=2\xi kTm_{b}\delta(t-t^{\prime})$.

To make the LD simulations comparable with the SRD simulations, the polymer model, potential parameters, simulation geometry, and simulation temperature were chosen the same. However, LD allowed us to use a larger time step $\delta t=0.001$ than SRD due to it being a more efficient thermostat. Only the friction parameter $\xi$ in LD does not have a direct mapping to the friction in SRD. In order to choose an appropriate value for $\xi$ in the LD capsid ejection model, we performed two straightforward simulations for different values of $\xi$.

Figure 2 shows the measured diffusion constant of a polymer of length $N_{0}=10$ in free solvent. Theoretically, the diffusion constant $D=kT/(N_{0}m_{b}\xi)$ Doi (1996), which is close to the value in our LD simulations. $D$ in LD and SRD were found to coincide when $\xi=1.63$. This is depicted by the dashed line in Fig 2. The double dashed line depicts the value of $\xi=1.58$ that we chose for the capsid ejection simulations. The slightly smaller value for $\xi$ was chosen because it resulted in a better correspondence of the total capsid ejection times $\tau$ from simulations using LD and SRD without hydrodynamics. In other words, for some reason $\tau$ were found to be larger for LD for the same values of friction parameter $\xi$. Even with the choice $\xi=1.58$ the ejection times are consistently larger in LD simulations.

To gain more confidence in the proper mapping of the correspondence between Langevin simulations and SRD without hydrodynamics, we performed simulations where a polymer of length $N_{0}=100$ is dragged by a constant force $f_{\rm drag}$ at the end bead. In the absence of hydrodynamics we expect the terminal velocity to follow

[TABLE]

which can be obtained from Eq. (8) by averaging over a long time, summing over all beads, and assuming that $\frac{{\rm d}\mathbf{p}_{i}}{{\rm d}t}(t)=0$ at terminal velocity. Figure 3 shows $v$ measured for different $f_{\rm drag}$. In these LD simulations $\xi=1.58$. By fitting we obtain the values $\xi_{\rm SRD,\;noHD}=1.65$ and $\xi_{\rm SRD,\;HD}\approx 0.50$ for SRD without and with hydrodynamics, respectively. Observe, however that $\xi_{\rm SRD,\;HD}$ is not well defined for SRD with hydrodynamics since Eq. (9) is not accurate when hydrodynamics is included. Nevertheless, it is a good estimate on the effect the hydrodynamics has on the effective viscosity. The diffusion constant and terminal velocity measurements show that Langevin and SRD without hydrodynamics are in reasonable accordance for these basic systems.

II.5 On polymer lengths and capsid volumes

In this section we comment on relating the length and timescales in the simulations to the corresponding real-world scales. The fairly generic model in the present study is appropriate for validating the SRD method for simulating capsid ejection model and for characterizing the effect of hydrodynamics.

Since we use the fully flexible chain model, the persistence length of the polymer is $\lambda_{p}=\frac{1}{2}b$, where $b$ is the polymer segment length. In our simulations the equilibrium distance between consequent beads is $0.97$, which we take as the (average) segment length. $\lambda_{p}$ is of the order of $4$ nm for ssDNA Tinland et al. (1997) and $50$ nm for dsDNA Manning (2006). Consequently, if we were to model ssDNA, one simulation unit would correspond to about $8$ nm. For dsDNA a simulation unit would correspond to about $100$ nm.

In our simulations, the capsid radii vary from $1.6$ ($N=25,$ $\rho_{0}=1.5$) to $5.1$ ($N=283,$ $\rho_{0}=0.5$). Thus, the persistence length is always an order of magnitude smaller than the capsid radius. For ssDNA the capsid radii would correspond to a range from 12.8 to 40.8 nm. For dsDNA the corresponding range is from 160 to 510 nm. A polymer of length $N_{0}=200$ would correspond to a ssDNA of length 1600 nm having about 4324 bases, since a single base is about 0.37 nm long Rechendorff et al. (2009). For dsDNA a polymer of length $N_{0}=200$ would correspond to a strand of length 20 000 nm with 59 000 base pairs (using 0.34 nm/bp).

If we expect for ssDNA $w\approx 1$ nm, then $\lambda_{p}\approx 8w$. dsDNA width $w$ is about $2$ nm. So, for dsDNA $\lambda_{p}\approx 25w$. In our simulations the repulsive LJ potential has an interaction distance of $\sim 1$, which can be taken as the approximate width of our flexible chain. Hence, $\lambda_{p}\approx 0.5w$, falling short even for the ssDNA that the flexible chain in principle models. For maximum correspondence of coarse-grained polymers modeling real-world polymers in confinements the ratio $\lambda_{p}/w$ has to be adjusted via bending potential. We will look into this more closely in a forthcoming paper.

The estimated volume fraction inside the bacteriophage lambda is

[TABLE]

This is in the same range as volume fractions in our simulations.

III Results

In what follows we refer to the SRD method with hydrodynamics as ’with hydrodynamics’ or ’HD’. The model where SRD is used but without hydrodynamic interactions is referred to as ’without hydrodynamics’ or ’noHD’. The LD method used as a reference does not include hydrodynamic interactions. The presented results are obtained by averaging over typically $50$ runs. For waiting time profiles $500$ ejections were simulated.

III.1 Ejection time

Polymer translocation processes are typically characterized by how the translocation time, here also called the ejection time, $\tau$ depends on the polymer length $N_{0}$. For the case of translocation through a nanometer-scale pore from one semi-infinite space to another it is established that $\tau\sim N_{0}^{\beta}$. Ejection time measurements would suggest that such a scaling relation would describe also polymers’ ejection from capsids. In our previous study using SRD without hydrodynamics we showed that actually $\beta\rightarrow 1$ as $N_{0}$ increases Piili and Linna (2015). As this contradicts the available theoretical treatments for the capsid ejection starting from moderate monomer densities $\rho_{0}$ Muthukumar (2001); Cacciuto and Luijten (2006a); Sakaue and Yoshinaga (2009), a verification using a more established method is called for. We make a close comparison of the polymer ejection models based on SRD and the well established LD. $\tau$ vs $N_{0}$ for the three different models are shown in Figs. 4 (a)-(c). The exponents extracted for the apparent relation $\tau\sim N_{0}^{\beta}$ are given in Table 1. LD and SRD without hydrodynamics are seen to give essentially identical scaling. Hydrodynamics is seen to reduce $\beta$ as has been found also for driven polymer translocation Lehtola et al. (2009). In both cases hydrodynamic interactions reduce the effective friction the polymer experiences outside the pore. Consequently, the effect of the pore friction, largely caused by the geometry, increases when hydrodynamics is included. Increasing this friction local to the pore with respect to the total friction takes the polymer translocation and ejection toward the linear dependence $\tau\sim N_{0}$, which explains the reduction of $\beta$ due to hydrodynamics.

Analogously to the case of driven translocation, where translocation time depends on the driving pore force $f_{d}$ as $\tau\sim f_{d}^{-\alpha}$, the ejection time decreases with increasing initial density as $\tau\sim\rho_{0}^{-\alpha}$, see Figs. 4 (d)-(f) and Table 2. For both LD and SRD $\alpha\to 1$ as $N_{0}$ increases. Inclusion of hydrodynamics decreases $\alpha$, again in analogy with driven translocation Lehtola et al. (2009). This is accounted for by the hydrodynamic interactions decreasing the effective length of the polymer due to increased correlation length along the polymer.

In summary, mere ejection time measurements would seem to confirm previous results on $\tau$ scaling with polymer length $N_{0}$. Also, the dependence of $\tau$ on the initial monomer density $\rho_{0}$ would seem to corroborate the scaling behavior. Results using SRD and LD are essentially identical. The theoretical arguments have typically been corroborated by ejection time measurements alone. However, inspection of the measured waiting time profiles changes the conclusions completely.

III.2 Waiting times

In this section we extract the waiting time profiles $t_{w}(s)$ for the different models. We obtain a more precise form for $t_{w}(s)$ than in our previous study Piili and Linna (2015) and find that it is universal for all models.

$t_{w}(s)$ is defined as the time it takes for the bead $s$ to eject the capsid after the ejection of the previous bead $s-1$.

[TABLE]

where $t(s)$ is the time when the bead $s$ exits the capsid for the last time, that is, the cumulative waiting time. $t_{w}(s)$ plotted as a function of $s\in[1,N_{0}-1]$ is the waiting time profile. For the process to be genuinely scale-invariant with respect to the polymer length, both $t(s)$ and $t_{w}(s)$ should scale with $s$. Figures 5 (a)-(c) show $t(s)$ obtained for the three models. For all the models the endpoints $t(N_{0}-1)$ scale with $N_{0}$ in accordance with the apparent scaling relation $\tau\sim N_{0}^{\beta}$. However, $t(s)$ do not scale with $s$.

Figures 5 (d)-(f) show the waiting time profiles $t_{w}(s)$ for the three different models. $t_{w}(s)$ is seen to be of the common form

[TABLE]

where $A=45,\ 42,$ and $32$ for LD, SRD without hydrodynamics, and SRD with hydrodynamics, respectively. Hence, for all the models polymer ejection slows down exponentially with the length of the ejected segment $s$. At a definite stage when approximately $63\%$ of the polymer has been ejected the ejection slows down more strongly with $s$. Presumably, the transition corresponds to $s=s_{0}$ when the monomer density inside the capsid is so low that internal pressure no longer exerts force on the ejecting polymer, see Section III.3. The exponential form $t_{w}\sim\exp(Cs)$ can only lead to linear dependence of the ejection time with $N_{0}$ for long polymers Piili and Linna (2015). This is true also for the sum of two exponential functions as in Eq. (12). Consequently, for sufficiently long polymers there is a scaling function $h$ such that $t(s)=N_{0}h(s/N_{0})$. In the present case

[TABLE]

It is seen that hydrodynamics only reduces the magnitude of waiting time profile $t_{w}(s)$ without changing its form. In other words, the waiting times are related via $t_{w}^{\rm HD}(s)=(\tau_{\rm HD}/\tau_{\rm noHD})t_{w}^{\rm noHD}(s)$, for the $N_{0}$ and $\rho_{0}$ selected. This is reminiscent of the driven polymer translocation where hydrodynamics speeds up translocation and scales down the length of the tensed segment $l$ on the cis side without changing the way the tension spreads on the polymer chain, that is, the form of $l(s)$ Moisio et al. (2016). Furthermore, $t_{w}(s)$ obtained by using LD aligns almost perfectly to that given by SRD without hydrodynamics within the precision of the mapping of the two models via the friction parameter, see Section II.4. Polymer ejects slightly faster in SRD without hydrodynamics than in LD simulations even though LD has a slightly smaller friction parameter $\xi$ in free solvent, as measured in Section II.4. This would indicate that for high monomer concentrations SRD has enhanced correlations between polymer beads residing in the same cell which would decrease the effective friction.

III.3 Force measured at the pore

Thanks to the implementation of the capsid geometry using constructive solid geometry no explicit forces at the pore are imposed. Hence, the ejection force results as far as possible from the pressure of the polymer segment confined inside the capsid. There remains contribution from local effects, such as the sharp edges of the pore restricting polymer movement at both openings, but these can be considered relevant also to real pores.

We characterize the dynamic force at the pore by measuring the force $f$ for fully and partly equilibrated polymer conformations for different reaction coordinate $s$. In SRD polymers of different lengths $N_{0}$ are packed to a random conformation inside a capsid of inner volume $V=4/3\pi R_{0}^{3}=N_{0}/\rho_{0}$ until the bead $s$ is at the pore entrance. Here $\rho_{0}$ is the initial monomer density of a corresponding capsid ejection.

In our previous study Piili and Linna (2015) the bead $s$ was attached to a point in the pore entrance via a FENE potential and the average force needed to keep the polymer in place was measured, see Fig. 6 (a). In response to the results from more precise measurements using LD, reported in what follows, we change here the measurement point to the middle of the pore, see Fig. 6 (b). In SRD, after attaching the bead at either the entrance or the middle of the pore we wait for a time $t_{\rm eq}=2.2\cdot 10^{4}$ before measuring the force over $t_{m}=2\cdot 10^{4}$ time steps, one measurement per step, and averaging over them. For a few $s$ we checked that setting the equilibration time $t_{\rm eq}=5\cdot 10^{4}$ did not change the measured average $f$. Hence, during the measurement the polymer is at or very close to an equilibrium conformation.

The measured force for equilibrated polymer conformations are not equal to the dynamic force during ejection. In order to be able to measure $f$ for conformations that are not fully equilibrated we do the measurement in a slightly more complicated way in the LD model. Here the polymer is initially packed inside the capsid and then freed for a single bead to eject. The appropriate bead is then pulled either to the pore opening or the middle of the pore and held fixed for a time $t_{\rm eq}$ after which the harmonic force that is needed to keep the bead fixed is measured for a time $t_{m}=t_{\rm eq}$. After this the polymer is again freed for a single bead to eject. This way force is measured for all $s$.

First we verify that the same equilibrium $f$ is obtained using SRD and LD. $f(s)$ for fully equilibrated conformations measured in the SRD model and in the LD model using $t_{\rm eq}=t_{m}=2000$ are shown in Fig. 7 (a). $f(s)$ for both models are seen to be identical. Force measured at the pore entrance $f_{\rm ent}(s)$ is seen to be larger than force measured in the middle of the pore $f_{\rm mid}(s)$; it does not decay to zero even at the end of the ejection. This is caused by the reduction of the degrees of freedom at the pore entrance due to the measured bead being held stationary there. This creates a bias toward the exit of the pore. The surprisingly strong bias created by an asymmetry in the pore was noticed already in Linna et al. (2014).

Eliminating the bias due to asymmetry imposed by the measurement affects the dependence of $f$ on $s$. Figure 7 (b) shows $f_{\rm mid}(s)$ for polymers of different lengths in the LD model together with $f_{\rm ent}(s)$ for $N_{0}=400$. Measurements are done for $t_{\rm eq}=t_{m}=2000$, so these $f$ are for conformations that are close to equilibrium. For the initially large monomer density $\rho$, that is, for relatively small $s$, all $f$ show close to exponential decay with $s$. This decay rate is very close to the rate of the exponential increase of the waiting time $t_{w}$ with $s$ for $s\leq 0.625N_{0}$, see the first exponential term in Eq. (12). The lines in Fig. 7 (b) show the function $f\propto\exp(-(2.8/N_{0})s)$ for the different $N_{0}$. For $s\leq 0.625N_{0}$ the ejection dynamics is thus seen to mainly result from the exponential decay of the pressure inside the capsid. Finally, as the pressure decreases more abruptly, there is a crossover to a stronger exponential increase of $t_{w}$, see Eq. (12).

The exponential decay of $f$ coinciding with the exponential increase of $t_{w}$ for $s\leq 0.625N_{0}$, that is, for the part of the process where the pressure resulting from packed monomers drives the ejecting polymer, is in keeping with our previous findings. For these densities and polymer lengths the number of beads per blob is very low, so the blob picture used in many theoretical approaches is not relevant but monomers interact individually. The exponential dependence of the resulting potential inside the capsid and hence the driving force has an exponential dependence on monomer density and so $s$. This can be obtained from the Flory-Huggins theory mixing free energy in the limit of high density Piili and Linna (2015). For large monomer densities and small $s$ we also find that $t_{w}(s)\sim 1/f(s)$. Observe that the small deviation from inverse proportionality to a power law relation $t_{w}(s)\sim 1/f(s)^{\gamma}$, where $\gamma\approx 0.95$, leads to the slightly deviating exponent as a form of $f(s)\sim\exp{\left(-\frac{2.8}{\gamma N_{0}}s\right)}$. This explains the slightly different exponent in the force curves in Fig. 7 (b).

The waiting time shows a stronger exponential for lower monomer densities and large $s$, see Eq. (12). Fig. 7 (c) shows force measured at the pore midpoint using LD and allowing for the polymers to equilibrate for different times $t_{\rm eq}=t_{\rm m}$, as explained above. It is seen that the further the polymer conformation is from the equilibrium, the more abruptly the measured $f$ falls off with increasing $s$ similarly to the corresponding stronger increase of $t_{w}$ with $s$ for $s>0.625N_{0}$. There are two potential reasons for this: First, on the cis side at the final stage tension may propagate in the remaining polymer segment, which increases friction and diminishes force measured at the pore. Second, on the trans side monomers may crowd thus possibly impeding the ejection of the polymer. We have shown that crowding plays no role in driven polymer translocation Suhonen et al. (2014). The force in the final ejection stage is much smaller than in the driven translocation, so here the effect cannot be ruled out offhand. The relevant non-equilibrium mechanism is determined in Section III.6.

The greater statistics due our LD model enables us to extract the dependence of the equilibrium force measured in the middle of the pore on the monomer density inside the capsid more precisely than using SRD. We find that the apparent dependence on the polymer length arises from the inevitable overlap of the repulsive monomer potentials with the capsid wall. This overlap is proportionally larger for small capsids and short polymers. We introduce the volume correction parameter $\epsilon$ to take this into account. The effective monomer density then becomes

[TABLE]

where $N$ is the number of monomers inside the capsid. The data for the measured $f$ vs $\rho_{e}$ for different $N_{0}$ falls onto the same curve when $\epsilon=0.3$.

Figure 8 shows the measured pore force in equilibrium as a function of the effective density inside the capsid $\rho_{e}$ in natural and logarithmic scales. For large and intermediate monomer densities the force grows exponentially with $\rho_{e}$ in keeping with the exponential decay of $f$ with $s$. However, for sufficiently small monomer densities the exponential relation does not hold precisely. This is due to the measured force decaying to zero and even going slightly negative for extremely small monomer densities.

In Fig. 8 we have fitted to the data the function $f(s)=C\left(\exp{\left(B\rho_{e}\right)}-1\right)$ which, unlike a pure exponential, has the property that it decays to zero when the monomer density is exactly zero. This function was chosen because we believe that the negative force is an artefact caused by the local pore geometry and the selected position for the force measurement. In measurements at the cis entrance the force is always positive while in the mid-pore measurement the minimum force is negative. This implies that there must be a measurement position where the minimum force is exactly zero. The relatively small offset in the force might seem like a minor detail, but as an ejecting polymer spends the majority of its time in the small force regime, the form of the force has a considerable effect on ejection time.

The form of the pore force leaves some room for speculation. For instance, it is also possible to describe the force in the small $\rho_{e}$ regime as a power law $f\sim\rho_{e}^{2}$ due to the shifted exponential and power law resembling each other in such a short range. Cacciuto and Luijten Cacciuto and Luijten (2006b) measured the scaling of the excess free energy with the number of monomers in the capsid as $\Delta F\sim N^{2.97}$ for $0.2\leq\phi<0.35$, where $\phi$ is the volume fraction. This would give $f=-\partial\Delta F/\partial N\sim N^{1.97}$, which is approximately the scaling obtained here. In Fig. 8 (b) we have plotted the scaling of this form on logarithmic scale alongside with force measurements. However, the scaling regime obtained here is for lower $\rho$ and cannot be of the same origin as assumed in Cacciuto and Luijten (2006b), namely the screening that steps in at higher $\rho$.

The main observation here is that in our model we do find a narrow interval at low $\rho$, where $f$ for equilibrium conformations may scale with $\rho$. However, the scaling exponent is not of the same magnitude that could be derived using the blob-scaling arguments. In addition, during the ejection the conformations are out of equilibrium.

III.4 First passage times derived from the pore force

In Ref. Piili and Linna (2015) we derived an analytical estimate for the ejection time of the polymer under the assumption of a purely exponential pore force. However, the presented treatment is not accurate at the very final stage of the ejection if the force is assumed to decay to zero in the end as presented in Fig. 8 (a). This is because the ejection stalls completely if the driving force vanishes. Thus, the ejection cannot complete without the help of diffusion. Also, as the ejection slows down considerably in the end, the final stage of ejection has a major effect on the total ejection time. We can estimate the first passage times $t_{\rm f}(s)$ based on the pore force solely by using the formula Gardiner (2009)

[TABLE]

where

[TABLE]

is the pore potential derived from the shifted exponential pore force $f(s)$ presented in Fig. 8 (a). Here, we use $\rho$ instead of the effective density for simplicity. The integral is not analytically solvable, but by numerical computation we obtain a good correspondence with simulations in the regime $N_{0}\leq 400$ as shown in Fig. 9 (a). The parameters $\xi m=254.5$ and $kT=5.23$ in Eq. (15) were chosen such that a good correspondence is obtained with all the presented curves, with the emphasis on the endpoints.

The parameter $\xi m$ describes the friction of both the solvent and of the pore while $kT$ describes the friction of the solvent multiplied by the temperature. Hence, we cannot expect to obtain a direct mapping of the parameters from the Langevin equation. The force parameter values $C=0.411$ and $B=2.70$ were obtained from a fit to the force measured for a polymer of length $N_{0}=200$ at the initial density $\rho_{0}=1.0$.

Figure 9 (b) shows the first passage time curves obtained from the numerical integration of Eq. (15) in the logarithmic scale. The endpoints of the curves reveal that the apparent scaling of $\tau$ vs. $N_{0}$ tends towards linear when $N_{0}$ grows extremely large. Therefore we can safely assume that the apparent scaling observed for $N_{0}\leq 400$ indeed is a finite size effect even for a strictly vanishing pore force.

When numerically integrating Eq. (15) we observed that the obtained first passage times are extremely sensitive to the choice of the force. A force that decays to negative values for very small densities leads to exponentially growing ejection times. If the force remains above zero throughout the ejection, including the endpoint $s/N_{0}=1$, the formula leads to linear scaling already with very small $N_{0}$. If in the end of the ejection there is a finitely long regime in $s/N_{0}$ where the force is exactly zero, this diffusive region gives a scaling with the exponent 2. The presented model where force decays to zero exactly at the end of the ejection presents a special case where it is not evident beforehand what happens to the scaling. It also seems to best describe the simulation results received and explain why we obtain such a nice apparent scaling while everything seems to imply linear scaling for very long polymers.

III.5 The radius of gyration

We have previously found that the radius of gyration of the polymer segment on the trans side grows roughly as $R_{g}\sim s^{0.6}$, which means that the ejecting polymer segment is not far from equilibrium Piili and Linna (2015), unlike in the case of driven translocation Lehtola et al. (2009). Here, using SRD we measure $R_{g}(s)$ for the ejected polymer segment with and without hydrodynamics and $R_{g,\rm eq}$ for polymers of different lengths $N=s$ at equilibrium, see Fig. 10. It is seen that throughout the ejection the ejected polymer segment is slightly compressed compared to the equilibrium conformation, that is, $R_{g}<R_{g,\rm eq}$. Hence, the ejected polymer segment remains slightly out of equilibrium throughout the ejection. $R_{g}$ is larger and closer to $R_{g,\rm eq}$ when hydrodynamics is included. So, hydrodynamics speeds up the relaxation of a polymer even more than it speeds up the ejection. The same observation was made concerning driven polymer translocation Moisio et al. (2016).

The radius of gyration of the trans side polymer segment manifests clearly one characteristic of the polymer ejection that makes it impossible for the ejection waiting time to scale with the reaction coordinate $s$. $R_{g}$ is seen to be smaller for long than short polymers. In other words, long polymers are driven farther out of equilibrium than short polymers starting from the same initial monomer density. The monomer density $\rho$ inside the capsid decays faster with increasing $s$ for a short polymer than for a long polymer. For example, for a polymer of $N_{0}=50$ $\rho=0$ at $s=50$, but not so for a polymer of $N_{0}=100$. Instead, $\rho$ decays with $s/N_{0}$ identically for all $N_{0}$. So, since polymers of different lengths are driven by different force at a same $s$, there can be no universal scaling of the waiting time $t_{w}$ with $s$. For the same reason, the longer the polymer, the farther out of equilibrium it is for any given $s>0$.

III.6 Modified models

From the above observations it is evident that the polymer ejection is inherently a non-equilibrium process whose non-equilibrium characteristics are more enhanced for long polymers. As stated in Section III.3, the two non-equilibrium mechanisms potentially affecting the ejection dynamics are tension propagation on the cis side and crowding on the trans side. In the case of driven polymer translocation we have shown that crowding has no effect, but the translocation dynamics is in practice determined by the dynamics of tension propagation Suhonen et al. (2014). For a worm-like-chain in capsid ejection it has been shown that crowding slows down the ejection Al Lawati et al. (2013). This result, however, cannot be generalized to the flexible chain polymer model.

In order to determine the dominating non-equilibrium effect in the polymer ejection from a capsid we remove beads from the trans side during ejection. We simulate ejection of polymers of $N_{0}=200$ such that the maximum number of beads on the trans side at any moment is $N_{\rm trans}=5$, $10$, or $160$. For $N_{\rm trans}=160$ the waiting time profile $t_{w}(s)$ is identical to that of the full polymer ejection, where no beads are removed. $t_{w}(s)$ for $N_{\rm trans}=5$ and $10$ deviate from $t_{w}(s)$ of full polymer, deviation being greater for $N_{\rm trans}=5$. Fig. 11 shows $t_{w}(s)$ for the full polymer and for $N_{\rm trans}=5$. $t_{w}(s)$ for $N_{\rm trans}=5$ is seen to be clearly smaller than $t_{w}(s)$ for the full polymer leading to the ejection time of 18% smaller with $N_{\rm trans}=5$ than with the full polymer. This can only be due to crowding not slowing down the ejection for $N_{\rm trans}=5$. This is in contrast to driven translation where the effect of crowding is far less significant than the effect the tension propagation has on the dynamics. In translocation with $N_{0}=200$, the ejection time with $N_{\rm trans}=5$ is 7.9% smaller with $f_{d}=2$ and 4.6% smaller with $f_{d}=8$ than ejection times with full polymer. In the ejection process tension can significantly propagate only at the final stages of the process when the monomers are less densely packed. Also the force driving the polymer has decreased at this point, so tension propagation will be mild. Under these circumstances crowding, although being weaker than in the driven polymer translocation, shows more prominently in the resulting polymer ejection dynamics.

IV Conclusion

We have studied the ejection of a fully flexible polymer chain from a spherical capsid modeled using simplistic boundary conditions. Hydrodynamic interactions were simulated using stochastic rotation dynamics (SRD) coupled with molecular dynamics. We also used SRD with hydrodynamic interactions switched off to better pin down effects due to hydrodynamics. The less established SRD method was compared to the thoroughly understood Langevin dynamics (LD) to verify its suitability for simulations of this kind. Results obtained using SRD without hydrodynamics were found to be in good agreement with those obtained using the LD model whose friction was matched. Computational efficiency of LD allowed us to perform more precise force measurements than before.

From measured ejection times $\tau$ for different polymer lengths $N_{0}$ the apparent scaling $\tau\sim N_{0}^{\beta}$ was obtained for $N_{0}\leq 400$. The differences in the fitted effective exponents $\beta$ obtained for different models make sense. Included hydrodynamic interactions not only reduce the ejection time but also reduce the $\beta$ exponent which is in line with results of forced translocation Moisio et al. (2016); Lehtola et al. (2009). This was addressed to hydrodynamic correlations reducing friction proportionately less in the pore region than outside it, which contributes towards linear scaling. We showed, however, that basing the analysis solely on ejection times, as in many previous studies, leads to an incorrect characterization of the process.

The waiting time $t_{w}(s)$ proved most valuable for understanding the polymer ejection dynamics. It describes how long it takes for the individual polymer bead $s$ to permanently exit the capsid after the final exit of the bead $s-1$. In the previous study we concluded that the waiting time $t_{w}(s)$ is of exponential form Piili and Linna (2015). The more precise measurements and analysis conducted here reveal that $t_{w}(s)$ is actually more accurately described by a sum of two exponentials. After about 63% of the polymer has ejected, the ejection starts to slow down even more considerably and the second exponential starts to dominate the waiting time. When hydrodynamics is included $t_{w}(s)$ is of almost identical form as without hydrodynamics, only by a constant factor smaller in magnitude.

Remarkably, waiting times $t_{w}(s)$ for polymers of different lengths $N_{0}$ starting their ejection at the same initial monomer density $\rho_{0}$ collapse when plotted as a function of the normalized reaction coordinate $s/N_{0}$. This implies, by definition, that the ejection times should scale linearly with $N_{0}$. Only the final retraction of $t_{w}$ makes the scaling superlinear for polymers of moderate length. However, the effect of the final retraction becomes proportionately smaller for longer polymers, ultimately vanishing for very long polymers. Hence, linear scaling is approached for extensively long polymers and the cumulative waiting time is given by $t(s)=N_{0}h(s/N_{0})$, where $h(s/N_{0})$ is a scaling function.

Our LD model allowed us to investigate the force exerted at the pore more carefully than in our previous study, due to its computational efficiency. We found that the position from where the force is measured has a large effect on the measurements. The force measured at the pore entrance is very accurately an exponential function of the monomer density inside the capsid. On the other hand, the force measured at the middle of the pore is exponential only for densities exceeding $\sim 0.4$. To capture the form of the force more precisely we need to reduce a constant factor from the exponential. This takes the force to zero for zero monomer density. This was observed to significantly affect the ejection times. We also found out that the force is a function of the effective density, where the effective radius of the inside of the capsid is $\sim 0.3$ larger than the real radius used in the simulations. This is due to the repulsive LJ potential overlapping with the capsid walls. This effect diminishes as the capsid volume increases.

By inserting the measured exponential form for the force, from which the offset term is reduced, into the first passage time formula for a random walk we are able to quite accurately reproduce the obtained first passage times from simulations. The numerical integration of the formula also reveals that this type of force leads to linear scaling for extremely large $N_{0}$.

The radius of gyration $R_{g}$ outside the capsid shows that the part of the polymer on the trans side is more compact than the corresponding polymer in equilibrium. When hydrodynamics is included the trans side is larger and therefore closer to equilibrium $R_{g}$ than without hydrodynamics. This occurs even though the faster ejection with hydrodynamics allows the polymer conformation less time to expand. In a modified model, where the number of monomers on the trans side are kept constant by continually removing them, a polymer ejects slightly faster than the corresponding full polymer. Hence, crowding has a small effect on ejection dynamics for the flexible polymer corroborating the findings from the measurements of $R_{g}$. Finally, hydrodynamics has a fairly weak effect on capsid ejection, including the crowding. Most importantly, it does not alter the universal form of the waiting time versus the number of ejected monomers. All evidence from our measurements goes to show that the apparent superlinear scaling of the ejection time with the polymer length $\tau\sim N_{0}^{\beta}$ tends to linear scaling for extremely long polymers.

Acknowledgements.

The computational resources of CSC-IT Center for Science, Finland, and Aalto Science-IT project are acknowledged. The work of Joonas Piili was supported by The Emil Aaltonen Foundation and Finnish Foundation for Technology Promotion. The work of Pauli Suhonen is supported by The Emil Aaltonen Foundation.