# The Spanning Tree Model and 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

Single-stranded (ss) RNA viruses self-assemble spontaneously in solutions that
contain the viral RNA genome molecules and the viral capsid proteins. The
self-assembly of empty capsids can be understood on the basis of free energy
minimization of rather simple models. However, during the self-assembly of
complete viral particles in the cytoplasm of an infected cell, the viral
genome molecules must be selected from a large pool of very similar host
messenger RNA molecules. It is known that the assembly process takes the form
of preferential heterogeneous nucleation of capsid proteins on viral RNA
molecules (“selective nucleation”). Recently, a simple mathematical model was
proposed for the selective nucleation of small ssRNA viruses [1]. In this
paper we present a statistical physics analysis of the thermal equilibrium and
kinetic properties of that model and show that it can account, at least
qualitatively, for numerous observations of the self-assembly of small ssRNA
viruses.

##  I Introduction

In 1956 Crick and Watson [2] noted that the total length of the genome of a
virus should scale with the interior volume of the protein shell (or “capsid”)
enclosing the genome. This observation was a central part of their argument
that sphere-like viruses should have icosahedral symmetry. In a similar way,
the length of a rod-like capsid should scale with the length of the genome.
This approach works well for small, single-stranded (ss) RNA viruses (such as
the polio and common cold viruses) as well as for double-stranded DNA viruses
(such as the herpes and bacteriophage viruses) that have close-packed genomes
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. while for the large ssRNA
retroviruses, which do not have close-packed genomes, the genome length scales
with the interior surface area of capsid [4].

Many ssRNA viruses self-assemble spontaneously in solutions that contain
concentrations of viral capsid proteins and RNA molecules [5, 6]. This is a
“passive” process, driven by free energy minimization. Early work by Aaron
Klug [7] proposed that the RNA genome molecules act as templates that direct
the viral assembly process. He proposed a physical model in which the
repulsive electrostatic interactions between positively charged groups of the
capsid proteins are just strong enough to overcome competing attractive
hydrophobic interactions, which just prevents the self-assembly of empty
capsids under physiological conditions. When viral RNA molecules are added to
the solution, the negative charges of the RNA nucleotides neutralize some of
the positive charges of the capsid proteins and this tilts the free energy
balance towards assembly 222For a quantitative treatment of this model, see
ref. [kegel2006].

The genome molecules of ssRNA 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 [9].
The “co-assembly” picture suggests that there could be connections between the
geometrical and topological features of the secondary structure and the capsid
geometry. The redundancy of the genetic code allows for the possibility of
“silent” (or synonymous) mutations that alter the secondary structure of a
viral RNA molecule without altering the structure of the proteins encoded by
the nucleotide sequence [10]. Viral RNA molecules indeed have specific
specific sections, known as packaging signals, with enhanced affinity for the
capsid proteins [11]. These packaging signals favor the encapsidation of viral
RNA molecules over host mRNA molecules with comparable length. Separately, the
topology of viral RNA molecules differs from that of generic single-stranded
RNA molecules of the same length: they are significantly more branched and
compact than generic RNA molecules of the same length [12].

Such observations have been understood in terms of the packaging signals
reducing the equilibrium free energy of fully assembled viral particles but
recent work has focused on the kinetics of viral assembly[13]. Kinetic studies
of the self-assembly of empty capsids [14, 15, 16] provided support for a
nucleation-and-growth scenario where assembly is triggered by the formation of
a “nucleation complex”, composed of a small number of capsid proteins. This
nucleation complex then extends or elongates by absorbing more proteins until
the capsid closes up. When the theory of nucleation-and-growth applied to
viral capsids [17], the nucleation complex is identified with the critical
nucleus, i.e., the state of an assembly of capsid proteins for which the free
energy has a maximum as a function of the number of capsid proteins. Formation
of the critical nucleus is an energetically “uphill” process while the
subsequent elongation is a “downhill” process. Under conditions of
thermodynamic equilibrium, the critical nucleus would be a half-formed shell.
However, experimentally measured values for the interaction energy between
capsid proteins are in the range of a few $k_{B}T$ and this leads to an
activation energy barrier for the assembly of capsids under conditions of
thermodynamic equilibrium in the range of hundreds of $k_{B}T$. The
corresponding assembly rates would be prohibitively low on laboratory time
scales. The size of the nucleation complex measured by light scattering [15]
is much smaller than a half-formed shell while assembly rates of empty are
generally in the range of $0.01-1.0Hz$. These observations indicate that
capsid assembly takes place under conditions that are well away from
thermodynamic equilibrium, more specifically under a relatively high level of
supersaturation. Finally, a recent study of the assembly kinetics of empty
capsids of the Hepatitis B virus indicated that the formation of the critical
nucleus takes place along well-defined assembly pathways [18], similar to the
distinct assembly pathways that govern the folding of proteins [19].

Does the nucleation-and-growth scenario describe the co-assembly of viruses
complete with encapsidated RNA genome molecules? Early time-resolved small-
angle X-ray scattering experiments on viral assembly in bulk solutions were
interpreted in terms of heterogeneous nucleation of capsid proteins on viral
RNA molecules with the resulting nucleo-protein complex attracting additional
capsid subunits via a “diffusive antenna” mechanism [20, 21]. Recent
observations on the assembly kinetics of individual MS2 viruses (a small ssRNA
bacteriophage virus) [22] also provided support for a heterogeneous
nucleation-and-growth scenario. An important observation was that the assembly
process was characterized by an extremely broad distribution of time scales.
Important information about the kinetics of viral co-assembly has been gleaned
as well from static structural studies. Until recently, reconstruction of
packaged genome molecules involved “icosahedral averaging”, which resulted in
RNA structures with imposed icosahedral symmetry [23]. Such studies showed
that the interior surface of the icosahedral capsids of certain viruses (e.g.,
the nodaviruses [24, 25]) 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 how individual ssRNA genome molecules are packaged inside spherical
capsids without icosahedral averaging (“asymmetric reconstruction” [26, 27]).
An important example is again the MS2 virus. It was found that sections of the
RNA genome rich in packaging signals reproducibly associated with roughly half
of the interior surface of the capsid [28, 29], as shown in Fig.1

Figure 1: CryoEM asymmetric reconstruction of the MS2 virion. The viral genome
molecules (in green and blue) associate reproducibly mostly with one half of
the capsid shell, which has been removed for clarity. From ref.[29]

This can be interpreted as evidence for a well-defined, reproducible assembly
pathway of a critical nucleus that is a compact nucleo-protein complex held
together by particular sections of the viral RNA molecule that are rich in
packaging sequences. The downhill part involves mostly generic electrostatic
RNA-protein attractions and involves non-reproducible RNA configurations. In
terms of the theory of nucleation and growth, the packaging signals reduce the
height of the assembly energy activation barrier for the packaging of viral
RNA molecules. Importantly, this provides a natural mechanism for the
selection of viral RNA molecules from a large pool of host messenger RNA
(“mRNA”) molecules .

Understanding how the nucleation-and-growth scenario applies to viral assembly
is important also for biological reasons. When the molecular components of a
single-stranded RNA virus assemble in the cytoplasm of an infected cell and
form virions, a limited number of viral RNA molecules compete with a large
pool of similar host messenger RNA (“mRNA”) molecules for packaging by viral
capsid proteins [30]. Reduction of the assembly activation energy barrier by
the packaging signals has the effect of selectively increasing the assembly
rate of viral particles with respect to the packaging of host mRNA molecules,
a mechanism that might be called “selective nucleation”333Selective nucleation
was proposed by I. Rouzina in the context of the assembly of retroviruses. In
this view, the action of packaging signals is similar to that of enzymes or
catalysts that increase the rate of a chemical reaction.

Nucleation-and-growth is not the only possible assembly scenario. The
packaging of linear genome molecules in small capsids has been numerically
simulated [32] and two different packaging scenarios were encountered. This
depended on the affinity ratio of the protein-genome and protein-protein
interactions. For low values of this affinity ratio, a nucleation-and-growth
assembly scenario was encountered of the diffusive antenna type. For larger
values of the affinity ratio, a so-called “en masse” scenario was encountered.
The first assembly step is, in that case, the formation of a disordered
protein/genome condensate that shrinks in size as more proteins are being
added. The condensate then undergoes an order-disorder transition in which the
spherical symmetry of the condensate is broken and an icosahedral particle
forms [33]. This form of assembly may have been reproduced in a two-stage in-
vitro assembly experiments of CCMV virus-like particles [34]. Micro-mechanical
experiments on the encapsidation of linear double-stranded genome molecules by
capsid proteins [35] have also been interpreted according to the en-massed
scenario.

