Fine Structure of Viral dsDNA Encapsidation

TL;DR

This paper models the liquid crystal ordering of viral dsDNA within bacteriophage capsids, using advanced numerical methods to analyze phase behavior, core disordering, and potential implications for phage design.

Contribution

It introduces a combined chromonic and liquid crystal model with a scalar order parameter and develops a finite element algorithm for equilibrium analysis and pressure calculation.

Findings

Identifies well-ordered and disordered DNA regions within capsids.

Provides phase space analysis of pressure and core radius.

Enables potential bacteriophage design through parameter exploration.

Abstract

In vivo configurations of dsDNA of bacteriophage viruses in a capsid are known to form hexagonal chromonic liquid crystal phases. This article studies the liquid crystal ordering of viral dsDNA in an icosahedral capsid, combining the chromonic model with that of liquid crystals with variable degree of orientation. The scalar order parameter of the latter allows us to distinguish regions of the capsid with well-ordered DNA from the disordered central core. We employ a state-of-the-art numerical algorithm based on the finite element method to find equilibrium states of the encapsidated DNA and calculate the corresponding pressure. With a data-oriented parameter selection strategy, the method yields phase spaces of the pressure and the radius of the disordered core, in terms of relevant dimensionless parameters, rendering the proposed algorithm into a preliminary bacteriophage designing tool. The presence of the order parameter also has the unique role of allowing for non-smooth capsid domains as well as accounting for knot locations of the DNA.

Figures

Tables

Table 1. Table 1: Physical measurements of four different bacteriophages.The symbol L p subscript 𝐿 𝑝 L_{p} denotes the persistence length of a DNA chain of length L 𝐿 L , effective diameter d 𝑑 d , in a sphere-like capsid of radius R c subscript 𝑅 𝑐 R_{c} with a measured radius r 0 subscript 𝑟 0 r_{0} of the disordered core. T4 [ 24 , 35 ] ; T5 [ 8 ] ; T7 [ 2 ] ; ϵ italic-ϵ \epsilon 15 [ 14 ] ; P22 [ 19 , 15 ] .

Virus	$L_p$ (nm)	$d$ nm	$L$ (nm)	$R_c$ (nm)	$r_0/R_c$	Pressure(atm)
T4	52.8	2.4	55047.6	40.0	0.5500	28.70
T5	58.38	2.94	39423.8	42.0	0.4286	30.17
T7	52.88	2.6	12932.0	26.05	0.5889	44.02
$ϵ$15	53.9	2.55	12846.0	28.37	0.5735	40.41

Equations

$$E=$$

$${\mathcal{W}}_1(\mathbf{n},\nabla \mathbf{n})=$$

$$+(k_2+k_4)[\mathrm{tr}([\nabla \mathbf{n}{]}^2)-(\mathrm{div}\mathbf{n}{)}^2]$$

$$\begin{array}{rl}{E}_{\mathrm{erk}}[s,\mathbf{n}]& :=\frac{1}{2}{\int }_{\Omega }{s}^2{\mathcal{W}}_1(\mathbf{n},\nabla \mathbf{n})+b|\nabla s{|}^{2}d\mathbf{x}+{\gamma }_0{\int }_{\Omega }w(s)d\mathbf{x},\end{array}$$

$$\begin{array}{rl}{E}_{\mathrm{a}}(s,\mathbf{n})& :=E_{\mathrm{a},\mathbf{n}}(s,\mathbf{n})+E_{\mathrm{a},s}(s),\\ E_{\mathrm{a},\mathbf{n}}(s,\mathbf{n})& :=\frac{1}{2}{\int }_{\Gamma }{s}^2\left[a_{\perp }(\mathbf{\nu }\cdot \mathbf{n}{)}^2+a_{\parallel }[1-(\mathbf{\nu }\cdot \mathbf{n}{)}^2]\right]dS(\mathbf{x}),\\ E_{\mathrm{a},s}(s)& :=\frac{1}{2}{\int }_{\Gamma }{a}_{\mathrm{ori}}(s-s_{\mathrm{a}}{)}^{2}dS(\mathbf{x}),\end{array}$$

