# The Spanning Tree Model for the Assembly Kinetics of RNA Viruses.

Inbal Mizrahi  Department of Physics and Astronomy, University of California,
Los Angeles, CA 90095  Robijn Bruinsma  Department of Physics and Astronomy,
University of California, Los Angeles, CA 90095  Department of Chemistry and
Biochemistry, University of California, Los Angeles, CA 90095  Joseph Rudnick
Department of Physics and Astronomy, University of California, Los Angeles, CA
90095

###### Abstract

We present a simple kinetic model for the assembly of small single-stranded
RNA viruses that can be used to carry out analytical packaging contests
between different types of RNA molecules. The RNA selection mechanism is
purely kinetic and based on small differences between the assembly energy
profiles. RNA molecules that win these packaging contests are characterized by
having a minimum Maximum Ladder Distance and a maximum Wrapping Number. The
former is a topological invariant that measures the “branchiness” of the
genome molecule while the latter measures the ability of the genome molecule
to maximally associate with the capsid proteins. The model can also be used
study the applicability of the theory of nucleation and growth to viral
assembly, which breaks down with increasing strength of the RNA-protein
interaction.

##  I Introduction

Most viruses have a rod-like or sphere-like shape. Crick and Watson [1] noted
that the volume enclosed by the protein shell (or “capsid”) of a sphere-like
virus determines the length of close-packed viral genome molecules. In a
similar way, the length of a rod-like capsid scales with the length of the
enclosed genome. They exploited the geometrical connection between capsid and
genome to relate the size of a virus to the number of encoded genes, and
predicted that spherical viruses should have icosahedral symmetry. This
approach works well for many smaller single-stranded (ss) RNA viruses (such as
the polio and common cold viruses) as well as for some larger double-stranded
DNA viruses (such as the herpes and bacteriophage viruses) 111For double-
stranded DNA viruses the genome molecules either are inserted into the
interior pre-fabricated spherical capsids by a molecular motor or the genome
is pre-condensed prior to assembly.. The genome molecules of the very large
retroviruses, like HIV-1, are not close-packed but the dimensions of the
capsid still scale with that of the enclosed genome because the genome
molecules adhere to the inner surface area of the capsid so the surface area
is proportional to the total length of the genome molecules [3].

Are there other relations between the geometry of that capsid and that of the
genome molecules? The co-assembly of small single-stranded RNA viruses
involves both the capsid proteins and the viral RNA molecules. Such viruses
often will self-assemble spontaneously in solutions that contain higher
concentrations of viral capsid proteins and RNA molecules [4, 5]. The process
is passive and driven by free energy minimization. Early work by Aaron Klug
[6] indicated that the RNA genome molecules act as templates that direct the
viral assembly process. He proposed a physical model in which repulsive
electrostatic interactions between capsid proteins are just strong enough to
prevent spontaneous capsid assembly driven by hydrophobic attractive
interactions between capsid proteins but if viral RNA molecules are present as
well then the negative charges of RNA molecules can neutralize some of the
positive capsid protein charges and thereby tilt the free energy balance
towards assembly 222For a quantitative treatment of this model, see ref.
[kegel2006].

The genome molecules of ss RNA viruses genome molecules have a tree-like
“secondary structure” produced by Watson-Crick base-pairing between
complementary RNA nucleotides of the primary sequence of RNA nucleotides [8].
The co-assembly picture suggests that there could be connections between the
geometric and topological features of the viral RNA secondary structure
directing assembly and the capsid geometry. The redundancy of the genetic code
indeed allows for the possibility of “silent” (or synonymous) mutations that
can alter the secondary structure of an ss RNA molecule without altering the
structure of the proteins encoded in the nucleotide sequence [9]. A systematic
comparison between the secondary structures of viral RNA molecules and
randomized versions of the same molecules indeed revealed that they are
significantly more branched and compact than generic RNA molecules [10].

The assembly of empty viral capsids is similar to that of micelles and other
amphipathic molecular systems that self-assemble spontaneously [11]. Empty
capsid assembly may be viewed as a chemical reaction where a certain number of
capsid proteins react to form the capsid. Application of the Law of Mass
Action (LMA) for reactions in chemical equilibrium produces a relation between
the concentrations of assembled capsids and of free capsid proteins with the
total protein concentration. This relationship has been confirmed
experimentally [12] (but see ref.[13]). An important feature of the LMA is the
presence of a critical aggregation concentration (or CAC). For the case of
empty capsid assembly, this means that capsids only can form if the total
capsid protein concentration exceeds this CAC. Another important feature is
that assembled capsids should spontaneously disassemble if the concentration
of capsid proteins in the surrounding solution is reduced to zero. This does
not occur in actuality—at least not on laboratory time scales—presumably
because of the presence of an activation energy barrier inhibiting disassembly
that is large compared to the thermal energy $k_{B}T$.

Kinetic studies of empty capsid assembly [14, 15, 16] report that there is a
delay or lag time before capsid assembly can take place in a solution of
capsid proteins that has been primed for assembly. Assembly starts with the
formation of a “nucleation complex” composed of a small number of capsid
proteins. The nucleation complex then extends or elongates by absorbing more
proteins until the capsid closes up. The theory of nucleation and growth [17]
has become an important tool for the interpretation of kinetic studies, with
the nucleation complex corresponding to the critical nucleus. The small size
of the critical nucleus means that in-vitro assembly assays are taking place
under conditions of a significant level of supersaturation and, based on
kinetic considerations, the same is probably true under in-vivo conditions. A
study of the assembly kinetics of empty capsids of the Hepatitis B virus
indicates that the assembly of such capsids is governed by distinct assembly
pathways [18], much like distinct assembly pathways govern the folding of
proteins [19]. While all-atoms simulations of capsids are possible [20],
coarse-grained models have been found to be useful to interpret equilibrium
and kinetic properties of empty capsid assembly[21, 22, 23, 24, 25, 26, 27,
13, 28, 29].

The co-assembly of capsid proteins with ss RNA molecules involves additional
thermodynamic parameters. One of these is the mixing ratio of the
concentration of RNA molecules to that of the capsid proteins. An in-vitro
study [30] of the co-assembly of CCMV (Cowpea Chlorotic Mottle Virus) with
viral RNA molecules reported that when the RNA-to-protein mixing ratio is low
then virus-like particles co-exist with excess proteins. By analogy with
formation reactions of a binary compound, one would expect that in the case of
large mixing ratios there is co-existence of virus-like particles with excess
RNA molecules. The border-line mixing ratio separating these two regimes
should correspond to the RNA-to-protein ratio of an assembled virus-like
particle, which can be thought of as a stoichiometric ratio. In actuality,
when the mixing ratio exceeded a certain threshold then a distribution of
disordered RNA-protein aggregates was observed [30]. Moreover, the border-line
mixing ratio was not the stoichiometric ratio. A second thermodynamic
parameter found to be important for co-assembly is the ratio of the RNA to
protein affinity with the protein to protein affinity. This affinity ratio can
be altered experimentally by changing the acidity and salinity of the solution
[31]. It was found that for increased values of this ratio, virus-like
particles were replaced by disordered RNA-protein aggregates.

Information about co-assembly also has been gleaned from structural studies.
Until recently, reconstruction of packaged genome molecules involved
“icosahedral averaging”, which resulted in RNA structures with imposed
icosahedral symmetry [32]. Such studies showed that the interior surface of
the icosahedral capsids of certain viruses (e.g., the nodaviruses [33]) is
decorated by paired RNA strands lining the edges of the “capsomers” (i.e.,
pentameric or hexameric groupings of capsid proteins)). Recent progress in
cryo-electron tomography has made it possible to image individual ss RNA
genome molecules packaged inside spherical capsids without icosahedral
averaging (“asymmetric reconstruction”[34, 35]). One example is the MS2 virus,
which is an ss RNA, quasi-icosahedral bacteriophage virus. It was found that a
subsection of the RNA genome reproducibly associated with a compact cluster of
capsid proteins [36]. This result was interpreted as evidence for a well-
defined assembly pathway where the energetically ”uphill” part of nucleation-
and-growth scenario produces a compact nucleation complex held together by
particular sections of the viral RNA molecule with enhanced affinity for the
capsid proteins (such sections are known as packaging signals [37]). In this
picture, the “downhill” part of the assembly process involves more generic
electrostatic RNA-protein attractions. According to this model, the nucleation
complex plays the important role of selecting the viral RNA molecules from
quite similar host messenger RNA molecules.

In this paper we propose a statistical mechanical model that can be used to
explore the question of the influence of geometric and topological features of
tree-shaped molecules on the promotion of packaging by capsid proteins and how
these features are connected to kinetically favored assembly pathways. It
allows testing of the theory of nucleation and growth. Other questions that
can be addressed are why the “chemical reaction picture” works well for empty
capsid assembly but not so well for co-assembly, whether packaging selectivity
is consistent with a high level of supersaturation.

The proposed model, the“spanning tree model” model, extends a simple model for
the assembly of empty dodecahedral capsids from pentamers in solution, due to
Zlotnick [21, 38], by allowing it to package branched genome molecules that
decorate the edges of the dodecahedral capsid. An important advantage of the
Zlotnick Model is that its assembly kinetics has been studied [25, 13] while
the packaging kinetics of linear genome molecules has been numerically
simulated [39]. Two different packaging scenarios were encountered. The first
scenario is similar to that of the assembly of empty capsids with assembly
intermediates in the form of compact pentamer clusters that grow in a
pentamer-by-pentamer fashion. In the second “en masse” scenario, the first
assembly step is the formation of a disordered pentamer/genome condensate
formed which then undergoes an ordering transition that can be described by
the Landau theory of symmetry-breaking on a spherical surface (see [40] and
references therein).

##  II The Spanning Tree Model.

We start with a brief review of the Zlotnick Model [21].

###  II.1 Zlotnick Model

In the Zlotnick Model, a capsid is modeled as a dodecahedral shell composed of
twelve pentamers. Assembly takes place in a reservoir of pentamers and is
driven by an attractive edge-edge interaction between pentamers. A minimum-
energy assembly pathway is defined as a pentamer-by-pentamer addition sequence
where each added pentamer is placed in a location on a partially assembled
dodecahedral shell that minimizes the total energy. An example of a minimum-
energy assembly pathway is shown in Fig. 1.

