# Markov State Model Approach to Simulate Self-Assembly

Anthony Trubiano<EMAIL_ADDRESS>Michael F. Hagan<EMAIL_ADDRESS>Martin Fisher School of Physics, Brandeis University, Waltham, Massachusetts
02454, USA

(August 31, 2024)

###### Abstract

Computational modeling of assembly is challenging for many systems because
their timescales vastly exceed those accessible to simulations. This article
describes the MultiMSM, which is a general framework that uses Markov state
models (MSMs) to enable simulating self-assembly and self-organization on
timescales that are orders of magnitude longer than those accessible to brute
force dynamics simulations. In contrast to previous MSM approaches to
simulating assembly, the framework describes simultaneous assembly of many
clusters and the consequent depletion of free subunits or other small
oligomers. The algorithm accounts for changes in transition rates as
concentrations of monomers and intermediates evolve over the course of the
reaction. Using two model systems, we show that the MultiMSM accurately
predicts the concentrations of the full ensemble of intermediates on the long
timescales required for reactions to reach equilibrium. Importantly, after
constructing a MultiMSM for one system concentration, a wide range of other
concentrations can be simulated without any further sampling. This capability
allows for orders of magnitude additional speed up. In addition, the method
enables highly efficient calculation of quantities such as free energy
profiles, nucleation timescales, flux along the ensemble of assembly pathways,
and entropy production rates. Identifying contributions of individual
transitions to entropy production rates reveals sources of kinetic traps. The
method is broadly applicable to systems with equilibrium or nonequilibrium
dynamics, and is trivially parallelizable and thus highly scalable.

##  I Introduction

The self-assembly of basic subunits into larger, more complex structures is
fundamental to life. Critical functions of cells and pathogens are performed
by assembled structures such as the outer shells (capsids) of viruses Caspar
and Klug (1962); Mateu (2013); Bruinsma and Klug (2015); Perlmutter and Hagan
(2015); Hagan and Zandi (2016); Twarock  _et al._ (2018); Zandi  _et al._
(2020); Hagan and Grason (2021) or bacterial microcompartments Kerfeld  _et
al._ (2010); Tanaka  _et al._ (2008); Rae  _et al._ (2013); Bobik  _et al._
(2015); Chowdhury  _et al._ (2014); Kerfeld and Melnicki (2016); Polka  _et
al._ (2016); Kerfeld  _et al._ (2018); Slininger Lee and Tullman-Ercek (2017),
cytoskeletal filaments Banerjee  _et al._ (2020); Cabeen and Jacobs-Wagner
(2010); Pilhofer and Jensen (2013), and ordered protein layers on bacteria
exteriors Whitelam (2010). Self-assembly is also transforming nanotechnology,
where designing synthetic building blocks that are preprogrammed to form
particular structures is enabling scalable bottom-up synthesis of materials
with desirable properties Garg T (2015); Beija  _et al._ (2012); Ebbens
(2016); Mallory  _et al._ (2018); Fan  _et al._ (2011); Huh  _et al._ (2020);
Ke  _et al._ (2012); Kraft  _et al._ (2009); Sacanna  _et al._ (2010, 2013);
Yi  _et al._ (2013); Wang  _et al._ (2014); He  _et al._ (2020, 2021); Chen
_et al._ (2011a, b); Zerrouki  _et al._ (2008); Yan  _et al._ (2013); Wolters
_et al._ (2015); Tikhomirov  _et al._ (2018); Oh  _et al._ (2019); Ben-Sasson
_et al._ (2021); Kahn  _et al._ (2022).

Since it is an inherently out-of-equilibrium process, understanding or
designing self-assembly processes requires detailed knowledge of assembly
intermediates and the dynamical transitions among them. Computational modeling
is an essential tool for revealing such assembly pathways, since most
intermediates are too transient to characterize in experiments. However,
simulating assembly at experimentally relevant conditions is intractable for
many models, since assembled structures are much larger than their components
and form on timescales that are orders of magnitude beyond computational
limitations. This article describes a Markov state model (MSM) framework that
overcomes this limitation for a broad array of self-assembly and self-
organization systems. The algorithm reduces computational times by orders of
magnitude while describing the time-dependent concentrations of subunits and
the complete ensemble of assembly intermediates and products. This capability
enables dynamical particle-based simulations of systems with unprecedented
size and complexity, at experimentally relevant conditions. The framework also
enables highly efficient analysis of the resulting simulation data.

MSMs are a powerful approach to simulate dynamics on long timescales; by
performing short simulations to estimate transition rates among system
configurations, one can construct an MSM that accurately describes dynamics on
timescales that are orders of magnitude longer than the individual simulations
Noé  _et al._ (2009); Bowman  _et al._ (2009); Pande  _et al._ (2010); Sarich
_et al._ (2010); Bowman  _et al._ (2011); Prinz  _et al._ (2011); Voelz  _et
al._ (2011, 2012); Noé and Nüske (2013); De Sancho  _et al._ (2013); Bowman
_et al._ (2014); Chodera and Noé (2014); Malmstrom  _et al._ (2014); Schwantes
_et al._ (2014); Hummer and Szabo (2015); Husic and Pande (2018); Zeng  _et
al._ (2018); Wu and Noé (2020); Weng  _et al._ (2020); Suárez  _et al._
(2021). MSMs also provide a means to coarse-grain complex dynamical processes
into reduced-order forms that facilitate identifying key slow degrees of
freedom and corresponding mechanisms. Furthermore, MSMs can enable designing
non-equilibrium assembly protocols that can accelerate assembly and increase
selectivity of a specific target state by orders of magnitude in comparison to
equilibrium processes Juárez and Bevan (2012); Tang  _et al._ (2016, 2017);
Grover  _et al._ (2019); Trubiano and Hagan (2022).

In contrast to previous MSM approaches to self-assembly that pre-assume the
state space and transition rates Jamalyaria  _et al._ (2005); Keef  _et al._
(2006); Hemberg  _et al._ (2006); Dykeman  _et al._ (2013); Sweeney  _et al._
(2008); Zhang and Schwartz (2006a); Misra  _et al._ (2008); Kumar and Schwartz
(2010); Xie  _et al._ (2012); Smith  _et al._ (2014), we seek a framework in
which the state space and transition rates are computed directly from
dynamical simulations, and the accuracy of the resulting MSM (including the
validity of the Markov assumption) can be directly tested against microscopic
dynamics. While several approaches have been developed to construct MSMs from
particle-based assembly simulations Kelley  _et al._ (2008); Perkett and Hagan
(2014); Dibak  _et al._ (2018); Olsson  _et al._ (2017); Sengupta  _et al._
(2019); Trubiano and Holmes-Cerfon (2021); Trubiano and Hagan (2022), these
algorithms are designed to track individual assembling clusters evolving under
constant conditions such as the concentration of free subunits. Thus, they
cannot describe a typical experiment in which a fixed total number of subunits
assemble into many structures. In this case, the concentrations of
intermediates and free subunits, and thus the transition rates, continuously
evolve over time. Moreover, some transitions involve association or
dissociation of oligomers or larger intermediates. Therefore, although
sophisticated approaches have been recently developed to design optimal
assemblies and compute free energy landscapes Miskin  _et al._ (2016);
Pietrucci (2017); Sidky and Whitmer (2018); Sherman  _et al._ (2020); Das and
Limmer (2021); Hénin  _et al._ (2022); Wang  _et al._ (2022); Das and Limmer
(2023); Goodrich  _et al._ (2021); Lieu and Yoshinaga (2022); Jhaveri  _et
al._ (2023); Curatolo  _et al._ (2023), to our knowledge, there is no existing
enhanced sampling method that can comprehensively model such self-assembly
experiments.

In this article we present the MultiMSM approach, which provides a complete
description of assembly reactions, accounting for changes in transition rates
as concentrations evolve, as well as association between intermediates. Using
two model systems, we show that that the MultiMSM algorithm accurately
predicts the depletion of monomers and the ensemble of resulting intermediate
and target assembly species, on the long timescales required for reactions to
reach equilibrium. While brute-force dynamics simulations with related models
have been limited to restricted parameter ranges, such as high subunit
concentrations Hagan and Chandler (2006); Perlmutter  _et al._ (2014);
Rapaport (2018); Baschek  _et al._ (2012); Gupta  _et al._ (2023), the
MultiMSM approach enables simulation over a broad range of experimentally
relevant parameter values. In particular, the algorithm reduces simulation
times by orders of magnitude for systems with large nucleation barriers, which
are typically required for productive assembly at experimental conditions
Zlotnick (2003); Ceres and Zlotnick (2002); Gartner  _et al._ (2020);
Schlicksup and Zlotnick (2020); Hagan and Grason (2021); Hagan and Mohajerani
(2023).

Crucially, once the MultiMSM has been constructed for one value of the total
subunit concentration, assembly dynamics can be simulated over a wide range of
lower concentrations without any additional sampling. This enables
representing a typical experiment in which assembly is performed over a range
of subunit concentrations, but with the computational cost of a single subunit
concentration. Further, the method provides a detailed analysis of assembly
mechanisms by computing quantities such as the free energy landscape,
nucleation timescales, committor probabilities and flux along different
assembly pathways, and entropy production rates. The latter quantify the
extent to which a reaction is out of equilibrium and identify sources of
kinetic trapping that impede productive assembly. These capabilities allow
analyzing data from particle-based assembly simulations in unprecedented ways.

We provide an open-source Python library Trubiano (2024a, b) which performs
all calculations required to construct MultiMSMs and the analysis described in
this work, from simulations performed with the open-source molecular dynamics
simulation package HOOMD-blue Anderson  _et al._ (2020). The library can be
readily generalized to other software packages.

Figure 1: Schematics of the subunits and their interactions for the two model
self-assembly systems. (a) (Left) The subunit is a rigid body with Attractors
at the vertices of a pentagon and a ‘Top’ and ‘Bottom’ pseudoatom above and
below the vertex plane. (Middle) Subunit-subunit interactions. Attractor-
Attractor interactions drive subunit association, with binding affinity
parameter $\epsilon_{11}$. Top-Top repulsions result in a preferred subunit-
subunit binding angle of $116.57^{\circ}$, making a dodecahedral capsid the
ground state. (Right) A snapshot of an assembled capsid from a simulation. (b)
(Left) The subunit is a rigid triangular body, motivated by recent DNA origami
experiments Sigl  _et al._ (2021); Wei  _et al._ (2024). The cyan pseudoatoms
enforce excluded volume, and complementary pairs of attractor beads (other
colors) on each side drive association between Side $1$ – Side $2$ and Side
$3$ – Side $3$ pairs. (Middle) These interactions, with binding affinity
parameters $\epsilon_{12}$ and $\epsilon_{33}$, respectively stabilize
pentamers and dimers. (Right) A snapshot of a complete $T=3$ capsid from a
simulation. (c) Snapshots of the simulation box during dodecahedron assembly
with $\epsilon_{11}=5.5$ at early times (left) and late times (right). (d)
Same as (c) for $T=3$ capsid assembly with $\epsilon_{12}=11$ and
$\epsilon_{33}=8$.

##  II Model Systems

We first describe the two model self-assembly systems that we use to test and
demonstrate our MultiMSM approach. Our first example, dodecahedral capsid
assembly from pentagonal subunits, is sufficiently tractable that brute-force
dynamics simulations can be performed on relatively long-time scales to
stringently test the MultiMSM results. Our second example, $T=3$ capsid
assembly from triangular subunits, is more complex and computationally
expensive, and shows that the MultiMSM method scales well for complicated
problems with a large state space.

For all models and results presented in this work, we give energies in units
of the thermal energy $k_{\text{B}}T$ and lengths and concentrations in units
of $l_{0}$ and $l_{0}^{-3}$ respectively, where $l_{0}$ is related to the
subunit size for each model (see supplement Section S1 sup ).

###  II.1 Dodecahedron Assembly

Our model subunit (Fig. 1a) is adapted from previous studies of dodecahedral
capsids, including assembly of empty capsids and assembly around RNA and
synthetic polyelectrolytes McMillan  _et al._ (2005); Fejer  _et al._ (2009);
Johnston  _et al._ (2010); Perlmutter  _et al._ (2013, 2014). The subunit is a
rigid body consisting of five attractor sites, placed at the vertices of a
regular pentagon, that have attractive interactions through a Morse potential
with well-depth $\epsilon_{11}$. Each subunit also has a ‘Top’ and ‘Bottom’
pseudoatom; Top-Top and Top-Bottom pairs on nearby subunits each have
repulsive Weeks-Chandler-Anderson (WCA) interactions Weeks  _et al._ (1971).
Top-Top interactions drive subunit-subunit binding angles consistent with a
dodecahedron ($116.57^{\circ}$), while Top-Bottom interactions suppress
subunit-subunit binding in inverted orientations McMillan  _et al._ (2005);
Perlmutter  _et al._ (2013, 2014). We perform simulations in three distinct
assembly regimes by setting $\epsilon_{11}\in\\{5.0,5.5,6.0\\}$ with total
subunit concentration $c_{0}=0.0156$, which spans from almost no assembly to
rapid assembly.

###  II.2 T=3 Capsid Assembly

Our second example model was previously developed as a simplified
representation of an experimental system in which DNA origami forms rigid
triangular subunits that assemble into $T=3$ icosahedral capsids Sigl  _et
al._ (2021); Wei  _et al._ (2024). The model builds on extensive previous
simulations of capsid assembly Rapaport  _et al._ (1999); Rapaport (2010,
2018); Nguyen  _et al._ (2007); Nguyen and Brooks (2008); Elrad and Hagan
(2010); Williamson  _et al._ (2011); Perlmutter  _et al._ (2013); Zhang and
Linse (2013, 2014); Zandi  _et al._ (2020); Baschek  _et al._ (2012);
Castelnovo  _et al._ (2014, 2013); Boettcher  _et al._ (2015); Mendoza and
Reguera (2020). The model subunit excluded volume shape is represented by
three layers of ‘excluder’ atoms arranged so that the edges have a bevel angle
consistent with an icosahedron ($156.72^{\circ}$). Each excluder interacts
with all pseudoatoms through a WCA potential. Subunit-subunit attractions in
the experimental system are driven by DNA blunt-end stacking and hybridization
of single-stranded DNA molecules on complementary subunit edges. In the
computational model, these short-ranged interactions are represented by
placing two ‘attractor’ atoms on each subunit edge, on the middle layer of
excluders. Complementary pairs of attractors interact through a Lennard-Jones
potential. To match the experimental $T=3$ system, attractors on Side 1 and
Side 2 of two interacting subunits are complementary with binding energy
(Lennard-Jones potential well-depth) $\epsilon_{12}$, and attractors on Side 3
are complementary with binding energy $\epsilon_{33}$. Attractors that are not
complementary interact through a repulsive WCA potential. Fig. 1b shows a
representation of the triangular subunit, the preferred intermediate for each
interaction, and an example of the fully assembled capsid in simulation.
Despite the simplicity of the model, Wei et al. 2024 Wei  _et al._ (2024)
found that the simulation results semi-quantitatively match experimental
observations of capsid assembly dynamics.

Simulations and experiments in Wei 2024 Wei  _et al._ (2024) showed that
sufficiently imbalanced values of $\epsilon_{12}$ and $\epsilon_{33}$ lead to
hierarchical assembly pathways, since $\epsilon_{12}$ and $\epsilon_{33}$
respectively stabilize intra-pentamer and intra-dimer interactions. Stronger
$\epsilon_{12}$ leads to pentamer-biased assembly pathways, in which subunits
rapidly form pentamers, which in turn undergo assembly into capsids; whereas
stronger $\epsilon_{33}$ leads to dimer-biased pathways, with rapid formation
of dimers and their subsequent assembly of capsids. In this work, we focus on
parameters that lead to pentamer-biased assembly pathways, $\epsilon_{12}=11$
and $\epsilon_{33}=8$, with total subunit concentration $c_{0}=1.7\times
10^{-4}$.

All dynamics simulations described in this work were performed with HOOMD-blue
Anderson  _et al._ (2020). Full simulation details can be found in supplement
Section S1.

##  III Methods

Here we describe the protocol to construct and use the MultiMSM. The procedure
is separated into four steps: selection of discrete states, data processing
and transition counting, monomer fraction discretization and transition matrix
construction, and model evaluation and prediction.

Our python libraries to construct MultiMSMs from HOOMD simulations and perform
all the calculations described below are available on Github Trubiano (2024a,
b).