$$\begin{array}{r}{E}_{\mathrm{chr}}=E_{\mathrm{erk}}[\mathbf{n},s]+E_{\mathrm{elas}}[\nabla \mathbf{u}],\end{array}$$

$$\nabla \cdot \mathbf{n}=0,\text{and}(\nabla \times \mathbf{n})\cdot \mathbf{n}=0.$$

$$\frac{d\omega }{d\mathbf{m}}=\frac{2\pi }{d},\frac{d\vartheta }{d\mathbf{p}}=\frac{2\pi }{d},s\ne 0.$$

$$\frac{d\mathbf{r}}{d\zeta }=\mathbf{n}(\mathbf{r}(\zeta )),\mathbf{r}(0)={\mathbf{r}}_0,$$

$$E_{\mathrm{elas}}=B|\nabla (\mathbf{m}+\mathbf{p}){|}^2+C|\nabla (\mathbf{m}-\mathbf{p}){|}^2,$$

$$E=$$

$$\mathbf{n}\cdot \mathbf{m}=0,\mathbf{m}\cdot \mathbf{p}=0,\mathbf{n}\cdot \mathbf{p}=0.|\mathbf{n}|=|\mathbf{m}|=|\mathbf{p}|=1,$$

$$k_1$$

$$b$$

$$a_0$$

$$\begin{array}{rl}w(s_{\mathrm{a}})& =\underset{0\le s\le 1}{\min }w(s)=0,w\text{ is quadratic and convex near }s=s_{\mathrm{a}},\\ w(0)& =1,w\text{ has a local max at }s=0,\end{array}$$

$$\begin{split}\mathcal{E}[s,\mathbf{n},\mathbf{m},\mathbf{p},{\Omega}]&=\frac{1}{2}\int_{\Omega}s^{2}\big{[}k_{1}(\mathrm{div}\,\mathbf{n})^{2}+k_{2}(\mathbf{n}\cdot\mathrm{curl}\,\mathbf{n})^{2}+k_{3}|\mathbf{n}\times\mathrm{curl}\,\mathbf{n}|^{2}\\ &\qquad\quad+k_{2}[\mathop{\rm tr}\nolimits([\nabla\mathbf{n}]^{2})-(\mathrm{div}\,\mathbf{n})^{2}]\big{]}d{\mathbf{x}}+\frac{1}{2}\int_{\Omega}b|\nabla s|^{2}d{\mathbf{x}}\\ &\quad+\gamma_{0}\int_{\Omega}w(s)\,d{\mathbf{x}}+\frac{1}{2}a_{0}\left[\int_{\Gamma}s^{2}(\bm{\nu}\cdot\mathbf{n})^{2}dS({\mathbf{x}})+\int_{\Gamma}(s-s_{\mathrm{a}})^{2}dS({\mathbf{x}})\right]\\ &\quad+\frac{1}{2}\int_{\Omega}\left[Bs^{2}|\nabla(\mathbf{m}+\mathbf{p})|^{2}+Cs^{2}|\nabla(\mathbf{m}-\mathbf{p})|^{2}\right]d{\mathbf{x}}.\end{split}$$

$${\tilde{k}}_3=k_3(1+\alpha ).$$

$$(\mathbf{n}\cdot \mathrm{curl}\mathbf{n}{)}^2+|\mathbf{n}\times \mathrm{curl}\mathbf{n}{|}^2+\mathrm{tr}([\nabla \mathbf{n}{]}^2)=|\nabla \mathbf{n}{|}^2,$$

$$\begin{array}{rl}\mathcal{E}[s,\mathbf{n},\Omega ]& =\frac{{\tilde{k}}_3}{2}{\int }_{\Omega }{s}^2{\tilde{\mathcal{W}}}_1(\mathbf{n},\nabla \mathbf{n})d\mathbf{x}+\frac{1}{2}{\int }_{\Omega }b|\nabla s{|}^{2}d\mathbf{x}+{\gamma }_0{\int }_{\Omega }w(s)d\mathbf{x}\\ +& \frac{1}{2}{a}_0\left[{\int }_{\Gamma }{s}^2(\mathbf{\nu }\cdot \mathbf{n}{)}^{2}dS(\mathbf{x})+{\int }_{\Gamma }(s-s_{\mathrm{a}}{)}^{2}dS(\mathbf{x})\right]\end{array}$$

