Spool-nematic ordering of dsDNA and dsRNA under confinement

Spatial organization of genomes in viral containers Sun et al. (2010) and its physical principles Zandi et al. (2020) is an outstanding problem with a long history Pollard (1953), resulting in refined models based on cryo-electron microscopy Luque and Castón (2020), X-ray scattering techniques Khaykelson and Raviv (2020), thermodynamic osmotic pressure methodology Gelbart and Knobler (2009) and single molecule studies Smith (2011). These experiments display a range of morphologies (for a recent review see Ref. Comas-Garcia and Rosales-Mendoza (2023)). In the case of bacteriophages, high-resolution 3D reconstructions indicate inhomogeneous ordering of the encapsulated genome, with the dsRNA/DNA partitioning into ordered Ilca et al. (2019); Leforestier (2013) or disordered Leforestier and Livolant (2010) nanodomains. Similarly, measurements of DNA conformational dynamics during packaging Berndsen et al. (2014), as well as intermittent ejection dynamics Chiaruttini et al. (2010) implicate multi-domain structures. However, so far the general assumption has been that the equilibrium ground state is a single inverted spool Purohit et al. (2003); Zandi et al. (2020), or multiple spools Curk et al. (2019).

Mesoscopic theories Klug et al. (2005); Šiber et al. (2008); Shin and Grason (2011); Podgornik et al. (2016); Liang et al. (2019) and coarse-grained simulations of double-stranded genome packaging Kindt et al. (2001); Petrov and Harvey (2008); Rapaport (2016); Marenduzzo et al. (2009); Curk et al. (2019); Marzinek et al. (2020) are usually performed using a semiflexible polymer model that assumes linear elasticity with bending modulus $B\approx 50\,k_{\textnormal{B}}T$nm, where $k_{\textnormal{B}}$ is the Boltzmann constant and $T$ the temperature. However, the realistic bending response of dsDNA and dsRNA is highly non-linear due to local melting, which enables the formation of kinks. For dsDNA, kinks form beyond a critical bending torque ${\tau_{\textnormal{c}}\approx 30}$ pN nm Yan and Marko (2004); Qu et al. (2011); Qu and Zocchi (2011), and typically comprise about three melted base pairs with the reversible work required to form a kink ${\mu\approx 12k_{\textnormal{B}}T}$ Qu et al. (2011). Such non-linear effects can alter optimal packing configurations; in particular, a spherically-confined polymer that is able to kink exhibits increased local nematic ordering near the sphere surface Myers and Pettitt (2017); Bores and Pettitt (2020).

A general theory of a semi-flexible kinkable chain postulates an elastic free-energy density per chain length

with $\kappa$ the local curvature, $\mu$ the free energy per kink and $\rho_{\textnormal{k}}$ the local density of kinks in the chain. Here, we first perform non-equilibrium simulations to understand the packing arrangement of a kinkable, semi-flexible chain actively pushed into a spherical enclosure, and then analyze the stability of different packings using analytical theory.

A coarse-grained model that accurately describes the properties of DNA or dsRNA at small curvatures is the worm-like-chain (WLC) model; a semi-flexible polymer of length $L$ and circular cross-section $\sigma$ with persistence length $\ell_{\textnormal{p}}=B/k_{\textnormal{B}}T$. Efficient description of nonlinear elasticity includes the kinks as regions with locally reduced persistence length $\ell^{\textnormal{m}}_{\textnormal{p}}$. In the two-state kinkable WLC (KWLC) model Wiggins et al. (2005), molten regions behave as freely-jointed chains ($\ell^{\textnormal{m}}_{\textnormal{p}}=0$), while the more accurate meltable WLC (MWLC) Yan and Marko (2004); Sivak and Geissler (2012) assigns a single-stranded persistence length $\ell^{\textnormal{m}}_{\textnormal{p}}\approx 2~{}\textnormal{nm}$ to melted sections.

