Spontaneous Domain Formation in Spherically-Confined Elastic Filaments

TL;DR

This study reveals that elastic filaments confined in spherical spaces naturally form multiple ordered domains, with morphologies depending on density, challenging the assumption of a single spool configuration.

Contribution

It demonstrates through analytical and simulation methods that the ground state of confined elastic filaments consists of multiple domains rather than a single spool.

Findings

Multiple homogeneously-ordered domains form instead of a single spool.

Different morphologies, including concentric spools and topological links, emerge at varying densities.

Results apply to viral DNA, metallic wires, and flexible polymers.

Abstract

Although the free energy of a genome packing into a virus is dominated by DNA-DNA interactions, ordering of the DNA inside the capsid is elasticity-driven, suggesting general solutions with DNA organized into spool-like domains. Using analytical calculations and computer simulations of a long elastic filament confined to a spherical container, we show that the ground state is not a single spool as assumed hitherto, but an ordering mosaic of multiple homogeneously-ordered domains. At low densities, we observe concentric spools, while at higher densities, other morphologies emerge, which resemble topological links. We discuss our results in the context of metallic wires, viral DNA, and flexible polymers.

Equations

$$E_{el}=\frac{2K_{c}}{\pi\sigma^{2}}{\int_{V}}t~{}\Big{(}\mathbf{\hat{t}}\times(\nabla\times\mathbf{\hat{t}})\Big{)}^{2}~{}d{\bf r}\;,$$

$$E_{el}^{1s}(\nu,R_{0})=\frac{\pi^{2}KR_{0}}{\sqrt{3}}\log\Biggl{[}-{\nu}^{1/3}+\textstyle{\frac{1}{2}}\ln\frac{1+\nu^{1/3}}{1-\nu^{1/3}}\Biggr{]}\;.$$

$$E_{el}^{2\mathrm{s}}({\nu }_1,{\nu }_2,R_0)=E_{el}^{1\mathrm{s}}({\nu }_1,R_0)+E_{el}^{1\mathrm{s}}({\nu }_2^{\mathrm{t}},R_1).$$

Full text

Introduction. Spatial organization of constrained elastic filaments underlies a range of packing problems from macroscopic wires in spherical cavities Stoop et al. (2011); Shaebani et al. (2017) to genomes in bacteriophages Knobler and Gelbart (2009); Keller et al. (2017). A theoretical interpretation of single molecule experiments on DNA packaging in viruses Purohit et al. (2003) suggests single-domain spool-like structures with axial symmetry as a prevailing morphology. However, recent experiments report a range of different structures with inhomogeneously-ordered genomes, defects, and phase transitions between the packings. Cryo-electron microscopy of bacteriophages Leforestier and Livolant (2010); Leforestier (2013) at various physicochemical conditions indicates a range of possible morphologies of DNA within the capsid: from isotropically disordered, cholesteric liquid crystalline ordering with a disordered core, to uniform, concentric spools with local hexagonal order, reminiscent of the nested spools first proposed in Hall and Schellman (1982). Measurements of DNA dynamics during packaging Berndsen et al. (2014), and intermittent ejection dynamics Chiaruttini et al. (2010) suggest multi-domain structures. Similarly, experiments and simulations of quasi-one-dimensional elastoplastic filaments forced into spherical confinement Stoop et al. (2011); Shaebani et al. (2017), as well as the theoretical description of elastic filaments within the continuum polymer nematic theory Shin and Grason (2011); Svenšek et al. (2010), indicate a multidomain structures with a mixture of local disorder and nematic order.

Coarse-grained simulations of confined semiflexible polymers mimicking DNA in viral shells have been performed at various levels of detail. DNA interactions and entropic contributions need to be considered to explain the energetics of confinement  Kindt et al. (2001); Chen et al. (2018), but the conformation of the confined DNA is dictated by elasticity Petrov and Harvey (2008). Packing into icosahedral capsids was shown to have a strong stochastic component Forrey and Muthukumar (2006), and structures without spool-like symmetry have been reported, including multi-domain spool packing, observed when DNA is pushed into the capsid through a portal  Rapaport (2016). The effect of elastic kinks on packing has been discussed in Myers and Pettitt (2017), while local liquid-crystalline order of compacted DNA was found to play a crucial role promoting the formation of knots Marenduzzo et al. (2009).

The slow dynamics of stiff chains in confinement implies that the configurations observed in the non-equilibrium simulations crucially depend on the applied protocol, e.g., on how the external force is applied in packaging-motor driven DNA encapsidation in bacteriophages Chen et al. (2018). It is therefore not surprising that, even in cases where the models are very similar, a wide variety of behaviours is reported. Here, we focus on the nature of the ground-state packing configurations. While the present analysis thus answers a different question from previous works, its results should nevertheless provide a well-defined background scenario with which non-equilibrium configurations can be compared and assessed.