In this article we examine the physical properties of a recently proposed
mathematical model – the “Spanning Tree Model” – for the selective nucleation
of small RNA viruses. The model is sufficiently simple so that its kinetics
can be determined by numerical integration of a Master Equation. The model is
used to determine thermodynamic parameters that optimize RNA selectivity, to
follow the effect of changing the affinity ratio on the assembly scenario and
to explore which topological and geometrical properties of RNA molecules
optimize selectivity. The model is a generalization of an early model for the
assembly of empty dodecahedral capsids in a solution of pentamers that was due
to Zlotnick [36, 37, 38]. The kinetics of the Zlotnick mode obeys a
nucleation-and-growth assembly scenario [39]. It was used to carry out
simulations of the packaging of linear genome molecules [32].

The spanning tree model is introduced in Section II together with the
topological and geometrical classification of the model genome molecules
followed by a discussion of structural transitions of partial assemblies as a
function of the affinity ratio. The thermodynamic properties of the model are
discussed in Section III. The Master Equation for the assembly kinetics is
introduced in Section IV. Numerical integration of the Master Equation is used
to determine the characteristic time scales of the kinetics. In Section V we
study the dependence of the packaging efficiency on the topological and
geometrical properties of the genome molecules for different values of the
affinity ratio. In the concluding Section VI we discuss the possible relevance
of these results for our understanding of viral assembly as well as
limitations and possible extensions of the model.

##  II The Spanning Tree Model.

###  II.1 Empty capsid assembly

The Zlotnick model treats the capsid as a dodecahedral shell composed of
twelve pentamers. The pentamers represent pentameric oligomers of capsid
proteins – known as “capsomers” – of a minimal (i.e., T=1) capsid composed of
sixty proteins 444The assembly of larger icosahedral capsids can be
represented by a version of the model that allows for pentameric and hexameric
capsomers.. Assembly is driven by an attractive edge-edge interaction between
pentamers. A minimum-energy assembly pathway will be defined as a pentamer-by-
pentamer addition sequence where each added pentamer is placed in a location
that minimizes the free energy of the partial shell. An example of one of the
many (of the order of $10^{5}$) degenerate minimum-energy assembly pathways is
shown in Fig.2.

Figure 2: Empty capsid assembly pathway. The figure shows a minimum-energy
pathway for the assembly of a dodecahedral shell composed of twelve pentamers
with adhesive edges. The edge-to-edge binding energy is $\epsilon$. The change
in the total energy of the cluster for each added pentamer is indicated. Note
that partially assembled capsids are compact.

The Helmholtz free energy $E(n)$ of a partial shell composed of $n$ pentamers
is expressed as $E(n)=n_{1}\epsilon-n\mu$ with $n_{1}$ the number of shared
pentamer edges of the cluster times and with $\mu$ the chemical potential of a
pentamer in solution. Assuming dilute solutions, we will set $\mu=\mu_{0}+\ln
c_{f}/\bar{c}$ with $c_{f}$ the concentration of free pentamers in solution
and with $\mu_{0}$ the chemical potential at the reference pentamer
concentration $\bar{c}$. The assembly free energy of a complete capsid then
equals $30\epsilon-12\mu$ for all minimum energy assembly pathways. Assembly
equilibrium is defined as the state where the chemical potential of a pentamer
in solution is the same as that of a pentamer that is part of a capsid. This
is the case if $E(12)=0$ so if $\mu$ equals $\mu^{*}=(5/2)\epsilon$.

Figure 3 (top) shows the free energy of the partial assemblies shown in Fig.2
for three different values of the chemical potential $\mu$ close to $\mu^{*}$
and $\epsilon=-1$ (energy will be expressed in units of $k_{B}T$).

Figure 3: Top: Free energy profiles of the minimum assembly energy pathways of
the Zlotnick model. Blue dots: The chemical potential $\mu$ is slightly below
$\mu^{*}=-2.5$, the value of the chemical potential for assembly equilibrium.
Orange squares: $\mu$ is equal to $\mu^{*}$. Green diamonds: $\mu$ is slightly
above $\mu^{*}$. Bottom: Assembly free energy profiles according to continuum
theory. Solid red dots: energy minima. Solid black dots: energy maxima. For
$\mu<\mu^{*}$, the absolute energy minimum is at $n=0$ while for $\mu>\mu^{*}$
the absolute minimum is at $n=12$, the assembled capsid. The assembly
activation energy barrier of a profile is the height of the maximum. The
location $n^{*}$ of the maximum is for a half-filled shell under equilibrium
conditions (around n=6), shifting to lower values as $\mu$ increases.

Figure 3 (bottom) shows the free energy of a spherical cap growing into a
spherical shell [17]. The initial rise of $E(n)$ with $n$ is due to dominance
of the line energy of the perimeter of the cap. The line energy decreases as a
function of $n$ for larger $n$. The curve has the standard nucleation-and-
growth form. Under assembly equilibrium conditions $\mu=\mu^{*}$ the critical
nucleus is the half-formed shell with the size of the critical nucleus
decreasing with increasing levels of supersaturation.

###  II.2 Spanning Trees

The second part of the definition of the model is the specification of the
section of the encapsidated RNA molecule that interacts with the capsid
proteins. The secondary structure of this outer section is represented by a
tree graph. Tree graphs are collections of nodes connected by links such that
there is one and only one path of links connecting any pair of nodes [41] (see
Fig.4, right). The nodes represent in our case specific interaction sites
between an encapsidated RNA molecule and the interior of the dodecahedral
capsid. 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 [42]. The outer section of the
packaged RNA molecule is, in the model, assumed to be a spanning tree of the
dodecahedron. Each spanning tree of the dodecahedron has the same number of
vertices (twenty) and the same number of links (nineteen). A spanning tree
thus covers only nineteen of the thirty edges of the dodecahedron. The
remaining edges of the dodecahedron represent parts of the capsid that
associate with the genome only through non-specific interactions. It should
however be emphasized that the spanning tree model is not supposed to be a
realistic representation of any particular virus.

Figure 4 (left) shows an example spanning tree graph of the dodecahedron. The
projection of this spanning tree on the plane is shown on the right (we will
refer to this example as “molecule (1)”).

Figure 4: Left: Spanning tree “molecule (1)” connecting the vertices of a
dodecahedron (solid lines). The dashed lines indicate edges of the
dodecahedron that are not part of the spanning tree. Six pentamers, shown in
blue, can be placed on the dodecahedron with each pentamer wrapped by four
links of the tree. Right: Planar graph of the same spanning tree.

Now place pentamers on the spanning tree. Each pentamer interacts with five
edges of the dodecahedron. From the geometry of the spanning tree it is easy
to see that a pentamer can only interact with four links of a spanning tree.
If interaction of pentamer edges with the links of the spanning tree is
energetically favorable – as we will assume – then a cluster of pentamers
minimizes the interaction free energy between pentamers and spanning tree by
maximizing the number of pentamers in contact with four links of the spanning
tree. The figure shows (in blue) that a maximum of six pentamers can be
positioned in this fashion on the spanning tree.

###  II.3 Classification of Spanning Tree Molecules.

The set of spanning trees that leave the dodecahedron invariant (modulo
symmetry operations of the dodecahedron) has about $10^{5}$ members. A
topological characteristic that has been previously applied to viral RNA in
terms of packaging efficiency is the Maximum Ladder Distance (or MLD) [12,
43]. This is the maximum number of paired RNA nucleotides separating any two
nucleotides. The MLD of a spanning tree graph is defined as 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 a tree graph [41]. The MLD of the tree molecule shown
in Fig.4 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. Minimum MLD spanning
trees resemble Cayley trees while maximum MLD spanning trees are Hamiltonian
Paths, i.e., walks without self-intersection that visit all vertices of a
polyhedron [44, 45]. It has been shown that, in the absence of interactions,
the solution radius of gyration of a branched polymeric molecule increases
with the MLD as a power law [46]. A systematic comparison between the genomic
RNA molecules of RNA viruses confirms that they have significantly lower MLDs
than randomized versions of the same molecules [12, 43].