$$\begin{array}{rl}{\tilde{\mathcal{W}}}_1(\mathbf{n},\nabla \mathbf{n})& =({\tilde{k}}_1/{\tilde{k}}_3)(\mathrm{div}\mathbf{n}{)}^2+({\tilde{k}}_2/{\tilde{k}}_3)(\mathbf{n}\cdot \mathrm{curl}\mathbf{n}{)}^2+|\mathbf{n}\times \mathrm{curl}\mathbf{n}{|}^2\\ & +({\tilde{k}}_2/{\tilde{k}}_3)[\mathrm{tr}([\nabla \mathbf{n}{]}^2)-(\mathrm{div}\mathbf{n}{)}^2].\end{array}$$

$$\Omega =L_0^3\hat{\Omega },\Gamma =L_0^2\hat{\Gamma },\partial {\Omega }_0=L_0^2\widehat{\partial {\Omega }_0},\nabla =\frac{1}{L_0}\hat{\nabla },$$

$$\begin{array}{rl}\mathcal{E}[s,\mathbf{n}]& ={\tilde{k}}_3{L}_0\hat{\mathcal{E}}[s,\mathbf{n}]\\ \hat{\mathcal{E}}[s,\mathbf{n}]& :=\frac{1}{2}{\int }_{\Omega }{s}^2{\tilde{\mathcal{W}}}_1(\mathbf{n},\nabla \mathbf{n})d\mathbf{x}+\frac{\hat{b}}{2}{\int }_{\Omega }|\nabla s{|}^{2}d\mathbf{x}+{\hat{\gamma }}_0{\int }_{\Omega }w(s)d\mathbf{x}\\ & +\frac{{\hat{a}}_0}{2}\left[{\int }_{\Gamma }{s}^2(\mathbf{\nu }\cdot \mathbf{n}{)}^{2}dS(\mathbf{x})+{\int }_{\Gamma }(s-s_{\mathrm{a}}{)}^{2}dS(\mathbf{x})\right],\end{array}$$

$$\begin{array}{r}\hat{b}=\frac{b}{{\tilde{k}}_3},{\hat{\gamma }}_0=\frac{L_0^2}{{\tilde{k}}_3}{\gamma }_0,{\hat{k}}_i=\frac{{\tilde{k}}_i}{{\tilde{k}}_3},i=1,2,{\hat{a}}_0=\frac{a_0{L}_0}{{\tilde{k}}_3}=\frac{100L_0}{R_c},\end{array}$$

$$k_3=K_{B}T{L}_p{m}_0.$$

$$\begin{array}{r}0.1\le \hat{b}\le 5.0,0.1\le {\hat{\gamma }}_0\le 20,\end{array}$$

$$\begin{array}{rl}{\mathbf{V}}_{\Gamma }& =\mathbf{\nu },\text{ on }\Gamma ,\text{ and smoothly decays to zero inside }\Omega .\end{array}$$

$$\begin{array}{rl}\text{outer pressure }& =\frac{1}{L_0^2|\Gamma |}{\tilde{k}}_3{\delta }_{\Omega }\hat{\mathcal{E}}[s,\mathbf{n},\Omega ;{\mathbf{V}}_{\Gamma }],\end{array}$$

Full text

The authors are very grateful to this support granted to the project by the National Science Foundation through grants NSF-DMS-1817156(Arsuaga and Vazquez), NSF-DMS-1816740 (Calderer) and NSF-DMS-1555222-CAREER (Walker)

A Introduction

The discovery of liquid crystal phases formed by DNA in free solution as well as inside viral capsids led to the chromonic denomination of a class of lyotropic liquid crystals. Different from liquid crystals consisting of elongated, rod-like molecules found, for instance, in display devices, chromonic liquid crystals consist of disk, plank-like molecules that form liquid crystal phases with varying concentration. The molecules stack to form cylindrical aggregates, and at even higher concentration, the cylinders cluster so their center axes align in a hexagonal lattice, forming what is known as the hexagonal columnar liquid crystal phase. Experimental and theoretical studies acquired over the last 30 years [22, 23, 25, 37, 39, 42, 36] have shown that encapsidated DNA molecules form a columnar hexagonal liquid crystal phase. In particular, liquid crystalline phases in bacteriophages were first proposed in [16], with an explicit reference to hexagonal packing made in [28] and since then, consistent data have been accumulating [21, 20, 30, 37].