Theoretical Model. We consider a spherical enclosure $\cal{C}$ with radius $R_{0}$ (volume $V_{0}$), and a hard semi-flexible cylinder with cross-section $\sigma$, length $L$, and a persistence length $l_{p}=K_{\rm c}/k_{\rm B}T$, where $K_{\rm c}$ is the bending rigidity. The chain cannot overlap with itself nor with the confining wall and it is torsionally relaxed; its elastic energy is described by the integral $E_{el}={\textstyle\frac{1}{2}}K_{\rm c}\int_{0}^{L}ds/R^{2}(s)$, with $R(s)$ the local radius of curvature.

We explore the limit of a rigid, long, thin chain, where the elastic energy dominates over the configurational entropy. In this limit, the chain can be described by the continuum theory of polymer nematics Svenšek et al. (2010); Shin and Grason (2011), where the polymer is described by a nematic-order vector field or local polymer “current” $\mathbf{t}(\mathbf{r})$ Svenšek et al. (2010) with $\mathbf{\hat{t}}=\mathbf{t}/t$ the unit tangent vector (director) to the chain. The total amount of material is $V_{p}=\int_{\cal{C}}t~{}d\mathbf{r}$, with $t\equiv|\mathbf{t}|$ the local polymer density. Due to the excluded volume of the chain, $t$ is bounded by the densest (hexagonal) packing of rigid disks in two dimensions: ($t^{m}=\frac{\pi\sqrt{3}}{6}\approx 0.907$) Podgornik et al. (2016). The ratio between the actual amount of polymer in the capsid and its theoretical maximum, $\nu\equiv V_{p}/(t_{m}V_{0})=\frac{3\sqrt{3}\sigma^{2}}{8\pi R_{0}^{3}}L$, is a key dimensionless parameter in the model.

The polymer current is only non-zero inside the capsid and has to be solenoidal ($\nabla\cdot\mathbf{t}=0$) except at the polymer end-points (the contribution of which vanishes in the long-chain limit). The lowest non-trivial order in the expansion of the elastic free energy of the divergence-free polymer current is given by Svenšek et al. (2010); Shin and Grason (2011)

[TABLE]

where twist deformations ($\mathbf{\hat{t}}\cdot(\nabla\times\mathbf{\hat{t}})$) are by assumption excluded. The problem is to find the solenoidal field $\mathbf{t}({\bf r})$ within $\cal{C}$ which minimizes $E_{el}$.

Assuming axial symmetry, the solution is an spool Podgornik et al. (2016). However, relaxing the symmetry constraints might lead to more complex structures such as those seen in experiments Leforestier (2013); Berndsen et al. (2014), theory Shin and Grason (2011); Svenšek et al. (2010), and simulations Myers and Pettitt (2017); Forrey and Muthukumar (2006); Petrov and Harvey (2008). Indeed, as we show later, a single-domain spool is almost never the lowest elastic energy configuration. For solutions containing more than one domain, $\mathbf{t}(\mathbf{r})$ is not continuous in space, but it can still be solenoidal if at discontinuity interfaces ${\cal S}$, the local field on both sides is perpendicular to the normal, $\mathbf{t}\cdot\mathbf{n}({\cal S})=0$. Multiple-domain structures should therefore be considered. Their advantage is flexibility in arranging the material so as to avoid highly curved regions close to the capsid center, which cannot be done in single spool structures.

With this in mind, we assume that the field is partitioned into $n_{d}$ domains, $\nu=\sum_{1}^{n_{d}}{\nu_{i}}$, with relative sizes $w_{i}\equiv\nu_{i}/\nu$. For $n_{d}\leq 3$, we construct candidate structures and search for the optimum partitioning $\{m_{i}\}$ at several packing fractions $0<\nu<1$. Inspired by observations Stoop et al. (2011); Shaebani et al. (2017); Svensek and Podgornik (2012); Rapaport (2016); Forrey and Muthukumar (2006); Petrov and Harvey (2008), we consider four classes of trial structures depicted in Fig. 1: single spools (a), double spools (b) comprising two nested tori with perpendicular axes, triply nested spools (c) with perpendicular axes, and—deviating from the nested-spool paradigm—Hopf links (d), with a symmetric partitioning into two equivalent, perpendicular, topologically-linked tori. We assume that the field density is $t_{m}$ within the defined domains and zero elsewhere and approximated (1) by a sum of integrals over independent domains. We neglect the short segments connecting the domains, arguing that this contribution is negligible in the thin, long chain limit.