###  III.1 Selection of Discrete States

We first need to define a state-space discretization, which maps system
configurations into MSM microstates. In this section we focus on
characterizing the state of an individual assemblage (cluster); we address
multiple clusters in section III.2. We have previously shown that a general
state decomposition for assembly is enabled by mapping an assemblage to an
undirected graph, with nodes and edges respectively corresponding to subunits
and ‘bonds’ (subunit-subunit interactions) Perkett and Hagan (2014); Trubiano
and Holmes-Cerfon (2021). Alternative approaches based on pairwise distances
between subunits and other structural properties have also been used Karpen
_et al._ (1993); de Groot  _et al._ (2001); Singhal  _et al._ (2004); Andrec
_et al._ (2005). However, these description can be simplified, and the size of
the state space significantly decreased, with a simplified state definition
that characterizes the number of subunits and bonds within an assemblage
Perkett and Hagan (2014); Perlmutter  _et al._ (2014); Trubiano and Hagan
(2022). We use the latter approach for both examples in this article; we
define a state as $\mathcal{S}=(\mathcal{N},\mathcal{B})$ where $\mathcal{N}$
is the number of subunits in the configuration and $\mathcal{B}$ is a count of
the number of each type of bond present in the configuration. For the examples
we consider in this work, $\mathcal{B}$ is a scalar for the dodecahedron
assembly since there is only one type of bond, but is a vector with two
components for the $T=3$ capsid assembly, since there are two interaction
types (the Side1-Side2 bond and the Side3-Side3 bond, see Fig. 1b). We find
that these coordinates are sufficient to accurately characterize the dynamics
of both systems studied here.

Since a bond refers to a pair of sufficiently strongly interacting subunits,
it must be defined based on a threshold. In this work, we use cutoff distances
between corresponding particle types to define a bond. See supplement Section
S1.3 for details on the bond definition for each system.

###  III.2 Processing Simulation Data and Counting Transitions

We seek to model self-assembly dynamics in the canonical (NVT) ensemble. Since
there may be many clusters undergoing different stages of nucleation and
growth at the same time, we must compute the time evolution of the joint
probability distribution of all cluster types $j$. A complete state
decomposition would classify the assembly configuration of every cluster at a
given time point. However, for a large system with many clusters the number of
such states would be intractable. Therefore, we use the independent Markov
decomposition (IMD) method Hempel  _et al._ (2021), in which each cluster is
considered as a quasi-independent local subsystem. However, note that the
clusters are not strictly independent since pairs of clusters can merge or
split during assembly. In this framework, the state probability distribution
at frame $i$ is given by Hempel  _et al._ (2021)

$\displaystyle\vec{p}_{\,i}=\vec{p}^{\,1}_{\,i}\otimes\vec{p}^{\,2}_{\,i}\otimes\ldots\vec{p}^{\,{n_{\text{types}}}}_{\,i}$
(1)

where $\vec{p}^{\,j}_{\,i}$ is the probability distribution for the state of
cluster $j$ at frame $i$, $n_{\text{types}}$ is the total number of cluster-
types, and $\otimes$ is the Kronecker product Hempel  _et al._ (2021); Satake
(1975). Here, each cluster is defined according to the state decomposition
described in section III.1.

For self-assembly, it is useful to cast the cluster probabilities as
concentrations

$\displaystyle c^{j}_{i}={\mathcal{N}}^{j}{n}^{j}_{i}/V$ (2)

where ${\mathcal{N}}^{j}$ is the number of subunits in cluster-type $j$,
${n}^{j}_{i}$ is the number of such clusters at a given frame, and $V$ is the
volume. It is important that we use mass-weighted concentrations to maintain
the constant total subunit concentration

$\displaystyle c_{0}=\sum_{j}{\mathcal{N}}^{j}c^{j}_{i}\quad\forall i.$ (3)

Using number-weighted concentrations would result in a probability that is
normalized to the cluster distribution, which depends on time and interaction
parameters, whereas the mass-weighted concentrations maintain a normalization
that depends only on the control parameter $c_{0}$.

The transition matrix elements can be estimated from the ensemble of short
simulations by recording the number of each cluster type $j$ at each
simulation window, and then computing the number of transitions between all
cluster types, ranging from monomer to dimer transitions to
association/dissociation of larger intermediates and complete capsids, as a
function of lag time. Importantly, to maintain the mass-weighted cluster
distribution, the transition counts need to be weighted by the number of
subunits involved in each transition. For example, if a cluster $j$ with
${\mathcal{N}}^{j}$ subunits transitions to a cluster $l$ with
${\mathcal{N}}^{l}>{\mathcal{N}}^{j}$ subunits, then all ${\mathcal{N}}^{j}$
of those subunits undergo the transition. That is, the transition
$j\rightarrow k$ occurs ${\mathcal{N}}^{j}$ times. Intuitively, this can be
thought of as viewing transitions from the perspective of individual subunits
rather than clusters (see supplement Fig. S2.1).

Procedure for counting transitions. For each frame in the simulation, we group
the subunits into clusters, with a cluster defined as any collection of more
than one subunit that is bonded together. On the first frame, all clusters are
given an ID. On subsequent frames, we check if any newly found clusters are
derived from existing clusters, either through merging or splitting of sub-
clusters, and update any matching existing cluster with the new configuration.
In the case of splitting, we record what the parent cluster was and form a new
cluster. Any cluster that was not derived from an existing one is assigned a
new ID.

We treat monomers separately, tracking their addition and removal as a
separate time series. For example, if an $10$-mer loses two subunits, but they
do not form a dimer after dissociating, we log a $-2$ in the monomer time
series. If they do form a dimer, we log nothing in the monomer time series,
since no monomers are involved in this transition, but record a transition to
the dimer. To track monomer-to-monomer transitions, we store the IDs of all
subunits that are not bonded in frame $i$ in a list $M_{i}$. To determine the
number of monomer-to-monomer transitions after a lag time $k$, we take the
cardinality of the intersection of these lists, $\lvert M_{i}\cap
M_{i+k}\rvert$. Finally, we also save the monomer fraction,
$f_{\text{1}}(i)=\lvert M_{i}\rvert/N_{\text{tot}}$, where $N_{\text{tot}}$ is
the total number of subunits, at every frame, and augment any transitions that
occur in that frame with this value.

We then construct the transition count matrix as follows. We loop over every
cluster’s time series of configurations and extract every pair of
configurations separated by a lag time $k$, incrementing the corresponding
entry of the count matrix by $1$. Then, to ensure mass-weighting, we multiply
that transition count by $\min({\mathcal{N}}_{t},{\mathcal{N}}_{t+k})$.

###  III.3 Monomer Fraction Discretization and Construction of MSMs

In this section we describe how to construct the component MSMs at the
different monomer concentrations which arise as subunits are depleted during
an assembly reaction. The component MSMs are then combined into the MultiMSM.

The first step is to discretize the monomer fraction in the interval $[0,1]$,
where the monomer fraction is defined as $f_{\text{1}}=c^{1}/c_{0}$ with
$c^{1}$ the concentration of monomers (free subunits). The discretization
contains $N+1$ intervals, $D_{N}=(0,d_{1},d_{2},\cdots,d_{N},1)$. We construct
an independent MSM on each of these intervals, following the approach
described in Section III.2. Although we will describe a smoothing procedure in
Section III.4.1 to interpolate transition matrices between interval edges, the
accuracy of the MultiMSM depends strongly on the choice of discretization.
While increasing the resolution leads to higher accuracy, it also increases
the amount of total sampling required for convergence of the MSMs. Thus, the
size of each interval must be chosen such that a sufficiently large number of
relevant transitions are sampled. Sampling efficiency can be improved using
adaptive sampling Bowman  _et al._ (2010); Doerr and De Fabritiis (2014);
Voelz  _et al._ (2014); Zimmerman and Bowman (2015); Shamsi  _et al._ (2017,
2018); Zimmerman  _et al._ (2018); Wan and Voelz (2020).

We have employed a heuristic, but adaptive discretization procedure, in which
intervals are defined to contain monomer fraction values that give rise to
similar dynamics. For example, in most cases monomer fractions between $0.95$
and $1$ correspond to the initial stage of a reaction when most transitions
correspond to monomer-dimer association or dissociation. In contrast, when the
monomer fraction approaches its infinite-time limit (e.g. the equilibrium
monomer concentration for reversible assembly), most transitions will involve
large intermediates. Ensuring that each interval separates different types of
dynamics improves the accuracy of the MultiMSM. We describe a systematic
method for refining the monomer fraction discretization in supplement Sections
S4 and S7.1.

We compute the MSM in each interval as follows. We initialize a sparse count
matrix to store the number of observed transitions between each pair of states
on each of these intervals. Using the output of our cluster analysis at a
given lag time (see section III.2), we identify the monomer fraction at which
each transition occurred and increment the count for the corresponding
transition and monomer fraction interval. Once all transitions have been
logged, we compute the probability transition matrix for each interval by
normalizing the rows of the count matrices to sum to $1$.

###  III.4 Solving the Forward Kolmogorov Equation

The forward equation for the MultiMSM is a straightforward extension of the
forward equation for a standard MSM, which is given by

$\vec{p}^{\,n+1}=\vec{p}^{\,n}{\bm{P}},$ (4)

with ${\bm{P}}$ as the transition matrix and $\vec{p}^{\,n}$ as a row vector
giving the probability distribution over discrete states at timepoint $n$.

For the MultiMSM, the equation has the same form, but involves a collection of
transition matrices. At each timepoint $n$ we use the transition matrix
corresponding to the current monomer fraction. By construction, we store the
monomer fraction in the first component of the probability distribution,
$p_{0}^{n}$. Formally, we can write

$\vec{p}^{\,n+1}=\vec{p}^{\,n}{\bm{P}_{m^{n}}},\quad
m^{n}=\text{index}(p_{0}^{n}),$ (5)

where the index operator converts a value in $[0,1]$ to its corresponding
interval in the discretization.

Note that some of the methods typically used to efficiently solve the forward
equation, such as computing the spectral decomposition or pre-computing large
powers of the transition matrix, cannot be directly used here since we do not
know a priori at what timepoints the monomer fraction will cross the
discretization boundaries and change the transition matrix. The most
straightforward approach is to solve Eq. (5) iteratively via vector-matrix
multiplication for each timepoint. However, it is also possible to pre-assume
the monomer concentration as a function of time, and then iteratively apply an
efficient approach such as spectral decomposition, in which the computed
monomer concentration is updated at each iteration.

####  III.4.1 Smoothing Solutions

Solving Eq. (5) following the above approach typically leads to solutions that
are well behaved within each interval of the discretization, but have ‘jumps’
(abrupt changes in slope) at time steps when a discretization boundary is
crossed. These jumps reflect the abrupt change in the transition dynamics due
to changing the transition matrix. The jumps introduce small errors in the
solution that can accumulate and reduce the accuracy of the MultiMSM
prediction at long times. To solve this problem, we describe a smoothing
procedure to continuously interpolate between the two transition matrices
across a discretization boundary.

Consider the two intervals closest to $1$, $I_{1}=(d_{N-1},d_{N})$ and
$I_{2}=(d_{N},1)$. Let $L_{1}=d_{N}-d_{N-1}$ and $L_{2}=1-d_{N}$ be the length
of each interval and let ${\bm{P}}_{1}$ and ${\bm{P}}_{2}$ be the transition
matrices on intervals 1 and 2, respectively. Since neighboring intervals can
have significantly different lengths, we define a smoothing region that is
agnostic of absolute interval sizes. Let $\chi$ be the fraction of each
interval that is used to smoothly interpolate between them. The transition
region will then begin at $a=d_{N}-\chi L_{1}$ and end at $b=d_{N}+\chi
L_{2}$. If the monomer fraction falls within $[a,b]$, we construct a linear
combination of each interval’s transition matrix to use at that value. We
choose the weights proportional to where in the region the monomer fraction
falls, with an even split if the monomer fraction is precisely $d_{N}$. In
general, we compute

$\alpha(f_{\text{1}})=\begin{cases}\frac{1}{2}\frac{f_{\text{1}}-a}{d_{N}-a}&\text{if
}a\leq f_{\text{1}}\leq d_{N},\\\
\frac{1}{2}+\frac{1}{2}\frac{f_{\text{1}}-d_{N}}{b-d_{N}}&\text{if }d_{N}\leq
f_{\text{1}}\leq b,\end{cases}$ (6)

where $f_{\text{1}}$ is the current monomer fraction. We then construct the
final transition matrix as

${\bm{P}}_{m}=(1-\alpha(f_{\text{1}})){\bm{P}}_{1}+\alpha(f_{\text{1}}){\bm{P}}_{2}.$
(7)

Since each individual transition matrix is normalized and their coefficients
sum to $1$, ${\bm{P}}_{m}$ is also normalized and thus a valid transition
matrix.

This smoothing procedure works remarkably well for $\chi\in[0.2,0.3]$. We
include $\chi$ as an optional parameter to our solvers, with a default value
of $0.25$. Setting a value of $0$ turns off all smoothing and solves Eq. (5)
as stated.

##  IV Results and testing of the MultiMSM

We constructed a MultiMSM using simulation data for each of the model systems
described in Section II. For each example, we solved Eq. (5) to predict the
time-dependent yields of each discrete state. See Section S6 in the supplement
for a detailed description of how many trajectories were used to build each
MSM, adaptive sampling strategies, and the values used for all parameters to
the models, such as the monomer fraction discretization.

###  IV.1 Dodecahedron Capsids

Figure 2: MultiMSM predictions for dodecahedron capsid assembly dynamics at
three values of the binding affinity parameter $\epsilon_{11}$, compared
against results from brute-force dynamics on accessible timescales. Note the
increasing timescale on the x-axis as binding affinity decreases. (a) Strong
affinity $\epsilon_{11}=6$. The blue, green, and purple curves denote the
MultiMSM predictions for mass-fraction of monomer, dodecahedron, and all size
12 structures (capsids and ‘danglers’), respectively, with representative
structures from simulations labeling each curve. The noisy curves (orange,
red, brown) show the same mass-fractions estimated from $50$ independent
brute-force dynamics trajectories. (b) Moderate affinity $\epsilon_{11}=5.5$.
MultiMSM predictions for the monomer (blue) and dodecahedron (green) mass
fractions; danglers do not form at this binding affinity. The noisy curves
(orange, red) show the same mass-fractions estimated from $20$ independent
brute-force dynamics trajectories. (c) Weaker affinity $\epsilon_{11}=5.0$.
MultiMSM predictions are shown for monomers and capsids. The dashed line is
$1-f_{\text{1}}$ with $f_{\text{1}}$ the monomer fraction. The total subunit
concentration for (a)-(c) is $c_{0}=0.0156$. Error bars are estimated for the
MultiMSM by bootstrapping with $1000$ resamplings (see Section IV.3). In this
work, all energies are given in units of the thermal energy $k_{\text{B}}T$
and all length scales in units of $l_{0}$ (see text).  Figure 3: Comparison of
MultiMSM intermediate size distributions as a function of time for
dodecahedron assembly with (a) $\epsilon_{11}=6.0$ and (b)
$\epsilon_{11}=5.0$. Times and yields are plotted on a log-log scale, over the
same time intervals for the corresponding plot in Fig. 2. (c), (d) Snapshots
of representative configurations at times corresponding to the red dashed
lines in (a), (b).

Fig. 2 shows example results of assembly dynamics predicted by the MultiMSM
for the dodecahedron system at several values of the binding energy
$\epsilon_{11}$, compared against brute-force dynamics. Fig. 2a shows results
for the strongest interactions, $\epsilon_{11}=6$, for which assembly is rapid
and thus the MultiMSM results can be directly compared against brute-force
dynamics simulations on all relevant timescales. Note that for this relatively
strong binding energy, it is common for the $12$-th subunit to bind to an
$11$-mer in the wrong orientation, and then become trapped for long times. We
refer to this off-target, metastable configuration as a ‘dangler’. The smooth
lines (blue, green, purple) in Fig. 2a show the MultiMSM prediction of the
mass-weighted yields of the monomer, dodecahedra, and the size-$12$ structures
(dodecahedron and dangler), respectively, while the noisy lines show estimates
from brute-force dynamics. The agreement is excellent, with the largest
differences being only a few percent at intermediate times, while the short-
and long-time behaviors show even closer agreement.