This article studies the liquid crystal ordering of viral dsDNA in a capsid with axial symmetry. We use advanced liquid crystal models and their state-of-the-art numerical implementation to determine the equilibrium structure of the DNA in the capsid and the forces that it sustains. The model is also aimed to be used as a tool for the design and prediction of viral encapsidated DNA in response to the renewed interest in bacteriophages in medicine and biotechnology. This goal guides our parameter selection strategy which is also a main component of the work. The proposed method delivers a parameter phase space for the pressure and the the size of the inner disordered core.

Icosahedral bacteriophages consist of a protein capsid whose assembly is followed by the packing, by means of a molecular motor, of a single naked dsDNA molecule [41]. The DNA molecule inside the viral capsid is found under extreme concentration, between $200$ and $800$ mg/ml [16], and pressure estimated in the range between 40 and 60 atmospheres [10, 13]. Three factors contribute to the excess pressure found inside the viral capsid: the decrease in entropy associated with the confinement imposed by the capsid, the high resistance of the DNA molecule to bending beyond its persistence length and the self-repulsion of the DNA molecule [38]. At the time of infection the DNA is released by a mechanism that suggests a phase transition, possibly into a ’liquid-like’ state [23, 27, 40]. Both processes, packing and releasing of the genome depend greatly on how the DNA molecule folds within the viral capsid; however our understanding of this folding remains very limited.

We develop an effective mechanical model of DNA encapsidation based on the theory of hexagonal chromonic liquid crystals, combined with that of liquid crystals with variable degree of orientation. We characterize the hexagonal columnar phase in terms of a triple of orthonormal vectors, $\mathbf{n}$, $\mathbf{m}$ and $\mathbf{p}$, the first giving the local direction of the columnar axis, with $\mathbf{m}$ and $\mathbf{p}$ denoting lattice vectors of the orthogonal cross sections. The unit vector $\mathbf{n}$ also represents the tangent vector to the DNA center curve at a point, with $\mathbf{m}$ and $\mathbf{p}$ representing the directions along which neighboring strands are found in the encapsidated arrangement. Moreover, we include the scalar order parameter $s\in(0,1)$ that describes the degree of orientation of filaments in the capsid [9]. The value $s=1$ corresponds to (ideal) perfect alignment and $s=0$ to the isotropic state, with no alignment ordering, so it naturally lends itself to quantitatively distinguishing between the region of the capsid with ordered DNA from the disordered central core. We point out that $\mathbf{n},\mathbf{m}$ and $\mathbf{p}$ are undefined at points where $s=0$. The order parameter in the model also yields a natural realization of the phase transition from the ordered to the liquid-like state. Furthermore, the model allows for the capsid domain to be of icosahedral shape, since $s=0$ at the edges where the vector fields are not defined. Ordering of the DNA molecule at the boundary as promoted by the capsid [7, 8, 25, 2] motivates our selection of the boundary conditions for the fields of the problem.

In addition to the genome length $L$, and the volume $V$ of the capsid, data at our disposal include cry-EM images of the bacteriophage yielding the size and shape of the inner disordered region, the DNA density graphs that dictate the number of concentric layers [4, 3, 8, 14] giving, in turn, the effective diameter $d$ of the DNA. We take the customary point of view that the DNA molecule is a semiflexible elastic polymer characterized by the persistence length $L_{p}$ which is about the size of the radius of the capsid. In the calculation of the bending modulus $k_{3}$, we also assume that the hydrated DNA fills the entire volume of the capsid [16], justifying our selection of the maximum linear density $m_{0}$.