Figure 5 is a plot of the number of spanning trees of the dodecahedron as a
function of the MLD.

Figure 5: 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 twelve. By comparison, the
configurational entropy of an annealed branched polymer composed of nineteen
monomers that is not constrained to be a spanning tree depends on the MLD as
$19-MLD^{2}/19$ [46] and has a maximum at the smallest possible MLD. It
follows that the demand that a tree molecule with a certain number of links is
also a spanning tree of a dodecahedron greatly constrains the branching
statistics.

Could the MLD fully characterize a spanning tree? Figure 6 shows a linear
molecule with nineteen links, so the MLD is equal to nineteen. In the form of
a nineteen step Hamiltonian walk, it also is a spanning tree:

Figure 6: A possible Hamiltonian path for a linear genome molecule (“molecule
(2)”). Only two of the pentamers can be maximally wrapped.

For this “molecule (2)”, only two pentamers can be placed on the dodecahedron.
After enumerating all possible Hamiltonian paths, one finds that the wrapping
numbers of a Hamiltonian path can be two, three, or four (see Fig.8) yet all
these spanning trees have the same MLD. This indicates that the MLD does not
suffice as a useful characteristic of a spanning tree.

A characteristic that, unlike the MLD, is sensitive not only to the topology
but also to the folding geometry of a spanning tree is the wrapping number
(denoted by $N_{P}$). The wrapping number of a spanning tree of the
dodecahedron is the maximum number of pentamers that can be placed on the
spanning tree with four of its edges wrapped by a link of the spanning tree.
For the spanning tree shown in Fig.4, $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.7:

Figure 7: 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$. A compact cluster
of pentamers that are maximally wrapped by the spanning tree would be a
candidate of a nucleation complex with a low energy activation barrier.
However, the wrapping number does not take compactness into account. In
particular, it does not distinguish whether maximally wrapped pentamers have
shared links (Fig.6 is an example of maximally wrapped pentamers with no
shared edges). The wrapping number also does not seem to be a sufficient
characteristic to predict whether a spanning tree has a low assembly energy
activation barrier.

Are the wrapping number and MLD characteristics independent? 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).

There is a correlation between the two characteristics. For the largest MLDs
where the wrapping number is restricted to be between two and four this
correlation is strong while it is weaker for the smallest MLDs when the
wrapping number can adopt nearly its full range.

###  II.4 Assembly Energy Profiles.

Next we construct the minimum free energy assembly pathways and energy
profiles. The initial state will be a spanning tree molecule folded over the
edges of a mathematical dodecahedron with no pentamers. Physically, this
starting state can be viewed as representing a folded or pre-condensed form of
the viral ssRNA genome molecule(s) 555In actuality, condensation of the RNA
genome molecules takes place during encapsidation.. It will be assumed that
different spanning trees have the same folding energy prior to assembly so the
spanning trees have the same a-priori probability. Next, pentamers are placed
on the dodecahedron, one after the other. The free energy of a cluster of $n$
pentamers associated with spanning tree $i$ is defined as:

$E^{i}(n)=n_{1}\epsilon_{1}+n_{2}\epsilon_{2}+n_{3}\epsilon_{3}+n_{4}\epsilon_{4}-n\mu$
(1)

Here, $n_{1}$ is the number of links of the spanning tree that lie along a
pentamer edge that is not shared with another pentamer; $n_{2}$ is the number
of spanning tree 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 not covered by a spanning tree link and finally $n_{4}$ is the number
of pentamer edges that is not shared with another pentamer and that also is
not associated with a spanning tree link. The associated energies are given as
$\epsilon_{i}$ with $i=1,2,3,4$. An energy scale will be used in which the
energy $\epsilon_{4}$ of a bare edge is equal to zero. Energies are expressed
in units of an overall energy scale $\epsilon_{0}$, itself expressed in units
of the thermal energy $k_{B}T$. In these units, the interaction energy
$\epsilon_{3}$ between two pentamer edges is equal to $-1$ while the physical
meaning of $-\epsilon_{1}$ is that of the ratio of the affinity of a spanning
tree link with a pentamer edge over the pentamer-pentamer affinity. Next,
interactions between edges and spanning tree links are assumed to be additive
with $\epsilon_{2}=-1+2\epsilon_{1}$. This leaves $\epsilon_{1}$, $\mu$, and
$\epsilon_{0}$ as remaining free parameters. The total assembly energy of a
complete particle is equal to $\Delta E=19\epsilon_{2}-11-12\mu$ for all
spanning trees. The model is thus constructed so there is no thermodynamic
bias for the different spanning tree so packaging selectivity is necessarily
kinetic in origin.

####  II.4.1 Small $|\epsilon_{1}|$

Minimum free energy assembly pathways are constructed in the same way as
before. Figure 9 compares the minimum free energy assembly profiles for the
$MLD=9$, $N_{P}=6$ tree of Fig.4 (from here on referred to as “molecule (1)”)
and for the $MLD=19$, $N_{P}=2$ Hamiltonian walk of Fig.6 (from here on
referred to as “molecule (2)”) for $\epsilon_{1}=-0.2$, so in the regime of
relatively low affinity of RNA for capsid proteins

Figure 9: Minimum free energy assembly pathways for the $MLD=9$, $N_{P}=6$
spanning tree of Fig.4 (left, molecule(1)) and the $MLD=19$, $N_{P}=2$
spanning tree of Fig.6 (right, molecule (2)) at the reference concentration.
Energy parameters are $\epsilon_{1}=-0.2$, $\epsilon_{0}=1$, and $\mu=-2.5$.
The total assembly energy is indicated by a dashed red line while the assembly
activation barrier by a solid red line. The red arrows indicate the activation
energy for assembly (left) and disassembly (right).

Both energy profiles are similar to that of an empty capsid (i.e., the case
$\epsilon_{1}=0.0$) and also resemble each other, though the two molecules
have very different structures. The total assembly energies (solid red lines)
are the same by construction but the heights of the assembly energy activation
barriers also are practically the same. The same is true for the energy
barrier for disassembly separating the n=12 and n=2 states (right red arrow,
about 10.0). The main difference is that the width of the energy barrier for
molecule (2) is somewhat larger than that of molecule (1). If the assembly
kinetics is viewed as a random walk over the energy landscapes of Fig.9 then
molecule (1) is expected to have somewhat faster assembly kinetics than
molecule (2).

####  II.4.2 Large $|\epsilon_{1}|$.

Figure 10 shows the energy profiles of molecules (1) and (2) for
$\epsilon_{1}=-1.0$, corresponding to a higher RNA/protein affinity. The value
of $\mu$ was decreased from $-2.5$ to $-5.2$ in order to keep the activation
energy barrier in the same range as for $\epsilon=-0.2$.

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

The assembly activation energy of molecule (1) is more than twice as small as
that of molecule (2) while the energy profile has developed a second minimum
at $n=10$. Metastable intermediate states are familiar from experimental
studies of viral assembly [48, 49, 50] as well as from numerical simulations
[51, 52, 53, 54]. Known as “kinetic traps”, they may retard or block assembly.
Intermediate metastable states become progressively more pronounced as the
chemical potential $\mu_{0}$ is reduced towards assembly equilibrium as
illustrated in Fig.11 for the case of $\epsilon_{1}=-1.2$ of molecule (1)

Figure 11: Assembly energy profile of molecule (1) for $\epsilon_{1}=-1.2$
(bottom) under conditions of assembly equilibrium.

The energy profile has acquired a deep metastable minimum at $n=6$, which is
the exact opposite of what is expected from a nucleation-and-growth scenario.
Moreover, the minimum energy state is no longer the $n=12$ assembled viral
particle but rather the defected $n=11$ state. In short, for larger ratios of
the pentamer-tree molecule and pentamer-pentamer affinity, the nucleation-and-
growth scenario breaks down.

####  II.4.3 Partial Assemblies and Structural Transitions.