We model dsDNA as a discrete MWLC: a bead-spring polymer with $N$ beads of diameter $\sigma=2\,\textnormal{nm}$. Consecutive beads are bonded by a two-body stretching term ${V_{\textnormal{b}}=K\left(r/\sigma-1\right)^{2}}$, where $r$ is the distance between two consecutive beads and the prefactor ${K=16k_{\textnormal{B}}T\ell_{\textnormal{p}}/\sigma}$ Curk et al. (2019), while excluded volume and short range hydration interactions at relevant packing densities Podgornik et al. (2016) are incorporated as a WCA repulsion Weeks et al. (1971) with strength ${\epsilon=k_{\textnormal{B}}T}$. To enable molecular dynamics (MD) simulations, we adapt the two-state MWLC model by writing it as a continuous, three–body bending interaction obtained as a canonical-ensemble superposition of a non-melted and a melted state. Assuming the states in each three-body segment are independent, and defining the Boltzmann factors ${q=\exp{\left[-\ell_{\textnormal{p}}(1-\cos{\theta})/\sigma\right]}}$ and ${q^{\textnormal{m}}=\exp{\left[-\beta\mu_{\textnormal{MWLC}}-\ell_{\textnormal{p}}^{\textnormal{m}}(1-\cos{\theta})/\sigma\right]}}$, the potential of mean force is

with $\theta$ the angle between consecutive bond vectors, and $\mu_{\textnormal{MWLC}}$ the melting penalty per bead (Fig. 1). We obtain the average kink density $\overline{\rho}_{\textnormal{k}}$ as the canonical average, $\overline{\rho}_{\textnormal{k}}=\left\langle q^{\textnormal{m}}/\left(q+q^{\textnormal{m}}\right)\right\rangle/\sigma$. By fitting to experimental data, we find the MWLC model captures the non-linear elasticity of DNA bending at parameters $\ell^{\textnormal{m}}_{\textnormal{p}}=\sigma$, $\mu_{\textnormal{MWLC}}=10k_{\textnormal{B}}T$ and reproduces the kink formation of dsDNA chain at critical torque $\tau_{\textnormal{c}}=32\,\textnormal{pN nm}$ (see Fig. S1 in Note (1)). Obtaining the kink formation free-energy at a specific base pair requires de-convolution of the coarse-grained model Schöpflin et al. (2012), but since $\mu\gg k_{\textnormal{B}}T$, the probability that more than one kink occurs within a single bead is negligible. Thus, the full ($180^{\circ}$) kink formation penalty at a specific base-pair is $\mu=\mu_{\textnormal{MWLC}}{(\sigma)}+2k_{\textnormal{B}}T+k_{\textnormal{B}}T\log(\sigma/d_{\textnormal{n}})\approx 13.8k_{\textnormal{B}}T$, where the last term re-scales the number of kink states at bead size $\sigma$ relative to length per nucleotide $d_{\textnormal{n}}\approx 0.33~{}\textnormal{nm}$, and we emphasize that $\mu_{\textnormal{MWLC}}$ is a function of the discretization length $\sigma$.

Using the MWLC model we perform MD simulations of DNA packing into a spherical capsid Plimpton (1995). Confinement is implemented as a WCA repulsion at $R=6\sigma$ with strength $\epsilon$ and diameter $\sigma$, giving an effective enclosure radius $R_{\textnormal{eff}}=5.5\sigma$. The maximum packing density of chains is determined by hexagonal close-packing Podgornik et al. (2016). Taking the volume of the chain relative to the maximum packed volume we define the packing fraction $\eta\equiv 3\sqrt{3}\sigma^{3}N_{\textnormal{in}}/(8\pi R_{\textnormal{eff}}^{3})$ where $N_{\textnormal{in}}$ the number of beads packed in the capsid.

We model the viral packing motor as a cylindrical portal of length $\lambda_{\textnormal{p}}=2\sigma$ and radius $R_{\textnormal{p}}=1.2\sigma$ (Fig. 1A inset). The motor force operates along the portal axis on beads inside the portal, implemented as a parabolically-modulated traveling wave

where $x$ is the bead position in the portal ($0\leq x\leq\lambda_{\textnormal{p}}$), $\omega$ the packing rate, $t$ the time, and $f_{\textnormal{s}}=50k_{\textnormal{B}}T\sigma^{-1}$ the stall force. Real viral motors pack in a discrete fashion, at a rate dependent on the availability of ATP Fuller et al. (2007a, b). Our periodic motor force mimics this behavior, and also allows us to tune the stall force and packing rate independently. To model diffusive dynamics in an aqueous solution we apply a Langevin thermostat with damping time ${\tau=\sqrt{m\sigma^{2}/(k_{\rm B}T)}}$ where $m$ is the mass of the beads, which via the Stokes-Einstein relation Note (2) results in a simulation timescale $\tau\approx 18$ ns. Since the torsional strain can relax on the packing timescale Rollins et al. (2008); Rapaport (2016), we omit torsional terms from our models.