At early times about half of the $12$-mers are danglers. These off-pathway
intermediates gradually anneal into the target dodecahedron structure.

Figs. 2b,c show the MultiMSM predictions for lower binding energy values
$\epsilon_{11}=5.5$ and $5.0$ respectively. Notably, the weaker subunit-
subunit attractions result in significantly longer assembly timescales (about
$20\times$ and $500\times$ respectively). Note that the dangler intermediates
do not occur for these weaker binding energies, and thus we focus on the most
common structures (monomers and dodecahedrons). Fig. 2b shows that the
MultiMSM predictions closely match the brute-force dynamics predictions, even
on the long timescales required for this system to approach equilibrium. We
provide further evidence that this model is accurate over longer timescales in
Section V.1 and Fig. 7b. In Fig. 2c we cannot directly compare the MultiMSM
predictions against brute-force dynamics across all timescales, as the
simulations would take about $2$ GPU-months per trajectory. We do perform a
comparison over accessible timescales ($T_{\text{f}}=2.5\times 10^{5}t_{0}$)
and see good agreement (see supplement Fig. S8.4).

To further test the accuracy of the MultiMSM prediction, we note that the
MultiMSM predictions (see Fig. 3), brute-force dynamics simulations, and
previous modeling results (e.g. Zlotnick (1994); Endres and Zlotnick (2002);
Hagan and Chandler (2006); Moisant  _et al._ (2010)) show that intermediates
are present at extremely low concentrations for weak binding energies such as
$\epsilon_{11}=5.0$. Thus, Fig. 2c also compares $1-f_{\text{1}}$ with the
dodecahedron fraction, showing that these two results are within a few percent
for all times as expected for small intermediate concentrations. Additionally,
Fig. 3 shows the MultiMSM results for the full cluster-size distribution as a
function of time for $\epsilon_{11}=6$ and $\epsilon_{11}=5$. In the case of
stronger binding, there is a broad distribution of detectable intermediate
sizes during the rapid assembly phase. In contrast, weak binding results in
approximately two-state kinetics — only dodecahedra and monomer occur at high
concentrations, with low concentrations of dimers and trimers, and trace
amounts of other transient intermediates. The snapshots in Fig. 3c,d show
representative system configurations during the rapid assembly phase for each
case. For $\epsilon_{11}=6$ (Fig. 3c) we see dodecahedra and monomers
coexisting with intermediates of various sizes, while for $\epsilon_{11}=5$
(Fig. 3d) we observe two dodecahedra and monomers along with a few transient
dimers.

These results demonstrate a powerful aspect of MSMs that is preserved in our
MultiMSM approach; the simulations used to construct the model are all of
length $0.2\times 10^{5}t_{0}$, which is orders of magnitude smaller than the
relevant assembly timescales.

###  IV.2 T=3 Capsids

Figure 4: MultiMSM predictions of $T=3$ capsid assembly dynamics with
$(\epsilon_{12},\epsilon_{33})=(11,8)$ and $c_{0}=1.7\times 10^{-4}$, compared
against results from brute-force dynamics on accessible timescales. (a)
MultiMSM predictions of mass-fractions (smooth curves) are shown for monomers
(blue), capsids (green), and ‘near-capsids’ (purple, $56-60$ subunits), with
example snapshots labeling each curve. These predictions are compared against
mass fractions estimated from $40$ brute-force dynamics trajectories performed
up to $12.5\times 10^{5}t_{0}$ (noisy curves). Error bars are estimated for
the MultiMSM by bootstrapping with $1000$ resamplings (see Section IV.3) (b) A
more detailed comparison of MultiMSM and brute-force results, also showing
mass fractions of pentamers (pink).

Fig. 4a shows MultiMSM predictions (smooth curves) for $T=3$ capsid assembly
dynamics, with results shown for monomers, complete capsids, and ‘near-
capsids’ which include any structure with $56$ or more subunits (including
complete capsids with 60 subunits). These predictions are compared against
results from brute-force dynamics (noisy curves) on accessible timescales (up
to $12.5\times 10^{5}t_{0}$). This comparison is shown in more detail in Fig.
4b, where we also include pentamers. We observe extremely close agreement,
within the statistical error of the brute-force simulations. In particular,
the MultiMSM captures the rapid conversion of monomers into pentamers at early
times, slow monomer depletion at late times, and the tendency of assembly
pathways to become trapped in near-capsid intermediates.

These behaviors arise because these simulations are performed at binding
affinity values $\epsilon_{12}=11$ and $\epsilon_{33}=8$ for which intra-
pentamer interactions are strong and intra-dimer interactions are relatively
weak Wei  _et al._ (2024). This imbalance leads to hierarchical assembly
pathways in which many subunits first form pentamers, which in turn assemble
into nearly complete capsids. However, due to a combination of steric effects,
monomer binding, and pentamer depletion, many assembly pathways become trapped
in near-capsid structures with $56-59$ subunits. While some of these
structures are converted into complete capsids by monomer additions, they
persist at $\approx 10\%$ mass fraction even at very long times.

###  IV.3 Error and Efficiency

| D12 | D12 | D12 | T3
---|---|---|---|---
| $\epsilon_{11}=6.0$ | $\epsilon_{11}=5.5$ | $\epsilon_{11}=5.0$ | $(\epsilon_{12},\epsilon_{33})=(11,8)$
Dynamics | $0.85(6)$ | $0.75(5)$ | $0.11(7)$ | $0.21(8)$
MultiMSM | $0.826(9)$ | $0.74(2)$ | $0.13(3)$ | $0.25(2)$
Table 1: Comparing the estimated capsid yield and error at $T_{\text{f}}$ for
each of the four examples with brute-force dynamics and the MultiMSM. Yields
and errors are estimated from brute-force dynamics by sample averages (see
Fig. 2 and 4), and from the MultiMSM by bootstrapping with $1000$ resamplings.
Final simulation times for the dodecahedron examples are
$T_{\text{f}}=0.5\times 10^{5}t_{0}$ for $\epsilon_{11}=6$,
$T_{\text{f}}=10\times 10^{5}t_{0}$ for $\epsilon_{11}=5.5$,
$T_{\text{f}}=2.5\times 10^{5}t_{0}$ for $\epsilon_{11}=5$; for the $T=3$
capsid $T_{\text{f}}=12.5\times 10^{5}t_{0}$.

Error. For a standard MSM, the uncertainty in the equilibrium distribution and
other quantities can be directly propagated from uncertainty estimates in
transition matrix entries Singhal and Pande (2005); Hinrichs and Pande (2007);
Bhattacharya and Chatterjee (2019); Kozlowski and Grubmüller (2023). However,
in the MultiMSM such propagation is complicated by the unknown switching times
between the component MSMs and the smoothing procedure. Therefore, we quantify
errors using bootstrapping Efron (1979, 1981) (see supplement Section S7 for
further details).

Table 1 shows a comparison of the estimated means and standard errors of the
capsid yield for each of the examples from the MultiMSM by bootstrapping with
$1000$ resamples. These results are compared against sample averages from the
brute-force dynamics simulations. In each example the comparison is shown for
the final simulation time point $T_{\text{f}}$. We see that the MultiMSM
yields are within the statistical error of the estimates from brute-force
dynamics, and that the statistical error for the MultiMSM is consistently
smaller than that from the dynamics.

We have used the same bootstrapping approach to compute the statistical error
of the MultiMSM yield predictions as a function of time for each of the
assembly examples (Fig. 2 and 4a). Where available, the sample-averaged yields
are typically within the error bars of the MultiMSM prediction. An exception
is dodecahedron assembly with $\epsilon_{11}=6$ (Fig. 2a), for which the
computed error bars are quite small at some times and sample averages lay
outside them. This is likely because this example used the fewest sample
trajectories to build the MultiMSM and we did not perform adaptive sampling,
so the sampling with replacement step performed for the bootstrap results in
very similar models. For the examples with longer assembly timescales (Fig. 2c
and 4a), the errors generally grow in time due to accumulated error from each
monomer fraction interval of the MultiMSM. While the error for dodecahedron
assembly with $\epsilon_{11}=5$ is particularly large at long time scales, our
analysis shows that this is because we had very limited sampling at the small
monomer fraction values that occur at long timescales. Importantly though,
this error could be significantly reduced with further sampling. The error
estimates provide a guide to refining the discretization of the monomer
fractions and performing additional sampling (see supplement Sections S5 and
S7.1).

As for a traditional MSM, the accuracy of the MultiMSM at long times is
limited by the sampling of the most relevant slow transitions. For example, in
the triangles system which involves strong intra-pentamer interactions, we
sampled pentamer-to-monomer transitions only $25$ times and dissociation from
complete capsids only $\sim 500$ times. As noted above, adaptive sampling
techniques focus on such transitions to improve statistics. However, estimates
of transition matrix elements that involve such rare events can be improved
much more efficiently by incorporating non-Boltzmann techniques Trendelkamp-
Schroer and Noé (2016); Noé  _et al._ (2019); Bonati  _et al._ (2021); Falkner
_et al._ (2023).

Figure 5: Estimate of the computational speedup provided by the MultiMSM. The
plot shows the ratio of the total simulation time for constructing a MultiMSM,
$t_{\text{MSM}}$, and the total simulation time to generate a given number of
brute-force dynamics trajectories, $t_{\text{MD}}$, for the four example
systems considered in this work: dodecahedron assembly with $\epsilon_{11}=6$
(green), $\epsilon_{11}=5.5$ (orange), $\epsilon_{11}=5$ (blue), and $T=3$
capsid assembly (red) with $\epsilon_{12}=11$ and $\epsilon_{33}=8$. We
compute $t_{\text{MD}}$ as the number of GPU-hours to simulate a system until
the capsid mass fraction reaches $99\%$ of its equilibrium yield (as
determined from the MultiMSM results), multiplied by the number of independent
trajectories (which determines the statistical relevance of the results). Note
that these estimates are approximately independent of system size, since both
the statistical relevance and computational costs per trajectory are
approximately linear in system size. The stars denote the point on each curve
for which the sampling error over brute-force trajectories is expected to be
comparable to the bootstrapping errors in Table 1. The dashed line denotes an
efficiency ratio of one, above which the MultiMSM is more efficient.

Efficiency. [This is confusing and I think underestimates the benefit of the
MSM because you are comparing different accuracy levels, as you allude to in
the next paragraph. Why not instead just compute how long it takes with each
method to simulate 99% of equilibrium yield at same level of accuracy?] [I
prefer having it in terms of number of trajectories so it is possible to
estimate the speedup for a given error tolerance, albeit indirectly as
$\approx\sigma_{0}/\sqrt{N}$. I added stars to the points denoting the same
level of accuracy as the bootstrap, which hopefully makes it less confusing.]
Next we compare the efficiency of constructing a MultiMSM compared to
performing an ensemble of straight-forward dynamics simulations, with the
estimate of the equilibrium capsid yield as a benchmark. Figure 6 shows the
efficiency ratio, comparing the number of GPU-hours needed to simulate $N$
brute-force dynamics trajectories until the average capsid mass fraction
reaches $99\%$ of its equilibrium value (with corresponding error $\sim
N^{-1/2}$), to the number of GPU-hours needed to build the final converged
MultiMSM with the accuracy indicated in Table 1). The time required to perform
the bootstrapping with the MultiMSM is negligible compared to the simulation
time, so we exclude this time from the calculation. For most of our examples
the efficiency ratio crosses $1$ for a few trajectories ($1<N\leq 10$), except
the dodecahedron example with rapid assembly, $\epsilon_{11}=6$ (yellow curve)
for which the MultiMSM becomes more efficient for $N\gtrsim 50$. The
difference arises from the timescale required to reach steady state — longer
relaxation timescales provide more opportunity for speedup with the MultiMSM.

We can also use this plot to determine error-weighted comparisons of the
efficiency. For example, the error bars from sample averaging are at least
double the error bars obtained from bootstrapping. If we wanted to improve the
sample averaged estimates to a comparable error, it would take roughly four
times the number of trajectories. The stars on each curve of Figure 6 denote
the point at which we expect these errors to be comparable. In all examples
but the rapidly assembling dodecahedron, we see speedups of $1$-$2$ orders of
magnitude over brute-force simulation to reach this level of accuracy.

Figure 6: Estimate of the computational speedup provided by the MultiMSM. The
plot shows the ratio of the total simulation time for constructing a MultiMSM,
$t_{\text{MSM}}$, and the total simulation time to estimate the equilibrium
capsid yield to a given error tolerance by sample averaging across brute-force
dynamics trajectories, $t_{\text{MD}}$, for the four example systems
considered in this work: dodecahedron assembly with $\epsilon_{11}=6$ (blue),
$\epsilon_{11}=5.5$ (orange), $\epsilon_{11}=5$ (green), and $T=3$ capsid
assembly (red) with $\epsilon_{12}=11$ and $\epsilon_{33}=8$. We compute
$t_{\text{MD}}$ as the number of GPU-hours to simulate a system until the
capsid mass fraction reaches $99\%$ of its equilibrium yield (as determined
from the MultiMSM results), multiplied by the expected number of independent
trajectories, $N$, needed to reach a given sample error tolerance (assumed to
scale as $\sigma_{0}/\sqrt{N}$, where $\sigma_{0}$ is fit for each example
using data in Table 1). The dashed line denotes an efficiency ratio of one,
above which the MultiMSM is more efficient.

Next we compare the efficiency of constructing a MultiMSM compared to
performing an ensemble of straight-forward dynamics simulations, with the
estimate of the equilibrium capsid yield as a benchmark. Fig. 6 shows the
efficiency as a function of sampling error for each example system, where we
have defined the efficiency as the ratio of computational time for brute-force
dynamics required to simulate dynamics until the average capsid mass fraction
reaches $99\%$ of its equilibrium value with the specified error, relative to
the time required to build a converged MultiMSM with the same accuracy. The
time required to perform the bootstrapping with the MultiMSM is negligible
compared to the simulation time, so we exclude this time from the calculation.
We obtain similar results for other metrics of comparison between the brute-
force dynamics and MultiMSM. The horizontal dashed line shows the threshold
above which the MultiMSM is more efficient.

We see that for all systems the MultiMSM is significantly more efficient than
brute-force dynamics for acceptable error levels, with better than $10^{4}$
speedup to achieve $1\%$ error for the dodecahedron system with
$\epsilon_{11}=5$. The speedup is significantly lower for the $T=3$ system —
this is because, to enable direct comparison between the MultiMSM predictions
and brute-force dynamics, we have run that system for a strong binding
affinity with correspondingly small nucleation barriers and rapid assembly. As
with any MSM algorithm, the MultiMSM speedup depends on the separation of
timescales. For self-assembly, the speedup will increase exponentially with
the height of nucleation barriers; i.e., with decreasing binding affinity or
subunit concentration.

##  V Applications

###  V.1 Concentration Sweeps

A powerful consequence of how the MultiMSM is constructed is that total
subunit concentrations $c_{0}$ below that used in the simulations to build the
MultiMSM can be analyzed with no additional sampling. Consider a
discretization of the monomer fraction, $D=[0,d_{1},...,d_{N},1]$, such that a
monomer fraction of $1$ corresponds to a maximum total subunit concentration
of $c_{0}^{\text{max}}$. For our primary analysis, we initialize the system
using the MSM defined on the interval $[d_{N},1]$ with a starting distribution
of all monomers, $p_{i}^{0}=\delta_{i0}$, and the solution to Eq. (5) gives
the system dynamics with concentration $c_{0}^{\text{max}}$.

Figure 7: Using the MultiMSM to estimate how assembly dynamics depends on
total subunit concentration, $c_{0}$, without requiring additional
simulations. (a) The MultiMSM prediction for the capsid mass fraction is shown
as a function of time and $c_{0}$ for the dodecahedron assembly with
$\epsilon_{11}=5.5$. Representative snapshots of the simulation box are shown
for four choices of time and total subunit concentration. (b) Comparison of
the MultiMSM predictions of monomer (dashed lines) and capsid (solid lines)
mass fractions against brute-force dynamics results (symbols) for total
subunit concentrations $0.7c_{0}^{\text{max}}$ (red), $0.35c_{0}^{\text{max}}$
(green), and $0.25c_{0}^{\text{max}}$ (blue). For the brute-force dynamics
simulations, the total number of subunits was fixed at $N_{\text{tot}}=125$,
the box size was increased to $L_{0}=22.525l_{0}$, $28.380l_{0}$, and
$31.750l_{0}$ respectively, and $20$ independent trajectories were performed
for each system. Results from the brute-force dynamics simulations at lower
concentrations were not used in construction of the MultiMSM.