Figure 1: Zlotnick Model. The figure shows a minimum-energy pathway for the
assembly of a dodecahedral shell composed of twelve pentamers with adhesive
edges. The edge-edge binding energy is $\epsilon$. The change in energy per
added pentamer is indicated.

The energy $E(n)$ of a cluster of $n$ pentamers is defined as the number of
shared pentamer edges times the binding energy $\epsilon$ per edge minus a
constant $\mu_{0}$ times the number of pentamers. Here, $\mu_{0}$ is the
chemical potential of a pentamer in solution at a certain reference pentamer
concentration $\bar{c}$ 333The assembly of the capsid is assumed here to take
place on a specific location. The quantity $\mu_{0}$ reflects the entropic
free energy cost of removing a pentamer from the solution to this location
plus that of any conformational change that is required for the pentamer prior
to joining a partial capsid. in units of $k_{B}T$. If the concentration of
free pentamers $c_{f}$ differs from $\bar{c}$ then $\mu_{0}$ must be replaced
by $\mu=\mu_{0}+\ln c_{f}/\bar{c}$ but in this section $c_{f}$ equals
$\bar{c}$. From here on, concentrations are expressed in terms of $\bar{c}$
and thus dimensionless. Since a dodecahedron has thirty edges, the assembly
energy of a complete capsid equals $30\epsilon-12\mu_{0}$. If the chemical
potential equals $(5/2)\epsilon$ then the assembly energy of a pentamer that
is part of a capsid is the same as that of a free pentamer in solution (i.e.,
zero). We will refer to $\mu^{*}=(5/2)\epsilon$ as the chemical potential for
assembly equilibrium. Below, we will use energy units in which $\epsilon=-1$
so with $\mu^{*}$ equal to -(5/2). An important feature of the Zlotnick model
is that the energy gain for the initiation of assembly is relatively low,
namely one factor of $\epsilon$, while the cost of removing a single pentamer
from an assembled capsid is relatively high: it requires breaking five bonds
with a total energy cost of minus $5\epsilon$. This allows for metastability
of assembled capsids in solutions with low pentamer concentrations.

Figure 2 (top) shows the minimum-energy assembly pathway corresponding to Fig.
1 for three different values of $\mu_{0}$ close to $\mu^{*}$ and for
$\epsilon=-1$.

Figure 2: Top: Assembly energy profile $E(n)$ for the minimum assembly energy
pathway of the Zlotnick Model shown in Fig. 1. Blue dots: $\mu_{0}$ is
slightly below the chemical potential $-(5/2)$ for assembly equilibrium.
Orange squares: $\mu_{0}$ is equal to $-(5/2)$. Green diamonds: $\mu_{0}$ is
slightly above $-(5/2)$. The edge-to-edge binding energy is equal to minus
one. Bottom: interpretation of the energy profile in terms of nucleation-and-
growth theory. Solid red dots: energy minima corresponding to the assembled
capsid. Solid black dots: energy maxima. For $\mu_{0}<\mu^{*}$, the absolute
energy minimum is at $n=0$ while for $\mu_{0}>\mu^{*}$ the absolute minimum is
at $n=12$, the assembled capsid. The height of the maxima corresponds to the
activation energy barrier. The location $n^{*}$ of the maximum is around $n=6$
in the assembly equilibrium state and shifts to smaller values as $\mu_{0}$
increases.

There is a large number of such minimum energy assembly pathways for the
Zlotnick model (of the order of $10^{5}$) with the energy profile of Fig. 2
(top). Figure 2(bottom) shows how these energy profiles can be interpreted in
terms of the theory of nucleation and growth [17]. The height of the energy
maximum plays the role of the energy activation barrier. The size of the
critical nucleus is the value $n^{*}$ for which $E(n^{*})$ has a maximum. If
$\mu_{0}\simeq\mu^{*}$, then $n^{*}$ is around 6 while $n^{*}$ decreases with
increasing levels of pentamer supersaturation. Note that the Zlotnick Model
effectively incorporates the surface or line tension that plays a key role in
the theory of nucleation and growth.

For the Zlotnick Model, as well as other simple models of capsid assembly, the
critical nucleus under assembly equilibrium conditions ($\mu_{0}=\mu^{*}$) is
a half-formed capsid. Experimentally measured values for the interaction
energies between capsid proteins are in the range of a few $k_{B}T$. If the
critical nuclei really were half-formed capsids then the activation energy
barrier for the assembly of actual capsids would be in the range of hundreds
of $k_{B}T$ under conditions of assembly equilibrium, which would mean
prohibitively slow kinetics. Actual measured values of the size of the
nucleation complex are much smaller [15]. As already noted, this indicates
that capsid assembly takes place under conditions of a high degree of
supersaturation.

###  II.2 Spanning Trees

The second part of the definition of the model is the specification of “toy”
ss RNA genome molecules. These are represented by tree graphs, i.e.,
collections of nodes connected by links such that there is one and only one
path of links connecting any pair of nodes [42] (see Fig. 3 top right). The
construction of a genomic tree graphs starts with a spanning tree graph of the
dodecahedron. A spanning tree graph of a polyhedron is defined as a tree graph
whose nodes are located on the vertices of the polyhedron with just enough
links to connect the nodes together in a tree structure without circuits [43].
Each member of the set of spanning trees of the dodecahedron has the same
number of vertices (twenty) and the same number of links (nineteen). Figure 3
shows three ways to represent the same spanning tree graph.

Figure 3: Top left: Spanning tree connecting the vertices of a dodecahedron
(solid lines). The dashed lines indicate edges of the dodecahedron that are
not part of the spanning tree of specific links. Six pentamers can be placed
on the dodecahedron with each one wrapped by the spanning tree with four links
per pentamer (numbered). Top right: Planar graph of the same tree. Bottom: The
spanning tree projected on a planar Schlegel graph of the dodecahedron (red).
The black arrows represent side branches that are added to the spanning tree
so that all edges of the dodecahedron are covered by a link of the tree.

The top left figure shows a three dimensional representation with the solid
black lines indicating the spanning tree. Note that the spanning tree covers
only 19 of the 30 edges of the dodecahedron. The top right figure is a two
dimensional representation of a tree graph with the dodecahedron removed. In
the bottom picture, the dodecahedron is represented in the form of a planar
Schlegel diagram [44] that is decorated by the spanning tree graph (solid red
lines). The final step of the construction of a genome molecule is the
extension of a spanning tree by adding side-branches so the tree covers the
remaining eleven edges. These are represented by black arrows in Fig. 3. The
completed genomic tree molecule covers all edges of the dodecahedron. The red
spanning tree links will represent the packaging signals that have a specific
binding affinity to pentamer edges that is enhanced with respect to generic
electrostatics while the black arrows represent side branches with only non-
specific generic affinity for pentamer edges.

###  II.3 Classification Indices for Spanning Tree Molecules.

The number of unique spanning trees 444The set of unique spanning trees
consists of all spanning trees that, when depicted as in the top left of Fig.
3, cannot be mapped into each other by rotations and/or reflections that leave
the dodecahedron invariant.  of the dodecahedron is on the order of $10^{5}$
and we need to classify them. A “global” characteristic that has been applied
to classify RNA secondary structures is the Maximum Ladder Distance (or MLD)
[10, 46]. This is the maximum number of paired RNA nucleotides separating any
two nucleotides. The MLD of a secondary structure is a global measure of its
size [10, 46]. Specifically, in the absence of interaction between nodes, the
solution radius of gyration of the molecule scales with the MLD as a power
law. The MLD is a topological invariant that does not change if the molecule
is folded up in different ways. A systematic comparison between the genomic
RNA molecules of RNA viruses revealed that they have significantly lower MLDs
than randomized versions of the same molecules [10, 46]. The analog of the MLD
for the toy genome molecules is the maximum number of links separating any
pair of nodes. In graph theory, the ladder distance between two nodes of a
tree graph is called “the” distance while the MLD is known as the “diameter”
of the tree graph [42]. The MLD of the tree molecule shown in Fig. 3 is nine.
It can be demonstrated that the smallest possible MLD for a spanning tree of
the dodecahedron is nine (as shown in Appendix A), while the largest possible
MLD of a spanning tree is nineteen. The former resemble a Bethe lattice while
the latter is a Hamiltonian Path. A walk that visits all vertices of a
polyhedron is said to trace a Hamiltonian Path [47, 48]. Figure 4 is a plot of
the number N of spanning trees of the dodecahedron as a function of the MLD.

Figure 4: The number of spanning trees on the dodecahedron as function of the
maximum ladder distance (MLD). Two spanning trees of the dodecahedron that are
related by a symmetry operation of the dodecahedron are treated as the same.

The plot has a pronounced maximum around MLD equal to 12. By comparison, the
configurational entropy of an annealed branched polymer composed of 19
monomers that is not constrained to be a spanning tree depends on the MLD as
$19-MLD^{2}/19$ [49] and thus has a maximum at the smallest possible MLD.
Demanding that a tree molecule with a certain number of links is also a
spanning tree of a dodecahedron constrains significantly the branching
statistics.