The repulsive energy between different DNA strands is not explicitly present in the model but is accounted for by the effective role of the parameters of the mechanical energy. We construct the latter in two different stages. We first consider the total energy of encapsidation consisting of the bending energy penalty of the DNA filament formulated in terms of the director field $\mathbf{n}$, together with a normalized potential $w(s)$ that attains its minimum value at the ordered configuration, labeled as $s=s_{a}\approx 1$, and the maximum at $s=0$, reflecting the preference of the DNA to organize inside the capsid. A penalty is also associated with changes in $s$ resulting in a transition region between the coiled state and the melted core. However, under parameter choices that provide good agreement between the values of experimental and the calculated disordered core sizes, we find that, the pressure is underpredicted by a factor of 100 with respect to experimentally measured values Finally, pressure measurements for different viruses are also available [5, 10, 11, 22, 43]. This leads us to include an additional term in the energy penalizing distortions of the cross section formulated as quadratic expressions of the gradients of the cross-sectional vectors $\nabla\mathbf{m}$ and $\nabla\mathbf{p}$. The choice follows from the elastic energy of hexagonal columnar phases proposed by de Gennes. It is the bulk and shear modulus, B and C, present in the transverse elastic energy, that effectively account for the cohesive effects postulated in earlier models of DNA encapsidation [17], [43] that, in turn, effectively account for the electrostatic interactions.

We find that the dimensionless energy restricted to spooling configurations involves three material constants, $\hat{b},\hat{\gamma}_{0}$ and $\alpha$, in addition to $k_{3}$ that can be directly calculated from DNA data (Table 1). We determine the orders of magnitude of $\hat{b}$ and $\hat{\gamma}_{0}$ from the size of the disordered core of the T4 virus, as measured by imaging. Likewise, the parameter $\alpha,$ used in scaling the elastic moduli $B$ and $C$ with respect to $k_{3}$, is estimated from experimental pressure values corresponding to the same virus. The numerical simulations deliver phase spaces for the radius of the inner disordered core and the pressure, in terms of the parameters $b$ and $\gamma_{0}$ (Figures 2-4).

Earlier mathematical results on Ericksen’s model guarantee existence of a global minimizer of the total energy [1]. In addition, numerical methods previously developed by Walker, et. al. [33, 34, 46] on such a model extend to the present case and provide a convenient algorithm to predict packing configurations and the pressure values that they sustain.

Sine cryoelectron microscopy (cryo-EM) images for most bacteriophages show multilayered spooling-like configurations on the outer layers of the packed genome [7, 2, 4, 3, 8, 22, 14], we present numerical simulations of minimizing configurations with spooling geometry.

In section A, we give a brief review of the nematic and chromonic models of liquid crystals, to explain and motivate our methods, with a special emphasis on Ericksen’s model that includes the variable degree of orientation.

In section B, we present the energy of packing, its scaling properties and the role of the spooling configurations. In section C, we present our approach for estimating parameters, concluding with numerical simulations that relate remaining phenomenological parameters to the radius of the disordered core. Section D uses shape optimization tools to calculate the pressures near the capsid surface.

This work is motivated by earlier studies of free boundary problems of chromonic liquid crystals and bacteriophages carried out in the Ph.D dissertation by Lindsey Hiltner [12].

A.1 Nematic liquid crystals

In the Oseen-Frank theory, equilibrium states of a nematic liquid crystal are minimizers of the Oseen-Frank energy given by

[TABLE]

The splay, twist and bend material constants, $k_{1}>0$, $k_{2}>0$ and $k_{3}>0$, respectively, are associated with the liquid crystal deformations that embody the corresponding geometries. The splay-bend term multiplying $k_{2}+k_{4}$ is a null Lagrangian and so, it depends on boundary values only. Defects in the Oseen-Frank theory correspond to locations, points, lines or surfaces, where $\mathbf{n}=0$, that is where the molecular order is lost, giving rise to locally disordered structures, that may have infinite energy.

In order to represent liquid crystal configurations where ordered and disordered states coexist as well as assign finite energy to defects, it is convenient to introduce the order parameter $s$ in the model, and the energy associated with changes of order. For this, we will adopt Ericksen’s model (1991) of liquid crystals with variable degree of orientation. It may be viewed as an extension of the Oseen-Frank model; alternatively, it may also be interpreted as a restriction of that of Landau-deGennes to uniaxial states. In this model, the equilibrium state of the liquid crystal material is described by the pair $(s,\mathbf{n})$ minimizing a bulk-energy functional

[TABLE]