We can then approximate the dynamics for lower subunit concentrations
$c_{0}=d_{k}c_{0}^{\text{max}}$ without any additional sampling, by
initializing the MSM defined on one of the inner intervals, $[d_{k-1,}d_{k}]$
for $k<N$, with the same initial distribution of all monomers. We define a new
monomer fraction discretization
$D_{k}=[0,d_{1}/d_{k},d_{2}/d_{k},...,d_{k-1}/d_{k},1]$, and assign
corresponding transition matrices to the same intervals. For example, if
${\bm{P}}_{1}$ originally was the transition matrix on the interval
$[0,d_{1}]$, in the reduced system it will be the transition matrix on the
interval $[0,d_{1}/d_{k}]$. A limitation is that we are constrained to the
finite set of concentrations corresponding to the bins of our discretization,
$d_{k}c_{0}$, and depending on the quality and amount of sampling initially
performed in the bins closer to zero, the error may increase as smaller
concentrations are probed. However, if needed, additional sampling can be
performed to refine the discretization based on error analysis as described in
Section IV.3 and supplement Section S7.1.

Fig. 7a shows the results of a concentration sweep performed using this method
for the dodecahedron system with $\epsilon_{11}=5.5$. The capsid yield curves
are computed for the initial concentration $c_{0}^{\text{max}}$ and smaller
values corresponding to the monomer fraction discretization bin cutoffs, and
then interpolated between these discrete values. The results show that the
final capsid yields decrease while assembly timescales increase as the total
subunit concentration is reduced, consistent with previous theory and
experiment Zlotnick (1994); Zlotnick  _et al._ (1999); Endres and Zlotnick
(2002); Zlotnick (2003); Hagan and Elrad (2010); Mizrahi  _et al._ (2022);
Zlotnick and Stray (2003); Mateu (2013); Perlmutter and Hagan (2015); Hagan
(2014); Zandi  _et al._ (2020); Prevelige  _et al._ (1993); Schwartz  _et al._
(1998); Sweeney  _et al._ (2008); Zhang and Schwartz (2006b); Tan  _et al._
(2021); Yu  _et al._ (2021); Panahandeh  _et al._ (2020); Rapaport (2008);
Whitelam  _et al._ (2009); Nguyen  _et al._ (2007); Wilber  _et al._ (2007,
2009); Hagan  _et al._ (2011); Cheng  _et al._ (2012); Rapaport (2004); Nguyen
_et al._ (2009); Johnston  _et al._ (2010); Rapaport (2010); Nguyen and Brooks
(2008); Bruinsma and Klug (2015); Hagan and Zandi (2016); Rapaport (2012);
Rapaport  _et al._ (1999); Mendoza and Reguera (2020); Wagner and Zandi
(2015); Castelnovo  _et al._ (2014, 2013); Boettcher  _et al._ (2015); van der
Schoot and Zandi (2007). We can also infer the critical assembly
concentration, as the MultiMSM initialized at concentration
$0.25c_{0}^{\text{max}}$ predicts a capsid yield of zero.

This procedure provides an excellent approximation to the dynamics when there
is a large separation of timescales between nucleation and growth timescales,
so that monomers are depleted slowly in comparison to the timescale for
transitions among larger intermediates. Such a separation of timescales is
consistent with the usual criteria for MSMs to provide effective computational
speed up. The dodecahedron assembly examples with lower binding energies
($\epsilon_{11}=5.0,5.5$) are good examples of this scenario. Cluster
nucleation is a rare event, and entire transition pathways from monomer to
capsid are frequently sampled in each of the monomer fraction discretization
bins. Consequently the approximation is highly accurate in these cases, as
shown in Fig. 7b, which compares the MultiMSM predictions made by
extrapolating to lower total concentrations against brute-force dynamics
trajectories. We performed simulations until $10\times 10^{5}t_{0}$, when the
capsid yield reaches $\approx 99\%$ of its equilibrium value according to the
MultiMSM for $c_{0}=0.7c_{0}^{\text{max}}=0.011$ and
$c_{0}=0.35c_{0}^{\text{max}}=0.0055$. For
$c_{0}=0.25c_{0}^{\text{max}}=0.004$, the MultiMSM predicts no assembly, and
we observe only a single capsid out of a possible $200$ in the brute-force
dynamics. These comparisons provide further evidence that the transition
matrix estimation is accurate in the lower monomer fraction bins of the
original MultiMSM. This reinforces the accuracy of the long-time predictions
at the original concentration $c_{0}^{\text{max}}$ in Fig. 2b, which depends
on the accuracy of the transition matrix on each interval.

For the example in Fig. 2c, approaching equilibrium with brute-force
simulations is computationally intractable, and thus we cannot directly test
the MultiMSM predictions. Importantly though, performing this concentration
sweep out to such timescales using the MultiMSM takes only a few minutes on a
CPU, demonstrating many orders of magnitude of computational speed-up.

Figure 8:  The $T=3$ capsid assembly example shows that the concentration
sweep approximation can become inaccurate for insufficient separation of
timescales. Comparison of the MultiMSM predictions of monomer (dashed lines)
and capsid (solid lines) mass fractions against brute-force dynamics results
(symbols) for total subunit concentrations $c_{0}^{\text{max}}$ (blue),
$0.63c_{0}^{\text{max}}$ (green), and $0.1c_{0}^{\text{max}}$ (purple). For
the brute-force dynamics simulations, the total number of subunits was fixed
at $N_{\text{tot}}=600$, the box size was set to $L_{0}=153l_{0}$,
$179.4l_{0}$, and $331.4l_{0}$ respectively, and $20$ independent trajectories
were performed for each system. Results from the brute-force dynamics
simulations at lower concentrations were not used in construction of the
MultiMSM.

The approximate concentration results are less accurate in cases with poor
separation of timescales. The triangle system at the simulated binding
affinity provides an example of this scenario, in which monomers deplete
quickly, but larger intermediates form over a longer timescale. Fig. 8 shows
the result of performing the concentration sweep for the $T=3$ capsid
assembly, comparing to brute-force dynamics simulations over the
computationally accessible timescale of $T_{\text{f}}=12.5\times 10^{5}t_{0}$.
The blue curves used the initial total subunit concentration,
$c_{0}=c_{0}^{\text{max}}$, as a reference, while the green and purple curves
used concentrations $c_{0}=0.63c_{0}^{\text{max}}=1.0\times 10^{-4}$ and
$c_{0}=0.1c_{0}^{\text{max}}=1.7\times 10^{-5}$, respectively. We see that for
the intermediate concentration (green), $c_{0}=1.0\times 10^{-4}$, the
MultiMSM predictions match well with the brute-force estimates except for a
slight overestimate of capsid yields at early times. Monomer depletion is
still accurately characterized by the MultiMSM in this case. However, for
lower concentrations such as $c_{0}=1.7\times 10^{-5}$ (purple), the method
breaks down; the MultiMSM predicts too rapid monomer depletion and a nonzero
capsid yield, even though the largest intermediate size is $10$ in brute-force
simulations on this timescale. This breakdown can be attributed to the
relatively high populations of intermediates; once the monomer fraction
becomes small enough ($<10\%$) in the initial system, most of the sampled
transitions are between monomers and larger intermediates. When we then use
these transition matrices to predict the dynamics of a system with all
monomers in the approximate concentration sweep, the prediction overestimates
the rate of monomers forming larger intermediates. This breakdown of the
concentration sweep approximation can be identified by comparing the yield of
intermediates to the total subunit concentration. When the triangle monomer
fraction reaches $0.63$, intermediates with sizes between $5$ and $55$
subunits account for $\sim 10\%$ of the yield, and the concentration sweep
remains accurate when initialized at this total subunit concentration. When
the triangle monomer fraction reaches $0.1$, the same set of intermediates
account for over $55\%$ of the yield and the concentration sweep approximation
breaks down. In our examples, an intermediate yield of about $20\%$ seems to
be an upper bound for getting reasonable results from the MultiMSM
concentration sweep.

###  V.2 Transition Path Theory

Transition path theory (TPT) is a powerful tool to study the statistics of
reactive trajectories, in which the system transitions between $2$
(meta)stable states Vanden-Eijnden (2006); Metzner  _et al._ (2009); Noé  _et
al._ (2009); E and Vanden-Eijnden (2010). For self-assembly, TPT provides
detailed information about kinetics, assembly pathways, and related
quantities. TPT is well established for equilibrium systems with stationary
dynamics, and has been previously applied to self-assembly systems Perkett and
Hagan (2014). However, the MultiMSM is inherently non-stationary and non-
equilibrium, so we must extend the TPT framework for self-assembly to this
context.

We take the approach of Helfmann to generalize TPT to non-stationary, finite-
time systems Helfmann  _et al._ (2020). The main quantity of interest in
transition path theory is referred to as the committor, $\vec{q}$. Given an
initial state (free monomers in our case) and a final state (e.g. a complete
capsid in our system), the _forward_ committor element $q_{i}^{+}(n)$ gives
the probability of reaching the target state before returning to the initial
state, given the system begins in an intermediate state, $i$, at time $n$. If
we let $A$ denote the set of initial states and $B$ denote the set of target
states, the forward committor can be computed as the solution of the iterative
equations:

$\displaystyle q_{i}^{+}(n)$
$\displaystyle=\sum_{j}P_{ij}(n)q_{j}^{+}(n+1),\quad i\notin A\cup B,$ (8)
$\displaystyle q_{i}^{+}(n)$ $\displaystyle=0,\quad i\in A,$ (9)
$\displaystyle q_{i}^{+}(n)$ $\displaystyle=1,\quad i\in B,$ (10)

with a final condition $q_{i}^{+}(N)=1_{B}(i)$ where $\vec{1}_{B}$ is the
indicator vector for the set $B$, and $P_{ij}(n)$ is the $(i,j)$ component of
the transition matrix at time $n$. This equation can be solved for all time by
starting at the final condition and backpropagating the solution to earlier
times via the transition matrix, while explicitly enforcing the boundary
conditions at every time step. Note that this calculation requires knowing the
transition matrix at every time step, so Eq. (5) must be solved first (see
supplement Section S3).

In addition to the forward committor, we can also compute a _backward_
committor, $q_{i}^{-}(n)$, which gives the probability that a trajectory being
in state $i$ at time $n$ came most recently from set $A$ rather than $B$. This
quantity is governed by the transition matrix of the time-reversed process,
which is given by

$P_{ij}^{-}(n)=\frac{p_{j}^{n-1}}{p_{i}^{n}}P_{ji}(n-1),$ (11)

where $p_{i}^{n}$ is the probability of being in state $i$ at time $n$
according to Eq. (5). Note that for a stationary, equilibrium process, this is
simply a statement of detailed balance. The backward committor can be computed
from a similar set of iterative equations:

$\displaystyle q_{i}^{-}(n)$
$\displaystyle=\sum_{j}P_{ij}^{-}(n)q_{j}^{-}(n-1),\quad i\notin A\cup B,$
(12) $\displaystyle q_{i}^{-}(n)$ $\displaystyle=1,\quad i\in A,$ (13)
$\displaystyle q_{i}^{-}(n)$ $\displaystyle=0,\quad i\in B,$ (14)

with initial condition $q_{i}^{-}(0)=1_{A}(i)$ with $\vec{1}_{A}$ as the
indicator vector for the set $A$. These equations can be solved simply by
forward propagation from the initial condition, which can be done concurrently
while solving Eq. (5) for $\vec{p}^{\,n}$. After computing both committors,
the current or flux of reactive trajectories can be defined as

$f_{ij}^{AB}(n)=q_{i}^{-}(n)p_{i}^{n}P_{ij}(n)q_{j}^{+}(n+1),$ (15)

with effective current
$f_{ij}^{+}(n)=\max\\{f_{ij}^{AB}(n)-f_{ji}^{AB}(n),0\\}$.

These TPT calculations are well-suited to the MultiMSM framework; Eqs. (5) and
(12) can be solved in the forward pass of the model, and then Eqs. (8) and
(15) can be solved with a backward pass. However, we must consider mass-
weighted transitions in the transition matrix, and focus the committer
calculation on the largest cluster involved in the transition. For example,
when an $n$-mer loses a subunit, we log a transition to both the $(n-1)$-mer
and the monomer. The committor would be greatly underestimated if all
intermediate states were able to transition directly back to the initial
state. We describe how to account for these complications in supplement
Section S8.1.

Figure 9: Transition path theory (TPT) on the MultiMSM for dodecahedron
capsids with $\epsilon_{11}=5.5$ and $c_{0}=0.0156$. The initial state is
entirely monomers and the target state is a complete capsid with
$(\mathcal{N},\mathcal{B})=(12,30)$. (a) The forward committor probability
$q_{i}$ as a function of intermediate size, $i$, for indicated values of the
monomer fraction $f_{\text{1}}$. Committor values are taken as a maximum over
all configurations of a given size at the first time step for which the
monomer fraction falls below the reported value. A dashed line is plotted at
probability $0.5$; the intersection points with the committor curves denote
the critical nucleus sizes. (b) The forward committor probabilities in the
full $(\mathcal{N},\mathcal{B})$ state space at the first time the monomer
fraction falls below $f_{\text{1}}=0.5$. The red, dotted line denotes the most
likely reactive path from the initial to target state.

Fig. 9 shows the results of performing these TPT calculations for the
dodecahedron system with $\epsilon_{11}=5.5$. In this example, the initial
state is entirely monomers, $(\mathcal{N},\mathcal{B})=(1,0)$, and the target
state is a complete capsid with $(\mathcal{N},\mathcal{B})=(12,30)$. Fig. 9a
shows the forward committor probability as a function of intermediate size,
with each curve corresponding to a different time point along the course of
the reaction, and thus a different free monomer fraction $f_{\text{1}}$. The
intermediate size at which the forward committor is $50\%$ (dashed line) gives
an estimate of the critical nucleus size, which thus can be identified as the
intersection between the dashed line and a given committor curve. Importantly,
we note that statistical tests would be required to establish whether the
intermediate size is a good reaction coordinate, and what, if any, additional
features would need to be added to the reaction coordinate to give a rigorous
description of the ensemble of critical nucleus structures Dellago  _et al._
(2002); Peters (2016). However, we expect the trend shown in Fig. 9a to be
reliable, and the results can be tested against the free energy calculations
shown in Section V.4. We find that as monomers are depleted, the committor
curves shift toward larger intermediate sizes and the critical nucleus size
increases from about $5$ subunits when $f_{\text{1}}=0.9$ to $7$ subunits when
$f_{\text{1}}=0.3$. This trend reflects the fact that while disassembly rates
are independent of $f_{\text{1}}$, association rates to a growing capsid are
proportional to $f_{\text{1}}$. Or, from a thermodynamic viewpoint, the
chemical potential cost of removing a free monomer from solution grows as
$k_{\text{B}}T\ln f_{\text{1}}$.

Fig. 9b shows the forward committor values for the dodecahedron system when
$f_{\text{1}}=0.5$, but in the full $(\mathcal{N},\mathcal{B})$ state space.
We also plot the most probable reaction pathway from a monomer to a capsid
(red line and points). This path is computed by determining the sequence of
states between monomer and capsid with the largest effective flux, given by
$f_{ij}^{+}$. The resulting path matches observations from previous
calculations of capsid assembly Endres and Zlotnick (2002); Hagan (2014);
Mizrahi  _et al._ (2022); the most probable assembly pathway proceeds through
intermediates with the maximum number of bonds between subunits, or those
intermediates that can easily rearrange to form new contacts and thus increase
the number of bonds without the addition of new subunits. We obtain similar
results for other monomer fraction values, except that committor probabilities
shift to slightly higher intermediate values with decreasing $f_{\text{1}}$.
The most probable path is not sensitive to the monomer fraction, although
sometimes the transitions that involve rearranging the bonds at fixed
intermediate size occur in one step rather than two.

###  V.3 Entropy Production Rates