For a single spool, $\mathbf{\hat{t}}^{1s}=t_{m}\left({-y},{x},0\right)/\sqrt{x^{2}+y^{2}}$, we have Purohit et al. (2003):

[TABLE]

The energy diverges for $\nu\to 1$, due to diverging curvature in the center of the capsid. The Hopf link comprises two perpendicular spools with a partitioning dictated by symmetry ($w_{1}=0.5$) and $E_{el}^{\rm hl}\left(\nu,R_{0}\right)\approx 2E_{el}^{1\rm s}\left(\nu/2,R_{0}\right)$. In the double spool, the outer domain ($\nu_{1}$) is the same as above, and the inner spool ($\nu_{2}$) with radius $R_{1}=R_{0}\sqrt{1-\nu_{1}^{2/3}}$ is composed of two small hemi-toroidal caps of combined length proportional to $\nu_{2}^{\rm t}$, connected by straight sections ($\nu_{2}^{\rm s}$): $\nu_{2}=\nu_{2}^{\rm t}+\nu_{2}^{\rm s}$. The elastic energy is localized in the hemispheres and is determined by $\nu_{2}^{\rm t}$ alone (see section SI(B)):

[TABLE]

Interestingly, at a given $\nu$, the ratio $E_{el}^{2\rm s}/E_{el}^{1\rm s}$ is completely determined by the partitioning of the material $w_{1}$. We calculate the optimal $w_{1}$ for several values of $\nu$ by minimizing (3). The calculation for triple spools is analogous: $\nu=\nu_{1}+\nu_{2}+\nu_{3}$, and the elastic energy a sum over two outer spools (3) and an innermost spool with an elaborate shape (see section SI(C)): $E_{el}^{3\rm s}\left(\nu_{1},\nu_{2},\nu_{3},R_{0}\right)=E_{el}^{2\rm s}\left(\nu_{1},\nu_{2},R_{0}\right)+E_{el}^{3^{\rm rd}\rm spool}\left(\nu_{1},\nu_{2},\nu_{3},R_{0}\right)$. The optimal structure is obtained by minimizing over $w_{1}$ and $w_{2}$.

Fig. 1e) shows the ratio of the elastic energies as a function of $\nu$ for optimal double- and triple-spool configurations. We note that both are always lower in energy than a single spool, and that the triple spool wins over the double spool. Adding a fourth spool in an analogous manner does not reduce the energy further (see online SI for details).

Our best Hopf link structure (Fig. S4) can have lower elastic energy than the single spool only at extremely low $\nu$ and even then by a tiny margin, but never wins over double or triple spools.

The robust conclusion within our analytical model is that in the long, thin, and hard cylinder limit, multiple domains are favored over a single-domain. In principle, this may differ if we consider spools with non-orthogonal axes, or any other, better arrangements not included in our theoretical considerations. Moving away from this limit (i.e., finite $\sigma$, repulsive interactions other than excluded volume), the volume accessible to the polymer will be modified, which can result in a domain-wall energy penalty, disfavoring multi-domain arrangements. At finite temperature, entropic terms due to conformational fluctuations of the chain also become important. We will address all of these issues with computer simulations.

Simulation Model. We consider a long elastic filament of diameter $\sigma$ confined to a sphere with radius $R_{0}=5\sigma$. We model the elastic filament as a sequence of $N$ soft, repulsive beads interacting with each other and the confining wall via the Weeks–Chandler–Anderson (WCA) repulsion Weeks et al. (1971). We use standard LJ units, expressing energies in $\epsilon$ and lengths in $\sigma$. The beads are bonded by a two-body stretching term ${V_{b}=16\,\Gamma\left(r-1\right)^{2}}$ and a three-body bending term ${V_{\theta}=\Gamma\left(1+\cos{\theta}\right)}$, where $r$ is the distance between two consecutive beads, and $\theta$ the angle between consecutive bond vectors. The ratio 16 follows the exact result for elastic rods Podgornik et al. (2000), and the prefactor $\Gamma$ determines the chain rigidity, $\Gamma=T^{*}l^{*}_{p}$, with $l^{*}_{p}$ the dimensionless persistence length and $T^{*}\equiv k_{B}T/\epsilon$ the reduced simulation temperature. Changing $T^{*}$ or $l^{*}_{p}$ has the same effect on the elastic energy; therefore, by fixing $\Gamma$, we can simulate polymers with different rigidity by changing the simulation temperature 111This is not entirely true, since the LJ term does not scale with $\Gamma$. However, in the rigid limit the effect of repulsive interactions on the persistence length is negligible.. In our simulations, we fix $\Gamma=1000$, so at low temperatures we are modeling fairly rigid chains (at $T^{*}=1$, $l_{p}/\sigma=1000$) and at $T^{*}\approx 40$ ($l_{p}/\sigma\approx 25$) the parameters correspond to DNA.