Using this setup we ran packing simulations of a 600-bead chain into a radius $R_{\textnormal{eff}}=5.5\sigma$ sphere at rates from $\omega=1\times 10^{-4}\sigma\tau^{-1}$ to $5\times 10^{-2}\sigma\tau^{-1}$, with 120 independent simulations per $\omega$ and evaluate the line density of spontaneously formed kinks (Fig. 1B). The initial peak corresponds to the appearance of the first kink during packing. The observed kink density is strongly rate-dependent; however, at slowest packing rates, $\omega\sim 1\times 10^{-4}\sigma/\tau\approx 40\,\textnormal{kb}/\textnormal{s}$, it converges to $\overline{\rho}_{\textnormal{k}}\approx 0.03/\sigma$. This rate is still an order of magnitude faster than in real viral motors $\omega\lesssim 2\textnormal{kb}/\textnormal{s}$ Fuller et al. (2007b, a), but the observed convergence suggests that reducing the rates further would not alter the packed configurations. Moreover, the equilibrium theory ($\omega\to 0$) shown below agrees well with the slow-packing simulations, suggesting that our predictions capture the biologically relevant packing rates.

Spherically-confined dsDNA is thought to arrange into spools Purohit et al. (2003); Curk et al. (2019), but kinking permits nematic ordering Myers and Pettitt (2017); Bores and Pettitt (2020), or a mixed phase of spools wrapped around a nematic core Ilca et al. (2019). To quantify spool and nematic order in our simulation configurations, we compute the order parameter tensor,

where $\hat{\upsilon}_{i\alpha}$ is the $\alpha^{\textnormal{th}}$ component of either the unit binormal vector, $\bm{\hat{b}}$, or the unit tangent vector, $\bm{\hat{\tau}}$, at chain segment $i$. The order parameter, $S$, is the principal eigenvalue of $\mathbf{Q}$.

Taking $\bm{\hat{\upsilon}}=\bm{\hat{b}}$, gives the spool order parameter, $S_{\textnormal{spl}}$, similar to that defined in Ref. Nikoubashman et al. (2017), while taking $\bm{\hat{\upsilon}}=\bm{\hat{\tau}}$ gives the nematic order parameter, $S_{\textnormal{nem}}$. $S_{\textnormal{nem}}$ is of limited practical use because the spool domain is always more massive than the nematic domain. A more useful mixed spool–nematic order parameter $S_{\textnormal{mix}}$ is defined by setting

where $r_{i}$ is the distance from the center-of-mass of the $i^{\textnormal{th}}$ chain segment to a director $\bm{\hat{z}}$, centered at the origin, and $r_{\textnormal{c}}$ defines the spool–nematic boundary. In this case, $S_{\textnormal{mix}}$ must be maximized with respect to both $r_{\textnormal{c}}$ and $\bm{\hat{z}}$. We denote the partial packing fractions of the spool and nematic domains $\eta_{\text{spl}}$ and $\eta_{\text{nem}}$, respectively.

At the slowest packing rates, we found that the chains initially arrange into a loose spool (Fig. 2A). At intermediate packing densities ($\eta\approx 0.3$), a distinct inner domain emerges exhibiting some nematic order (Fig. 2B), which is reflected in the growth of $S_{\textnormal{mix}}$ while $S_{\textnormal{spl}}$ remains constant (Fig. 2D). Upon further packing, $S_{\textnormal{spl}}$ begins to grow faster than $S_{\textnormal{mix}}$, indicating that the growth of the nematic core has stopped ($\eta_{\text{nem}}$ levels off at $\approx 0.16$), while chain continues to be packed into the outer spool, and the core shrinks in size from $\approx 2.8\sigma$ to $\approx 2.1\sigma$ (see Figs. S3 and S4 in Note (1)). These trends, and the snapshot in Fig. 2C, confirm our main observation: for dense packings, the genome adopts a mixed spool–nematic configuration.

Fast and slow packings are indistinguishable at low densities, but the growth of $S_{\textnormal{mix}}$ is dramatically reduced at faster packing rates, even falling at high densities (Fig. 2D) indicating kinetically-arrested disorder. Moreover, for the fastest packing rate $r_{\textnormal{pack}}=0.05\sigma\tau^{-1}$, the size of the core region $r_{\textnormal{c}}$ remains roughly constant $r_{\textnormal{c}}\approx 2.8\sigma$ during the time the kink density doubles from $0.04$ to $0.08$, indicating that the chain has insufficient time to relax and grow the spool domain as packing proceeds. At slower packing rates the structures are more ordered and approach the spool–twisted-nematic configuration, while the values of all three order parameters and the kink density are well-converged among the three slowest packing rates.