The entropy production rate provides a means to quantitatively measure how far
from equilibrium a process is, thus elucidating the irreversibility of a
process (which produces entropy), the heat produced or work done by the
system, or the efficiency of a process Bag (2002); Vaikuntanathan and
Jarzynski (2009); Mandal  _et al._ (2017); Landi and Paternostro (2021). Much
recent work has developed optimal control algorithms for non-equilibrium
systems to minimize entropy production Gomez-Marin  _et al._ (2008); Sivak and
Crooks (2016); Large and Sivak (2019); Tafoya  _et al._ (2019); Gingrich  _et
al._ (2016); Nakazato and Ito (2021); Dechant (2022); Blaber and Sivak (2023);
Brown and Sivak (2020); Abiuso  _et al._ (2022); Louwerse and Sivak (2022);
Whitelam and Schmit (2022); Zhong and DeWeese (2022); Gupta  _et al._ (2023);
Davis  _et al._ (2024). However, entropy production is frequently difficult to
measure in experiments or simulations of complex models, due to the large
amount of data needed to reliably estimate probability distributions and
currents, although approaches based on machine learning Kim  _et al._ (2020);
Nir  _et al._ (2020); Otsubo  _et al._ (2020, 2022); Boffi and Vanden-Eijnden
(2023a, b); Whitelam (2023) and automatic-differentiation Engel  _et al._
(2023) can help. Fortunately, the MSMs enable computationally efficient
computation of entropy production.

The time-dependent entropy production rate for a Markov chain at step $n$ is
given by the expression

$\epsilon_{\text{p}}^{n}=\frac{1}{2}\sum_{i,j}\left(p_{i}^{n}P_{ij}-p_{j}^{n}P_{ji}\right)\log\frac{p_{i}^{n}P_{ij}}{p_{j}^{n}P_{ji}},$
(16)

where $p_{i}^{n}$ is the probability of being in state $i$ at time $n$, and
$P_{ij}$ is the transition probability from state $i$ to state $j$. This
extends naturally to the MultiMSM framework, with the transition matrix
replaced by the time-dependent transition matrix across the monomer fraction
bins

$\epsilon_{\text{p}}^{n}=\frac{1}{2}\sum_{i,j}\left(p_{i}^{n}P_{ij}(n)-p_{j}^{n}P_{ji}(n)\right)\log\frac{p_{i}^{n}P_{ij}(n)}{p_{j}^{n}P_{ji}(n)},$
(17)

where ${\bm{P}}(n)$ is the transition matrix used to update the system at time
step $n$. This quantity can be computed while solving the forward Kolmogorov
equation as described in Section III.4.

Figure 10: Entropy production rates as a function of time as computed by the
MultiMSM for (a) dodecahedron assembly with $\epsilon_{11}=5.5$ and
$c_{0}=0.0156$ and (b) T=3 capsid assembly with
$(\epsilon_{12},\epsilon_{33})=(11,8)$ and $c_{0}=1.7\times 10^{-4}$. Selected
peaks and other areas of interest are labeled with the transition or
transitions that contribute most to the total entropy at that time.

Figs. 10a,b show the result of computing the entropy production rate as a
function of time using the MultiMSM for dodecahedron assembly at the
intermediate binding energy, $\epsilon_{11}=5.5$, and the T=3 model. For the
dodecahedron case, we observe an initial spike in the entropy production rate
at early times, corresponding to a non-equilibrium flux of monomers into
larger intermediates. The entropy production rate then decreases rapidly and
remains small throughout the remainder of the simulation. Importantly, note
the difference in timescales for the entropy production rate decreasing toward
zero and the yields reaching a steady state in Fig. 3b. This suggests the
assembly is near-equilibrium once the monomer concentration reaches roughly
half its starting value. This is consistent with the accepted notion that
productive assembly occurs when the system is near equilibrium Zlotnick
(2007); Hagan (2014); Zandi  _et al._ (2020).

We can also track the dominant contributions to the sum in Eq. (17) to better
understand the assembly process. The two initial peaks are labeled with a
representation of the transition contributing the most to the entropy at that
time. The first peak corresponds to the transition between
$(\mathcal{N},\mathcal{B})=(7,12)$ to $(\mathcal{N},\mathcal{B})=(8,15)$,
accounting for approximately $7\%$ of the entropy production rate, while the
second peak corresponds to the transition between
$(\mathcal{N},\mathcal{B})=(9,18)$ to $(\mathcal{N},\mathcal{B})=(10,21)$,
accounting for approximately $11\%$ of the entropy production rate. These two
transitions are dominant considering that the next largest contributions
account for less than $2\%$. This observation suggests these are the dominant
nearly irreversible transitions in the assembly process, and indeed we can see
these transitions occur on the most likely assembly pathway identified in Fig.
9b. The enhanced lack of reversibility for these transitions can be understood
because both transitions add three bonds to the configuration, compared to
just one or two for most of the other early transitions along the pathway.
This result demonstrates that the entropy production can provide important
insights into a reaction pathway by identifying the key transitions that
stabilize intermediates or products.

The $T=3$ model (Fig. 10b) exhibits a similar large spike in entropy
production at early times due to rapid formation of dimers and small
intermediates, followed by a decay over time. The initial peak corresponds to
the monomer to dimer transition, accounting for approximately half the entropy
production at this time. The next peak is dominated by transitions from a
dimer, trimer, or tetramer to a pentamer, with each transition contributing
approximately $2.5\%$ for a total of $7.5\%$. These are the expected dominant
contributions, as the strong intra-pentamer interactions make pentamerization
a nearly irreversible process under these parameter values. There are no
intermediate peaks that have a dominant contribution to the total entropy
production, but rather many transitions with roughly the same contribution. An
interesting observation for this example is that the entropy production rate
does not tend to zero at large times, but rather decays to a small positive
number, approximately $0.006$. It does not decay further even if we increase
the final time by up to two orders of magnitude. This entropy production is
almost entirely ($>99\%$) due to configurations with $59$ subunits
transitioning to the T=3 capsid, as well as configurations with $60$ subunits
but the wrong bond configuration transitioning to the T=3 capsid. At this
point in the assembly, the free monomer concentration has already decayed to
less than $1\%$ of its initial value, which indicates that the system has
become trapped in metastable states. This result demonstrates that entropy
production provides a useful measure of the extent to which a system is
kinetically trapped.

###  V.4 Free Energy

Figure 11: MultiMSM free energy calculations for dodecahedron assembly. (a)
The per-subunit Helmholtz free energy $f_{n}$ computed from each component MSM
as a function of intermediate size $n$ for $\epsilon_{11}=5.5$ (solid lines),
showing that the results are independent of monomer concentration. The free
energy averaged over the component MSMs is shown as a dashed line. (b) The
averaged Helmholtz free energy profiles as a function of binding energy,
$\epsilon_{11}$. The inset compares the free energy for the dimer computed
from the MultiMSM against the estimated value and standard error obtained
directly from dynamical simulations (without using an MSM) (see supplement
Section S8.3.3). (c) The grand free energy, computed analytically from the
Helmholtz free energy using Eq. (18) and evaluated at the average monomer
fraction within each discretization bin (solid lines). Dashed lines show the
grand free energy computed directly from the MSM equilibrium distribution, for
monomer fractions and intermediate sizes that are sufficiently close to
equilibrium to enable comparison.

Another useful application of MSMs is calculating the equilibrium free energy
of each microstate, which can then be projected onto reaction coordinate(s)
for mechanistic insight. Notably, the MultiMSM allows computing the
equilibrium free energy profile independent of free monomer concentration,
from nonequilibrium simulations at arbitrary monomer concentration.

While a typical MSM enables computing the Helmholtz free energy from the
equilibrium distribution of microstates $\vec{\pi}$ as
$F_{i}=-k_{\text{B}}T\log(\pi_{i})$, in the MultiMSM framework the transition
matrix depends on the free subunit concentration and thus the corresponding
chemical potential. Therefore, the MultiMSM gives the grand free energy for
the $j$-th component MSM as

$\Omega_{i}^{j}=-k_{\text{B}}T\log(\pi_{i}^{j})=F_{i}-\mu_{j}(n_{i}-1),$ (18)

where $\Omega_{i}^{j}$ is the grand free energy for state $i$ in the $j$-th
MSM, $\mu_{j}=k_{\text{B}}T\log(c_{j}/c_{\text{ss}})$ is the chemical
potential for the $j$-th MSM with $c_{j}$ as the average monomer concentration
within the bin and $c_{\text{ss}}$ as the standard state concentration, and
$n_{i}$ is the number of subunits in microstate $i$.

Importantly, since the Helmholtz free energy $F_{i}$ in Eq. (18) is
independent of the chemical potential (and thus the free monomer
concentration), it should be the same for each component MSM of the MultiMSM.
Therefore, its statistics can be improved by averaging over each of the
component MSMs.

To simplify the following presentation, we project the free energy onto a
reaction coordinate where there is one state per intermediate size $n$, but
the approach readily generalizes to multiple states per size. To compute
$F_{i}$ within a particular MSM, we compute the equilibrium distribution from
the transition matrix and then the set of equilibrium constants for the
formation of each $n$-mer,

$K_{n}=\left[n\right]/\left[1\right]^{n},$ (19)

where $\left[n\right]$ denotes the equilibrium concentration of $n$-mers. Note
that any absorbing states (states or groups of states for which there were an
insufficient number of exit transitions to estimate an outward transition
rate) should be eliminated from the transition matrix before computing
equilibrium quantities. As noted above, outward transition rates could be
computed for such states by combining free energy calculations with the
dynamical simulations used to estimate transition rates Trendelkamp-Schroer
and Noé (2016), but we have not implemented this approach for the present
work. Then, the concentration of subunits in absorbing states should be
subtracted from the total subunit concentration to ensure proper normalization
of the grand free energy (see supplement Sections S8.3.1 and S8.3.2 for
details).

The Helmholtz free energy can be computed as

$F_{n}=-k_{\text{B}}T\log(c_{\text{ss}}^{n-1}K_{n}).$ (20)

While the equilibrium constant corresponds to the true observable and is thus
independent of standard state, the free energy values necessarily depend on
$c_{\text{ss}}$. For the results presented here, we choose $c_{\text{ss}}$ in
simulation units such that there is one subunit per circumscribed volume
occupied by the subunit: $c_{\text{ss}}=0.66/l_{0}^{3}$ for the pentagonal
subunit and $c_{\text{ss}}=0.17/l_{0}^{3}$ for the triangular subunits.

Figure 12: MultiMSM free energy calculations for $T=3$ capsids. (a) The per-
subunit Helmholtz free energy for each component MSM (solid lines) and the
averaged profile (dashed) line, showing that the results are independent of
monomer concentration. The inset compares the free energy estimates from the
MultiMSM against results from dynamical simulations up to pentamers (see
supplement Section S8.3.4). (b) The per-subunit grand free energy, computed
analytically from the Helmholtz free energy using Eq. (18) and evaluated for
several free monomer fractions.

Figs. 11a and 12a show the Helmholtz free energy per subunit $f_{n}=F_{n}/n$
as a function of intermediate size, computed from the MultiMSM using Eq. (20),
across the monomer fraction bins for dodecahedron assembly with
$\epsilon_{11}=5.5$ and $T=3$ capsid assembly, respectively. We see that the
curves from different component MSMs collapse, showing that the Helmholtz free
energy is independent of monomer concentration as it should be. We then
average this data to construct the best statistical estimate. As expected, the
barrier height decreases and the stability of each intermediate increases with
increasing binding affinity. As an independent test, we compare the dimer
Helmholtz free energy against the free energy computed directly from
simulations (without using MSMs) in the inset (see supplement Sections S8.3.3
and S8.3.4 for details). For the $T=3$ system we performed the same procedure
and compared against direct simulations for all intermediates up to a
pentamer, which also show strong agreement in Fig. 12a.

The Helmholtz free energy enables computing the grand free energy at any free
monomer concentration by adding the chemical potential term as in Eq. (18).
For example, Figs. 11c and 12b show the grand free energy for dodecahedra with
$\epsilon_{11}=5.5$ and $T=3$ capsids, respectively, at several representative
subunit concentrations. For the dodecahedron system, we evaluate the grand
free energy at monomer fractions corresponding to the average monomer fraction
within each discretization window (solid lines). This allows for a direct
comparison of our computed grand free energy with the result of computing it
directly from the component MSM equilibrium distribution (dashed lines). Note
that the nucleation barrier increases, and the stability of intermediates
decrease, as free monomer concentrations decrease due to the increased monomer
chemical potential. Note that these curves should coincide only when the
assembly is at or near equilibrium. Fig. 10a shows that the entropy production
vanishes once the monomer fraction is below roughly $0.5$, and we see
excellent agreement between the two calculations at lower monomer fractions.
For larger monomer fractions, we may expect structures smaller than the
critical nucleus to be in quasi-equilibrium Hagan (2014), so we make
comparisons up to an intermediate size of $5$. In this case, we find loose
agreement between the curves that improves as the monomer fraction decreases.

For $T=3$ capsids, we do not make comparisons with the equilibrium MSM as we
have previously determined that this example is kinetically trapped and thus
not in equilibrium. We instead evaluate this for the range of monomer
fractions that we observe at the final simulation time of our long
trajectories. Interestingly, these profiles exhibit local minima at multiples
of $5$, reflective of the pentamer biased assembly parameters and consistent
with previous computational and experimental results Wei  _et al._ (2024).

##  VI Conclusions

We have described the MultiMSM, a general framework to construct MSMs for
systems in which many clusters assemble simultaneously and concentrations (and
potentially other parameters) change over the course of the reaction. Using
two model systems, we show that MultiMSMs can accurately describe assembly
dynamics over the long timescales required to approach equilibrium, even when
constructed from trajectories that are orders of magnitude shorter. This
capability enables particle-based simulations with complex models to be
simulated at experimentally relevant concentrations. The degree of speed up
enabled by the MultiMSM increases exponentially for systems with large
nucleation barriers, and the method is well-suited for systems that assemble
by diverse pathways. Moreover, the method is trivially parallelizable and thus
highly scalable.

In addition to extending on previous work by allowing for multiple clusters
and the depletion of free subunits, the MultiMSM approach allows for a number
of further applications. Notably, a MultiMSM constructed at one total subunit
concentration can be used to perform an approximate parameter sweep over a
wide range of concentrations with no additional sampling. This capability
corresponds to orders of magnitude additional speed up in comparison to brute-
force simulations. The results are highly accurate for conditions leading to
productive assembly, under which concentrations of intermediates remain
relatively low, and qualitatively accurate for more aggressive assembly
conditions that lead to a buildup of intermediates. The MultiMSM framework
computes transition path theory quantities, such as the committor probability
describing the extent of progress along a reaction coordinate and the relative
flux along different assembly pathways. It can also be used to estimate the
Helmholtz free energy of a system, as well as compute the grand free energy as
a function of the monomer concentration. Further, the method allows for
efficient calculation of entropy production rates, a quantity which has been
difficult to compute from particle-based simulation trajectories. We find that
the entropy production rate provides a useful quantification of how far an
assembly reaction is from equilibrium and whether it is susceptible to kinetic
traps. More interestingly, by analyzing which transitions contribute most to
the entropy production rate, one can identify factors that stabilize critical
nuclei or other key intermediates, and whether these transitions lead to
productive assembly or engender kinetic traps. These insights can form the
basis for rational design of synthetic assembly systems, or, in the case of
biomedically relevant assembly systems such as viruses, help to identify
targets for antiviral molecules that interfere with assembly Kim  _et al._
(2021); Kondylis  _et al._ (2018); Schlicksup  _et al._ (2018); Pavlova  _et
al._ (2022); Wang  _et al._ (2023); Kra  _et al._ (2023).

###### Acknowledgements.

This work was supported by the NSF through DMR 2309635 and the Brandeis Center
for Bioinspired Soft Materials, an NSF MRSEC (DMR-2011846). We also
acknowledge computational support from NSF XSEDE computing resources
allocation TG-MCB090163 (NCSA Delta GPU) and the Brandeis HPCC which is
partially supported by the NSF through DMR-MRSEC 2011846 and OAC-1920147.

## Conflicts of Interest

The authors declare no conflicts of interest.

## Data Availability Statement