A second characteristic, complementary to the MLD, is the wrapping number (or
$``N_{P}"$). The wrapping number of a spanning tree of the dodecahedron is the
maximum number of pentamers that can be placed on the dodecahedron such that
all pentamers have four edges covered by a link of the spanning tree (four—not
five—is the maximum number of specific links that can be associated with a
pentamer). For the spanning tree shown in Fig. 3, $N_{P}$ equals six. The
maximum $N_{P}$ for a spanning tree of the dodecahedron is eight while the
minimum is two. The distribution of wrapping numbers is shown in Fig. 5:

Figure 5: The number of spanning trees on the dodecahedron as a function of
the wrapping number. Two spanning trees of the dodecahedron that are related
by a symmetry operation of the dodecahedron are treated as the same.

The wrapping number distribution has a maximum at $N_{P}=5$.

Unlike the MLD, the wrapping number of a tree molecule is not an invariant of
the tree topology. Instead, it depends on the spatial configuration of a tree
molecule distributed over the edges of a dodecahedron. An example is a linear
spanning tree as shown in Fig. 6.

Figure 6: A possible Hamiltonian path for a linear genome molecule. Only two
of the pentamers are maximally wrapped. They have no shared edge.

The Hamiltonian path of Fig. 6 has a wrapping number of two. After enumerating
all possible Hamiltonian paths, one finds that their wrapping number can be
two, three, and four. In general, a spanning tree molecule can be placed in
different ways on the edges of a dodecahedron and these different
configurations may have different wrapping numbers. Figure 7 shows two
spanning tree configurations with different wrapping numbers that belong to
the same tree molecule.

Figure 7: Top: An MLD=13 tree molecule with a two-fold symmetry site (marked).
Bottom: Two dodecahedral spanning tree configurations of this molecule. Both
configurations retain two-fold symmetry (marked). The configuration on the
left has a wrapping number $N_{P}=2$ while the one on the right has a wrapping
number $N_{P}=6$

Tree structures can have multiple wrapping numbers. The wrapping number is a
geometrical characteristic of the different ways to distribute nineteen
specific links over the edges of a dodecahedron. A network with circuits that
would visit all nodes of a dodecahedron would not have a well-defined MLD but
it would still have a wrapping number.

The wrapping number and MLD are correlated. For example, a spanning tree with
a large wrapping number is expected to be have many branches in order to
maximize the number of pentamers that can be accommodated. Hence it is
expected to have a small MLD (and vice versa). Figure 8 is a plot of the range
of allowed wrapping numbers for given MLD:

Figure 8: Plot of the range of wrapping numbers ($N_{P}$) of the spanning tree
molecules of the dodecahedron as a function of the maximum ladder distance
(MLD).

The largest possible wrapping number is $N_{P}=8$. In that case, the MLD can
only be nine, ten or eleven. For the smallest possible wrapping number
$N_{P}=2$ the MLD ranges from eleven to nineteen.

###  II.4 Assembly Energy Profiles.

We now construct the minimum energy assembly energy profiles for the spanning
tree model. The starting state is a tree molecule, composed of a spanning tree
(the specific links) plus the additional eleven side branches (the non-
specific links), that decorates all edges of a mathematical dodecahedron.
Physically, the starting state can be viewed as representing a folded or pre-
condensed form of the viral ss RNA genome molecule(s) prior to its
encapsidation by pentamers 555In actuality, condensation of the RNA genome
molecules takes place during encapsidation. It is driven by positively charged
polypeptide chains associated with the capsid proteins. Next, place $n$
pentamers on the dodecahedron. The energy $E(n)$ of the assembly is defined to
be
$E(n)=n_{1}\epsilon_{1}+n_{2}\epsilon_{2}+n_{3}\epsilon_{3}+n_{4}\epsilon_{4}-k_{B}T\mu_{0}n$.
Here, $n_{1}$ is the number of specific links of the tree that lie along a
pentamer edge that is not shared with another pentamer, $n_{2}$ is the number
of specific links that lie along a pentamer edge that is shared with another
pentamer, $n_{3}$ is the number of edges shared between two pentamers that are
associated with a non-specific link and $n_{4}=11$ is the number of non-
specific links that lie along a pentamer edge that is not shared with another
pentamer. The associated binding energies are given as $\epsilon_{i}$ with
$i=1,2,3,4$. We will use an energy scale in which the binding energy
$\epsilon_{4}$ of a non-specific link is equal to zero. We also will assume
that the interactions between edges and links are additive. This means that
$\epsilon_{2}$ is given by $\epsilon_{2}=\epsilon_{3}+2\epsilon_{1}$. Finally,
$\mu_{0}$ is again the reference pentamer chemical potential. The assembly
energy of a completed particle is equal to
$19\epsilon_{2}+11\epsilon_{3}-12\mu_{0}$. All spanning trees thus have the
same assembly energy. If $\epsilon_{1}$ is zero so $\epsilon_{2}=\epsilon_{3}$
then the assembly energy profiles are the same as that of the Zlotnick Model
with $\epsilon=\epsilon_{2}=\epsilon_{3}$.

The energy parameters enter in the physics of assembly discussed in the next
sections always in the form of $\beta\epsilon_{i}$ with $\beta=1/k_{B}T$. We
will use dimensionless energy parameters
$\epsilon^{\prime}_{i}=\epsilon_{i}/|\epsilon_{3}|$. In these units
$\epsilon^{\prime}_{3}=-1$ and $\beta^{\prime}=\beta/|\epsilon_{3}|$. That
leaves only $\epsilon^{\prime}_{1}$, $\beta^{\prime}$ and $\mu_{0}$ as three
free dimensionless energy parameters. Physically, $|\epsilon^{\prime}_{1}|$
represents the ratio between the affinities of a specific link with a pentamer
edge and that of pentamer edges with each other $1/\beta^{\prime}$ is the
dimensionless attractive interaction strength between two pentamers in units
of the thermal energy. Below we will drop the primes.

Minimum energy assembly pathways can now be constructed in the same way as
before. The degeneracy of the energy profile Zlotnick Model is greatly
reduced. As an example, Fig. 9 shows a minimum energy assembly profiles for
the $MLD=9$, $N_{P}=6$ tree of Fig. 3 (from here on referred to as “molecule
(1)”) and for the $MLD=19$, $N_{P}=2$ tree of Fig. 6 (from here on referred to
as “molecule (2)”).

Figure 9: Assembly energy profiles for minimum energy assembly pathways for
the $MLD=9$, $N_{P}=5$ spanning tree of Fig. 3 (left) and the $MLD=19$,
$N_{P}=2$ spanning tree of Fig. 6. Energy parameters are $\epsilon_{1}=-0.2$,
$\epsilon_{3}=-1$, and $\mu_{0}=-2.5$. The total assembly energy is indicated
by a dashed red line and the assembly activation barrier by a solid red line.

These two energy profiles are consistent with what is expected from nucleation
and growth theory for a state of elevated supersaturation. Note how similar
they are: the total assembly energies (solid red lines) are of course the same
by construction (about 7.5 in our units) but the heights of the assembly
energy activation barriers also are the same (about 2.5 in our units). The
main difference is that the width of the energy barrier for molecule (2) is
somewhat larger than that of molecule (1). Mathematically, the assembly
process can be viewed as a random walk over energy landscapes of the form
shown in Fig. 9, which suggests that molecule (1) will have somewhat faster
assembly kinetics.

Next, we increased minus $|\epsilon_{1}|$ from $0.2$ to $1.0$ so the
attractive interaction between a tree link and a pentamer edge is equally
strong as that between two pentamer edges. We also decreased the reference
chemical potential from $-2.5$ to $-5.2$ in order to keep the activation
energies roughly in the same range as for $\epsilon=-0.2$. The assembly energy
profiles are shown in Fig. 10.

Figure 10: Assembly energy profiles for molecule (1) (left) and molecule (2)
(right) for energy parameters $\epsilon_{1}=-1.0$, $\epsilon=-1$, and
$\mu_{0}=-5.2$

The total assembly energy is still the same (which is by construction) but the
two profiles now have a quite different appearance. The activation energy
barriers also are quite different The reason is that the minimum energy
placement for the second pentamer for molecule (2) is, for large values of
minus $\epsilon_{1}$ determined by the requirement that it is maximally
wrapped. The two maximally wrapped pentamers of molecule (2) do not have a
shared edge (see Fig. 6). The $n=2$ state of molecule (2) has a higher energy
than that of molecule (1) because it lost the attractive energy between the
edge shared between the two first pentamers of molecule (1). Examining the
assembly energy profiles of molecules with different MLD and $N_{P}$ we found
that, in most cases, genome molecules with the same MLD and $N_{P}$ have the
same, or closely similar, minimum energy assembly profiles and pathways. On
the other hand, genome molecules with the same MLD but different $N_{P}$ can
have quite different profiles and pathways.

Finally, one also can define disassembly energy pathways using the same method
except that now pentamers are removed successively, each time with the minimum
energy disassembly cost. By enumeration we found that, in nearly all cases,
the disassembly path simply reverses the assembly path.

###  II.5 Mechanical Rigidity.

Figure 3 (top left) shows a six pentamer structure placed on the dodecahedron.
Each of these pentamers has the maximum of four contacts with a specific (red)
link per pentamer. For any (negative) value of $\epsilon_{1}$, the energy
$E(6)$ of this state has the lowest value any six-pentamer structure could
have when placed on the dodecahedron for any spanning tree. In other words,
for an $n^{*}=6$ nucleation complex (with $\mu_{0}=\mu^{*}$) the activation
energy barrier could not be any lower. One thus might expect the tree
structure of Fig. 3 to be a particularly strong contender in a packaging
competition experiment. This structure has however another interesting
feature. If one were to (i) relax the constraints that keep the links of the
genome molecule associated with the edges of a mathematical dodecahedron, (ii)
allow tree links to rotate freely around the nodes of the tree and (iii) allow
pentamers to swivel around shared edges then the pentamers still could not
move with respect to each other without breaking pentamer-pentamer bonds. Only
rigid-body translations and rotations would be possible. We will say that
pentamer assemblies with this property are mechanically rigid. The pentamer
structures of the minimum energy assembly pathway of the Zlotnick Model (see
Fig.1) are all rigid past $n=2$.

As a counterexample, the linear genome molecule with 19 edges of Fig. 6 with
wrapping number two can accommodate only two pentamers with a maximum of four
favorable contacts. Placing a third and fourth pentamer on the dodecahedron
along a minimum energy assembly path generates the linear four-pentamer
structure shown in Fig. 11.

Figure 11: Intermediate $n=4$ assembly state of molecule (2) for large values
of $|\epsilon_{1}|$. The structure is mechanically flaccid.

If this structure were released from the dodecahedron and allowed to freely
fluctuate then the four pentamers could freely swivel along shared edges. We
will call such an unstable assembly mechanically flaccid. Assembly
intermediates that are flaccid are expected to be characterized by strong
conformational thermal fluctuations. Whether or not an assembly intermediates
are flaccid is controlled by the parameter $\epsilon_{1}$. We find that in the
limit that $\epsilon_{1}$ is small compared to the other energy parameters,
all assembly intermediates are rigid and basically reproduce the scenario of
the Zlotnick Model. The assembly intermediates of some—though not all—pathways
become flaccid when $|\epsilon_{1}|$ is increased past a critical value of the
order of one. For molecule (2) the transition from rigid to flaccid
intermediates takes place exactly at $\epsilon_{1}=-1$ (it is easy to show
that the critical value of $\epsilon_{1}$ is independent of $\mu_{0}$ and
$\beta$). Another example is the tree in Fig. 7 with $N_{P}=6$ which shows the
flaccid minimum energy state when minus $\epsilon_{1}$ exceeds a critical
value. Transitions from rigid to flaccid as a function of $|\epsilon_{1}|$
become more common for increasing MLD.

##  III Boltzmann Distribution.

The minimum energy assembly pathways can be used to define “low-temperature”
Boltzmann Distributions for the thirteen probabilities $P_{n}$ that a tree
molecule is associated with $n$ pentamers. The low-temperature limit applies
if assembly and disassembly processes are entirely restricted to the minimum
energy pathways discussed earlier. Under those conditions, the Boltzmann
distribution equals

$P_{n}\propto\exp\left(-\beta E(n)+n\ln(c_{f})\right)$ (1)

The last term in the argument of the exponential is obtained by replacing the
reference chemical potential $\mu_{0}$ with the actual chemical potential
$\mu=\mu_{0}+\ln(c_{f})$ where $c_{f}$ is the concentration of free pentamers
in solution (in units of the reference concentration). The validity condition
of the low-temperature approximation is that the quantities
$\beta|\epsilon_{i}|$ are larger than one.

It can be checked that it follows from the Boltzmann distribution that

$\frac{c_{f}^{12}r_{0}}{c_{12}}=K$ (2)

where $K=\exp\beta E(12)$, $r_{0}$ the concentration of free RNA molecules and
where the $c_{n}$ are the concentrations of RNA molecules associated with $n$
pentamers. If $r_{t}$ is the total solution concentration of genome molecules
then $P_{n}=c_{n}/r_{t}$. Equation 2 has the form of the Law of Mass Action
(LMA) for a chemical reaction in which twelve free pentamers in solution
“react” with a genome molecule not associated with pentamers, forming a viral
particle with twelve pentamers. The low temperature Boltzmann distribution Eq.
(1) is thus consistent with a chemical reaction picture. In the context of the
LMA, the quantity $K$ is known as the dissociation constant of the reaction.
This dissociation constant depends on the energy parameters only through the
total assembly energy so it is independent of the MLD and $N_{P}$ numbers.

###  III.1 Conservation Laws

The probability distribution $P_{n}$ is related to conservation laws for the
number of pentamers and genome molecules. Conservation of genome molecules
requires that $r_{0}+\sum_{n=1}^{12}c_{n}=r_{t}$. This is assured if
$\sum_{n=0}^{12}P_{n}^{eq}=1$. The proportionality factor in the definition of
$P_{n}$ is determined by this condition. Next, conservation of pentamer
molecules is assured if $c_{f}=c_{0}\Gamma\thinspace([P_{n}])$ where

$\Gamma\thinspace([P_{n}])\equiv 1-(D/12)\sum_{n=1}^{12}nP_{n}$ (3)

with $D=12r_{t}/c_{0}$, $c_{f}$ the concentration of free pentamers and
$c_{0}$ the total pentamer concentration.

The quantity $D$ is the mixing ratio mentioned in Section I as an important
thermodynamic control parameter. If $D=1$ then there are exactly enough
pentamers to encapsidate all genome molecules so $D=1$ corresponds to the
stoichiometric ratio. Because $c_{f}=c_{0}\Gamma\thinspace([P_{n}])$ depends
on the set of all $P_{n}$, the thirteen Boltzmann equations are in fact
coupled and need to be solved self-consistently.

###  III.2 Self-Assembly

A standard diagnostic plot for self-assembling systems is that of the
concentration of free building blocks and that of the assembled particles as a
function of the total concentration of building blocks [51]. Such a plot was
computed from the self-consistent solution of the Boltzmann distribution Eq.1
for example molecule (1) shown in 3. The energy parameters are the same as for
Fig. 9. Figure 12 shows a plot of the concentrations of free pentamers in
solution and that of pentamers that are part of an assembled particle as a
function of the total pentamer concentration $c_{0}$.

Figure 12: Equilibrium self-assembly diagram for molecule (1) with the
parameter values of Fig. 9. Horizontal axis: total pentamer concentration
$c_{0}$. Vertical axis: either the free pentamer concentration $c_{f}$ (blue)
or the concentration $c_{0}-c_{f}$ of pentamers that are associated with a
tree molecule (ochre). Solid lines: solution of Eq.8

For pentamer concentrations that are small compared to one (in our units),
nearly all pentamers are free in solution so $c_{f}\simeq c_{0}$. As the total
concentration increases, the concentration of free pentamers saturates when
particle assembly starts. In this regime the concentration of pentamers that
are part of an assembled particle increases roughly proportional to the total
concentration. The resulting diagram is typical of that of self-assembling
systems [51].

The occupation probability $P_{12}$ for completed particles and the
probability $P_{0}$ for pentamer-free genome molecules are, for these
parameter values, large compared to the intermediate occupation probabilities
with $1\leq n\leq 11$. If one can neglect these assembly intermediates then
the conservation law for genome molecules reduce to $P_{0}\simeq(1-P_{12})$
and that of pentamers to $c_{f}\simeq c_{0}(1-DP_{12})$. Inserting these two
relations into the LMA equation Eq.2 produces an (approximate) equation for
the concentration $c_{f}$ of free pentamers:

$\left(\frac{c_{f}}{c_{0}}\right)^{12}\left(\frac{D-1+\frac{c_{f}}{c_{0}}}{1-\frac{c_{f}}{c_{0}}}\right)\simeq\left(\frac{K}{c_{0}^{12}}\right)$
(4)

For $\left(\frac{c_{0}^{12}}{K}\right)$ small compared to one, this equation
has a solution with $c_{f}$ close to $c_{0}$:

$\frac{c_{f}}{c_{0}}\simeq 1-\frac{Dc_{0}^{12}}{K}$ (5)

For $\left(\frac{c_{0}^{12}}{K}\right)$ large compared to one and $D$ larger
than one, the equation has a solution with $c_{f}$ independent of $c_{0}$:

${c_{f}}\simeq\left(\frac{K}{D-1}\right)^{1/12}$ (6)

Finally, for $\left(\frac{c_{0}^{12}}{K}\right)$ large compared to one but $D$
less than one, the equation has a solution with $c_{f}$ independent of
$c_{0}$:

$\frac{c_{f}}{c_{0}}\simeq 1-D+\frac{K}{c_{0}^{12}}\frac{D}{(1-D)^{1/2}}$ (7)

There is thus a change in regimes when the depletion factor is equal to one.
For the special case that $D=1$, the LMA equation reduces to

$\left(\frac{c_{f}}{c_{0}}\right)^{13}\left(\frac{1}{1-\frac{c_{f}}{c_{0}}}\right)\simeq\left(\frac{K}{c_{0}^{12}}\right)$
(8)

The solid lines in Fig. 12 for $K$ equal to $\exp\beta\Delta E$ with $\Delta
E=13.6$ are solutions of this equation. There is reasonable agreement with the
concentrations computed from the full Boltzmann equation.

###  III.3 Mixing Ratio.

We can use the Boltzmann distribution to address the question how the mixing
ratio $D$ affects assembly under equilibrium conditions. This is done by
computing contours in the $c_{0}-D$ plane along which the packaging fraction
$P_{12}$ is fixed. The dots in Fig. 13 were obtained were computed numerically
for two different values of $P_{12}$ for the same tree molecule and energy
parameters as Fig. 9:

Figure 13: Phase coexistence of assembled particles with excess free monomers
and excess free tree molecules for $\epsilon_{1}=-0.2$, $\epsilon_{3}=-1$, and
$\mu_{0}=-2$. Horizontal axis: Depletion factor $D$. Vertical axis: Pentamer
concentration $c_{0}$. Blue dots: points where 95 percent of the genome
molecules are packaged. Most of the remaining pentamers are free. Ochre dots:
points where 35 percent of the tree molecules are packaged. The green sector
is below the critical aggregation concentration at $c_{0}\simeq 0.2$ so there
are no capsids in the green sectors. Solid blue and ochre lines: contours of
fixed packaging fraction as computed from Eq.9 under the neglect of assembly
intermediates. Deviations between the solid lines and the dots are a measure
of the importance of assembly intermediates. The red star is a region of the
coexistence diagram where there is strong kinetic selection of genome
molecules.

Contours of fixed genome packing fraction can be computed analytically if one
again neglects assembly intermediates. This approximation produces the family
of hyperbolae

$c_{0}(D)\simeq\frac{1}{(1-D\thinspace
P_{12})}\left(\frac{KP_{12}}{1-P_{12}}\right)^{1/12}$ (9)

Note that $c_{0}(D)$ diverges at $D=1/P_{12}$. The agreement between the
approximate equilibrium theory and the actual values is reasonable for the
$P_{12}=0.95$ contour (solid blue line vs. blue dots). The blue curve in Fig.
13 can be viewed as an equilibrium phase boundary that separates two forms of
phase coexistence. To the left, nearly all tree molecules are encapsidated
with assembled particles in coexistence with excess free monomers while to the
right most pentamers are part of assembled particles in coexistence with
excess free tree molecules. This is consistent with the intuitive chemical
reaction picture discussed in the introduction section. It seems surprising
that the boundary line is shifted to values of $D$ below the stoichiometric
ratio $D=1$ when $c_{0}$ is reduced. The reason for this shift becomes evident
if one recalls there can be no assembly for $c_{0}$ less than the CAC for
empty capsid assembly. According to Fig. 13, this CAC is around $0.2$ (green
sector). The threshold value of $D$ below which the state of assembled
particles dominates thus necessarily has to go to zero as $c_{0}$ approaches
the CAC.

The agreement of the $P_{12}=0.35$ contour (ochre line) with the actual values
(ochre dots) is poor. This can be traced to the neglect of assembly
intermediates. Note also that the $P_{12}=0.35$ contour diverges at larger
values of $D$ than the $P_{12}=0.95$ contour. Figure 14 (Top) shows that for
larger values of $D$ the relative contribution of assembly intermediates
indeed becomes comparable to the concentration of assembled capsids, which
invalidates the assumption used to obtain Eq.9.

Figure 14: Equilibrium occupation numbers as a function of the depletion
factor $D$ for $c_{0}=1.0$. Top: same parameters as Fig. 10. The
stoichiometric ratio D=1 is indicated. Bottom: same parameters as Fig.10
except that the ratio of the attractive interactions between the RNA-capsid
and the capsid-capsid interactions has increased from $0.2$ to $1.2$. The
reference chemical potential $\mu_{0}$ was reduced to $-5$ so the total
assembly energy remains approximately the same.

When the ratio between the genome/pentamer and pentamer/pentamer interaction
strengths is increased from $-\epsilon_{1}=0.2$ to $-\epsilon_{1}=1.2$ then
the contribution from the partially assembled particles overwhelms that of
fully assembled capsids for $D$ larger than one, as shown in Fig. 14 (bottom).
This is due to the fact that when $D$ increases more and more pentamer binding
sites become available as there are more tree molecules. Breaking up assembled
particles and distributing the pentamers over the additional binding sites
increases the entropy of the system while larger values of the binding energy
of individual pentamers to specific edges of the tree molecules make up for
the loss of pentamer-pentamer adhesion. As a result a high-entropy state with
a distribution of partial shells has a lower free energy than the state of
coexistence of viral particles with excess genome molecules. This transition
can be viewed as a form of an order-disorder transition. While such a
transition would seem to conflict with the earlier chemical reaction picture,
in actuality the LMA itself remains valid for Fig. 14(Bottom) provided all
assembly intermediates.

##  IV Kinetics.

In this section we complete the definition of the model by specifying the
kinetics. We first consider the case where there is only one kind of spanning
tree molecule in solution.

###  IV.1 Master Equation.

The assembly kinetics of the model is defined in terms of thirteen occupation
probabilities $P_{n}(t)$ that now depend on time. The assembly and disassembly
dynamics for a given tree molecule with a particular assembly pathway is
assumed to obey Markov chain statistics [52] for which the occupation
probabilities evolve in time according to the master equation [53]:

$\frac{dP_{n}(t)}{dt}=\sum_{m=n\pm
1}\left(W_{m,n}P_{m}(t)-W_{n,m}P_{n}(t)\right)$ (10)

Here, $W_{m,n}$ is a thirteen-by-thirteen matrix of transition rates from
state $m$ to state $n$. We only include transitions with $m=n\pm 1$ so with
only one pentamer gained or lost at a time. Note that the diagonal entries of
$W_{m,n}$ are not specified at this point. The matrix of transition rates
$W_{m,n}$ will be defined in terms of simple diffusion-limited chemical
kinetics (see Eq. 8.35 of ref.[54]) where the addition of a pentamer to a
genome molecule associated with a pentamer cluster of size $n$ is treated as a
bimolecular reaction with a rate $k_{n,n+1}{r_{n}}{c_{f}}$ where $k_{n,n+1}$
is the on-rate defined as

$k_{n,n+1}=\lambda\begin{cases}e^{-\beta(E(n+1)-E(n))}\quad\quad\thinspace\thinspace\text{if}\thinspace\quad
E(n+1)>E(n)\\\ 1\qquad\qquad\qquad\quad\qquad\thinspace\text{if}\quad
E(n+1)<E(n)\end{cases}$ (11)

Here, $\lambda$ is a base rate that depends on quantities like the diffusion
coefficient but that is independent of concentration. If adding a pentamer
reduces the energy then the rate is equal to the base rate. If there is an
energy cost to adding a pentamers, then the base rate is reduced by an
Arrhenius factor, similar to the Metropolis algorithm of Monte-Carlo
simulations.

The on-rates are related to the rate matrix by $W_{n,n+1}=k_{n,n+1}{c_{f}}$.
The entries of the rate matrix that correspond to adding a pentamer are then

$W_{n,n+1}=\lambda c_{0}\Gamma\thinspace([P_{i}])\begin{cases}e^{-\beta\Delta
E_{n,n+1}}\qquad\text{ if }\Delta E_{n,n+1}>0\\\
1\qquad\qquad\qquad\thickspace\thinspace\text{ if }\Delta
E_{n,n+1}<0\end{cases}$ (12)

where $\Delta E_{n,n+1}\equiv(E(n+1)-E(n))$. The off-rate entries $W_{n+1,n}$
are determined by the on-rates through the condition of detailed balance:

$\frac{W_{n+1,n}}{W_{n,n+1}}=\frac{P_{n}}{P_{n+1}}$ (13)

where on the right hand side the Boltzmann Distribution Eq. (1) must be
inserted. This results in

$\begin{split}&W_{n,n+1}=c_{0}\Gamma\thinspace([P_{i}])\begin{cases}e^{-\beta\Delta
E_{n,n+1}}\qquad\Delta E_{n,n+1}>0\\\ 1\qquad\qquad\qquad\thickspace\Delta
E_{n,n+1}<0\end{cases}\\\
&W_{n+1,n}=\qquad\qquad\thickspace\begin{cases}1\qquad\qquad\qquad\thickspace\Delta
E_{n,n+1}>0\\\ e^{\beta\Delta E_{n,n+1}}\qquad\thickspace\thickspace\Delta
E_{n,n+1}<0\end{cases}\end{split}$ (14)

where the base rate $\lambda$ has been absorbed in a redefinition of time.
Note that the master equation is nonlinear because of the dependence of the
factor $\Gamma(|P_{i}|)$, as defined in Eq. (3), on the occupation
probabilities.

Some examples of numerical solutions of the master equation are shown in Fig.
15 for molecule (1) under reference conditions $c_{0}=D=\beta=1$.

Figure 15: Numerical solution of the master equation for molecule (1). Top:
parameter values are those of Fig. 9 with $\epsilon_{1}=-0.2$. Bottom: same
except that $\epsilon_{1}=-1.0$. Shown are the occupation probabilities
$P_{n}(t)$. The color code is the same as that of Fig.14. Time is in
dimensionless units of $1/\lambda$.

The top figure shows the case of $\epsilon_{1}=-0.2$. Other energy parameters
are those of Fig. 12), which shows the energy profile that corresponds to the
figure. The red curve is the probability $P_{12}(t)$ for a tree to be
encapsidated, the black curve the probability $P_{0}(t)$ for a tree to be free
of pentamers. Eventually, about forty percent of the genome molecules are
encapsidated by complete capsids with $n=12$ pentamers while about five
percent are encapsidated by partial capsids with $n$ less than 12. We checked
the assembly pathway and it is consistent with the Zlotnick Model shown in
Fig. 1. The occupation probabilities exponentially approach constant values at
late times that agree with the Boltzmann distribution. Closely similar results
are obtained if one replaces molecule (1) by molecule (2). At the this point,
the packaging kinetics does not seem to be able to really distinguish between
the two molecules.