For small values of $|\epsilon_{1}|$, the structure of a partial assembly for
given $n$ is the same as that of the empty capsid, as shown shown in Fig.2.
However, for larger $|\epsilon_{1}|$ the first $N_{p}$ pentamers will be
placed on maximum wrapping sites. It follows that partial assemblies are
expected to undergo a structural transition as a function of $|\epsilon_{1}|$.
For molecule (1) this structural transition has secondary effects. For small
$|\epsilon_{1}|$, a six-pentamer clusters has five-fold symmetry with one
central pentamer sharing its five edges with five other pentamers that share
three edges with their neighbors (see Fig.2). For large $|\epsilon_{1}|$ the
minimum energy $n=6$ cluster has its six pentamers placed on the six available
maximum wrapping sites of an $N_{P}=6$ spanning tree shown in Fig.4. By moving
only one pentamer, this structure can be transformed into the small
$|\epsilon_{1}|$ $n=6$ structure with five-fold symmetry. The transition takes
place at $\epsilon_{1}=-1$. However, for molecule (2) with $N_{p}=2$, the
transition is more dramatic. The $n=4$ pentamer cluster of an empty caspid has
a two-fold symmetry axis (see Fig.2). Only two pentamers can be maximally
wrapped (see Fig.6). The full $n=4$ minimum energy structure is shown in
Fig.12

Figure 12: The minimum energy $n=4$ assembly state of molecule (2) for
$\epsilon_{1}=-1.2$.

This structural transition has an interesting feature. If one would relax the
constraints that keep the links of the genome molecule associated with the
edges of a mathematical dodecahedron, allow tree links to rotate freely around
the nodes of the tree, and finally allow pentamers to swivel around shared
edges then the pentamers of the empty-capsid partial assemblies of Fig.2 would
– for $n>2$ – not be able to move with respect to each other without breaking
pentamer-pentamer bonds. The empty-capsid partial assemblies can be said to be
mechanically rigid. The same is true as well for the $n=6$ structure of Fig.4.
However, this is not the case for the four-pentamer structure shown in Fig.12.
If this structure were released from the dodecahedron, and allowed to
fluctuate freely, then the four pentamers could freely swivel along the three
shared edges. We will call such an assembly flaccid. Transitions from rigid to
flaccid as a function of $\epsilon_{1}$ become are common for larger values of
the MLD. We would expect that actual flaccid RNA/capsid protein assembly
intermediates free in solution would be characterized by strong conformational
thermal fluctuations.

##  III Boltzmann Distribution.

In the configuration space of minimum energy partial assemblies, the Boltzmann
distribution for the pentamer occupation probabilities $P^{i}_{n}$ of spanning
tree $i$ is

$P^{i}_{n}\propto\exp-E^{i}(n)$ (2)

By substitution, it can be checked that

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

with $c_{f}$ the concentration of free pentamers, $r_{0}^{i}$ the
concentration of pentamer-free spanning trees of type $i$ and
$c_{12}^{i}=r^{i}P_{12}$ the concentration of assembled particles. Here,
$r^{i}$ is the total concentration of this type of tree molecules while
$r_{0}^{i}=r^{i}P_{0}$. Concentrations are all expressed in units of the
reference concentration $\bar{c}$. Finally, $K=\exp E_{0}(12)$ where
$E_{0}(12)$ is the particle assembly free energy for $c_{f}=1$. The
superscript $i$ was dropped here because $E_{0}(12)$ is independent of $i$.
This relation has the form of the Law of Mass Action (LMA) of physical
chemistry where $K$ is known as the dissociation constant.

###  III.1 Conservation Laws

The probability distribution $P^{i}_{n}$ is related to conservation laws for
the number of pentamers and of tree molecules. Conservation of tree molecules
requires that $r_{0}^{i}+\sum_{n=1}^{12}c_{n}^{i}=r^{i}$ where
$c_{n}^{i}=r^{i}P_{n}$ is the concentration of partial assemblies with $n$
pentamers. This conservation law is assured if the probabilities sum to one
for each type of tree molecules, i.e., if $\sum_{n=0}^{12}P^{i}_{n}=1$. Next,
conservation of pentamer molecules is assured if
$c_{f}=c_{0}\Gamma\thinspace([P])$ with $c_{0}$ the total pentamer
concentration and with

$\Gamma\thinspace([P])\equiv 1-(D/12)\sum_{i}\sum_{n=1}^{12}nP^{i}_{n}$ (4)

the probability that a pentamer is free in solution. The quantity
$D=12r_{0}/c_{0}$ plays the role of the mixing ratio, i.e., the ratio of the
total number of RNA molecules over the total number of capsid proteins. Here,
$r_{0}=\sum_{i}r^{i}$ is the total concentration of tree molecules. If $D=1$
then there are exactly enough pentamers to encapsidate all tree molecules. At
first sight it would seem that $D=1$ should awlays be the best mixing ratio if
one wanted to minimize the number of unpackaged RNA molecules and free
pentamers at the end of the assembly process. D=1 is similar to the
stoichiometric ratio of physical chemistry. Because the concentration of free
pentamers $c_{f}=c_{0}\Gamma\thinspace([P])$ depends on the Boltzmann
distributions of all types of tree molecules, the different distribution are
coupled.

###  III.2 Self-Assembly Diagram.

A standard diagnostic plot for self-assembly processes are plots of the
concentration of free building blocks and of assembled particles as a function
of the total concentration of building blocks [55]. Such a plot was computed
from the Boltzmann distribution Eq.2, assuming that only one type of tree
molecule as in solution, (molecule (1), with the energy parameters of Fig.9.
The result is shown in Fig.13. The dots show the concentrations of free
pentamers in solution and of pentamers that are part of an assembled particle,
both as a function of the total pentamer concentration $c_{0}$.

Figure 13: Equilibrium self-assembly diagram for molecule (1) for
$\epsilon_{1}=-0.2$, $\epsilon_{0}=1$, $\mu=-2.5$, and $D=1$. 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.5

For low pentamer concentrations, nearly all pentamers are free in solution so
the concentration of free pentamers is close to $c_{0}$. As $c_{0}$ increases,
the concentration of free pentamers stops increasing and then saturates. The
concentration of pentamers that are part of an assembled particle now starts
to increase proportional to $c_{0}$. The transition point, which is around
$c_{0}=0.2$ between these two regimes is known in the soft-matter physics
literature as the critical aggregation concentration (or CAC) [55].

The intermediate occupation probabilities $P_{n}$ with $2\leq n\leq 11$ had
low values for this parameter setting. If one entirely neglects these
intermediates then the conservation law for genome molecules reduces 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.3
produces an closed form expression 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)=\left(\frac{K}{c_{0}^{12}}\right)$
(5)

Importantly, like the LMA, this relation is independent of the type of tree
molecules (Explicit forms in certain limiting regimes are given in Appendix
B). Figure 13 compares Eq.5 (solid lines) with the numerical solution of the
thirteen coupled Boltzmann equations (dots) that does allow for assembly
intermediates.

###  III.3 Phase-Diagram.

A second way to display self-assembly measurements under equilibrium
conditions is in the form of a quasi phase-diagram that shows the dominant
type of assembly as a function of thermodynamic parameters 666Since a virus is
a system of limited size true phase transitions are not possible. For viral
assembly, the protein and RNA concentrations are a natural choice for such a
phase-diagram. For the case of the spanning-tree model, we will use the
pentamer concentration $c_{0}$ and the mixing ratio $D$ as thermodynamic
parameters. The blue dots in Fig.14 show points in a $c_{0}$ vs $D$ diagram
where 95 percent of the spanning trees of molecule (1) are fully encapsidated.

Figure 14: Quasi phase-diagram 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 have been packaged according to the Boltzmann distribution. Solid
blue line: same as computed from Eq.6. In the green sector there is
practically no capsid assembly. The red star marks a possible operating point
for viral assembly inside infected cells.

To the right of the blue dots, most tree molecules are encapsidated and
coexist with excess free pentamers. To the left of the blue dots most
pentamers are part of assembled particles and coexist with excess free tree
molecules. The blue dots can be viewed as optimal mixing states that leave
minimal excess free pentamers of tree molecules. For high pentamer
concentrations, the line of blue dots approaches $D=1$, which was the expected
optimal mixing state. This is however definitely not the case for smaller
values of $c_{0}$.