Our Python library to post-process the results of a HOOMD simulation Trubiano
(2024a) and our Python library to construct and perform calculations on
MultiMSMs Trubiano (2024b) is freely available on Github. Our simulation
scripts, a subset of post-processed trajectory data, analysis scripts, and
figure generating scripts are hosted on the Open Science Framework OSFHome
(https://osf.io/ak49n/).

## References

  * Caspar and Klug (1962) D. L. D. Caspar and A. Klug, Cold Spring Harbor symposia on quantitative biology 27, 1 (1962).
  * Mateu (2013) M. G. Mateu, Arch. Biochem. Biophys. 531, 65 (2013).
  * Bruinsma and Klug (2015) R. F. Bruinsma and W. S. Klug, Annu. Rev. Condens. Matter Phys. 6, 245 (2015).
  * Perlmutter and Hagan (2015) J. D. Perlmutter and M. F. Hagan, Annu. Rev. Phys. Chem. 66, 217 (2015), http://dx.doi.org/10.1146/annurev-physchem-040214-121637 .
  * Hagan and Zandi (2016) M. F. Hagan and R. Zandi, Curr. Opin. Vir. 18, 36 (2016).
  * Twarock  _et al._ (2018) R. Twarock, R. J. Bingham, E. C. Dykeman, and P. G. Stockley, Curr. Opin. Vir. 31, 74 (2018).
  * Zandi  _et al._ (2020) R. Zandi, B. Dragnea, A. Travesset, and R. Podgornik, Physics Reports 847, 1 (2020), on virus growth and form.
  * Hagan and Grason (2021) M. F. Hagan and G. M. Grason, Rev. Mod. Phys. 93, 025008 (2021).
  * Kerfeld  _et al._ (2010) C. A. Kerfeld, S. Heinhorst, and G. C. Cannon, Microbiology+ 64, 391 (2010).
  * Tanaka  _et al._ (2008) S. Tanaka, C. A. Kerfeld, M. R. Sawaya, F. Cai, S. Heinhorst, G. C. Cannon, and T. O. Yeates, Science 319, 1083 (2008), http://www.sciencemag.org/content/319/5866/1083.full.pdf .
  * Rae  _et al._ (2013) B. D. Rae, B. M. Long, M. R. Badger, and G. D. Price, Microbiol. Mol. Biol. Rev. 77, 357 (2013), http://mmbr.asm.org/content/77/3/357.full.pdf+html .
  * Bobik  _et al._ (2015) T. A. Bobik, B. P. Lehman, and T. O. Yeates, Mol. Microbiol. 98, 193 (2015).
  * Chowdhury  _et al._ (2014) C. Chowdhury, S. Sinha, S. Chun, T. O. Yeates, and T. A. Bobik, Microbiol. Mol. Biol. Rev. 78, 438 (2014).
  * Kerfeld and Melnicki (2016) C. A. Kerfeld and M. R. Melnicki, Curr. Opin. Plant Biol. 31, 66 (2016).
  * Polka  _et al._ (2016) J. K. Polka, S. G. Hays, and P. A. Silver, Cold Spring Harb Perspect Biol 8 (2016), 10.1101/cshperspect.a024018.
  * Kerfeld  _et al._ (2018) C. A. Kerfeld, C. Aussignargues, J. Zarzycki, F. Cai, and M. Sutter, Nat. Rev. Microbiol. 16, 277 (2018).
  * Slininger Lee and Tullman-Ercek (2017) M. Slininger Lee and D. Tullman-Ercek, Current Opinion in Systems Biology (2017), 10.1016/j.coisb.2017.05.017.
  * Banerjee  _et al._ (2020) S. Banerjee, M. L. Gardel, and U. S. Schwarz, Annual Review of Condensed Matter Physics 11, 421 (2020).
  * Cabeen and Jacobs-Wagner (2010) M. T. Cabeen and C. Jacobs-Wagner, Annual Review of Genetics 44, 365 (2010).
  * Pilhofer and Jensen (2013) M. Pilhofer and G. J. Jensen, Current Opinion in Cell Biology 25, 125 (2013), cell architecture.
  * Whitelam (2010) S. Whitelam, Phys. Rev. Lett. 105, 088102 (2010).
  * Garg T (2015) G. A. Garg T, Rath G, Critical Reviews in Therapeutic Drug Carrier Systems 32(2), 89 (2015).
  * Beija  _et al._ (2012) M. Beija, R. Salvayre, N. Lauth-de Viguerie, and J.-D. Marty, Trends in Biotechnology 30, 485 (2012).
  * Ebbens (2016) S. Ebbens, Current Opinion in Colloid & Interface Science 21, 14 (2016).
  * Mallory  _et al._ (2018) S. A. Mallory, C. Valeriani, and A. Cacciuto, Annual Review of Physical Chemistry 69, 59 (2018), pMID: 29106809.
  * Fan  _et al._ (2011) J. A. Fan, Y. He, K. Bao, C. Wu, J. Bao, N. B. Schade, V. N. Manoharan, G. Shvets, P. Nordlander, D. R. Liu, and F. Capasso, Nano Letters 11, 4859 (2011), pMID: 22007607.
  * Huh  _et al._ (2020) J.-H. Huh, K. Kim, E. Im, J. Lee, Y. Cho, and S. Lee, Advanced Materials 32, 2001806 (2020).
  * Ke  _et al._ (2012) Y. Ke, L. L. Ong, W. M. Shih, and P. Yin, Science 338, 1177 (2012).
  * Kraft  _et al._ (2009) D. J. Kraft, J. Groenewold, and W. K. Kegel, Soft Matter 5, 3823 (2009).
  * Sacanna  _et al._ (2010) S. Sacanna, W. T. M. Irvine, P. M. Chaikin, and D. J. Pine, Nature 464, 575 (2010).
  * Sacanna  _et al._ (2013) S. Sacanna, M. Korpics, K. Rodriguez, L. Colón-Meléndez, S.-H. Kim, D. J. Pine, and G.-R. Yi, Nature Communications 4, 1 (2013).
  * Yi  _et al._ (2013) G.-R. Yi, D. J. Pine, and S. Sacanna, Journal of Physics: Condensed Matter 25, 193101 (2013).
  * Wang  _et al._ (2014) Y. Wang, Y. Wang, X. Zheng, G.-R. Yi, S. Sacanna, D. J. Pine, and M. Weck, Journal of the American Chemical Society 136, 6866 (2014).
  * He  _et al._ (2020) M. He, J. P. Gales, É. Ducrot, Z. Gong, G.-R. Yi, S. Sacanna, and D. J. Pine, Nature 585, 524 (2020).
  * He  _et al._ (2021) M. He, J. P. Gales, X. Shen, M. J. Kim, and D. J. Pine, Langmuir (2021).
  * Chen  _et al._ (2011a) Q. Chen, S. C. Bae, and S. Granick, Nature 469, 381 (2011a).
  * Chen  _et al._ (2011b) Q. Chen, J. K. Whitmer, S. Jiang, S. C. Bae, E. Luijten, and S. Granick, Science 331, 199 (2011b).
  * Zerrouki  _et al._ (2008) D. Zerrouki, J. Baudry, D. Pine, P. Chaikin, and J. Bibette, Nature 455, 380 (2008).
  * Yan  _et al._ (2013) J. Yan, K. Chaudhary, S. C. Bae, J. A. Lewis, and S. Granick, Nature Communications 4, 1 (2013).
  * Wolters  _et al._ (2015) J. R. Wolters, G. Avvisati, F. Hagemans, T. Vissers, D. J. Kraft, M. Dijkstra, and W. K. Kegel, Soft Matter 11, 1067 (2015).
  * Tikhomirov  _et al._ (2018) G. Tikhomirov, P. Petersen, and L. Qian, Journal of the American Chemical Society 140, 17361 (2018).
  * Oh  _et al._ (2019) J. S. Oh, S. Lee, S. C. Glotzer, G.-R. Yi, and D. J. Pine, Nature Communications 10, 1 (2019).
  * Ben-Sasson  _et al._ (2021) A. J. Ben-Sasson, J. L. Watson, W. Sheffler, M. C. Johnson, A. Bittleston, L. Somasundaram, J. Decarreau, F. Jiao, J. Chen, I. Mela, _et al._ , Nature 589, 468 (2021).
  * Kahn  _et al._ (2022) J. S. Kahn, B. Minevich, A. Michelson, H. Emamy, K. Kisslinger, S. Xiang, S. K. Kumar, and O. Gang, ChemRxiv https://doi.org/10.26434/chemrxiv-2022-xwbst (2022).
  * Noé  _et al._ (2009) F. Noé, C. Schütte, E. Vanden-Eijnden, L. Reich, and T. R. Weikl, Proceedings of the National Academy of Sciences 106, 19011 (2009), https://www.pnas.org/doi/pdf/10.1073/pnas.0905466106 .
  * Bowman  _et al._ (2009) G. R. Bowman, K. A. Beauchamp, G. Boxer, and V. S. Pande, The Journal of Chemical Physics 131, 124101 (2009).
  * Pande  _et al._ (2010) V. S. Pande, K. Beauchamp, and G. R. Bowman, Methods 52, 99 (2010), protein Folding.
  * Sarich  _et al._ (2010) M. Sarich, F. Noé, and C. Schütte, Multiscale Modeling & Simulation 8, 1154 (2010).
  * Bowman  _et al._ (2011) G. R. Bowman, V. A. Voelz, and V. S. Pande, Current Opinion in Structural Biology 21, 4 (2011).
  * Prinz  _et al._ (2011) J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera, C. Schütte, and F. Noé’, The Journal of Chemical Physics 134, 174105 (2011), https://doi.org/10.1063/1.3565032 .
  * Voelz  _et al._ (2011) V. A. Voelz, M. Jäger, L. Zhu, S. Yao, O. Bakajin, S. Weiss, L. J. Lapidus, and V. S. Pande, Biophysical Journal 100, 515a (2011).
  * Voelz  _et al._ (2012) V. A. Voelz, M. JÃ¤ger, S. Yao, Y. Chen, L. Zhu, S. A. Waldauer, G. R. Bowman, M. Friedrichs, O. Bakajin, L. J. Lapidus, S. Weiss, and V. S. Pande, Journal of the American Chemical Society 134, 12565 (2012).
  * Noé and Nüske (2013) F. Noé and F. Nüske, Multiscale Modeling & Simulation 11, 635 (2013).
  * De Sancho  _et al._ (2013) D. De Sancho, J. Mittal, and R. B. Best, Journal of Chemical Theory and Computation 9, 1743 (2013).
  * Bowman  _et al._ (2014) G. R. Bowman, V. S. Pande, and F. Noé,  _An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation_ (Springer Netherlands, 2014).
  * Chodera and Noé (2014) J. D. Chodera and F. Noé, Current Opinion in Structural Biology 25, 135 (2014), theory and simulation / Macromolecular machines.
  * Malmstrom  _et al._ (2014) R. D. Malmstrom, C. T. Lee, A. T. Van Wart, and R. E. Amaro, Journal of Chemical Theory and Computation 10, 2648 (2014).
  * Schwantes  _et al._ (2014) C. R. Schwantes, R. T. McGibbon, and V. S. Pande, The Journal of Chemical Physics 141, 090901 (2014), https://doi.org/10.1063/1.4895044 .
  * Hummer and Szabo (2015) G. Hummer and A. Szabo, The Journal of Physical Chemistry B 119, 9029 (2015).
  * Husic and Pande (2018) B. E. Husic and V. S. Pande, Journal of the American Chemical Society 140, 2386 (2018), pMID: 29323881.
  * Zeng  _et al._ (2018) X. Zeng, L. Zhu, X. Zheng, M. Cecchini, and X. Huang, Phys. Chem. Chem. Phys. 20, 6767 (2018).
  * Wu and Noé (2020) H. Wu and F. Noé, Multiscale Modeling & Simulation 30, 23 (2020).
  * Weng  _et al._ (2020) J. Weng, M. Yang, W. Wang, X. Xu, and Z. Tian, Journal of the American Chemical Society 142, 21344 (2020).
  * Suárez  _et al._ (2021) E. Suárez, R. P. Wiewiora, C. Wehmeyer, F. Noé’, J. D. Chodera, and D. M. Zuckerman, Journal of Chemical Theory and Computation 17, 3119 (2021).
  * Juárez and Bevan (2012) J. J. Juárez and M. A. Bevan, Advanced Functional Materials 22, 3833 (2012), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adfm.201200400 .
  * Tang  _et al._ (2016) X. Tang, B. Rupp, Y. Yang, T. D. Edwards, M. A. Grover, and M. A. Bevan, ACS Nano 10, 6791 (2016).
  * Tang  _et al._ (2017) X. Tang, J. Zhang, M. A. Bevan, and M. A. Grover, Journal of Process Control 60, 141 (2017), dYCOPS-CAB 2016.
  * Grover  _et al._ (2019) M. A. Grover, D. J. Griffin, and X. Tang, IFAC-PapersOnLine 52, 1 (2019), 12th IFAC Symposium on Dynamics and Control of Process Systems, including Biosystems DYCOPS 2019.
  * Trubiano and Hagan (2022) A. Trubiano and M. F. Hagan, The Journal of Chemical Physics 157, 244901 (2022), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0130407/16648222/244901_1_online.pdf .
  * Jamalyaria  _et al._ (2005) F. Jamalyaria, R. Rohlfs, and R. Schwartz, J. Comput. Phys. 204, 100 (2005).
  * Keef  _et al._ (2006) T. Keef, C. Micheletti, and R. Twarock, J. Theor. Biol. 242, 713 (2006).
  * Hemberg  _et al._ (2006) M. Hemberg, S. Yaliraki, and M. Barahona, Biophys. J. 90, 3029 (2006).
  * Dykeman  _et al._ (2013) E. C. Dykeman, P. G. Stockley, and R. Twarock, Phys Rev E 87, 022717 (2013).
  * Sweeney  _et al._ (2008) B. Sweeney, T. Zhang, and R. Schwartz, Biophys. J. 94, 772 (2008).
  * Zhang and Schwartz (2006a) T. Q. Zhang and R. Schwartz, Biophys. J. 90, 57 (2006a).
  * Misra  _et al._ (2008) N. Misra, D. Lees, T. Q. Zhang, and R. Schwartz, Comput. Math. Method Med. 9, 277 (2008).
  * Kumar and Schwartz (2010) M. S. Kumar and R. Schwartz, Phys. Biol. 7, 045005 (2010).
  * Xie  _et al._ (2012) L. Xie, G. R. Smith, X. Feng, and R. Schwartz, Biophys. J. 103, 1545 (2012).
  * Smith  _et al._ (2014) G. R. Smith, L. Xie, B. Lee, and R. Schwartz, Biophys. J. 106, 310 (2014).
  * Kelley  _et al._ (2008) N. W. Kelley, V. Vishal, G. A. Krafft, and V. S. Pande, The Journal of Chemical Physics 129, 214707 (2008), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3010881/13346829/214707_1_online.pdf .
  * Perkett and Hagan (2014) M. R. Perkett and M. F. Hagan, The Journal of Chemical Physics 140, 214101 (2014).
  * Dibak  _et al._ (2018) M. Dibak, M. J. del Razo, D. De Sancho, C. Schütte, and F. Noé, The Journal of Chemical Physics 148, 214107 (2018).
  * Olsson  _et al._ (2017) S. Olsson, H. Wu, F. Paul, C. Clementi, and F. Noé’, Proceedings of the National Academy of Sciences 114, 8265 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1704803114 .
  * Sengupta  _et al._ (2019) U. Sengupta, M. Carballo-Pacheco, and B. Strodel, The Journal of Chemical Physics 150, 115101 (2019), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5083915/16723630/115101_1_online.pdf .
  * Trubiano and Holmes-Cerfon (2021) A. Trubiano and M. Holmes-Cerfon, Soft Matter 17, 6797 (2021).
  * Miskin  _et al._ (2016) M. Z. Miskin, G. Khaira, J. J. de Pablo, and H. M. Jaeger, Proceedings of the National Academy of Sciences 113, 34 (2016).
  * Pietrucci (2017) F. Pietrucci, Reviews in Physics 2, 32 (2017).
  * Sidky and Whitmer (2018) H. Sidky and J. K. Whitmer, The Journal of Chemical Physics 148, 104111 (2018), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5018708/15543653/104111_1_online.pdf .
  * Sherman  _et al._ (2020) Z. M. Sherman, M. P. Howard, B. A. Lindquist, R. B. Jadrich, and T. M. Truskett, The Journal of Chemical Physics 152, 140902 (2020), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5145177/13478051/140902_1_online.pdf .
  * Das and Limmer (2021) A. Das and D. T. Limmer, The Journal of Chemical Physics , 1 (2021).
  * Hénin  _et al._ (2022) J. Hénin, T. Lelieèvre, M. R. Shirts, O. Valsson, and L. Delemotte, Living Journal of Computational Molecular Science 4 (2022), 10.33011/livecoms.4.1.1583.
  * Wang  _et al._ (2022) D. Wang, Y. Wang, J. Chang, L. Zhang, H. Wang, and W. E., Nature Computational Science 2, 20 (2022).
  * Das and Limmer (2023) A. Das and D. T. Limmer, Proceedings of the National Academy of Sciences 120, e2217242120 (2023), https://www.pnas.org/doi/pdf/10.1073/pnas.2217242120 .
  * Goodrich  _et al._ (2021) C. P. Goodrich, E. M. King, S. S. Schoenholz, E. D. Cubuk, and M. P. Brenner, Proceedings of the National Academy of Sciences 118, e2024083118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2024083118 .
  * Lieu and Yoshinaga (2022) U. T. Lieu and N. Yoshinaga, The Journal of Chemical Physics 156, 054901 (2022), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0072234/16770292/054901_1_online.pdf .
  * Jhaveri  _et al._ (2023) A. Jhaveri, S. Loggia, Y. Qian, and M. E. Johnson, bioRxiv (2023), 10.1101/2023.08.30.555551, https://www.biorxiv.org/content/early/2023/09/01/2023.08.30.555551.full.pdf .
  * Curatolo  _et al._ (2023) A. I. Curatolo, O. Kimchi, C. P. Goodrich, R. K. Krueger, and M. P. Brenner, Nature Communications 14, 8328 (2023).
  * Hagan and Chandler (2006) M. F. Hagan and D. Chandler, Biophysical Journal 1, 42 (2006).
  * Perlmutter  _et al._ (2014) J. D. Perlmutter, M. R. Perkett, and M. F. Hagan, Journal of Molecular Biology 426, 3148 (2014).
  * Rapaport (2018) D. C. Rapaport, Journal of Biological Physics 44, 147 (2018).
  * Baschek  _et al._ (2012) J. E. Baschek, H. C. R. Klein, and U. S. Schwarz, Bmc Biophysics 5 (2012), 10.1186/2046-1682-5-22.
  * Gupta  _et al._ (2023) M. Gupta, A. J. Pak, and G. A. Voth, Sci Adv 9, eadd7434 (2023).
  * Zlotnick (2003) A. Zlotnick, Virology 315, 269 (2003).
  * Ceres and Zlotnick (2002) P. Ceres and A. Zlotnick, Biochemistry 41, 11525 (2002).
  * Gartner  _et al._ (2020) F. M. Gartner, I. R. Graf, P. Wilke, P. M. Geiger, and E. Frey, eLife 9, e51020 (2020).
  * Schlicksup and Zlotnick (2020) C. J. Schlicksup and A. Zlotnick, Current Opinion in Virology 45, 43 (2020), 45. Virus structure & Expression (Dec 2020).
  * Hagan and Mohajerani (2023) M. F. Hagan and F. Mohajerani, PLoS Comput. Biol. 19, e1010652 (2023).
  * Trubiano (2024a) A. Trubiano, “Self-Assembly Analysis Suite for HOOMD,” https://github.com/onehalfatsquared/SAASH (2024a).
  * Trubiano (2024b) A. Trubiano, “MultiMSM,” https://github.com/onehalfatsquared/MultiMSM (2024b).
  * Anderson  _et al._ (2020) J. A. Anderson, J. Glaser, and S. C. Glotzer, Computational Materials Science 173, 109363 (2020).
  * Sigl  _et al._ (2021) C. Sigl, E. Willner, W. Engelen, J. Kretzmann, K. Sachenbacher, A. Liedl, F. Kolbe, F. Wilsch, S. Aghvami, U. Protzer, M. Hagan, S. Fraden, and H. Dietz, Nature Materials (2021), 10.1038/s41563-021-01020-4.
  * Wei  _et al._ (2024) W.-S. Wei, A. Trubiano, C. Sigl, S. Paquay, H. Dietz, M. F. Hagan, and S. Fraden, Proceedings of the National Academy of Sciences 121, e2312775121 (2024), https://www.pnas.org/doi/pdf/10.1073/pnas.2312775121 .
  * (113) See Supplemental Material at URL-will-be-inserted-by-publisher for additional details on the performed simulations and implementation details outside the scope of the main text.
  * McMillan  _et al._ (2005) P. F. McMillan, D. C. Clary, and D. J. Wales, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 363, 357 (2005).
  * Fejer  _et al._ (2009) S. N. Fejer, T. R. James, J. HernÃ¡ndez-Rojas, and D. J. Wales, Phys. Chem. Chem. Phys. 11, 2098 (2009).
  * Johnston  _et al._ (2010) I. G. Johnston, A. A. Louis, and J. P. K. Doye, Journal of Physics: Condensed Matter 22, 104101 (2010).
  * Perlmutter  _et al._ (2013) J. D. Perlmutter, C. Qiao, and M. F. Hagan, eLife 2, e00632 (2013).
  * Weeks  _et al._ (1971) J. D. Weeks, D. W. Chandler, and H. C. Andersen, Journal of Chemical Physics 54, 5237 (1971).
  * Rapaport  _et al._ (1999) D. Rapaport, J. Johnson, and J. Skolnick, Computer Physics Communications 121-122, 231 (1999), proceedings of the Europhysics Conference on Computational Physics CCP 1998.
  * Rapaport (2010) D. C. Rapaport, Journal of Physics: Condensed Matter 22, 104115 (2010).
  * Nguyen  _et al._ (2007) H. D. Nguyen, V. S. Reddy, and C. L. Brooks, Nano Letters 7, 338 (2007).
  * Nguyen and Brooks (2008) H. D. Nguyen and C. L. I. Brooks, Nano Letters 8, 4574 (2008).
  * Elrad and Hagan (2010) O. M. Elrad and M. F. Hagan, Physical Biology 7, 045003 (2010).
  * Williamson  _et al._ (2011) A. J. Williamson, A. W. Wilber, J. P. K. Doye, and A. A. Louis, Soft Matter 7, 3423 (2011).
  * Zhang and Linse (2013) R. Zhang and P. Linse, The Journal of Chemical Physics 138, 154901 (2013), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4799243/13247395/154901_1_online.pdf .
  * Zhang and Linse (2014) R. Zhang and P. Linse, The Journal of Chemical Physics 140, 244903 (2014), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4883056/14085490/244903_1_online.pdf .
  * Castelnovo  _et al._ (2014) M. Castelnovo, T. Verdier, and L. Foret, Epl-europhys Lett 105 (2014), 10.1209/0295-5075/105/28006.
  * Castelnovo  _et al._ (2013) M. Castelnovo, D. Muriaux, and C. Faivre-Moskalenko, New J. Phys. 15 (2013), 10.1088/1367-2630/15/3/035028, castelnovo, M. Muriaux, D. Faivre-Moskalenko, C.
  * Boettcher  _et al._ (2015) M. A. Boettcher, H. C. R. Klein, and U. S. Schwarz, Phys. Biol. 12 (2015).
  * Mendoza and Reguera (2020) C. I. Mendoza and D. Reguera, eLife 9, e52525 (2020).
  * Karpen  _et al._ (1993) M. E. Karpen, D. J. Tobias, and C. L. I. Brooks, Biochemistry 32, 412 (1993).
  * de Groot  _et al._ (2001) B. L. de Groot, X. Daura, A. E. Mark, and H. GrubmÃ¼ller, Journal of Molecular Biology 309, 299 (2001).
  * Singhal  _et al._ (2004) N. Singhal, C. D. Snow, and V. S. Pande, The Journal of Chemical Physics 121, 415 (2004), https://pubs.aip.org/aip/jcp/article-pdf/121/1/415/19156582/415_1_online.pdf .
  * Andrec  _et al._ (2005) M. Andrec, A. K. Felts, E. Gallicchio, and R. M. Levy, Proceedings of the National Academy of Sciences 102, 6801 (2005), https://www.pnas.org/doi/pdf/10.1073/pnas.0408970102 .
  * Hempel  _et al._ (2021) T. Hempel, M. J. del Razo, C. T. Lee, B. C. Taylor, R. E. Amaro, and F. Noé, Proceedings of the National Academy of Sciences 118, e2105230118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2105230118 .
  * Satake (1975) I. Satake,  _Linear Algebra_, Pure and applied mathematics (M. Dekker, 1975).
  * Bowman  _et al._ (2010) G. R. Bowman, D. L. Ensign, and V. S. Pande, Journal of Chemical Theory and Computation 6, 787 (2010).
  * Doerr and De Fabritiis (2014) S. Doerr and G. De Fabritiis, Journal of Chemical Theory and Computation 10, 2064 (2014).
  * Voelz  _et al._ (2014) V. A. Voelz, B. Elman, A. M. Razavi, and G. Zhou, Journal of Chemical Theory and Computation 10, 5716 (2014).
  * Zimmerman and Bowman (2015) M. I. Zimmerman and G. R. Bowman, Journal of Chemical Theory and Computation 11, 5747 (2015).
  * Shamsi  _et al._ (2017) Z. Shamsi, A. S. Moffett, and D. Shukla, Sci Rep 7, 12700 (2017).
  * Shamsi  _et al._ (2018) Z. Shamsi, K. J. Cheng, and D. Shukla, The Journal of Physical Chemistry B 122, 8386 (2018).
  * Zimmerman  _et al._ (2018) M. I. Zimmerman, J. R. Porter, X. Sun, R. R. Silva, and G. R. Bowman, J Chem Theory Compute 14, 5459 (2018).
  * Wan and Voelz (2020) H. Wan and V. A. Voelz, J Chem Phys 152, 024103 (2020).
  * Zlotnick (1994) A. Zlotnick, Journal of Molecular Biology 241, 59 (1994).
  * Endres and Zlotnick (2002) D. Endres and A. Zlotnick, Biophysical Journal 83, 1217 (2002).
  * Moisant  _et al._ (2010) P. Moisant, H. Neeman, and A. Zlotnick, Biophys. J. 99, 1350 (2010).
  * Singhal and Pande (2005) N. Singhal and V. S. Pande, The Journal of Chemical Physics 123, 204909 (2005), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2116947/15379877/204909_1_online.pdf .
  * Hinrichs and Pande (2007) N. S. Hinrichs and V. S. Pande, The Journal of Chemical Physics 126, 244101 (2007), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2740261/15400057/244101_1_online.pdf .
  * Bhattacharya and Chatterjee (2019) S. Bhattacharya and A. Chatterjee, The Journal of Chemical Physics 150, 044106 (2019), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5066278/15556619/044106_1_online.pdf .
  * Kozlowski and Grubmüller (2023) N. Kozlowski and H. Grubmüller, Journal of Chemical Theory and Computation 19, 5516 (2023), pMID: 37540193.
  * Efron (1979) B. Efron, The Annals of Statistics 7, 1 (1979).
  * Efron (1981) B. Efron, Biometrika 68, 589 (1981), https://academic.oup.com/biomet/article-pdf/68/3/589/581658/68-3-589.pdf .
  * Trendelkamp-Schroer and Noé (2016) B. Trendelkamp-Schroer and F. Noé, Phys. Rev. X 6, 011009 (2016).
  * Noé  _et al._ (2019) F. Noé, S. Olsson, J. Köhler, and H. Wu, Science 365, eaaw1147 (2019), https://www.science.org/doi/pdf/10.1126/science.aaw1147 .
  * Bonati  _et al._ (2021) L. Bonati, G. Piccini, and M. Parrinello, Proceedings of the National Academy of Sciences 118, e2113533118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2113533118 .
  * Falkner  _et al._ (2023) S. Falkner, A. Coretti, S. Romano, P. L. Geissler, and C. Dellago, Machine Learning: Science and Technology 4, 035050 (2023).
  * Zlotnick  _et al._ (1999) A. Zlotnick, J. M. Johnson, P. W. Wingfield, S. J. Stahl, and D. Endres, Biochemistry 38, 14644 (1999).
  * Hagan and Elrad (2010) M. F. Hagan and O. M. Elrad, Biophysical Journal 98, 1065 (2010).
  * Mizrahi  _et al._ (2022) I. Mizrahi, R. Bruinsma, and J. Rudnick, Phys. Rev. E 106, 044405 (2022).
  * Zlotnick and Stray (2003) A. Zlotnick and S. J. Stray, Trends Biotechnol. 21, 536 (2003).
  * Hagan (2014) M. F. Hagan, “Modeling viral capsid assembly,” in  _Advances in Chemical Physics: Volume 155_ (John Wiley & Sons, Ltd, 2014) Chap. 1, pp. 1–68, https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118755815.ch01 .
  * Prevelige  _et al._ (1993) P. E. Prevelige, D. Thomas, and J. King, Biophys. J. 64, 824 (1993).
  * Schwartz  _et al._ (1998) R. Schwartz, P. W. Shor, P. E. Prevelige, and B. Berger, Biophys. J. 75, 2626 (1998).
  * Zhang and Schwartz (2006b) T. Q. Zhang and R. Schwartz, Biophys. J. 90, 57 (2006b).
  * Tan  _et al._ (2021) A. Tan, A. J. Pak, D. R. Morado, G. A. Voth, and J. A. G. Briggs, Proc Natl Acad Sci U S A 118 (2021), 10.1073/pnas.2020054118.
  * Yu  _et al._ (2021) A. Yu, A. J. Pak, P. He, V. Monje-Galvan, L. Casalino, Z. Gaieb, A. C. Dommer, R. E. Amaro, and G. A. Voth, Biophysical journal 120, 1097 (2021).
  * Panahandeh  _et al._ (2020) S. Panahandeh, S. Li, L. Marichal, R. Leite Rubim, G. Tresset, and R. Zandi, ACS Nano 14, 3170 (2020).
  * Rapaport (2008) D. Rapaport, Phys. Rev. Lett. 101, 186101 (2008).
  * Whitelam  _et al._ (2009) S. Whitelam, C. Rogers, A. Pasqua, C. Paavola, J. Trent, and P. L. Geissler, Nano Lett. 9, 292 (2009).
  * Wilber  _et al._ (2007) A. W. Wilber, J. P. K. Doye, A. A. Louis, E. G. Noya, M. A. Miller, and P. Wong, The Journal of Chemical Physics 127, 085106 (2007).
  * Wilber  _et al._ (2009) A. W. Wilber, J. P. K. Doye, A. A. Louis, and A. C. F. Lewis, J. Chem. Phys. 131, 175102 (2009).
  * Hagan  _et al._ (2011) M. F. Hagan, O. M. Elrad, and R. L. Jack, J. Chem. Phys. 135, 104115 (2011).
  * Cheng  _et al._ (2012) S. Cheng, A. Aggarwal, and M. J. Stevens, Soft Matter 8, 5666 (2012).
  * Rapaport (2004) D. Rapaport, Phys. Rev. E. 70, 051905 (2004).
  * Nguyen  _et al._ (2009) H. D. Nguyen, V. S. Reddy, and C. L. Brooks, J. Am. Chem. Soc. 131, 2606 (2009).
  * Rapaport (2012) D. C. Rapaport, Phys. Rev. E 86, 051917 (2012).
  * Wagner and Zandi (2015) J. Wagner and R. Zandi, Biophys. J. 109, 956 (2015).
  * van der Schoot and Zandi (2007) P. van der Schoot and R. Zandi, Phys. Biol. 4, 296 (2007).
  * Vanden-Eijnden (2006) E. Vanden-Eijnden, “Transition path theory,” in  _Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1_, edited by M. Ferrario, G. Ciccotti, and K. Binder (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 453–493.
  * Metzner  _et al._ (2009) P. Metzner, C. Schütte, and E. Vanden-Eijnden, Multiscale Modeling & Simulation 7, 1192 (2009).
  * E and Vanden-Eijnden (2010) W. E and E. Vanden-Eijnden, Annual Review of Physical Chemistry 61, 391 (2010).
  * Helfmann  _et al._ (2020) L. Helfmann, E. Ribera-Borel, C. Schütte, and P. Koltai, The Journal of Nonlinear Science 30, 3321 (2020).
  * Dellago  _et al._ (2002) C. Dellago, P. G. Bolhuis, and P. L. Geissler, Adv. Chem. Phys. 123, 1 (2002).
  * Peters (2016) B. Peters, Annu. Rev. Phys. Chem. 67, 669 (2016).
  * Bag (2002) B. C. Bag, Phys. Rev. E 66, 026122 (2002).
  * Vaikuntanathan and Jarzynski (2009) S. Vaikuntanathan and C. Jarzynski, Europhysics Letters 87, 60005 (2009).
  * Mandal  _et al._ (2017) D. Mandal, K. Klymko, and M. R. DeWeese, Phys. Rev. Lett. 119, 258001 (2017).
  * Landi and Paternostro (2021) G. T. Landi and M. Paternostro, Rev. Mod. Phys. 93, 035008 (2021).
  * Gomez-Marin  _et al._ (2008) A. Gomez-Marin, T. Schmiedl, and U. Seifert, The Journal of Chemical Physics 129, 024114 (2008), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2948948/6770253/024114_1_online.pdf .
  * Sivak and Crooks (2016) D. A. Sivak and G. E. Crooks, Phys. Rev. E 94, 052106 (2016).
  * Large and Sivak (2019) S. J. Large and D. A. Sivak, Journal of Statistical Mechanics: Theory and Experiment 2019, 083212 (2019).
  * Tafoya  _et al._ (2019) S. Tafoya, S. J. Large, S. Liu, C. Bustamante, and D. A. Sivak, Proceedings of the National Academy of Sciences 116, 5920 (2019), https://www.pnas.org/doi/pdf/10.1073/pnas.1817778116 .
  * Gingrich  _et al._ (2016) T. R. Gingrich, G. M. Rotskoff, G. E. Crooks, and P. L. Geissler, Proceedings of the National Academy of Sciences 113, 10263 (2016), https://www.pnas.org/doi/pdf/10.1073/pnas.1606273113 .
  * Nakazato and Ito (2021) M. Nakazato and S. Ito, Phys. Rev. Res. 3, 043093 (2021).
  * Dechant (2022) A. Dechant, Journal of Physics A: Mathematical and Theoretical 55, 094001 (2022).
  * Blaber and Sivak (2023) S. Blaber and D. A. Sivak, Journal of Physics Communications 7, 033001 (2023).
  * Brown and Sivak (2020) A. I. Brown and D. A. Sivak, Chemical Reviews 120, 434 (2020).
  * Abiuso  _et al._ (2022) P. Abiuso, V. Holubec, J. Anders, Z. Ye, F. Cerisola, and M. Perarnau-Llobet, Journal of Physics Communications 6, 063001 (2022).
  * Louwerse and Sivak (2022) M. D. Louwerse and D. A. Sivak, The Journal of Chemical Physics 156, 194108 (2022), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0086079/16541262/194108_1_online.pdf .
  * Whitelam and Schmit (2022) S. Whitelam and J. D. Schmit, Journal of Crystal Growth 600, 126912 (2022).
  * Zhong and DeWeese (2022) A. Zhong and M. R. DeWeese, Phys. Rev. E 106, 044135 (2022).
  * Davis  _et al._ (2024) L. K. Davis, K. Proesmans, and E. Fodor, Phys. Rev. X 14, 011012 (2024).
  * Kim  _et al._ (2020) D.-K. Kim, Y. Bae, S. Lee, and H. Jeong, Phys. Rev. Lett. 125, 140604 (2020).
  * Nir  _et al._ (2020) A. Nir, E. Sela, R. Beck, and Y. Bar-Sinai, Proceedings of the National Academy of Sciences 117, 30234 (2020), https://www.pnas.org/doi/pdf/10.1073/pnas.2017042117 .
  * Otsubo  _et al._ (2020) S. Otsubo, S. Ito, A. Dechant, and T. Sagawa, Phys. Rev. E 101, 062106 (2020).
  * Otsubo  _et al._ (2022) S. Otsubo, S. K. Manikandan, T. Sagawa, and S. Krishnamurthy, Communications Physics 5, 11 (2022).
  * Boffi and Vanden-Eijnden (2023a) N. M. Boffi and E. Vanden-Eijnden, Machine Learning: Science and Technology 4, 035012 (2023a).
  * Boffi and Vanden-Eijnden (2023b) N. M. Boffi and E. Vanden-Eijnden, “Deep learning probability flows and entropy production rates in active matter,” (2023b), arXiv:2309.12991 [cond-mat.stat-mech] .
  * Whitelam (2023) S. Whitelam, Phys. Rev. X 13, 021005 (2023).
  * Engel  _et al._ (2023) M. C. Engel, J. A. Smith, and M. P. Brenner, Phys. Rev. X 13, 041032 (2023).
  * Zlotnick (2007) A. Zlotnick, J. Mol. Biol. 366, 14 (2007).
  * Kim  _et al._ (2021) C. Kim, C. Schlicksup, L. Barnes, M. Jarrold, A. Patterson, B. Bothner, and A. Zlotnick, Microscopy and Microanalysis 27, 1130 (2021).
  * Kondylis  _et al._ (2018) P. Kondylis, C. J. Schlicksup, N. E. Brunk, J. Zhou, A. Zlotnick, and S. C. Jacobson, Journal of the American Chemical Society 141, 1251 (2018).
  * Schlicksup  _et al._ (2018) C. J. Schlicksup, J. C.-Y. Wang, S. Francis, B. Venkatakrishnan, W. W. Turner, M. VanNieuwenhze, and A. Zlotnick, Elife 7, e31473 (2018).
  * Pavlova  _et al._ (2022) A. Pavlova, L. Bassit, B. D. Cox, M. Korablyov, C. Chipot, D. Patel, D. L. Lynch, F. Amblard, R. F. Schinazi, and J. C. Gumbart, J. Med. Chem. 65, 4854 (2022).
  * Wang  _et al._ (2023) M. Wang, J. Zhang, Y. Dou, M. Liang, Y. Xie, P. Xue, L. Liu, C. Li, Y. Wang, F. Tao, X. Zhang, H. Hu, K. Feng, L. Zhang, Z. Wu, Y. Chen, P. Zhan, and H. Jia, J. Med. Chem. 66, 13968 (2023).
  * Kra  _et al._ (2023) K. Kra, S. Li, L. Gargowitsch, J. Degrouard, J. Pérez, R. Zandi, S. Bressanelli, and G. Tresset, ACS Nano 17, 12723 (2023).
  * Wu  _et al._ (2018) H. Wu, A. Mardt, L. Pasquali, and F. Noé, in  _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS 18 (Curran Associates Inc., Red Hook, NY, USA, 2018) pp. 3979–3988.
  * Swenson  _et al._ (2019) D. W. H. Swenson, J.-H. Prinz, F. Noé, J. D. Chodera, and P. G. Bolhuis, Journal of Chemical Theory and Computation 15, 813 (2019).
  * Bolhuis and Swenson (2021) P. G. Bolhuis and D. W. H. Swenson, Advanced Theory and Simulations 4, 2000237 (2021), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adts.202000237 .

Markov State Model Approach to Study Self-Assembly - Supplement

Anthony Trubiano and Michael Hagan
---

Martin Fisher School of Physics, Brandeis University, Waltham, Massachusetts
02454, USA

##  S1 Simulation Details

###  S1.1 Dodecahedron system

The subunits are rigid bodies composed of several kinds of pseudoatoms. There
are attractor pseudoatoms (‘A’) at the vertices of a regular pentagon, which
facilitate subunit assembly via an attractive Morse potential with an
equilibrium length of $L_{0}=0.2$, a range parameter $\alpha=2.5/L_{0}$, a
cutoff distance of $2$, and a well-depth (subunit-subunit binding strength)
$\epsilon_{11}$ that can be varied. In this work, we use values
$\epsilon_{11}\in\\{5.0,5.5,6.0\\}$. We also include a top pseudoatom (‘T’)
and a bottom pseudoatom (‘B’), located at positions $z=\pm 0.5$ with respect
to the center of the regular pentagon in the xy plane. The ‘T’ pseudoatoms
interact with other ‘T’ pseudoatoms via a repulsive Lennard-Jones potential
with $\sigma_{\text{TT}}=2.1$, cutoff distance equal to $\sigma_{\text{TT}}$,
and well-depth $\epsilon_{\text{TT}}=\epsilon_{11}/4$. These values favor a
subunit-subunit binding angle consistent with that of a dodecahedron. The ‘B’
pseudoatoms have a similar repulsive interaction with ‘T’ pseudoatoms, with
$\sigma_{\text{TB}}=1.8$, a cutoff distance equal to $\sigma_{\text{TB}}$, and
well-depth $\epsilon_{\text{TB}}=\epsilon_{11}/4$. This interaction helps to
prevent upside down assembly, i.e. ensuring that the top atom is in the
direction of the outward normal vector. Finally, we add edge (‘E’) pseudoatoms
at the midpoint between each adjacent vertex. These have no interactions and
thus do not affect the simulation, but they are used to track assembly
progress more easily compared to vertex pseudoatoms.

Results are reported in units for which the unit length $l_{0}$ corresponds to
the edge length of the pentagonal subunit and energies are measured in units
of $k_{\text{B}}T$. Simulations are initialized with $125$ pentagonal
subunits, enough to form $10$ dodecahedral capsids ($12$ subunits each), on an
equally spaced lattice. Subunit positions and orientations are then
equilibrated with a purely repulsive potential for $2\times 10^{5}$ time steps
before writing any output. The simulation box is a cube with periodic boundary
conditions and side lengths $20l_{0}$, giving a total subunit concentration of
$c_{0}=0.0156/l_{0}^{3}$. The time step is $0.001t_{0}$, and the base
simulations to construct MSMs are run for $5\times 10^{6}$ time steps unless
otherwise specified. The simulations use the HOOMD-blue Anderson  _et al._
(2020) version 3.9.0 Langevin integrator, with an inverse temperature
$\beta=1$. The configurations are logged every $\Delta t=t_{0}$ units of
simulation time.

###  S1.2 $T=3$ capsid system

The subunits are rigid triangles with each edge consisting of three stacked
layers of six overlapping ‘excluder’ pseudoatoms at a specified bevel angle.
These excluder pseudoatoms interact with all other pseudoatoms in the
simulation through a Weeks-Chandler-Anderson (WCA) potential to enforce
excluded volume. Embedded in the middle row of each edge of the triangle are
two attractor pseudoatoms, which have attractive Lennard-Jones interactions
with the pseudoatoms on complementary edges to facilitate edge-edge binding of
the subunits. Each pseudoatom has the same diameter, $\sigma$, which we set as
the unit distance $l_{0}$. To match the dimensions of the experimental subunit
Wei  _et al._ (2024), we can set $\sigma=18$nm in real units.

As described so far, the triangular subunits can be designed to form a broad
variety of target structures by tuning the interaction strengths, side
lengths, and bevel angles. For a $T=3$ capsid target, we make two sides of the
triangle equivalent. Sides 1 and 2 have an edge length of $3\sigma=54$nm and
have complementary pseudoatoms (Side 1 has pseudoatoms ‘4’ and ‘5’, which bind
with pseudoatoms ‘7’ and ‘6’ on Side 2, respectively) that attract with a
binding energy of $\epsilon_{12}$. These two sides do not interact at all with
Side 3, which has a slightly longer edge length of $3.35\sigma=60.3$nm, with
attractive pseudoatoms (pseudoatoms ‘2’ and ‘3’) that are self-complementary
with a binding well-depth of $\epsilon_{33}$. The bevel angle of each side is
the same, approximately $11.64^{\circ}$. This design can produce hierarchical
assembly Wei  _et al._ (2024). When $\epsilon_{33}\gg\epsilon_{12}$, dimers
form rapidly via the Side 3 – Side 3 interaction, and then the dimers more
slowly assemble into larger structures via the Side 1 – Side 2 interaction.
Conversely, when $\epsilon_{33}\ll\epsilon_{12}$, pentamers form rapidly and
then subsequently assemble into larger structures. In this work, we consider
an example of pentamer-biased assembly, with $\epsilon_{12}=11$ and
$\epsilon_{33}=8$.

The simulations contain $N=600$ triangular subunits, enough to form $10$ $T=3$
capsids ($60$ subunits each). The simulation domain is an $L\times L\times L$
box with periodic boundary conditions, whose side lengths determine the total
subunit concentration, $C=N/(N_{\text{A}}L^{3})$, where $N_{\text{A}}$ is
Avogadro’s number. We set $C=50$nM to be on the order of experimental
conditions, for which the corresponding box side length is $L=2.71$ microns.
In simulation units, the box side lengths are $153l_{0}$, giving a total
subunit concentration of $c_{0}=1.7\times 10^{-4}/l_{0}^{3}$.

We initialize subunit positions on an equally spaced, truncated lattice.
Subunit positions and orientation are then equilibrated with a purely
repulsive potential for $8\times 10^{4}$ time steps, after which attractive
interactions are turned on. We use a time step of $0.0025t_{0}$, and the base
simulations to construct MSMs are run for $5\times 10^{8}$ time steps unless
stated otherwise.

###  S1.3 Bond Definitions and Discrete States

For the dodecahedron system, the pentagonal subunits have edge (‘E’)
pseudoatoms that align when the adjacent vertex attractors bind with another
subunit. We use a cutoff distance of $0.3l_{0}$ between these edge pseudoatoms
to define a bond. For the $T=3$ capsid assembly, there are a pair of
complementary pseudoatoms on Sides 1 and 2 (‘4’ and ‘7’ pseudoatoms), and Side
3 is self-complementary (‘2’ and ‘3’ pseudoatoms). We use a cutoff distance of
$1.3l_{0}$ for both pairs of complementary pseudoatoms in this case. These
values were chosen by selecting the minimum distance that would correctly
identify fully formed capsids across many different realizations of the
assembly. Our results are insensitive to increasing these cutoffs up to the
next-nearest neighbor bond distances, which are approximately $0.8l_{0}$ for
the pentagonal subunit and $2.4l_{0}$ for the triangular subunits. Decreasing
these cutoffs results in frequent oscillations in the number of bonds, even
for stable configurations like the fully formed dodecahedron. Such behavior is
undesirable when building a Markov state model, as it overestimates the rate
of bond breakage in the model, so we ensure our choice of cutoff correctly
identifies the target structure over long timescales at the strongest binding
energies.

Using this definition of a bond, we construct our discrete state space as all
observed combinations of $(\mathcal{N},\mathcal{B})$, where $\mathcal{N}$ is
the number of subunits in a cluster and $\mathcal{B}$ is the distribution of
the number of bonds in a cluster. Defined in this way, the number of unique
discrete states we observe is system dependent. There are $180$ states for
dodecahedral capsid assembly with binding energies $\epsilon_{11}=5$ and
$\epsilon_{11}=5.5$. For $\epsilon_{11}=6$, the interactions are strong enough
to result in transient intermediates that are larger than $12$ subunits,
resulting in a state space size of $447$. Despite the enhanced size, most of
these states are extremely low probability and could likely be excluded
without affecting the model predictions, though we have not tested this. For
$T=3$ capsid assembly, the number of states is much larger due to a larger
capsid size and having two distinguishable bond types. For the pentamer-biased
assembly set that we consider, the state space size is $1662$. Since the
hierarchical assembly pathways likely decreases the number of intermediates,
the state space might be larger for other parameter sets.

##  S2 Transition Counting Example

Figure S2.1: A simple assembly example showing two equivalent approaches to
count transitions. (a) A simple example time series of four disk shaped
particles with attractive interactions undergoing clusterization. (b) A table
detailing how to count transitions from both a subunit and cluster perspective
for this example. In the subunit perspective, a time series of the cluster
size containing each subunit is computed, and a count matrix is built from the
number of times each transition is observed. In the cluster perspective, only
merging and splitting of clusters is tracked. The count matrix is built by
tallying the number of times a cluster of size $i$ points to a cluster of size
$j$, and then multiplied by $W_{ij}=\min(i,j)$ to get the counts in the
subunit perspective.

As stated in the main text, the construction of the MSM count matrix from the
set of cluster transitions must maintain mass weighting. There are two
equivalent ways to count transitions such that the mass weighting is
preserved. The first is to take a subunit perspective; for each individual
subunit, compute the time series of its cluster size and count every
transition that occurs. Since this procedure loops over each subunit in a
given cluster, it automatically accounts for the mass-weighting. This is at
the cost of efficiency however, since each subunit in a cluster stores the
same transition information. Instead, we take a cluster perspective, which is
described in Section III.2. Here, each cluster transition is only counted a
single time and then weighted in the count matrix by the number of subunits in
the smaller cluster. Since each entry of the count matrix will be weighted by