To investigate whether our packing is under thermodynamic or kinetic control, we consider an analytical theory where a chain confined to a radius $R_{0}$ sphere adopts a helical configuration parametrized as ${u\mapsto\left(r\cos{u},r\sin{u},bu\right),u\in\left[-h/b,h/b\right]}$ with radius $r$, pitch length $2\pi b$, and helix height $h=\sqrt{R_{0}^{2}-r^{2}}$. The slope $p=b/r$ corresponds to the helix angle $\theta=\arctan{\left(1/p\right)}$. The limit $p\to\infty$ corresponds to a nematic, and $p\to 0$ to a spool phase. To obtain the optimal slope, we minimize the free energy [Eq.  (1)] assuming that each radial shell is independent of its neighbors. The local curvature of a helix is ${\kappa(r,b)=|r|/\left(r^{2}+b^{2}\right)}$, and the contour length of a helix is $L{\left(r,b\right)}=\left(2h/b\right)\sqrt{r^{2}+b^{2}}$. Assuming a $180^{\circ}$ kink forms whenever the chain hits the confining capsid, the kink density at radius $r$ is $\rho_{\textnormal{k}}{\left(r,b\right)}=1/L{\left(r,b\right)}$. The optimal non-zero slope $p^{\textnormal{opt}}$ is then the real positive root of

where ${A=\left(r^{*}\right)^{2}/B^{*}\sqrt{1-(r^{*})^{2}}}$, with ${B^{*}=B/(\mu R_{0})}$ the ratio of bending stiffness to kink energy and ${r^{*}=r/R_{0}}$. Surprisingly, we find that the optimal slope is a discontinuous function of the helix radius (Fig. 3): for small $r$, i.e., $A<A_{\textnormal{c}}\equiv\sqrt{32/27}$, the optimal slope is given by

whereas $p^{\textnormal{opt}}=0$ for $A>A_{\textnormal{c}}$. At $A=A_{\textnormal{c}}$ the optimal slope jumps from zero to $p^{\textnormal{opt}}=\sqrt{2}$, which, interestingly, describes a helix with an angle complementary to the magic angle $\theta_{m}=\arctan{\sqrt{2}}$. This observation implies that the outer shells form a spool ($p=0$), while the inner shells form a twisted nematic domain with a large slope (Fig. 3A).

Whether the spool or the nematic core forms at low packing density depends on the ratio $B^{*}$. $B^{*}>1$ favors the spool phase, whereas at $B^{*}<1$ the nematic core is more stable (Fig. 3A, inset). For a typical virus ($2R_{0}\approx 50$ nm) and DNA parameters ($B=50k_{\textnormal{B}}T\,\textnormal{nm},\,\mu=13.8k_{\textnormal{B}}T$) we find $B^{*}\approx 0.14$ and so predict that a spool forms initially with a transition to a nematic core at radius $r_{\textnormal{c}}\approx 0.38R_{0}$. Hence, the majority of the DNA should form a spool, enclosing a twisted-nematic core.

This result explains the experimental observation of a helical core structure in a dsRNA bacteriophage Ilca et al. (2019); taking $R_{0}=20.8\,\textnormal{nm}$ and $B=63k_{\textnormal{B}}T\,\textnormal{nm}$ Abels et al. (2005); Hyeon et al. (2006), and assuming $\mu=13.8k_{\textnormal{B}}T$ (the value for dsDNA), we find excellent agreement with experiment on both the critical radius $r_{\textnormal{c}}$ and the helix angle at the transition (Fig. 3A). The analytical prediction of the helix angle also agrees with simulation results at low packing rates ($\omega\leq 0.0005\sigma\tau^{-1}$, Fig. 3B), suggesting that, in the slow packing regime, our simulations are under thermodynamic control, and the packing morphology is determined by equilibrium free-energy minimization.

To obtain the full phase diagram of packing configurations we use the continuum theory of polymer nematics Shin and Grason (2011); Svenšek et al. (2010); Curk et al. (2019) where the polymer is described by a vector field $\mathbf{t}(\mathbf{r})$. Assuming that kinks form only at the enclosure surface, and disregarding entropic contributions due to conformational fluctuations, the lowest-order non-trivial terms of the elastic free-energy functional are given by,