with $\mathcal{W}_{1}$ as in (2). The material constants $b>0$ and $\gamma_{0}>0$ are of the same of the general one, that contains coupling terms of $\nabla\mathbf{n}$ and $\nabla s$. (A full description of the model can be found in [29] and its analysis can be found in [34, 33, 46]).

The double well function $w(s)$ is of class $C^{2}$, defined on $-1/2<s<1$, that penalizes the phase transition between the nematic and the isotropic phases. It is chosen so that it has a minimum at $s=s_{a}$, $s_{a}$ being the preferred order parameter value at a given temperature, or concentration; it has a maximum at $s=0$ indicating the ordering preference of the DNA. We list its properties as follows [9, 1, 26]:

$\lim_{s\rightarrow 1}w(s)=\lim_{s\rightarrow-1/2}w(s)=\infty$, 2. 2.

$w(0)>w(s_{\mathrm{a}})=\min_{s\in[-1/2,1]}w(s)=0$ for some $s_{\mathrm{a}}\in(0,1)$, 3. 3.

$w^{\prime}(0)=0$.

Next, we discuss the approach to imposing boundary conditions. Rather than prescribing pointwise values of the fields $(\mathbf{n},s)$ on the boundary, that is, imposing strong anchoring boundary conditions, we appeal to weak anchoring, which consists in penalizing the energy of departure of the boundary values from the ideal ones [44]:

[TABLE]

where $\bm{\nu}$ is the oriented unit normal vector of $\Gamma:=\partial{\Omega}$. The first term $a_{\perp}$ tends to align the minimizing director field $\mathbf{n}$ perpendicular to $\bm{\nu}$, and the second term $a_{\parallel}$ favors a parallel alignment. The last contribution to the energy is needed to ensure that $s$ does not trivially vanish on the interface, and so cause $E_{\mathrm{a},\mathbf{n}}(s,\mathbf{n})$ to vanish as well [34, 6, 32]. The parameters $a_{\perp},a_{\parallel},a_{\mathrm{ori}}:\Gamma\to[0,\infty)$ may vary (in general), but here we take them to be fixed constants.

In determining a configuration of a liquid crystal, one minimizes the total energy consisting of (3) added to (4). A thorough study of existence of minimizer together with the numerical implementation of the minimizing algorithm can be found in [34, 33, 46]. In particular, the theory allows for $\Omega$ to be of icosahedral shape, since $s=0$ at the edges where the vector fields are not defined.

A.2 Chromonic Liquid Crystals

We are concerned with hexagonal columnar phases of chromonic liquid crystals. The energy consists of the Oseen-Frank contribution (or, the Ericksen’s model analog) together with an elastic energy of the cross-section perpendicular to the columnar director field $\mathbf{n}$. The latter also effectively includes the electrostatic repulsive energy of the negatively charged DNA strands. We then represent the energy as

[TABLE]

where ${\mathbf{u}}$ represents the displacement vector associated with deformation of the cross section. The minimization of the chromonic energy is done subject to the constraints

[TABLE]

The splay of a liquid crystal is zero whenever dislocations do not occur, that is, the same number of filaments that enter a unit area exit that cross-section. In the case of the hexagonal columnar phase, nonzero splay would allow for deviations from the lattice structure. Twist is prohibited because of its incompatibility with the two-dimensional lattice order in planes perpendicular to the director.

We conclude this section with the observation that the previously models described correspond to a continuum and do not reflect the discrete nature of the system made of filaments. These do not overlap except, perhaps, in possible knot locations. In order to associate a minimizing set $\{\mathbf{n},\mathbf{m},\mathbf{p},s\}$ to a DNA configuration, given the inter-strang DNA spacing taken as its effective diameter, $d>0$, we further introduce the scalar fields $\omega,\vartheta$ as solutions of the differential equations

[TABLE]

Then, the intersections of level surfaces $\omega=\text{constant}$, $\vartheta=\text{constant}$ give the location of the DNA strands. Furthermore, the piecewise smooth curve $\mathbf{r}=\mathbf{r}(\zeta),\zeta\in[0,L]$, such that

[TABLE]

describes the axis of the DNA. Here $L>0$ represents the length of the genome and $\mathbf{r}_{0}$ its entry point in the capsid. Restricting the elastic energy in (5) to discrete deformations of the level curves sustaining the DNA filaments, we arrive at the expression