If one neglects assembly intermediates then it can be shown from Eq.5 that the
relation $c_{0}(D)$ for 95 percent occupancy is a hyperbola in the $c_{0}-D$
plane:

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

with $P_{(}12)=0.95$. The hyperbola diverges at $D=1/P_{12}$, which is close
to one for a 95 percent packaging fraction. It shifts to smaller values of $D$
as the pentamer concentration $c_{0}$ is reduced with $c_{0}(D)$ always larger
than $K^{1/12}$, which is of the order of the CAC. The entropy-dominated green
sector in Fig.14 is below the CAC.

Figure 15(Top) shows how the occupation probabilities change as one crosses
the boundary between the pentamer-rich and genome-rich states for fixed
$c_{0}$

Figure 15: Equilibrium occupation numbers as a function of the depletion
factor $D$ for fixed $c_{0}=1.0$. Top: same parameters as Fig.14. Bottom:
$-\epsilon_{1}$ is increased from $0.2$ to $1.2$ and the reference chemical
potential $\mu_{0}$ is reduced from $-2.5$ to $-5$ so the total assembly
energy remains approximately the same.

The transition region between the two states is quite broad. Note that for
large $D$ the occupation probabilities of assembly intermediates start to
become comparable to the occupation probability of assembled capsids. This is
due to the fact that as $D$ increases, more and more pentamer binding sites
become available. Breaking up assembled particles and distributing the
resulting pentamers over the additional binding sites increases the system
entropy. When $-\epsilon_{1}$ is increased to $1.2$ (see Fig.15(Bottom)), this
entropic effect becomes dominant and causes partially assembled particles to
overwhelm fully assembled particles for larger $D$. The structural transition
around $\epsilon_{1}=-1$ of assembly intermediates discussed in the previous
section here shows up as an order-disorder transition for mixing ratios larger
than one.

##  IV Assembly and Disassembly Kinetics.

###  IV.1 Master Equation.

The kinetics of the model is defined in terms of the thirteen time-dependent
occupation probabilities $P_{n}(t)$. We will again limit ourselves to the case
of a single spanning tree and drop the subscript $i$ 777The case of competing
spanning trees is discussed in the paper, submitted to Plos Biocomputation..
The assembly kinetics is assumed to obey Markov chain statistics [58], in
which case the occupation probabilities evolve in time according to the Master
Equation [59]:

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

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 only
with one pentamer gained or lost at a time. The matrix $W_{m,n}$ is defined in
terms of a simplified form of the diffusion-limited chemical kinetics of Eq.
8.35 of ref.[60]. The addition of a pentamer to a genome molecule associated
with a pentamer cluster of size $n$ is treated as a bi-molecular reaction with
a reaction rate $k_{n,n+1}{c_{n}}{c_{f}}$. The on-rate $k_{n,n+1}$ is given by

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

Here, $\lambda$ is a base rate that depends on molecular quantities such as
diffusion coefficients and reaction radii but not on the various
concentrations. If the addition of a pentamer reduces the free energy of a
partial assembly then the on-rate is equal to this base rate. If there is an
energy cost, then the base rate is reduced by an Arrhenius factor.
Importantly, the assembly energies $E(n)$ the on-rates should be computed here
with $c_{f}=c_{0}$.

In terms of the on-rates, the rate matrix is given by
$W_{n,n+1}=k_{n,n+1}{c_{f}}$ after canceling a common factor $r_{t}$. The
entries of the rate matrix that correspond to adding a pentamer are then

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

where $\Delta E_{n,n+1}\equiv E(n+1)-E(n)$ again with $c_{f}=c_{0}$. The
prefactor $\lambda$ has been absorbed in a redefinition of the time scale. As
for the equilibrium case

$\Gamma\thinspace([P])=1-(D/12)\sum_{n=1}^{12}nP_{n}(t)$ (10)

is the probability for a pentamer to be free in solution. 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}}$ (11)

where the right hand side the time-independent equilibrium Boltzmann
Distribution Eq. (2) must be inserted.

Typical examples of numerical solutions of the Master Equation are shown in
Fig.16 for molecule (1) for reference conditions $c_{0}=D=\beta=1$.

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

The top figure shows the case where the energy parameters are those of
Fig.13), so $\epsilon_{1}=-0.2$. The red curve is the probability $P_{12}(t)$
for a tree molecule to be encapsidated while the black curve is the
probability $P_{0}(t)$ for a tree molecules to be free of pentamers. The
probabilities for partial assemblies are small compared to $P_{12}(t)$ and
$P_{0}(t)$. Eventually, about forty percent of the tree molecules are
encapsidated in assembled particles while about five percent are associated
with partial assemblies. The dominant assembly pathway is the same as the one
shown in Fig.2. The occupation probabilities exponentially approach constant
values at late times that agree with the Boltzmann distribution. The bottom
figure shows what happens when the strength of the interaction between
pentamers and tree molecules is increased, with $\epsilon_{1}=-1.0$. The
reference chemical potential is reduced to $\mu_{0}=-5.2$ in order to keep the
activation energy roughly the same. The role of partial assemblies becomes
significantly larger.

###  IV.2 Time Scales for Assembly and Disassembly.

The next step is to determine the characteristic time scales. The occupation
probabilities display a rather complex kinetics, as can seen by amplifying the
scale, which is done in Fig.17.

Figure 17: Assembly shock-wave for molecule (1) for $\epsilon=-0.2$ in terms
of the occupation probabilities $P_{n}(t)$. 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 color coding is the same as that of Fig.17.

The functions $P_{n}(t)$ in Fig.17 have, for $n$ less than 12, a maximum as a
function of $n$. The peak times of the different distributions increase with
$n$. This describes a type of assembly shock-wave propagating in configuration
space from small to large $n$. Similar assembly shock waves have been reported
for the dynamical versions of the Zlotnick model discussed in Section II [37,
39]. A characteristic assembly delay time $t_{d}$ can be defined from the
early time behavior of $P_{12}(t)$ as the intercept of the maximum tangent of
$P_{12}(t)$ with the time axis. For the case of molecule (1), this gives
$t_{d}\simeq 11$. The delay time is roughly an average of the various peak
times. Assembly delay times are a familiar feature of the assembly kinetics of
empty capsids [15] and of aggregation phenomena in general [61]. Figure 18
shows that the assembly delay time is not strongly dependent on the pentamer
concentration.

Figure 18 shows the effect of varying the total pentamer concentration $c_{0}$
on the assembly delay time.

Figure 18: 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 arrow shows the assembly delay time.

The assembly delay time depends only weakly on the pentamer concentration.
Note that assembly stops when the pentamer concentration drops below a
threshold of about $0.2$, consistent with the CAC of equilibrium
thermodynamics (see Figs 13 and 14). The values of $P_{12}(t)$ at late times
also are consistent with the Boltzmann distribution 888Similar plots are
obtained by time-dependent light scattering studies of the assembly of empty
capsids with varying protein concentrations. [15]..

It is clear from Fig.18 that final thermal equilibration takes places on time
scales that are much longer than $t_{d}$ and that do depend strongly on
$c_{0}$. To compute the equilibration time, complete the definition of the
transition matrix by introducing diagonal entries
$W^{\prime}_{n,n}=-\sum_{m\neq n}W(m,n)$. The column elements of the resulting
matrix $W^{\prime}_{m,n}$ add to zero. The Master Equation can be rewritten in
terms of $W^{\prime}_{m,n}$ as $\frac{d\bf{P}}{dt}=\bf{W^{\prime}P}$. The
eigenvalues of the completed transition matrix $W^{\prime}_{m,n}$ are then
(minus) the decay rates of the different modes of relaxation. These modes
correspond to the eigenvectors of $W^{\prime}_{m,n}$. The smallest eigenvalue
of $W^{\prime}_{m,n}$ must correspond to the Boltzmann distribution and has to
be zero. The rate of approach to final thermal equilibrium is the negative of
the next smallest eigenvalue of $W^{\prime}_{m,n}$.