The first term is the bending term Svenšek et al. (2010); Shin and Grason (2011) with $\hat{\mathbf{t}}=\mathbf{t}/|\mathbf{t}|$ the unit tangent vector to the chain contour, while the second describes the kink density at the enclosure surface with normal vector $\mathbf{n}$. This expression is an integral form of Eq. (1) with prefactors $K=2B/{(\pi\sigma^{2})}$ and $\varepsilon_{\textnormal{kink}}=2\mu/(\pi\sigma^{2})$.

The optimal structure of a confined meltable filament results from competition between the elastic and melting terms in Eq. (8). For stiff polymers, kinking is preferable to bending, so they will tend to arrange into nematically-aligned straight segments with kinks at the boundaries. In contrast, flexible polymers are expected to form spool structures without kinks. To explore these competing mechanisms, we assume that the chains form a two-domain structure: an outer spool and a nematic core. The spool contains no kinks ($\mathbf{n}\cdot\mathbf{t}=0$), while the nematic core does not contribute to bending ($\nabla\times\hat{\mathbf{t}}=0$). Assuming phases are close packed, the spool free-energy is Curk et al. (2019); Purohit et al. (2003),

while the nematic free energy is

Since the spool and the nematic core occupy disjoint volumes, and connections between domains can be neglected in the long-chain limit, the total free energy is the sum, $E(\eta,\eta_{\textnormal{spl}},R_{0})=E_{\textnormal{spl}}(\eta_{\textnormal{spl}},R_{0})+E_{\textnormal{nem}}(\eta-\eta_{\textnormal{spl}},R_{0})$, with $\eta_{\textnormal{spl}}$ the spool volume. Minimizing $E(\eta,\eta_{\textnormal{spl}},R_{0})$ we find a phase diagram (Fig. 4) that is a function of two dimensionless parameters: $B^{*}=B/(\mu R_{0})$ and $\eta$. The boundary delineating the pure nematic and the mixed phase is $B^{*}=(1-\eta)^{-1/3}$, while the boundary between the pure spool and the mixed phase is $B^{*}=1-\eta^{2/3}$.

The coexistence region emerges at $B^{*}\sim 1$, growing with $\eta$ and dominating in the high packing fraction limit. For a typical $50\,\textnormal{nm}$ virus, $B^{*}\approx 0.14$, implying a spool phase at low packing fractions, enclosing a nematic core at ${\eta>(1-B^{*})^{3/2}\approx 0.8}$. This theory assumes coupled shells with no nematic twist in the core ($p\to\infty$); considering the opposite limit of free shells slightly stabilizes the nematic phase by allowing non-zero twist [Eq. (7)], but the shift in the phase boundaries is less than $0.1B^{*}$ (dashed line in Fig. 4).

Experimental data can be mapped onto the phase diagram by assuming that $\mu$ for RNA and DNA are the same, and the range of $\eta$ is bounded by hexagonal and cubic packings (both observed experimentally). With such a mapping, the observed spool–twisted-nematic structure packings in $\phi 6$ bacteriophage Ilca et al. (2019) and herpes simplex virus Chen et al. (2023) are consistent with our theoretical predictions (Fig. 4).

In summary, our theoretical and simulation results imply that dsDNA/dsRNA packed into a sphere adopts a spool–nematic configuration. Coexistence between an outer spool and a nematic core arises spontaneously both in packing simulations and equilibrium theory based on the minimization of elastic energy. We explored the strong and weak limits of orientational coupling, finding nearly-identical phase diagrams for a meltable chain, with the main difference being whether the nematic core exhibits non-zero twist. We expect that packing of dsDNA and dsRNA in viruses falls between these two limits, forming a large outer spool with a (twisted-)nematic core emerging at high packing densities ($\eta>0.8$ for a typical $2R_{0}=50$ nm virus). These results explain the experimental data for dsRNA packing in a bacteriophage virus Ilca et al. (2019).

We have considered packing into a pre-formed spherical capsid, but our findings are based on global free-energy minimization and thus expected to hold qualitatively for other packing situations such as DNA condensation or ordering in bacterial nucleomes. For example, kink formation and nematic ordering may explain why plasmid DNA assembles into rod-like structures Jiang et al. (2013) rather than tori Hoang et al. (2014). Moreover, our phase diagram shows that spool–nematic ordering is expected for any strongly confined kinkable polymer. Non-linear bending and kinking behavior might be a common characteristic of helical polymers, as their structure is usually determined by reversible supramolecular interactions.