[TABLE]

where $B$ and $C$ correspond to bulk and shear moduli, respectively.

B DNA encapsidation

We assume that the capsid corresponds to a bounded, axisymmetric region $\Omega\in\mathbb{R}^{3}$ with $\partial\Omega$ representing the piecewise smooth, faceted, viral capsid. We let $l_{0}>0$ denote the length of the axis that connects the location of the connector with its antipodal site in the viral capsid. Assuming the DNA is tightly packed inside the capsid, we adopt a continuum approach for representing the configuration of the DNA [17, 18]. To this end, we let $\mathbf{n}$ represent the unit tangent vector of an “ensemble” of packed DNA strands. Moreover, we use the $s$ variable (degree-of-orientation) to represent how well ordered the DNA strands are, e.g. $s=0$ indicates that the DNA strands are oriented in all directions equally. Furthermore, the unit vectors $\{\mathbf{m},\mathbf{p}\}$ span the cross sections perpendicular to the DNA axis, with its associated level curves representing strand locations.

B.1 Total energy

The model that we propose consists of the total free energy together with a list of constraints, expressed as follows:

[TABLE]

We make the following choices for the parameters of the model:

[TABLE]

where $\alpha>0$ in (13) is a scaling parameter latter determined using pressure measurements. The choices for $\{k_{i}\}_{i=1}^{4}$ indicate that $\mathrm{div}\,\mathbf{n}$ and $\mathbf{n}\cdot\mathrm{curl}\,\mathbf{n}$ are almost negligible for the DNA coiling configurations that we consider here. Dimensional analysis as well as information derived from Onsager’s theory of rigid rods will further assist us in specifying relations (13). The selection of the anchoring constants (14) enforces tangential (planar) boundary conditions for $\mathbf{n}$, as this is consistent with DNA coiling inside capsids.

The bulk potential function $w(s):[0,1]\to[0,\infty)$ is modified to satisfy:

[TABLE]

The optimal degree of orientation is chosen to be $s_{\mathrm{a}}:=0.8$.

The total energy of the DNA coiling system is

[TABLE]

B.2 Model simplification: spooling configurations

From now on, we consider axisymmetric capsids and spooling configurations. Adopting cylindrical coordinates $(r,\theta,z)$, with the $z$-axis along the capsid, spooling DNA arrangements are characterized by $\mathbf{p}=\mathbf{e}_{z}$. Furthermore, for such configurations, it is natural to set $\mathbf{\nu}={\bm{e}}_{r}$ and $\mathbf{n}={\bm{e}}_{\theta}$ and $s=s_{a}$ on the boundary. Since $\nabla\mathbf{p}=0$ and $\mathbf{m}=\mathbf{n}\times\mathbf{p}$, so $|\nabla\mathbf{m}|=|\nabla\mathbf{n}|$. We now define ${\tilde{k}}_{i}=k_{i}+A,\,\,i=1,2,3,$ with $A$ as in (13). In particular, note that

[TABLE]

Finally, using the identity:

[TABLE]

we obtain

[TABLE]

with

[TABLE]

B.3 Dimensional Analysis

Introducing a length scale $L_{0}$ to be of the same order as $R_{c}$, and noting the following scalings:

[TABLE]

where $~{}\widehat{}~{}$ denotes non-dimensional quantities, (18) becomes

[TABLE]

where

[TABLE]

and we have dropped $~{}\widehat{}~{}$ on all variables for convenience (the domains are now non-dimensional).

C Parameter estimates and numerical simulations

As for semiflexible polymers and its application to the DNA molecule, following Tzlil et al. [43] and [17], we take the bending constant $k_{3}$ as

[TABLE]

The DNA cross-sectional length density, $m_{0}$, is estimated by taking the ratio of the length of a DNA strand to the volume of the container in which it is packed. Thus, for a perfectly packed DNA strand (i.e. no empty space in the volume), the maximum value is $m_{0}=1/(\pi(d/2)^{2})=1.85131\times 10^{17}$ m*-2*, where we used an average DNA diameter of $d=2.6225$ nanometers from Table 1. More realistically, for an average bacteriophage genome (length $30.0624$ micrometers) and average capsid head (radius $34.11$ nanometers), we get an average $m_{0}=1.809173\times 10^{17}$ m*-2*. Furthermore, taking an average persistence length of $L_{p}=54.49$ nanometers, our effective bending constant is $k_{3}=4.058027\times 10^{-11}$ J/m (assuming room temperature).