We performed $NVT$ replica-exchange molecular dynamics simulations Hansmann (1997); Sugita and Okamoto (1999) of chains of up to $N=250$ beads using the LAMMPS package Plimpton (1995). We simulated 128 replicas with geometrically-distributed temperatures $T^{*}(n)=1.045^{n}$ regulated by a Nosé–Hoover thermostat for up to $2\times 10^{9}$ time steps, attempting exchanges between replicas every $10^{3}$ steps, and dumping snapshots every ${10^{6}}$ steps. We employ a variable time step which is dynamically determined such that no bead moves by more than $0.005\sigma$ per time step 222The time step is thus smaller than in standard LJ simulations because we need to resolve the time scale of bond vibrations.. The polymer packing fraction is $\nu=3N(R_{0}-0.5)^{-3}/16t_{m}$, where the accessible capsid volume is rescaled due to finite chain thickness. Thus, for $100\leq N\leq 250$, $0.23\leq\nu\leq 0.57$. For larger packing fractions (we tried up to $N=350$, $\nu=0.8$), slow kinetics prevents equilibration.

We a posteriori minimized the energy of each snapshot using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (l-BFGS) algorithm  Nocedal (1980); Liu and Nocedal (1989) as implemented in the pele package Stevenson et al. , converging the root-mean-squared gradient to a tolerance of $1\times 10^{-3}\frac{\epsilon}{\sigma}$. We further optimized the energy of selected structures with a combination of random displacement, reptation, and l-BFGS minimization cycles (see algorithm S1).

Our theoretical analysis suggests that the elastic energy of a polymer packed into a sphere can be minimized by segregating into up to three distinct spools. Chain segments that belong to the same spool lie in almost parallel planes. As such, we can characterize the degree to which two beads belong to the same spool domain by comparing the binormal vectors at their positions. The binormal to the chain is given by $\mathbf{\hat{b}}=\mathbf{\hat{t}}\times\mathbf{\hat{n}}=\frac{\mathbf{r^{{}^{\prime}}}\left(s\right)\times\mathbf{r^{{}^{\prime\prime}}}\left(s\right)}{\left|\mathbf{r^{{}^{\prime}}}\left(s\right)\times\mathbf{r^{{}^{\prime\prime}}}\left(s\right)\right|}$, where $\mathbf{r}$ is the position of the particle, $s$ is the arc length along the chain, $\mathbf{\hat{t}}$ is the unit tangent vector, and $\mathbf{\hat{n}}$ is the unit normal vector. A simple way to visualize the structure of a compacted chain is to plot a correlation plot of the angles between binormal vectors at each pair of particle centers (e.g. Fig. 2). We identify the principal domains by searching for cliques using the networkx package Hagberg et al. (2008).

The size of a domain $D_{i}$ is denoted as $m_{i}=w_{i}N$. The largest domain, $D_{1}$, is the largest set of beads such that the angle between the binormal to the chain for every pair of beads is less than $30^{\circ}$ (see Fig. S5). The $n^{\mathrm{th}}$ largest domain, $D_{n}$ is then the largest such set disjoint with the union of all larger domains, $\bigcup_{i=1}^{n-1}D_{i}$.

Results. Analysis of the partitioning at different conditions reveals the following picture. In the lowest-temperature replicas we observe spool morphologies for all chain lengths (see Fig. 2). The ground state of a short chain ($N=100$) is a single spool whose ends curve inward toward the sphere center (Fig. 2a). Simulations of 150- and 200-bead chains show the emergence of double-spool structures (Fig. 2b-c), while the variety of low-energy morphologies sharply increases at $N=200$. This trend is illustrated by histograms of $m_{1}$ (see Fig. S6).

A double spool is the putative ground state structure for 200 beads, but other morphologies, resembling topological links, are also observed (Fig. 3a-b) 333We stress that these structures only resemble topological links, in a superficial sense; the simulated chains are not closed, nor, in general, are the chain segments that constitute the ‘components’ of the links contiguous.. The spool-like structures are consistent with our theoretical assumptions (Fig.1a-c). We also observe a hierarchical mass distribution $m_{1}\gg m_{2}\gg m_{3}$ characteristic of the optimal analytic solutions, distinguishing spools from other morphologies, though quantitative comparison of $w_{i}$ between theory and simulations is difficult due to the sensitivity of the partitioning to the model parameters.