The bottom figure shows the effect of increasing the relative strength of the
interaction between pentamers and genome molecules to $\epsilon_{1}=-1.0$,
while the reference chemical potential is reduced to $\mu_{0}=-5.2$ for the
reason discussed above Fig. 12. The final equilibrium state is now a
polydisperse mixture of aggregates of various sizes, consistent with Fig:14.
While for $\epsilon_{1}=-1.0$, the assembly pathway still follows that of the
Zlotnick Model for both molecule (1) and (2), for $\epsilon_{1}=-1.2$ the
assembly pathway of molecule (2) is completely changed. Now flaccid assembly
intermediates appear, such as the shown in Fig. 11.

###  IV.2 Time Scales

The occupation probabilities in Fig. 15 are governed by multiple scales. The
first important time scale is defined in Fig. 16, which shows the early stages
of the assembly of molecule (1) for $\epsilon=-0.2$ but now in magnified form
as compared with Fig.15(top).

Figure 16: Assembly shock-wave for molecule (1) for $\epsilon=-0.2$. The
intersection of the maximum tangent of the $P_{12}(t)$ curve with the time
axis (solid black line) defines the assembly lag time $t_{d}$.

The delay time $t_{d}$ is defined as the intercept of the maximum tangent of
$P_{12}(t)$ with the time axis. This gives $t_{d}\simeq 11$ for the case (of
molecule (1) with $\epsilon_{1}=-0.2$. As mentioned in Section I, assembly
delay times are a familiar feature of the assembly kinetics of empty capsids
[15] and of aggregation phenomena in general [55]. A second important time
scale refers to the late-time relaxation of the occupation probabilities
towards the Boltzmann distribution. Suppose one completes the definition of
the transition matrix by introducing the diagonal entries
$W_{n,n}=-\sum_{m\neq n}W(m,n)$. The resulting matrix $W_{m,n}$ now has column
elements adding to zero. Using this completed transition matrix, the master
equation can be rewritten in the form of a matrix equation
$\frac{d\bf{P}}{dt}=\bf{WP}$. It immediately follows that eigenvalues of the
completed transition matrix are the decay rates of the various modes that
correspond to the eigenvectors. The relaxation rate determining the approach
to final equilibrium is the smallest eigenvalue of $W_{m,n}$. The relaxation
time $t_{r}$ is defined as the inverse of the smallest eigenvalue so
$P_{n}(t)-P_{0}(\infty)\propto\exp-t/t_{r}$ in the late time limit for any
$n$. Figure 15 shows that this time-scale is much longer than the delay time
$t_{d}$, specifically $t_{r}\simeq 644$ so about two orders of magnitude
longer than $t_{d}$. The different functions $P_{n}(t)$ in Fig. 16 display a
maximum as a function of $n$. The corresponding peak times increase with $n$,
thus describing a type of assembly shock-wave propagating in configuration
space from small to large $n$. Similar assembly shock waves have been reported
for dynamical versions of the Zlotnick Model [38, 13]. The dependence of the
two time scales on the energy scale is also quite different. For $\beta=3$,
$t_{d}\simeq 11.7$ while $t_{r}\simeq 1.7\times 10^{5}$ producing a
significant separation in time-scales. Relaxation to equilibrium is evidently
an activated process, which is not surprising given the form of the assembly
energy profiles. In addition, increasing $\beta$ from one to three also has
the effect of suppressing the partially assembled capsids in Fig. 15 (top).
This remains the case if the mixing ratio $D$ is increased from one to two,
which now produces simple phase coexistence of fully assembled particles and
bare genome molecules.

Figure 17 shows the effect on the kinetics of varying the total pentamer
concentration $c_{0}$.

Figure 17: Plots of the fraction $P_{12}(t)$ of genome molecules that are
fully encapsidated for increasing total concentration for different total
pentamer concentration $c_{0}$ ranging from 1.0 to 0.2. Particles do not form
for $c_{0}$ less than 0.2.

The plot qualitatively reproduces the time-dependent light scattering studies
of the assembly of empty capsids [15]. Note that there is no assembly if the
pentamer concentration drops below about $0.2$, which is close to the CAC
shown in Figs 12 and 13. Fixing the concentration and varying the depletion
factor produces similar-looking plots.

###  IV.3 Disassembly Kinetics.

An important question concerns the fate of assembled particles when the
solution concentration $c_{0}$ of pentamers is reduced. This can be probed by
using the outcome of an assembly run with the total pentamer concentration
$c_{0}=1$ as the initial condition for a second run with a ten-fold reduced
pentamer concentration, so with $c_{0}=1$ reduced to $c_{0}=0.1$. The results
are shown in Fig. 18, still for molecule (1).

Figure 18: Top: The assembly run of Fig. 14 at $c_{0}=1$ and $\beta=1$
followed by a disassembly run at $c_{0}=0.1$. The color coding is that of
Fig.12. The assembled particles evaporate. Bottom: Same as the top figure
except that the temperature is reduced to $\beta=3$.

For $\beta=1$, the assembled particles rapidly disintegrate when the pentamer
concentration is reduced. The time-scale for the disintegration is of the
order of the relaxation time $t_{r}\simeq 31.6$, which is reduced compared to
the case of $c_{0}=1$. The energy activation barrier is, in this case, unable
to “protect” the assembled particle from disintegration by thermal activation.
Recall we found that when $\beta$ was increased from $1$ to $3$, then the
delay time $t_{d}$ barely increased while the relaxation time $t_{r}$
increased by three orders of magnitude. Repeating the assembly-disassembly run
for $\beta=3$ (bottom figure), produced the result that the concentration of
assembled particles still decreased with time during the disassembly run but
now on time-scales of the order of $10^{6}$. In contrast, during the assembly
run about fifty percent of the genome molecule had been encapsidated by a time
$t\simeq 10^{4}$. This assembly time scale is significantly longer than the
delay time but still is two orders of magnitude shorter than the
characteristic time scale for disassembly. We attribute this to the fact the
activation energy barrier for assembly is significantly lower than for
disassembly for a supersaturated system. Such a separation in assembly and
disassembly time scales would seem to be an essential condition for a
functioning viral particle. Below we will focus on the case of $\beta=3$.

Next, we carried out assembly-disassembly runs for $\epsilon_{1}=-1.2$ and
$\mu_{0}=-5.6$. The assembled particles and the various intermediate
aggregates rapidly “evaporated” during the disassembly run even for $\beta=3$.
The reason for the loss of stability can be seen by comparing the minimum
energy assembly profiles for $\epsilon_{1}=-0.2$ and $\epsilon_{1}=-1.2$ if
the reference chemical potential is set equal to $\mu^{*}$.

Figure 19: Assembly energy profiles for $\epsilon_{1}=-0.2$ (top) and
$\epsilon_{1}=-1.2$ (bottom) under conditions of assembly equilibrium with
$\mu_{0}=\mu^{*}$.

For $\epsilon_{1}=-0.2$, the energy profile is quite similar to that of the
Zlotnick Model (see Fig. 2). A standard nucleation and growth scenario is
expected for the assembly process. On the other hand, $E(n)$ has a deep
minimum at $n=6$ for $\epsilon_{1}=-1.2$ under assembly equilibrium
conditions. This corresponds to half-assembled particles, which is the exact
opposite of what is predicted by the nucleation and growth scenario. Recall
that in that case $n=6$ is expected to be close to an energy maximum.
Increasing the magnitude of $\epsilon_{1}$ indeed progressively destroys the
activation energy barrier that “protects” the assembled state.

##  V Packaging Competition

We are now in a position to carry out packaging competitions between two
genome molecules that either have a different tree topology or that have the
same tree topology but different dodecahedral wrapping configurations. We did
this for a system containing a concentration $c_{0}$ of pentamers as well as
equal concentrations of molecules (1) and (2) that competed with each other
for packaging for the case $\epsilon_{1}=-0.2$. Let $P_{n}^{(1)}(t)$ and
$P_{n}^{(2)}(t)$ be the corresponding sets of occupation probabilities.
Conservation of both types of tree molecules requires that the occupation
probabilities $P_{n}^{(1)}(t)$ and $P_{n}^{(2)}(t)$ separately sum to one.
Conservation of pentamers is satisfied if
$(c_{f}/c_{0})=\Gamma\thinspace([P_{n}])$ obeys the condition

$\Gamma\thinspace([P_{n}])=1-(D/24)\sum_{n=1}^{12}n(P_{n}^{(1)}+P_{n}^{(2)})$
(15)

The two occupation probabilities $P_{n}^{(1)}(t)$ and $P_{n}^{(2)}(t)$ obey
separate master equations of the form of Eq. (10). The only difference is that
now the factor $\Gamma\thinspace([P_{n}])$ in the expression for the
transition matrix Eq.15 couples the two master equations. Figure 20 shows an
example of the solution of the two coupled master equations for two different
parameter sets.

Figure 20: Time-dependence of the fraction of packaged tree molecules of
example molecules (1) and (2). (A): $\beta=3$, $c_{0}=D=1$ and $\mu_{0}=-2.5$.
The relaxation times of the two molecules are indicated by the two arrows.
(B): Same as (A) except that $c_{0}=0.22$ and $D=0.3$.

The relaxation times are indicated by arrows. Based on the discussion in
Section 2, it is expected that the highly branched molecule (1) will
outcompete linear molecule (2) since it can accommodate more optimally wrapped
pentamers during assembly. Recall however that in actuality the assembly
energy profiles of the two molecules are not very different (see Fig.9). The
maximum selectivity in Fig. 20(A) of molecule (1) over molecule (2) for
parameter values $\beta=3$, $c_{0}=D=1$ and $\mu_{0}=-2.5$. It is achieved
very early in the assembly process. In order to be quantitative, define the
relative packaging selectivity $S(1|2)$ of molecule (1) with respect to
molecule (2) to be

$S(1|2)\equiv(P_{12}^{(1)}-P_{12}^{(2)})/(P_{12}^{(1)}+P_{12}^{(2)})$ (16)

which ranges from one to minus one. This fraction depends on the time instant
when probabilities are being computed. The maximum value of the relative
packaging selectivity will be denoted by $S_{m}(1|2)$. For 20(A), the maximum
selectivity $S_{m}(1|2)\simeq 0.43$. The selectivity $S(1|2)(t)$ then slowly
decays to zero on a time scale of the order of the relaxation time
$t_{r}^{(2)}$. The fact that the selectivity goes to zero in the long time
limit can be understood from the fact that the total assembly energies of
molecules (1) and (2) are the same, along with the fact that there are
practically no intermediate assemblies. Selectivity is a purely kinetic
effect.

To further increase the selectivity, we varied the total pentamer
concentration $c_{0}$. For given $c_{0}$, the mixing ratio $D$ was given the
maximum value for which at least 90 percent of the type (1) genome molecules
were still encapsidated. The maximum selectivity was obtained when $c_{0}$ was
reduced to a value close to the lower bound of $0.2$ (the CAC). Figure 20(B)
shows the case $c_{0}\simeq 0.22$ with $D=0.3$. The maximum selectivity is
$S_{m}(1|2)\simeq 0.8$. This is the highest selectivity that we were able to
achieve. This increase in selectivity near the CAC is not the result of
increasing relaxation times (see Fig. 20(B)). If $c_{0}$ drops below the CAC
then the selectivity of the particles that assemble remains about $0.8$ but
only a very small fraction of the genome molecules are encapsidated.

The fact that the kinetic selectivity is largest in the limit of small $c_{0}$
is intriguing. In an Arrhenius description the nucleation rate for particle
assembly is the product of an attempt frequency and a Boltzmann factor
$\exp-\beta\Delta E$ with $\Delta E$ the assembly activation barrier. The
assembly energy profiles do not depend on the pentamer concentration $c_{0}$
so this increase of the selectivity when $c_{0}$ is reduced must be due to a
dependence of the attempt frequency on $c_{0}$ in a manner that differentiates
between the two molecules. According to Fig. 9, the first two assembly steps
of molecules (1) and (2) that bring the process to the energy maximum at $n=2$
are the same. Once past $n=2$, the negative slope of the energy profile for
molecule (1) is significantly steeper than that of molecule (2). Only after
three additional assembly steps does the energy profile of molecule (2)
acquire a similarly large negative slope. This distinction is a consequence of
the fact that the wrapping number of molecule (1) is six while it is two for
molecule (2). The energy of capsid assemblies with $n$ greater than two thus
decreases more rapidly for molecule (1). Moreover, the two maximally wrapped
pentamers of molecule (2) are not adjacent. If the pentamer concentration is
low then the packaging of molecule (2) will be significantly delayed as
compared to molecule (1). The system performs a weakly directed random walk
along the flattish top of the energy barrier with a significant probability
for disassembly of the assembly nucleus. This happens if the random walk
returns to n=2. As a result, the effective attempt frequency of molecule (2)
is significantly reduced. The reason why in the limit of small $c_{0}$ the
mixing ratio $D$ must be reduced to well below stoichiometric ratio $D=1$ in
order to achieve a 90 percent packaging probability is a consequence of the
equilibrium thermodynamics discussed earlier. Note that the point $c_{0}(D=0)$
on the contour in the equilibrium phase-coexistence diagram Fig. 13 for fixed
90 percent packaging probability is also the point with the lowest value for
$c_{0}$.

We constructed a table of packaging selectivities of pairs of trees with
different MLD and wrapping numbers (Fig. 21). This was done for $\beta=3$ and
$c_{0}=0.3$ instead of $0.2$, which produced an increase of the mixing ratio
for 90 percent packaging from $D=0.3$ to a more reasonable $D=0.5$.

Figure 21: Plot of wrapping numbers ($N_{P}$) and maximum ladder distances
(MLD showing the outcomes of packaging contests indicated by arrows. Each
arrow indicates the outcome of a contest, with the arrow pointing from the
loser to the the winner. The width of the arrow is a measure of the maximum
selectivity $S_{m}$. A dashed arrow indicates a weak selectivity. The energy
parameters were $\beta=3$, $c_{0}=0.3$, $D=0.5$, $\epsilon_{1}=-0.2$,
$\epsilon_{3}=-1$, and $\mu_{0}=-2.5$

The result is clear: when two trees compete with the same wrapping number but
different MLD, then the tree with the smaller MLD outcompetes the tree with
the larger MLD. When two trees with the same MLD but different wrapping number
compete, then the tree with the larger wrapping number outcompetes the tree
with the smaller wrapping number.

It would seem that one should be able to achieve an even higher selectivity by
increasing $|\epsilon_{1}|$, the ratio of strength of the genome-pentamer
attraction and the pentamer-pentamer attraction. The outcome of a packaging
contest between molecules (1) and (2) is shown in Fig. 22 for the case
$|\epsilon_{1}|=1$. The reference chemical potential was adjusted so the
energy activation barrier was the same as for our standard value
$|\epsilon_{1}|=0.2$.

Figure 22: Time-dependence of the packaging selectivity $S(1|2)(t)$ of
molecule (1) with respect to molecule (2). The energy parameters are
$\epsilon_{1}=-1.0$, $\epsilon=-1$, and $\mu_{0}=-5.2$

Figure 22 shows the time-dependence of the selectivity $S(1|2)$. During the
first stage of the assembly—with $t$ less than $10^{4}$—$S(1|2)$ rises very
quickly to (nearly) one, which agrees with the expectation that increasing
$|\epsilon_{1}|$ should enhance the packaging selectivity. However, as time
progresses, fully packaged type (1) molecules particles start to disassemble
while additional molecule (2) assemblies start to appear. The packaging
selectivity drops, goes to zero and then changes sign. This “inversion” is in
fact a generic feature of packaging competitions for large values of
$|\epsilon_{1}|$. Eventually, about 70 percent of the type (2) molecules are
packaged in complete capsids. A large fraction of type (1) molecule are
associated with a polydisperse distribution of incomplete particles. Recall
from Section IV that increasing $|\epsilon_{1}|$ produces a large fraction of
incomplete particles. Recall also that for large $|\epsilon_{1}|$ assembled
particles easily disintegrate when the pentamer concentration is reduced.

Since increasing $|\epsilon_{1}|$ leads to thermodynamically unstable
particles, how small can one make $|\epsilon_{1}|$ and still have a reasonable
selectivity? Figure 23 shows how the maximum packaging selectivity
$S(1|2)_{m}$ depends on $|\epsilon_{1}|$ in the regime of small
$|\epsilon_{1}|$.

Figure 23: Maximum packaging selectivity $S(1|2)_{m}$ of molecule (1) with
respect to molecule (2) as a function of $|\epsilon_{1}|$. The reference
chemical potential $\mu_{0}$ was adjusted to maintain the activation energy
barrier at an approximately constant value.

The maximum packaging selectivity rises very steeply from zero for increasing
$|\epsilon_{1}|$ and reaches it maximum value around $|\epsilon_{1}|\simeq
0.4$. Polydispersity and selectivity inversion only appears for larger values
of $|\epsilon_{1}|$.

##  VI Conclusion

We have presented a simple, quasi-analytical model for the study of the
packaging of tree-like genome molecules inside dodecahedral capsids composed
of twelve pentamers. The model is sufficiently simple that the kinetics can be
obtained from the numerical solution of a set of coupled master equations. The
genome molecules are represented by a “core” of specific links with enhanced
affinity for the edges and that visits all vertices of the dodecahedron but
covers only nineteen of the edges. The spanning tree is complemented by eleven
non-specific links with reduced edge affinity. The minimum energy assembly
pathways are characterized by two numbers: the maximum ladder distance (MLD)
and the wrapping number ($N_{P}$.) The MLD of a spanning tree is a topological
invariant of the tree that is a global measure of the “branchiness” of a tree
while $N_{P}$ is a geometrical measure of the surface distribution of specific
links over the dodecahedron in terms of the number of pentamers associated
with a maximum of four specific links that can be accommodated.

The assembly kinetics is characterized by a delay time $t_{d}$ for the onset
of particle production and a thermodynamic relaxation time $t_{r}$. When the
solution concentration $c_{0}$ of pentamers is lowered, assembled particles
are stable against disassembly by thermal fluctuations on time scales shorter
than $t_{r}$. If the energy scale of the binding energy is large compared to
the thermal energy and if the ratio $|\epsilon_{1}|$ between specific and non-
specific affinity of links of the genome molecules for the pentamer edges is
small compared to one then there is a pronounced separation in time scale
between particle assembly at high $c_{0}$ and particle disassembly at low
$c_{0}$. We carried out packaging selection contests between molecules with
different MLD and $N_{P}$. The selectivity is a purely kinetic effect, based
on minor but systematic differences between the minimum energy assembly
pathways of the different genome molecules. For small values of
$|\epsilon_{1}|$ molecules with small MLD and large $N_{P}$ outcompete
molecules with large MLD and small $N_{P}$. For increasing $|\epsilon_{1}|$,
genome selectivity increases as well but if $|\epsilon_{1}|$ becomes
comparable to one then assembly produces a polydisperse solution of partially
assembled particles that readily disintegrate when the pentamer solution is
reduced. In addition, for large $|\epsilon_{1}|$ the assembly energy profiles
become increasingly complex, leading to the failure of the nucleation-and-
growth scenario. Another interesting consequence of increasing
$|\epsilon_{1}|$ is a change from mechanically rigid to mechanically flaccid
assembly intermediates.

In the Introduction we mentioned a number of experimental observations that
motivated the construction of the model. The first was an in-vitro study of
the co-assembly of CCMV with non-CCMV RNA molecules [30]. When the RNA-to-
protein mixing ratio was low, virus-like particles (VLPs) formed with excess
proteins, in agreement with simple arguments based on phase-coexistence. When
the mixing ratio exceeeded a threshold, the virus-like particles were replaced
by disordered RNA-protein aggregates instead of the expected coexistence of
VLPs with excess RNA. Moreover, the threshold point separating the two regimes
was well below the stoichiometric ratio. Can this be understood in the light
of our results? When the mixing ratio (i.e., the depletion factor $D$) is
increased at fixed total protein concentration $c_{0}$ in Fig. 13 then for
smaller$D$ practically all tree molecules are encapsidated (i.e., to the left
of the blue line in Fig. 13) while larger values of $D$ an increasing fraction
of the tree molecules are not encapsidated. The transition point separation
the two regimes drops below the stoichiometric ratio when $c_{0}$ approaches
the CAC. The authors of the experimental study [30] relate this displacement
away from the stoichiometric ratio to electrostatic effects, which are not
included in our model. The model shows that there is a separate mechanism that
could produce the same effect. An experimental study in which the capsid
protein concentration is reduced towards the CAC, perhaps accompanied by
changes in the salinity, may be able to distinguish between the two
mechanisms. Next, Fig. 14 shows that with increasing $D$, partially assembled
capsids start to appear as observed experimentally. This would seem to be
consistent with the proposed model but this polydispersity only becomes a
dominant effect if $|\epsilon_{1}|$ approaches one (see Fig. 14 bottom). In
that case assembled particles should be thermodynamically unstable, as
discussed in section IV-C. In other words, they could not be true VLPs. This
certainly would be an interesting result if it were confirmed experimentally.
A second interesting observation was the fact that for increasing values of
the ratio between the RNA-protein and the protein-protein affinities, virus-
like particles were replaced by disordered aggregates [31]. In our model, this
ratio is represented by $|\epsilon_{1}|$. For increasing $|\epsilon_{1}|$
fully assembled shells indeed are replaced by partially assembled shells.
Depending on the MLD and $N_{P}$ numbers, the assemblies can be mechanically
unstable. A third observation concerned the fact that asymmetric
reconstruction of the MS2 virus showed that a subsection of the RNA genome
reproducibly associated with a compact cluster of capsid proteins [36]. Our
model reproduces this observation for genome molecules with small MLD and
large $N_{P}$. The initial assembly is a compact cluster composed of maximally
wrapped pentamers (see Fig.3, top left).

While the model qualitatively reproduces these observations there is an
important quantitative difference. As compared with the experiments of ref.
[30] on CCMV, the model appears to underestimate the ability of increased RNA-
to-protein mixing ratios to suppress the assembly of complete capsids. While
this could be something specific for CCMV, we believe that this due to the
fact that the model underestimates the entropy of a partial assembly. This is
particularly obvious for flaccid intermediates that are expected to undergo
strong thermal conformational fluctuations that are not represented in the
model. While it may be possible to include conformational fluctuations in the
model, a simpler route may be to carry out Brownian Dynamics simulations such
as those of Ref. [39] but then for a system with pentamers that interact with
tree molecules. Another natural extension of the model would be to larger
viruses. The Zlotnick Model can be viewed as a representation of the smallest
capsids, known as $T=1$ capsids, in which all capsid proteins have the same
local packing organization. The CCMV and MS2 viruses that are an important
testing ground for physical theories of viral assemblies are T=3 viruses that
have a more complex architecture.

###### Acknowledgements.

We would like to thank Alexander Grosberg for drawing our attention to
spanning trees in the context of virus structure and assembly. We benefitted
from discussions with Justin Little, Chen Lin, Zach Gvildys and William Vong.
RB would like to thank the NSF-DMR for continued support under CMMT Grant
No.1836404.

## References

  * Crick and Watson [1956] F. H. C. Crick and J. D. Watson, Nature 177, 473 (1956).
  * Note [1] For double-stranded DNA viruses the genome molecules either are inserted into the interior pre-fabricated spherical capsids by a molecular motor or the genome is pre-condensed prior to assembly.
  * Ganser-Pornillos  _et al._ [2008] B. K. Ganser-Pornillos, M. Yeager, and W. I. Sundquist, Curr. Opin. Struct. Biol. 18, 203 (2008).
  * Fraenkel-Conrat and Williams [1955] H. Fraenkel-Conrat and R. C. Williams, Proc. Natl. Acad. Sci. U. S. A. 41, 690 (1955).
  * Butler and Klug [1978] P. J. G. Butler and A. Klug, Sci. Am. 239, 62 (1978).
  * Klug [1999] A. Klug, Philos. Trans. R. Soc. Lond. B. Biol. Sci. 354, 531 (1999).
  * Note [2] For a quantitative treatment of this model, see ref. [kegel2006].
  * Mathews  _et al._ [1999] D. H. Mathews, J. Sabina, M. Zuker, and D. H. Turner, Journal of molecular biology 288, 911 (1999).
  * Tubiana  _et al._ [2015] L. Tubiana, A. L. Božič, C. Micheletti, and R. Podgornik, Biophysical journal 108, 194 (2015).
  * Yoffe  _et al._ [2008] A. M. Yoffe, P. Prinsen, A. Gopal, C. M. Knobler, W. M. Gelbart, and A. Ben-Shaul, Proceedings of the National Academy of Sciences 105, 16153 (2008).
  * Safran [1994] S. Safran,  _Statistical Thermodynamics of Surfaces, Interfaces, and Membranes_ (Addison-Wesley Pub., 1994).
  * Ceres and Zlotnick [2002] P. Ceres and A. Zlotnick, Biochemistry 41, 11525 (2002).
  * Morozov  _et al._ [2009] A. Y. Morozov, R. F. Bruinsma, and J. Rudnick, J. Chem. Phys. 131, 155101 (2009).
  * Prevelige  _et al._ [1993] P. E. Prevelige, D. Thomas, and J. King, Biophys. J. 64, 824 (1993).
  * Casini  _et al._ [2004] G. L. Casini, D. Graham, D. Heine, R. L. Garcea, and D. T. Wu, Virology 325, 320 (2004).
  * Medrano  _et al._ [2016] M. Medrano, M. Á. Fuertes, A. Valbuena, P. J. Carrillo, A. Rodríguez-Huete, and M. G. Mateu, Journal of the American Chemical Society 138, 15385 (2016).
  * Zandi  _et al._ [2006] R. Zandi, P. van der Schoot, D. Reguera, W. Kegel, and H. Reiss, Biophys. J. 90, 1939 (2006).
  * Asor  _et al._ [2019] R. Asor, L. Selzer, C. J. Schlicksup, Z. Zhao, A. Zlotnick, and U. Raviv, ACS nano 13, 7610 (2019).
  * Bryngelson  _et al._ [1995] J. D. Bryngelson, J. N. Onuchic, N. D. Socci, and P. G. Wolynes, Proteins: Structure, Function, and Bioinformatics 21, 167 (1995).
  * Freddolino  _et al._ [2006] P. L. Freddolino, A. S. Arkhipov, S. B. Larson, A. McPherson, and K. Schulten, Structure 14, 437 (2006).
  * Zlotnick [1994] A. Zlotnick, J. Mol. Biol. 241, 59 (1994).
  * Bruinsma  _et al._ [2003] R. F. Bruinsma, W. M. Gelbart, D. Reguera, J. Rudnick, and R. Zandi, Phys. Rev. Lett. 90, 248101 (2003).
  * Lidmar  _et al._ [2003] J. Lidmar, L. Mirny, and D. R. Nelson, Phys. Rev. E 68 (2003).
  * Zandi  _et al._ [2004] R. Zandi, D. Reguera, R. F. Bruinsma, W. M. Gelbart, and J. Rudnick, Proc. Natl. Acad. Sci. U. S. A. 101, 15556 (2004).
  * Zlotnick [2007] A. Zlotnick, J. Mol. Biol. 366, 14 (2007).
  * Rapaport [2008] D. Rapaport, Phys. Rev. Lett. 101, 186101 (2008).
  * Arkhipov  _et al._ [2006] A. Arkhipov, P. L. Freddolino, and K. Schulten, Structure 14, 1767 (2006).
  * Mannige and Brooks [2009] R. V. Mannige and C. L. Brooks, Proc. Natl. Acad. Sci. U. S. A. 106, 8531 (2009).
  * Kaplan  _et al._ [2014] R. Kaplan, J. Klobusicky, S. Pandey, D. H. Gracias, and G. Menon, Artif. Life 20, 409 (2014).
  * Comas-Garcia  _et al._ [2012] M. Comas-Garcia, R. D. Cadena-Nava, A. L. N. Rao, C. M. Knobler, and W. M. Gelbart, Journal of Virology 86, 12271 (2012).
  * Garmann  _et al._ [2013] R. F. Garmann, M. Comas-Garcia, A. Gopal, C. M. Knobler, and W. M. Gelbart, J. Mol. Biol. 32, ??? (2013).
  * Baker  _et al._ [1999] T. S. Baker, N. H. Olson, and S. D. Fuller, Microbiol. Mol. Biol. Rev. 63, 862 (1999).
  * Tihova  _et al._ [2004] M. Tihova, K. A. Dryden, T. V. L. Le, S. C. Harvey, J. E. Johnson, M. Yeager, and A. Schneemann, J. Virol. 78, 2897 (2004).
  * Koning  _et al._ [2016] R. I. Koning, J. Gomez-Blanco, I. Akopjana, J. Vargas, A. Kazaks, K. Tars, J. M. Carazo, and A. J. Koster, Nature communications 7, 1 (2016).
  * Beren  _et al._ [2020] C. Beren, Y. Cui, A. Chakravarty, X. Yang, A. Rao, C. M. Knobler, Z. H. Zhou, and W. M. Gelbart, Proceedings of the National Academy of Sciences 117, 10673 (2020).
  * Dykeman  _et al._ [2011] E. C. Dykeman, N. E. Grayson, K. Toropova, N. A. Ranson, P. G. Stockley, and R. Twarock, J. Mol. Biol. 408, 399 (2011).
  * Patel  _et al._ [2015] N. Patel, E. C. Dykeman, R. H. A. Coutts, G. P. Lomonossoff, D. J. Rowlands, S. E. V. Phillips, N. Ranson, R. Twarock, R. Tuma, and P. G. Stockley, Proc. Natl. Acad. Sci. U.S.A. 10.1073/pnas.1420812112 (2015).
  * Endres and Zlotnick [2002] D. Endres and A. Zlotnick, Biophys. J. 83, 1217 (2002).
  * Perlmutter  _et al._ [2014] J. D. Perlmutter, M. R. Perkett, and M. F. Hagan, J. Mol. Biol. 10.1016/j.jmb.2014.07.004 (2014).
  * Rudnick and Bruinsma [2019] J. Rudnick and R. Bruinsma, Physical Review E 100, 012145 (2019).
  * Note [3] The assembly of the capsid is assumed here to take place on a specific location. The quantity $\mu_{0}$ reflects the entropic free energy cost of removing a pentamer from the solution to this location plus that of any conformational change that is required for the pentamer prior to joining a partial capsid.
  * Bollobás [2013] B. Bollobás,  _Modern graph theory_ , Vol. 184 (Springer Science & Business Media, 2013).
  * Graham and Hell [1985] R. L. Graham and P. Hell, Annals of the History of Computing 7, 43 (1985).
  * Loeb [1991] A. L. Loeb, in  _Space Structures_ (Springer, 1991) pp. 45–50.
  * Note [4] The set of unique spanning trees consists of all spanning trees that, when depicted as in the top left of Fig. 3, cannot be mapped into each other by rotations and/or reflections that leave the dodecahedron invariant.
  * Fang  _et al._ [2011] L. T. Fang, W. M. Gelbart, and A. Ben-Shaul, The Journal of Chemical Physics 135, 10B616 (2011).
  * Rudnick and Bruinsma [2005] J. Rudnick and R. Bruinsma, Phys. Rev. Lett. 94, 038101 (2005).
  * Dykeman  _et al._ [2013] E. C. Dykeman, P. G. Stockley, and R. Twarock, J. Mol. Biol. 425, 3235 (2013).
  * A.M.Gutin  _et al._ [1993] A.M.Gutin, A.Y.Grosberg, and E.I.Shakhnovich, Macromolecules 26, 1293 (1993).
  * Note [5] In actuality, condensation of the RNA genome molecules takes place during encapsidation. It is driven by positively charged polypeptide chains associated with the capsid proteins.
  * Safran [2018] S. Safran,  _Statistical thermodynamics of surfaces, interfaces, and membranes_ (CRC Press, 2018).
  * Perkett and Hagan [2014] M. R. Perkett and M. F. Hagan, J. Chem. Phys. 140, 214101 (2014).
  * Van Kampen [1992] N. G. Van Kampen,  _Stochastic processes in physics and chemistry_ , Vol. 1 (Elsevier, 1992).
  * Schulten [1999] K. Schulten,  _Non-equilibrium Statistical Mechanics_ (UIUC, 1999).
  * Wu [1992] D. T. Wu, J. Chem. Phys. 97, 2644 (1992).

##  Appendix A Demonstration hat the smallest MLD for spanning trees on the
dodecahedron is nine

We begin by noting that for every vertex on the dodecahedron there is a vertex
on the opposite side of the polyhedron that is a ladder distance five away.
That is, getting from one of the two vertices to the other requires traversing
at least five edges. Figure 24 shows such a path.

Figure 24: Two maximally separated vertices on the dodecahedron and one of the
12 shortest paths consisting of five edges that join them.

For each such pairs of vertices there are 12 minimal paths.

Now, assume that there is a spanning tree with MLD 8. In such a case, we can
pick out a path of ladder distance eight in that tree. All other elements of
the tree will consist of trees that branch out from that path. Figure 25 is a
figurative depiction of the path along with the longest allowed branch
sprouting off each vertex on that path. The likelihood of branching off those
“side branch” paths is ignored; such branching does not alter the argument
below.

Figure 25: A ladder distance 8 path in the hypothetical MLD 8 spanning tree on
the dodecahedron. The path is shown as a thick line, and the nine vertices are
labeled for easy reference. The thinner vertical lines represent longest
allowed paths branching off the ladder distance 8 path.

Consider first the central vertex on the ladder distance eight path, labeled 5
in Fig. 25. The side path with ladder distance four is the longest that can
attach to it. A longer path increases the MLD of the tree. Clearly, there is
no possibility of reaching a point a ladder distance five from vertex 5 along
any path with ladder distance four, so the path shown cannot connect the
central vertex to the vertex a distance five away from it. Next, consider the
two sites flanking the central vertices, labeled 4 and 6. Attached to each is
the longest possible path branching out from them, Such a path has ladder
distance three. If either of these paths reached to the vertex a ladder
distance five away from the central vertex, then there would be a ladder
distance four (or less) path from that vertex through one of the flanking
vertices to the maximally separated vertex, and we know that no such path
exists. We can continue this argument to encompass all allowed paths sprouting
from vertices on the chosen path. Thus, there is a vertex on the dodecahedron
that cannot be a part of the MLD 8 tree containing this path. Consequently no
tree with MLD 8 can be a spanning tree on the dodecahedron. The argument above
can clearly be applied to the possibility of a spanning tree with MLD less
than eight. That there is a spanning tree with MLD 9 is readily established by
construction.