The largest eleven eigenvalues of $W^{\prime}_{m,n}$ are all of the order of
$1/t_{d}$. The twelfth eigenvalue $1/t_{r}$ is however much smaller, with
$t_{r}\simeq 644$ while the thirteenth eigenvalue is indeed consistent with
zero within our numerical precision. We will refer to $t_{r}$ as “the”
relaxation time. The dependence of these two key time scales on the overall
energy scale $\epsilon_{0}$ is quite different. For $\epsilon_{0}=3$, $t_{d}$
barely changes ($t_{d}\simeq 11.7$) but the relaxation time increases
exponentially to $t_{r}\simeq 1.7\times 10^{5}$. Figure 18 shows that the
relaxation time strongly increases as the pentamer concentration decreases.

In order to investigate how assembled particles “fall apart” when the pentamer
concentration is lowered below the CAC, we used the outcome of an assembly
calculation with $c_{0}$ above the CAC as the initial condition for a second
calculation with $c_{0}$ reduced below the CAC. The results are shown in Fig.
19, still for molecule (1).

Figure 19: Top: The assembly calculation of Fig. 16 at $c_{0}=1$ is followed
by a disassembly calculation at $c_{0}=0.1$. The color coding is that of
Fig.15. Bottom: Same except that the overall energy scale is increased to
$\epsilon_{0}=3$. The arrows indicate the relaxation time $t_{r}$.

The assembled particles disintegrate on time scales of the order of the
relaxation time $t_{r}\simeq 31.6$. Repeating the calculation for
$\epsilon_{0}=3$ (bottom figure), disassembly is now much slower. The
characteristic disassembly time scale remains of the order of the relaxation
time, now with $t_{r}\simeq 5\times 10^{5}$. Assembled viral particles must be
at least metastable in solutions where the concentration of capsid proteins is
well below the CAC where they are subject to an entropic thermodynamic drive
towards disassembly. In the following, we will focus on time scales large
compared to the assembly delay time but short compared to the relaxation time,
which governs the disassembly process. This is possible only if there is a
large separation of time scales, so for values of $\epsilon_{0}$ that
significantly larger than one. We will investigate how the topology and
geometry of the spanning trees affects assembly under these conditions.

##  V Single-Stage and Two-Stage Assembly.

We first examined normal pentamer-by-pentamer (“single-stage”) assembly.

###  V.1 Single-Stage Assembly.

A full survey was performed for every spanning tree of the probability
$P_{12}(\tau)$ of a tree molecule to be fully encapsidate at a time
$\tau=5\times 10^{5}$. This was done for the parameter setting
$\epsilon_{0}=3$, $\epsilon_{1}=0.2$, $\epsilon_{3}=-1.0$, $\mu=-2.5$,
$\epsilon_{0}=3$, $c_{0}=0.22$ and $D=0.3$ since for this setting, (i) partial
assemblies are compact and rigid, (ii) the nucleation-and-growth scenario
applies, (iii) the system is close to the CAC (as is probably the biologically
relevant regime, see Conclusion) with few assembly intermediates, and (iv)
there is large separation between the assembly time scale $t_{d}$ and the
disassembly time scale $t_{r}$. For these values, $\tau$ is about
$0.25\thinspace t_{r}$. Finally, for every pair of MLD and $N_{p}$ numbers,
the average of $P_{12}(\tau)$ and the standard deviation were computed over
all spanning trees that belonged to that category. These averages were then
plotted in the form of a three-dimensional bar graph with the MLD and $N_{P}$
numbers as the horizontal axes, as shown in Fig.20.

Figure 20: Plot of the probability $P_{12}(t=5\times 10^{5})$ (vertical axis)
that a genome molecule is fully packaged. The parameters are
$\epsilon_{1}=0.2$, $\epsilon_{3}=-1.0$, $\mu=-2.5$, $\epsilon_{0}=3$,
$c_{0}=0.22$ and $D=0.3$.

The largest value of $P_{12}(\tau)$ that we encountered was for MLD=9 and
$N_{p}=8$ with $P_{12}(\tau)=0.155\pm 0.016$. There were practically no
assembly intermediates so the remaining tree molecules were not associated
with pentamers, as was the case for other MLD and $N_{p}$ pairs. The lowest
value of $P_{12}(\tau)$ was $0.033$ for MLD=16 and $N_{p}=2$. In general, the
value of $P_{12}(\tau)$ systematically decreased when $N_{p}$ was decreased
for given MLD. For given $N_{p}$ larger than two, the packaging probability
decreased with increasing MLD, at least for MLD less than about thirteen. For
MLD larger than thirteen the packaging probability was not much dependent on
the MLD. This was also the case for $N_{p}=2$. The outliers that did not
follow these trends, such as $MLD=11$ with $N_{p}=4$, were consistent with the
statistical error bars. Single-stage assembly selects reasonably well for
spanning-trees with large Wrapping Number and small MLD.

###  V.2 Two-Stage Assembly.

In the introduction it was pointed out that a second mode of assembly – “en-
masse” – is encountered for higher values of the affinity between the RNA
molecules and the capsid proteins. In the first stage, a disordered protein-
RNA condensate forms. This condensate undergoes an ordering transition, either
after a certain period of time has elapsed or after an increase in the
protein-protein affinity induced by a change in the thermodynamic parameters
of the system. In order to mimic the en-masse assembly scenario, we carried
out assembly calculations in two steps. During the first step, pentamers had a
significant affinity for the tree molecule while the interaction between
pentamers was made weakly repulsive. The outcome of this first calculation was
then used as the starting state of a second calculation now with standard
parameter settings.

We again used molecule (1). During the first assembly stage $\epsilon_{3}$ was
equal to $+0.5$, which means that the pentamer-pentamer interaction was
repulsive. The reference chemical potential $\mu_{0}$, which includes a
contribution associated with conformational changes required to allow for
bonding between proteins, was set to zero. During the second assembly stage,
the pentamer-pentamer affinity and the reference chemical potential were reset
to the values of Fig.20. The assembly energy profile and occupation
probabilities of the first stage are shown in Fig.21.

Figure 21: Top: Assembly energy profile for the first assembly stage of a
molecule (1) for $\epsilon_{1}=-0.2$, $\epsilon_{3}=0.5$, $\epsilon_{0}=3$,
$D=1$, $c_{0}=1$ and $\mu=0$. Bottom: Minimum energy assembly path.

The top figure shows the minimum energy assembly profile at the reference
concentration. The minimum at $n=4$ can be understood from the bottom figure,
which shows the placing sequence of the first four pentamers. The first two
pentamers can be placed on maximum wrapping sites without being forced to
share any edges, which would be energetically unfavorable. Even though
additional maximum wrapping sites remain available, the third pentamer is
placed on a site where it makes contact with only three links of the tree
molecules in order to avoid pentamer-pentamer repulsion. Even that option is
no longer possible for the fourth pentamer, which is forced to share an edge
with the first pentamer. The assembly energy no longer decreases significantly
when the fourth pentamer is added. Note that the $n=4$ structure is flaccid,
as is typical for first-stage structures.

The occupation probabilities for the first stage, computed from this energy
profile, are shown in Fig.22 for times less than $t=20.000$. The parameter
settings are the same as Fig.21, except that the concentration was reduced
from the reference concentration to $c_{0}=0.22$, which is at the border of
the CAC. This increases the free energy of the $n=4$ structure with respect to
the $n=3$ structure by an amount $-\ln 0.22$ and by $-2\ln 0.22$ for the $n=2$
structure. As a result, the occupation probabilities of the $n=2$ and $n=3$
structures are increased with respect to that of the $n=4$ structures, as can
be seen in Fig.21. The occupation probabilities are time-independent,
indicating that the system has equilibrated, which is expected since the delay
time $t_{d}$ and the relaxation time $t_{r}$ are only of the order of $10$.

Next, at $t=20.000$, the parameters were returned to same values as for
single-stage assembly (so with $\epsilon_{3}=-1$ and with $\mu=-2,5$). The
parameter reset changes all minimum-energy assembly and disassembly pathways
so for every one of the thirteen starting states produced by the first
assembly stage at $t=20.000$, a different minimum-energy assembly/disassembly
pathway must be generated for the second pathway. The outcome of such a two-
stage assembly calculation is shown in Fig.22