Upon increasing the polymer length to 250 beads, a triple spool (Fig. 2d) becomes more favourable than double spools, but non-spool-like morphologies are more commonly observed. The structure with the lowest elastic energy ($V_{b}+V_{\theta}$) is a Hopf link-like morphology (Fig. 3c), some $36\epsilon$ lower in total energy than the triple spool. The putative ground state, another $\approx 59\epsilon$ lower in total energy, resembles three linked rings (L6n1 in Dowker-Thistlethwaite notation Doll and Hoste (1991); Figs. 3d, S7). For the dependence of the lowest energies of these motifs on chain length, see Fig. S8. Twisted packings have been proposed as models of toroidal aggregates of long chiral molecules such as DNA Kulić et al. (2004); Charvolin and Sadoc (2008); even in the absence of intrinsic twist, such packings were shown in Shin and Grason (2011) to be preferred at high densities, as they can fill the sphere without disclinations or voids. Space-filling considerations may also drive the transition to twisted, link motifs we see here.

The degree of binormal ordering decreases with temperature (Fig. 4). The distribution of largest domain sizes $m_{1}$ for a chain of 250 beads is visualized as a heat map for temperatures $1\leq T^{*}\leq 270$. spool structures ($m_{1}\gtrsim 100$) are rarely observed beyond $T^{*}=10$. For $10<T^{*}<40$, there is less order, but there are clear indications of link motifs. At high temperatures ($T^{*}>100$), while patches of order are retained near the boundary wall—due to local hexagonal packing—large ordered domains seldom emerge.

对手机会员入

Discussion & Conclusions. According to our results, at low temperature, spherical confinement is sufficient to order the majority of the filament into tight spools. We predict a variety of multi-domain morphologies such as nested spools and links of two or three components competing in the low temperature regime. Their energetic ordering may be sensitive to the interplay between elastic and repulsive forces and other microscopic details of the system. Our low-temperature results can be compared with previous experiments and simulations of packing metallic wires into spheres Stoop et al. (2011); Shaebani et al. (2017), where ring-like coiling was observed, reminiscent of our nested spools.

However, those systems are in the athermal limit, and the structures adopted by wires that are pushed into the system with large forces are protocol-dependent. In our work, we avoid this protocol dependence by employing the replica exchange method; as such, our low-temperature structures may describe kinetically-inaccessible ground states of their systems. Furthermore, in the case of metallic wires, the confinement size was varied, with a focus on relatively large containers. In our work, we fixed the enclosure size to $R_{0}=5\sigma$, which is at the very lowest end of their parameter range, where little data is available about the observed structures. Further work is needed to explore the packing as a function of the capsid size.

On the other hand, the chosen capsid size is appropriate for modeling DNA packing into viral capsids. For DNA, $\sigma\approx 2$ nm, which means that the viral capsid corresponding to our model would be about 20 nm in diameter, which is at the small end of the range for viral capsids in nature (20 - 200 nm). Given a persistence length of 50 nm, the DNA regime is described by simulation temperatures $T^{*}\approx 40$. The structures we observe in this regime are only partially ordered, retaining only some features of low-temperature motifs. Approaching the high-temperature limit (the flexible-chain regime), we eventually observe little structure, in agreement with previous studies Cacciuto and Luijten (2006).

Projecting our model on real systems, we acknowledge that system-specific details may influence packing motifs. Specific DNA–capsid interactions enhance the ordering of DNA in viral capsids via built-in folding pathways Twarock et al. (2018). Structures incompatible with the $I_{h}$ symmetry of many viral capsids (e.g., the orthogonal Hopf link with $D_{2d}$ symmetry) might be disfavoured. Moreover, the kinetic pathway of motor-driven insertion likely selects spool-like morphologies Rapaport (2016) even if they are metastable. Similarly, capsid ejection dynamics might depend crucially on the topology of the DNA structure, with knotted Marenduzzo et al. (2013) or linked morphologies possibly obstructing the ejection.

Acknowledgements. We acknowledge enlightening discussions with Arman Boromand, Daan Frenkel, Erika Eiser, Erik Luijten, and William Gelbart. The work was supported by the EU’s Horizon 2020 Program through grants ETN 674979-NANOTRANS and FET-OPEN 766972-NANOPHLOW, by the ARRS grant Z1-9170, by the 1000-Talents Program of the Chinese Foreign Experts Bureau, and by the Chinese National Science Foundation through grants 11874398, 11850410443, and 21850410459.