It is convenient to distinguish three types of parameters in the energy. First of all, $\tilde{k}_{3}$ appears as a multiplicative constant whose role will be adjusting the pressure. The ratios $\tilde{k}_{1}/\tilde{k}_{3},\tilde{k}_{2}/\tilde{k}_{3}$ should be interpreted as having large values to penalize the constraints (6). The ratio $\hat{a}_{0}=a_{0}L_{0}/\tilde{k}_{3}$ is used to enforce the planar boundary conditions, so it also has a large penalty value.

For computations, we take $\tilde{k}_{1}=\tilde{k}_{2}=10\tilde{k}_{3}$, where $\tilde{k}_{3}=k_{3}(1+\alpha)$ (recall 13) and $\alpha>0$ is determined by the pressure measurements (see Section D). The length scale $L_{0}=R_{c}/0.4$, where $0.4$ is the non-dimensional radius of the computational capsid domain; thus, $\hat{a}_{0}=250$. Moreover we explore the parameter range for $\hat{b}$ and $\hat{\gamma}_{0}$ given by:

[TABLE]

and compare the predicted cores sizes with the available virus measurements shown in Table 1. The outcome from the numerical simulations is summarized in Figure 2. Figure 3 shows how the internal “core radius” was estimated for a particular choice for $b$ and $\gamma_{0}$. Figure 4 shows how the capsid pressure, and the core radius, varies with $b$ and $\gamma_{0}$.

D Shape optimization

The minimum energy state of the packing problem is the one that minimizes $\mathcal{E}$ in (20) with respect to $s$ and $\mathbf{n}$. With this, one can compute pressures in the capsid by perturbing $\mathcal{E}$ with respect to the domain ${\Omega}$, i.e. take shape derivatives [45].

More specifically, given a (non-dimensional) perturbation field $\mathbf{V}$ that deforms ${\Omega}$, the shape variation $\delta_{{\Omega}}\widehat{\mathcal{E}}[s,\mathbf{n},{\Omega};\mathbf{V}]$ gives the (non-dimensional) change in energy due to a change in the volume. For example, let $\mathbf{V}_{\Gamma}$ be a smooth vector field such that:

[TABLE]

Then, the dimensional change in energy due to expanding the outer boundary $\Gamma$ is $\tilde{k}_{3}\delta_{{\Omega}}\widehat{\mathcal{E}}[s,\mathbf{n},{\Omega};\mathbf{V}_{\Gamma}]$ which has units of force: $J/m=N$. Thus, the average (dimensional) pressure being exerted on the boundary of the capsid $\Gamma$ is

[TABLE]

where we note that $|\Gamma|$ is the non-dimensional surface area of the capsid. Note that $\mathbf{V}_{\Gamma}$ is normalized so that $|\mathbf{V}_{\Gamma}|=1$ on $\Gamma$. The particular form of the extension in (24) does not matter because ultimately the shape derivative only depends on what happens at the boundary $\Gamma$ by the Hadamard-Zolésio structure theorem [45]. The graphs of the pressure in terms of the parameters $b$ and $\gamma_{0}$ are shown in Figure 4.

E Conclusions

We constructed a model of encapsidation of viral DNA, that combines that of hexagonal chromonic liquid crystals with the theory of liquid crystals with variable degree of orientation. The model is accompanied by a rigorous, finite-element based, algorithm of energy minimization, with design and prediction capabilities that may prove to be a real asset in applications. Due to the presence of the order parameter, the algorithm is able to deal with very general faceted capsid domains, including the geometry of the packing motor features, and also allow for the presence of knots in the DNA. The outcome of the model is a phase space for the radius of the disordered core and the pressure near the boundary of the capsid, with the potential to aid in the design of viral DNA to meet desired specifications. The model also allows for the exploration of large sets of parameters, reiterating its value as a design tool. Follow up work addresses the presence of ions and the dynamics of DNA ejection from the capsid.