Figure 22: Two-stage assembly of molecule (1). The parameter settings of the
first stage for times less than $t=20.000$ are those of Fig.21 except that
$c_{0}=0.22$. At $t=20,000$, the parameters of the second stage were set to
those of Fig.20.

Fully assembled particles appear practically immediately at $t=20.000$. The
occupation probabilities remain constant for times larger than $t=20.000$.
This is rather surprising because the relaxation time now is longer than
$20.000$ (about $10^{5}$). The kinetics of two-stage assembly somehow does not
seem to “see” the relaxation time. While the solution was poly-disperse for
times less than $t=20.000$, for times greater than $t=20.000$ the tree
molecules are either fully encapsidated (red line, about 5 percent) or
pentamer free (black line, about 95 percent). There were no assembly
intermediates. Next, the calculation was repeated except that the reference
chemical potential was reset from $\mu_{0}=-2.5$ to the assembly equilibrium
condition with $\mu_{0}=-3.14$. Curiously, this did not change the results for
the second stage in Fig.22.

These results can be understood by comparing the second-stage minimum energy
assembly profiles.

Figure 23: Second-stage assembly energy profiles of the thirteen different
states produced by the first stage. The index $N_{s}$ gives the occupancy at
the start of the second stage. Only the profiles for $N_{s}=0,1,2,3,4$ are
shown. The energy parameters of the second stage are $\epsilon_{1}=-0.2$,
$\epsilon_{3}=-1.0$, $\epsilon_{0}=3$, with $\mu_{0}=-2.5$ (top) and
$\mu_{0}=-3.1$ (bottom). The starting points for the second stage are
indicated by circles. The vertical energy scale is in units of $\epsilon_{0}$.

These plots show the minimum energy assembly profiles for different values of
$N_{s}$ (color indexed), which is the number of pentamers that are initially
present on the tree molecule as a result of the first assembly stage. The
probabilities for the different $N_{s}$ states are simply the occupation
probabilities of the first stage prior to the parameter switch. Second stage
kinetics thus always starts from the state $n=N_{s}$ (circled). For example,
in the top figure the initial state $n=4$ for the $N_{s}=4$ energy profile,
shown in purple, is circled. Note that $n=4$ is the maximum of the energy
profile (“transition state”). Imagine the state of the system as a particle
performing a random walk along this energy profile, starting at $n=4$. At the
first time step, the particle can hop either to $n=3$ or to $n=5$. Both are
energetically downhill. For our kinetic scheme, the partition ratio
$\frac{W_{4,5}}{W_{4},3}$ between the rates to hop forwards – from $n=4$ to
$n=5$ – or backwards – from $n=4$ to $n=3$ – is fixed by detailed balance to
equal the ratio $(P_{4}/P_{3})\exp\Delta E_{3,4}$, which in turn is equal to
$c_{f}$. The left/right partition ratio for the $n=2$ and $n=3$ initial state,
both of which are transition states, is also $c_{f}$. The energy cost of
uphill steps is about $5\epsilon_{3}\simeq 15k_{B}T$. If one neglects such up-
hill steps, then the particle will end up either at $n=0$ (i.e., free tree
molecules) or at $n=12$ (i.e., encapsidated tree molecules) with a partition
ratio that is equal to the concentration $c_{f}\simeq c_{0}$ of free pentamers
in the equilibrium state (in units of the reference concentration).

The bottom plot shows the case where the chemical potential is tuned to the
state of assembly equilibrium. For single-stage assembly, starting at $n=0$,
the energy activation barrier would be a prohibitive $30\epsilon_{3}\simeq
100k_{B}T$. However, because assembly now starts (mostly) at $N_{s}=3$ or
$N_{s}=2$, there is a reasonable probability that the particle will reach the
$n=12$ state. The partition ratio is not exactly the same as for the previous
case because the $n=2$ state and the $n=3$ state no longer are transition
states. Adding a pentamer is now in both cases an energetically uphill step.
However, the energy cost of the uphill step is rather modest in this case
(about $2\epsilon_{3}$). This should lead to only a modest reduction in the
fraction of assembled viral particles, provided $c_{f}$ remains the same. This
last condition is consistent with the observation that the partition ratio did
not change significantly with respect to the previous case.The physical
mechanism that allows for this amplification effect is that when a few
pentamers are already present on the tree molecule prior to the parameter
switch, because of the first stage “priming” process, then these pentamers
will act as nucleation sites for the completion of the capsid. The production
yield of assembled particles can be much larger in the two-stage process.

###  V.3 Two-Stage Selectivity.

Finally, we investigated the selectivity of two-stage assembly in terms of the
Wrapping Number and the MLD. A full survey of tens of thousands of spanning
trees was not practical because two-stage assembly is computationally quite
intensive. Instead, single representatives were picked from each pair of
Wrapping Number and MLD. The results are shown in Fig.20

Figure 24: Plot of the probability $P_{12}(t_{r})$ (vertical axis) that a tree
molecule is fully packaged for two-stage assembly. The parameters are the same
as for Fig.20

For low values of $N_{p}$, there appears to be no selectivity with respect to
the MLD. Similarly, for low values of the MLD, there appears to be no
selectivity with respect to $N_{p}$ except that $N_{p}=3$ and $N_{p}=8$ have
lower packaging probability. While there is no clear selectivity, the
packaging probabilities are on average higher for $N_{p}$ between four and
seven and MLD less than twelve.

##  VI Conclusion

We analyzed the thermodynamics and kinetic properties of the spanning-tree
model for dodecahedral particles using numerical solutions of a set of coupled
master equations.

The dependence on the time-dependent occupation probabilities of partial and
complete assemblies on the topology and geometry of the spanning-tree
molecules – which represent the outer part of condensed viral RNA – was shown
to be largely determined by two correlated criteria: the maximum ladder
distance (MLD), which is a topological measure of the degree of branching of
the spanning tree, and the Wrapping Number ($N_{P}$), which is geometrical
measure of the spanning tree that counts the number of sites that can
accommodate the pentamers – which represent capsid protein oligomers – with
maximum binding energy to the spanning tree.

The assembly kinetics is characterized by two important time scales: the delay
time $t_{d}$ for the onset of production of fully assembled particles and
$t_{r}$, the relaxation time for thermal equilibration. While delay time
increases modestly as a function of the energy scale $\epsilon_{0}$ of the
interactions, the relaxation increases exponentially. When the solution
concentration $c_{0}$ of pentamers is lowered below the critical aggregation
concentration of the equilibrium phase diagram, assembled particles remain
stable against spontaneous disassembly on time scales shorter than $t_{r}$.
The dependence of the probability for a spanning tree to be completely
packaged on the $N_{P}$ and MLD numbers is a purely kinetic effect that
disappears on time-scales large compared to the relaxation time. This
selectivity is due to small differences between the activation energy barriers
spanning trees with different $N_{P}$ and MLD (“selective nucleation”).
Selective nucleation only “works” if the assembly activation energy barrier is
large compared to the thermal energy.

As a function of the ratio $|\epsilon_{1}|$ between protein/RNA and
protein/protein affinity, partial assemblies undergo a structural transition
around $|\epsilon_{1}|=1$ from predominantly rigid to predominantly flaccid
structures, which has the character of an order-disorder transition in the
equilibrium phase diagram. For $|\epsilon_{1}|$ small compared to one,
assembly kinetics has the character of the classical nucleation-and-growth
scenario under conditions of supersaturation. This scenario breaks down for
$|\epsilon_{1}|$ of the order of one or larger. At the same point, assembled
particles start to readily disintegrate when the pentamer concentration is
reduced after assembly. The minimum-energy assembly pathway is then
characterized by pronounced kinetic traps.

Finally, we carried out an exhaustive survey of the assembly kinetics of all
spanning trees and found that, for direct assembly, the probability for the
appearance of fully assembled particles progressively increases with
increasing Wrapping Number and decreasing MLD. When direct single-stage
assembly was replaced by two-stage assembly – as a model for the en-masse
scenario – then the yield of fully assembled particles sharply increased, even
under conditions of assembly equilibrium where single-stage assembly is
blocked. However, the dependence of two-stage assembly on the $N_{P}$ and MLD
numbers is weak and non-systematic as compared to single-stage assembly.

We mentioned in the introduction 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 [63], which reported that
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 exceeded a threshold, the virus-like
particles were replaced by disordered RNA-protein aggregates instead of the
expected coexistence of VLPs with excess RNA while the threshold point
separating the two regimes was well below the stoichiometric ratio. In the
spanning tree model, the mixing ratio is represented by $D$. If $D$ is
increased at fixed total protein concentration $c_{0}$ in the quasi phase-
diagram Fig.14 then for smaller $D$ practically all tree molecules are
encapsidated (to the left of the blue line in Fig.14) while for larger values
of $D$ an increasing fraction of the tree molecules are not encapsidated. As
$|\epsilon_{1}|$ increases, a significant fraction are partially encapsidated.
The transition point separating the two regimes drops below the stoichiometric
ratio when $c_{0}$ approached the critical aggregation concentration. These
results are consistent with the experimental observations.

Another important motivation was the asymmetric reconstruction of the MS2
virus, which provided evidence that a subsection of the RNA genome that is
rich in packaging signals, reproducibly associates with a compact cluster of
capsid proteins [28]. Our model reproduces this observation for genome
molecules with small MLD and large $N_{P}$, in which case the intermediate
states that regulate the assembly rate are compact and rigid for smaller
values of $\epsilon_{1}$. In the spanning-tree model, the role of packaging
signals is played by the $N_{p}$ maximum wrapping sites. This reinforces the
notion that packaging selection is carried not by individual RNA stem-loops
but rather by spatially correlated patterns of stem-loops, as proposed in
ref.[28].

An important aspect of the model that should be improved is the fact that it
does not account for conformational fluctuations of the genome molecules,
prior to assembly, as well as thermal fluctuations of flaccid intermediates.
For the same reason, the model also does not include the “antenna” effect for
the diffusive influx of capsid proteins to partial assemblies [64]. Another
aspect of the model that should be improved in the context of asymmetric
reconstructions is a more realistic treatment of the the condensed RNA
molecule. Such a generalization would start from a determination of the key
interaction sites between capsid proteins and the surface of the enclosed RNA.
These sites would span a polyhedron with many more sites than the vertices of
a dodecahedron. The spanning tree should reproduce the outer density of the
RNA molecule, such as Fig.1. Classifying the set of all spanning trees would
be significantly more challenging and require more extensive numerical work.

###### Acknowledgements.

We would like to thank Alexander Grosberg for drawing our attention to
spanning trees and Ioulia Rouzina for introducing us to the concept of
selective nucleation. We would like to thank Reidun Twarock, Charles Knobler
and William Gelbart for reading a first draft and commenting on it. We also
benefitted from discussions with Chen Lin, Zach Gvildys and William Vong. RB
would like to thank the NSF-DMR for continued support under CMMT Grant
No.1836404.

## References

  * [1] I. Mizrahi, R. Bruinsma and J. Rudnick, submitted to Plos Biocomputation.
  * 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).
  * 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).
  * 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).
  * Bruinsma  _et al._ [2021] R. F. Bruinsma, G. J. Wuite, and W. H. Roos, Nature Reviews Physics , 1 (2021).
  * 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).
  * Hu and Shklovskii [2007] T. Hu and B. I. Shklovskii, Phys. Rev. E 75, 051901 (2007).
  * Kler  _et al._ [2012] S. Kler, R. Asor, C. Li, A. Ginsburg, D. Harries, A. Oppenheim, A. Zlotnick, and U. Raviv, Journal of the American Chemical Society 134, 8823 (2012).
  * Garmann  _et al._ [2019] R. F. Garmann, A. M. Goldfain, and V. N. Manoharan, Proceedings of the National Academy of Sciences 116, 22485 (2019).
  * 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).
  * Johnson  _et al._ [2004] J. M. Johnson, D. A. Willits, M. J. Young, and A. Zlotnick, J. Mol. Biol. 335, 455 (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).
  * Dai  _et al._ [2017] X. Dai, Z. Li, M. Lai, S. Shu, Y. Du, Z. H. Zhou, and R. Sun, Nature 541, 112 (2017).
  * Dimmock  _et al._ [2001] N. Dimmock, A. Easton, and K. Leppard,  _Introduction to modern virology_ (Blackwell Publishing, Malden, MA, 2001).
  * Note [3] Selective nucleation was proposed by I. Rouzina in the context of the assembly of retroviruses.
  * 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).
  * Garmann  _et al._ [2013] R. F. Garmann, M. Comas-Garcia, A. Gopal, C. M. Knobler, and W. M. Gelbart, J. Mol. Biol. 32, ??? (2013).
  * Van Rosmalen  _et al._ [2020] M. G. Van Rosmalen, D. Kamsma, A. S. Biebricher, C. Li, A. Zlotnick, W. H. Roos, and G. J. Wuite, Science advances 6, eaaz1639 (2020).
  * Zlotnick [1994] A. Zlotnick, J. Mol. Biol. 241, 59 (1994).
  * Endres and Zlotnick [2002] D. Endres and A. Zlotnick, Biophys. J. 83, 1217 (2002).
  * Zlotnick [2007] A. Zlotnick, J. Mol. Biol. 366, 14 (2007).
  * Morozov  _et al._ [2009] A. Y. Morozov, R. F. Bruinsma, and J. Rudnick, J. Chem. Phys. 131, 155101 (2009).
  * Note [4] The assembly of larger icosahedral capsids can be represented by a version of the model that allows for pentameric and hexameric capsomers.
  * 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).
  * 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.
  * Parent  _et al._ [2006] K. N. Parent, A. Zlotnick, and C. M. Teschke, J. Mol. Biol. 359, 1097 (2006).
  * Tuma  _et al._ [2008] R. Tuma, H. Tsuruta, K. H. French, and P. E. Prevelige, J. Mol. Biol. 381, 1395 (2008).
  * Basnak  _et al._ [2010] G. Basnak, V. L. Morton, O. Rolfsson, N. J. Stonehouse, A. E. Ashcroft, and P. G. Stockley, J. Mol. Biol. 395, 924 (2010).
  * Johnston  _et al._ [2010] I. G. Johnston, A. A. Louis, and J. P. K. Doye, J. Phys.: Condens. Matter 22, 104101 (2010).
  * Hagan  _et al._ [2011] M. F. Hagan, O. M. Elrad, and R. L. Jack, J. Chem. Phys. 135, 104115 (2011).
  * Baschek  _et al._ [2012] J. E. Baschek, H. C. R. Klein, and U. S. Schwarz, Bmc Biophysics 5, 10.1186/2046-1682-5-22 (2012), baschek, Johanna E. Klein, Heinrich C. R. Schwarz, Ulrich S. Schwarz, Ulrich/K-4111-2014 Schwarz, Ulrich/0000-0003-1483-640X.
  * Perlmutter and Hagan [2015] J. D. Perlmutter and M. F. Hagan, J. Mol. Biol. 10.1016/j.jmb.2015.05.008 (2015).
  * Safran [2018] S. Safran,  _Statistical thermodynamics of surfaces, interfaces, and membranes_ (CRC Press, 2018).
  * Note [6] Since a virus is a system of limited size true phase transitions are not possible.
  * Note [7] The case of competing spanning trees is discussed in the paper, submitted to Plos Biocomputation.
  * 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).
  * Note [8] Similar plots are obtained by time-dependent light scattering studies of the assembly of empty capsids with varying protein concentrations. [15].
  * 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).
  * Hu  _et al._ [2006] T. Hu, A. Y. Grosberg, and B. I. Shklovskii, Biophys. J. 90, 2731 (2006).

##  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 25 shows such a path.

Figure 25: 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 26 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 26: 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.26. 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.

##  Appendix B Limiting solutions of Eq.5

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

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

For $\left(\frac{c_{0}^{12}}{K}\right)$ large compared to one and mixing ratio
$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}$ (13)

Finally, for $\left(\frac{c_{0}^{12}}{K}\right)$ large compared to one but the
mixing ratio $D$ less than one, the equation has a different 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}}$ (14)

There is thus a change in regimes near the point where the mixing ratio 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)$
(15)