Robust nonequilibrium pathways to microcompartment assembly
Grant M. Rotskoff Phillip L. Geissler




Courant Institute of Mathematical Sciences, New York University, New York, NY 10002
Department of Chemistry, University of California, Berkeley, CA 94720
 
Department of Chemistry, University of California, Berkeley, CA 94720




  Cyanobacteria sequester photosynthetic enzymes into microcompartments which facilitate the conversion of carbon dioxide into sugars. 
  Geometric similarities between these structures and self-assembling viral capsids have inspired models that posit microcompartments as stable equilibrium arrangements of the constituent proteins. 
  Here we describe a different mechanism for microcompartment assembly, one that is fundamentally nonequilibrium and yet highly reliable. 
  This pathway is revealed by simulations of a molecular model resolving the size and shape of a cargo droplet, and the extent and topography of an elastic shell. 
  The resulting metastable microcompartment structures closely resemble those of carboxysomes, with a narrow size distribution and faceted shells. 
  The essence of their assembly dynamics can be understood from a simpler mathematical model that combines elements of classical nucleation theory with continuum elasticity. 
  These results highlight important control variables for achieving nanoscale encapsulation in general, and for modulating the size and shape of carboxysomes in particular.



Robust nonequilibrium pathways to microcompartment assembly
    Phillip L. Geissler
    December 30, 2023
===========================================================



Spatial segregation is an ubiquitous strategy in biology for organizing the crowded, active viscera of the cell <cit.>.
Viral capsids exemplify this organization at very small scales, sequestering genetic material from the cytosol and recapturing it for delivery to new hosts.  
Extensive work has explored the structure, stability, and assembly dynamics of viruses, highlighting generic design principles and physical origins of the spontaneous assembly process <cit.>.  
Overall capsid structure typically follows from the arrangements of neighboring proteins that are preferred by their noncovalent interactions.  
Strong preferences yield regular and highly stable structures, but at the same time impair kinetic accessibility by producing deep kinetic traps <cit.>.

Bacterial microcompartments serve a very different biomolecular purpose from viruses but have striking structural similarities, namely, a quasi-icosahedral protein shell that assembles around a fluctuating cargo <cit.>.
This comparison raises the question: Do the same assembly principles, based on a balance between equilibrium stability and kinetic accessibility, apply to microcompartments as well?
Here we focus on a paradigmatic example of a microcompartment, the carboxysome.
Carboxysomes play an essential role in the carbon fixation pathway of photosynthetic cyanobacteria <cit.>. 
The shell proteins of the carboxysome, which have hexameric and pentameric crystallographic structures, assemble to encapsulate a condensed globule of protein including the enzymes RuBisCO and carbonic anhydrase <cit.>.  
These vital, organelle-like structures regulate a microscopic environment to enhance catalytic efficiency, which has made them an attractive target for bioengineering applications <cit.>.





Like some viral capsids, carboxysomes feature pronounced facets and vertices, and crystal structures of shell components suggest highly symmetric local protein arrangements <cit.>.  
Unlike in viral capsids, however, these apparent preferences do not directly indicate the assembled structure, nor do they clearly point to a characteristic microcompartment size <cit.>.  
Indeed, experimental evidence has not constrained a particular mechanism for assembly: live cell images support the notion that the cargo assembles first and is subsequently coated by the shell <cit.>, while direct observations of partially formed carboxysomes in electron micrographs point to a cooperative mechanism for assembly <cit.>.
Identifying essential variables for controlling or emulating carboxysome assembly awaits a clearer understanding of its dynamical pathways and underlying driving forces.

Theoretical and computational models for nanoshell assembly have
primarily followed approaches inspired by viral capsid assembly.
These approaches have emphasized the role of preferred shell curvature and
shell-cargo interactions as a template for the final
structure <cit.>.
For example, irreversible growth of shells comprising monomers
that prefer a bent binding geometry has been shown in simulations
to successfully assemble empty shells
 <cit.>.
The assumption of such a preferred local curvature, however, is at odds
with structural models of the carboxysome based on crystallography of
shell proteins, which feature tiling of the shell by hexameric proteins
that appear to bind optimally with zero curvature <cit.>
Indeed, a number of experimental measurements bolster the view that
bacterial microcompartments are bounded by protein sheets that are
essentially flat away from localized regions of high curvature.  The
predominant constituents of diverse microcompartments (e.g.,
β-carboxysomes, α-carboxsyomes, eut
microcompartments, and pdu microcompartments) crystallize as
layers of flat sheets comprising hexamers whose lateral contacts are
genetically conserved
 <cit.>.  Recent atomic force
microscopy measurements have demonstrated
that these constituents also form flat monolayers in physiological buffer
conditions <cit.>.
Furthermore, in vivo tomography
data strongly suggest that biological carboxysomes have extended flat
faces with curvature sharply localized at the joint of a
facet <cit.>.

Here we introduce and explore a model based on a mechanistic
perspective that is fundamentally different from the one commonly
applied to virus assembly and that does not require assuming an
innately preferred curvature for contacts between shell proteins.  The
basic components are a cargo species that is prone to aggregation, a
shell species that spontaneously forms flat,
hexagonally symmetric
elastic sheets, and an
attractive interaction between the inside of the shell and the cargo,
depicted in Figs. 1A, S1-S3.  These ingredients
appear to be
the
essential
constituents for carboxysome biogenesis in vitro—mutagenesis
experiments have shown that pentameric proteins
sometimes presumed to stabilize shell curvature are in fact
not necessary for the formation of faceted shell
structures <cit.>.
For such a basic model, thermodynamic
considerations imply that finite encapsulated structures
have negligible weight at equilibrium (Sec. S3).  Nevertheless,
we show that regularly sized microcompartments can emerge reliably in
the course of natural dynamics.



§ METHODS


We regard the protein shell of a growing carboxysome as
a thin elastic sheet, whose energy 
can be
discretized and expressed as a sum over the edges ij of a
triangulated surface <cit.>.
Its Hamiltonian is

    ℋ_ shell = ϵ/2∑_ij (l_ij-l_0)^2 + κ/2∑_ij 1-cos(θ_ij) + ℋ_ steric,

where l_ij is the length of an edge, whose deviations from the
preferred length l_0 experience a restoring force with stiffness
ϵ. The bending rigidity κ sets the energy scale for
developing a nonzero angle θ_ij between normal vectors of
triangles sharing the edge ij <cit.>.

We associate each triangular face in the discretized sheet with a
protein monomer.  The term ℋ_ steric in (<ref>)
imposes steric constraints that prevent overlap of these monomers, by
placing a hard sphere of diameter 0.45 l_0  at the center of each triangle
in the sheet.  
Because the deformations of individual shell proteins
are expected to be small relative to the bending fluctuations between
adjacent monomers, we work in a limit where ϵ l_0^2 ≫κ. The cohesive interactions that bind adjacent shell proteins
in the sheet are accounted for separately, with an energy
ℋ_ bind that tightly associates the vertices of contacting
monomers. In the limit that bound vertices coincide exactly, the
thermodynamic influence of ℋ_ bind depends only on
the connectivity of the sheet, contributing a binding affinity factor
K for each constrained vertex, as described in SI.

We represent aggregating cargo as an Ising lattice gas on an FCC lattice
with chemical potential μ_c, restricted to configurations
with a single connected cargo droplet.
Nearest neighbors are coupled through an interaction energy ϵ_c,

    ℋ_ cargo = -μ_c∑_i=1^N_Lσ_i - ϵ_c∑_ijσ_i σ_j,

where σ_i = {0, 1} for occupied and unoccupied lattice sites,
respectively, ij indicates that the sum is taken over
nearest neighbors, and N_L is the number of lattice sites.  
The lattice spacing used throughout is l_0, the
average length of a shell monomer edge.
The primary cargo species in a carboxysome, the protein RuBisCO, has a
diameter L_0 that is larger than l_0; adopting a finer lattice
spacing amounts to averaging approximately over configurational
fluctuations of RuBisCO about a close-packed lattice.

Interactions between cargo and shell species in our model mimic a
short-ranged directional attraction suggested by structural
data <cit.>, and the steric repulsion intrinsic to
compact macromolecules. For a particular lattice site i and shell
monomer j, 

    ℋ_int(σ_i, n̂_j) = 
    
      -γ_in   if  -n̂_j ∈𝒱(i)  and σ_i = 1 
    
        +γ_out   if n̂_j ∈𝒱(i)  and σ_i = 1 
    
        0    otherwise,

where n̂_j
is
the outward facing normal vector of shell monomer j.  
The notation ∈𝒱(i) means that the point specified by the vector lies within the volume 𝒱(i) occupied by lattice cell i.
The full interaction is the sum of H_int over all i and j.  
This interaction, which contributes a favorable energy if the 
inward-facing normal lies within a cargo-occupied lattice site
and an unfavorable energy if the outward-facing normal vector
does so, is depicted schematically in Fig. S5.

The thermal relaxation of a shell with a fixed number of vertices and
fixed connectivity, subject to the energetics described above,
can be simulated with a straightforward Monte Carlo
dynamics.  In each move, a single vertex
is chosen
at random.  As depicted in Fig. S2 (A), a random perturbation in
three-dimensional space is made to the selected vertex,
attempting to change its position from a to a'.
The resulting energy difference,
Δℋ_ shell + Δℋ_ int,
is used in a standard Metropolis acceptance
criterion,

    acc(a→ a') = min[ 1, e^ -β
          (Δℋ_ shell + Δℋ_ int)].

This procedure ensures that the configurations sampled in the Markov
chain are consistent with a Boltzmann distribution.

A Monte Carlo dynamics for growth of the shell is significantly more
involved.  In order that the procedure satisfy detailed balance, care
must be taken
to ensure that every step is reversible and
that algorithmic asymmetries of forward and
reverse moves are fully accounted for.
Our goal is to draw samples from a grand
canonical distribution,

    p(X, N_s, {σ}) = Ξ^-1 z_s^N_s e^-β (ℋ_ shell
        +ℋ_ bind + ℋ_ cargo + ℋ_ int),

where Ξ denotes the grand canonical partition function, X gives
the coordinates of the shell protein configuration, N_s is the number
of shell monomers, and z_s∝ e^βμ_s is an activity set by
the concentration of free shell monomers in solution.
We implement two types of Monte Carlo moves
that allow the structure to grow,
increasing the number of monomers by adding new vertices at the edge
of the shell.  Fusion of pre-existing vertices can also occur,
representing the binding between monomers in the shell that are nearby
in space but not necessarily in connectivy.  Details of these moves
and their reverse counterparts, along with acceptance criteria that
ensure detailed balance with respect to p(X, N_s, {σ}), are described in
Sec. SI 1.  In microcompartment growth simulations, they are
performed in tandem with standard trial moves that change the
occupation state of the cargo lattice, acting to grow or shrink the
cargo droplet.

The routes of Monte Carlo trajectories propagated in this way are
influenced by energetic parameters like ϵ, κ, and
μ, and by the shell binding affinity K.
The pathways are also shaped by the relative
frequencies of proposing each type of move. A basic time scale τ
is set by the duration of a “sweep” in which each vertex experiences
a single attempted spatial displacement (on average). For every n_c
sweeps we attempt roughly N_ boundary moves that add or remove
material at the surface of the cargo droplet, where N_ boundary
is the number of lattice sites defining the droplet surface. This procedure
establishes an effective rate k_c^0 = (2 n_cτ)^-1 at which
cargo monomers arrive from solution at a given surface
site. Similarly, a basic rate k_s^0 = (2 n_sτ)^-1 for shell
monomer arrival at a given site on the shell's perimeter is
established by attempting N_perim shell monomer
addition/removal moves for every n_s sweeps, where N_perim
is the number of edges at the shell boundary. In experiments these
arrival rates are approximately proportional to the solution
concentrations of the respective monomers. 
As a result,
microcompartment mass changes slowly on the time scale of elastic
fluctuations.  Dynamics of vertex fusion and fission also proceed much more
rapidly than growth.
In essence, all
processes other than monomer addition and removal follow along
adiabatically.
For the fate of assembly, the key dynamical parameter in our model is
therefore the ratio k_c^0/k_s^0 of arrival rates for cargo and shell
species.



§ MOLECULAR SIMULATIONS


Fig. <ref> depicts an example assembly pathway of this molecular model.  
Trajectories are initiated with a small droplet of cargo
(comprising a few hundred cargo monomers)
and a handful of proximate shell monomers, as described in Sec. S5.  
Under conditions favorable for assembly,
such a droplet faces no thermodynamic barrier to growth; absent interactions with the shell, its radius R would increase at a constant average rate as cargo material arrives from the droplet's surroundings.  
Growth of the shell, on the other hand, is impeded by the energetic cost of wrapping an elastic sheet around a highly curved object.  
In early stages of this trajectory, shown in Fig. <ref> (C), the net cost of encapsulation is considerable, and the population of shell monomers remains small as a result.

As the droplet grows,
this elastic penalty
is eventually overwhelmed by
attractions between shell and cargo, similar to the mechanism of
curvature generation by nanoparticles adsorbed on membrane
surfaces <cit.>.  At a
characteristic droplet size R^*, encapsulation becomes
thermodynamically favorable.  If the arrival rate of shell monomers is
much higher than that of cargo (e.g., due to a higher concentration in
solution), then a nearly complete shell will quickly develop,
Fig. <ref> (D), hampering the incorporation of additional
cargo.  The nascent shell that results, while sufficient to block
cargo arrival, is highly defective and far from closed.  A slow
healing process ensues, dominated by relaxation of grain boundaries
between growth faces, Fig. <ref> (E).  This annealing process
is, in part, achieved by relocating topological defects in the
structure, a phenomenon that has been previously encountered in ground
state calculations for closed elastic
shells <cit.>.


In the vast majority of the hundreds of assembly trajectories we have
generated, healing leads ultimately to a completely closed structure
with exactly 12 five-fold defects.  Placement of these defects is
often irregular and unlikely to yield minimum elastic energy, but
further evolution of the structure is extraordinarily slow.
Subsequent defect dynamics could produce more ideal shell structures
as in Ref. <cit.>, and transient shell opening
could allow additional cargo growth. But these relaxation processes
require the removal of shell monomers that are bound to several others
and that interact strongly with enclosed cargo. Under the conditions
of interest, the time scale for such a removal is vastly longer than
the assembly trajectories we propagate. 
Our model microcompartments are
equilibrated with respect to neither shell geometry nor droplet size,
yet they are profoundly metastable, requiring extremely rare events to
advance towards the true equilibrium state.

In addition to generating molecular trajectories of microcompartment
assembly, we have performed umbrella sampling simulations to compute
the equilibrium free energy of the molecular model as a function of
N_s and the number N_c of cargo
monomers. Results, presented in Sec. S4, underscore the mechanistic
features described above: a small critical nucleus size for cargo
condensation, and a size-dependent thermodynamic bias on encapsulation
that becomes sharply favorable at a characteristic length scale.  In
the next section we examine the origins and consequences of these
features in a more idealized context.





§ MINIMALIST MODEL

The scenario described above for our molecular model can be cast in a
simpler light by considering as dynamical variables only the amounts
of cargo and shell material in an idealized geometry.  In this
minimalist approach we focus on the radius R of a spherical droplet
and the polar angle θ subtended by a contacting spherical cap
of shell material, as depicted in Fig. <ref> (A). Taking the
cargo and shell species to have uniform densities ρ and ν
within the droplet and cap, respectively, the geometric parameters R
and θ can be simply related to the monomer populations
N_c = (4πρ/3)R^3 and N_s = 2πν R^2
(1-cosθ).


In the
spirit of classical nucleation theory, we estimate a free energy
landscape in the space of R and θ through considerations of
surface tension and bulk thermodynamics.  For our system these
contributions include free energy of the condensed cargo phase,
surface tension of a bare cargo droplet, free energy of a macroscopic
elastic shell, and line tension of a finite shell, together with the
energetics of shell bending and shell-cargo attraction.  As shown in
the Sec. S3, this free energy F(R,θ) can be written in the form



    F(R, θ) = a_c[ ( R/R_c)^2 - ( R/R_c)^3 ] + a_s[ {(R^*/R_s)^2 - (R/R_s)^2 }(1-cosθ) + R/R_ssinθ],



where the energy scales a_c and a_s and the length scales R_c, R^*, and R_s (all of which can be related to parameters of the molecular model) transparently reflect the thermodynamic biases governing growth and encapsulation.
The function F(R,θ) can be minimum only at R=0 and R=∞.  As in the molecular model, finite microcompartments can appear and persist only as kinetically trapped nonequilibrium states.

Contours of this schematic free energy surface, plotted in
Fig. <ref> (A), are consistent with basic thermodynamic trends
observed for the molecular model. Most importantly, they manifest the
prohibitive cost of encapsulating small droplets.  Only for radii
R>R^* is shell growth thermodynamically favorable.  This
characteristic length scale is determined by a straightforward balance
among the shell's bending rigidity, the macroscopic driving force for
flat shell growth, and the attraction between shell and cargo, as
detailed in Sec. S3.  Whether encapsulation occurs immediately once
the droplet exceeds the critical size R^* depends on kinetic factors
that are only partly governed by F(R,θ).  Line tension of the
shell, represented by the final term in Eq. <ref>, effects a
barrier to encapsulation.  As the droplet size increases, the height
of this barrier declines, and encapsulation becomes facile only when
it is sufficiently low. As with our molecular model, the subsequent
dynamics may deviate significantly from equilibrium expecations based
on the free energy surface. To explore these nonequilibrium outcomes,
we must supplement the free energy surface F(R,θ) with a
correspondingly minimal set of dynamical rules.





We take the rate of shell monomer addition to be a product of the
per-edge rate of monomer arrival k_s^0 and the number of
edges that can bind new monomers,

    k_s^+(N_s, N_c)   = k_s^0 2π
        R/l_0sinθ

Similarly accounting for the number of sites
where cargo can be added, we take the rate of cargo monomer
addition to be

    k_c^+(N_s, N_c) = k_c^0 4πρ L_0 R^2 (1+cosθ),

The ratio of addition and removal rates is completely specified by the constraint of detailed balance.
Defining the shell monomer and cargo monomer free energy differences

    Δ F_s   = F(N_s, N_c) - F(N_s-1, N_c), 
    Δ F_c   = F(N_s, N_c) - F(N_s, N_c-1),

we require that the rates of monomer removal satisfy

    k_s^-(N_s, N_c) = k_s^+(N_s-1,N_c) e^βΔ F_s,

and

    k_c^-(N_s, N_c) = k_s^+(N_s,N_c-1) e^βΔ F_c.


Stochastic trajectories for this minimalist description were generated
by advancing the populations N_s and N_c in time with a kinetic
Monte Carlo algorithm based on the transition rates k_s^+, k_c^+,
k_s^-, and k_c^-. By construction, these rates are consistent with
the equilibrium statistics implied by (<ref>), but the geometry
dependence of
addition and removal rates
ensures that typical trajectories do not simply
follow a gradient descent on the free energy surface F(R,θ).
The sampled paths shown in Fig. <ref> (A) (initiated with
N_s=1 and N_c=10, and terminated when N_s≥ 2 N_c^2/3 as
described in Sec. S4) indeed exhibit strong geometric biases.

Most notably, θ increases rapidly at late times in these
successful assembly trajectories, despite equilibrium forces that
encourage droplet growth at the expense of reduced shell coverage.
This route reflects an imbalance in the base addition rates
k_s^0 and k_s^0 of cargo and shell material,
which suppresses droplet growth while the shell develops.  If the
addition of shell material does not outpace that of cargo, then the
angle θ decreases in time—the microcompartment grows but
never advances towards closure, Figs. S5 and S6.  Since the addition
rate of cargo scales with the droplet's surface area, while the rate
of adding shell material scales only with the shell's perimeter,
successful encapsulation generally requires that the imbalance in
addition rates be substantial, which could be straightforwardly
achieved through a large excess of shell monomers in solution relative
to cargo.  When the shell is able to envelop the cargo with only
modest changes in droplet size, the distribution of sizes can be quite
narrow, as shown in Fig. <ref> (B).






§ DISCUSSION

The minimalist model defined by Eqs. 7-13 explains the general course
of our more detailed molecular assembly trajectories.  But its
assumption of ideally spherical geometry is a rough approximation, as
evidenced by asymmetric droplet shapes in the molecular trajectory of
Fig. 1.  Fluctuations in cargo droplet shape are discouraged by its
surface tension, but they are inevitable on small length scales.  In
extreme cases they lead to highly aspherical microcompartments akin to
distended structures that have been observed in
experiments <cit.>.  More common shape fluctuations have
mechanistic influences that are more subtle but nonetheless
important.  These act to reduce local curvature and thus facilitate
growth of a nearly flat shell patch over considerable areas even in
the early stages of assembly.  Enveloping the entire cargo globule
still awaits sufficient cargo growth to reduce the average curvature
below a threshold of roughly 1/R^*.  The necessary curvature that
ultimately develops in the shell is also spatially heterogeneous.
Protein sheets, which generally resist stretching more strongly than
bending, tend to concentrate curvature at topological defects and
along lines connecting them <cit.>.  As a kinetic
consequence, bending is advanced primarily by discrete, localized
events, in particular the formation of fivefold defects during shell
growth.  Indeed, the dynamics of encapsulation are noticeably
punctuated by the appearance of topological defects, as seen in Movie
S1.

Systems out of equilibrium often support rich spatial and temporal patterns, but in many cases they are highly variable and challenging to harness or control.  
By contrast, the encapsulation process we have
described is quite reliable. 
Fig. <ref> (A) shows closed structures resulting from assembly
trajectories of our molecular model.  
They are not perfectly icosahedral (Figs. S7 (A) and S7 (B)), echoing the discernibly nonideal geometries observed in electron micrographs of carboxysomes <cit.>.  
Our assemblies are nonetheless roughly isotropic and distinctly faceted, as quantified by the distributions of sphericity and defect angles shown in Fig. <ref> (B) and Fig. S7.  
Gently curved regions do sometimes appear, as they do in the laboratory.
In our model this smoothness reflects the presence of paired defects (one fivefold and one sevenfold) that sometimes form in the process of shell closure, screening the elastic interactions that give rise to faceting <cit.>.

In published <cit.> and concurrent work, Hagan et
al. have computationally examined assembly dynamics of
microcompartments stabilized by spontaneous curvature of the protein
shell.  Experimental data suggest that key shell proteins of the
carboxysome tend, by contrast, to bind in an unbent geometry. But
existing measurements do not definitively exclude a role for
mechanically preferred curvature in carboxysome assembly. Quantitative
characterization of the elastic properties of shell protein sheets in
the absence of cargo would help clarify this issue.  Systematic
measurement of microcompartment geometry as a function of protein
concentrations would provide a more stringent test of our
nonequilibrium mechanistic hypothesis, which can be achieved with reconstitution experiments.

The nonequilibrium assembly mechanism revealed by our model has a
particularly spare set of essential requirements: One needs only a
monomer that assembles planar, elastic sheets and has a directional
affinity for cargo that is prone to aggregation.  Tuning the
concentrations of these components provides a direct route to
generating monodisperse, nanoscale, self-organizing structures.  The
simplicity of the process presents opportunities to investigate the
dynamics of shell assembly in highly controlled experimental settings,
for example using colloidal particles. Such a platform has already
been demonstrated for a related example of kinetically defined
microsctructure assembly, namely bicontinuous gels stabilized by
colloidal interfaces <cit.>.  Our results suggest that 
kinetic control can be used for precise, nanoscale assembly in a surprisingly general context. 



We thank Michael Hagan, Christoph Dellago, David Savage, Laura
Armstrong, Avi Flamholz, and Luke Oltrogge for helpful conversations.
GMR acknowledges support from National Science Foundation Graduate
Research Fellowship and the James S. McDonnell Foundation. 
PLG acknowledges support from the Erwin Schrödinger International
Institute for Mathematics and Physics. This work
was supported by the US Department of Energy, Office of Basic Energy
Sciences, through the Chemical Sciences Division (CSD) of the Lawrence
Berkeley National Laboratory (LBNL), under Contract
DE-AC02-05CH11231(PLG).













§ SUPPLEMENTARY MATERIALS











§ MONTE CARLO DYNAMICS OF THE MOLECULAR MODEL



The Monte Carlo dynamics of our molecular model for
microcompartment assembly involves several kinds of trial moves, each
constructed to be reversible and to preserve the grand canonical
probability distribution p(X,N_s, {σ}) in Eq. [5]. Here we detail these
moves and their associated acceptance probabilities. For the cargo we
require only a cargo addition/removal move.  For the shell we
employ: vertex displacement moves, simple monomer
  insertion/deletion moves, wedge insertion/deletion moves, and
vertex fusion/fission moves.



 §.§ Vertex distribution function


Representing the protein shell as a triangulated sheet introduces a
subtlety to the grand canonical distribution function p(X,N).
Triangles in the sheet each correspond to a molecular monomer with
3× 3 degrees of freedom, corresponding to the positions of each
of its constituent vertices. Physically, a pair of such monomers would
possess 2× 3× 3 degrees of freedom even when
bound. Represented as a “sheet” comprising two juxtaposed triangles
(with a total of four vertices), however, a bound dimer possesses only
4× 3 degrees of freedom. The implied constraint mimics strong
and short-ranged forces of cohesion between bound monomers, which
limit the fluctuations of a bound vertex to a very small volume v_a
about its binding partner. For a favorable binding energy
ϵ_ bind that is constant within this volume, the binding
equilibrium for a pair of monomers is parameterized by the resulting
vertex affinity K=v_a e^-βϵ_ bind  <cit.>.
In setting v_a=0, we implicitly take
ϵ_ bind→ -∞ such that K remains
nonzero and finite.

The triangulated sheet representation does not resolve fluctuations
of bound vertices on the irrelevantly small length scale
v_a^1/3. These fluctuations have effectively been integrated out,
with each bound pair contributing a factor of K to the probability
distribution for sheet vertices:

    p̅(X̅,N_ vertex,{σ}) = K^N_ bindz_s^N_se^-β( H_ el +
         H_ cargo + H_ int )/Ξ,

where N_ bind = 3 N_s - N_ vertex, and X̅ denotes the configuration of the
set of N_ vertex distinct vertices comprising the triangulated
sheet. The activity of shell monomers,

    z_s = e^βμ_s/λ^9

has units of (length)^-9 reflecting the 9 degrees of freedom for
an unbound monomer. As written in Eq. <ref>, λ is an
arbitrary length scale that sets the compensatingly arbitrary zero of
chemical potential μ_s.



 §.§ Metropolis acceptance probability

Our simulations are designed to sample the probability distribution
p̅(X̅,N_ vertex,{σ}). In a Metropolis Monte
Carlo scheme, detailed balance with respect to this distribution is
ensured by accepting proposed moves with probability

    acc(Γ→Γ') = 
      min[ 1, gen(Γ' →Γ) p̅(Γ')/gen(Γ→Γ') p̅(Γ')],

where Γ is shorthand for the variables (X̅,N_
  vertex,{σ}) that define the system's microscopic state,
and gen(Γ→Γ') is the probability of proposing
Γ' as a trial state from the existing state Γ.

Details of our Monte Carlo moves, and the resulting changes in
equilibrium probability p̅(Γ) that result, are described
below. For those that change the number of shell monomers, acceptance
probabilties involve factors of z_s and K in addition to the
energetics described by
H_ shell and
H_ int.  Several of these trial moves have asymmetric
generation probabilities, gen(Γ' →Γ) ≠gen(Γ→Γ'), which must be accounted for as well.





 §.§ Vertex displacement






Elastic relaxation of the shell is accomplished by conventional Monte
Carlo displacement moves. Such a move selects a vertex at random and
attempts to displace it by an amount Δ a that is chosen at random from a symmetric distribution, as depicted in Fig. S2.
The standard acceptance criterion is given in Eq. [4].

For a state Γ with N_ vertex(Γ) shell monomers, a
Monte Carlo sweep is defined as a sequence of N_ vertex(Γ)
displacement moves, each acting on a randomly chosen vertex.



 §.§ Simple monomer insertion/deletion








Our shell growth moves are depicted in Fig. <ref>.
They are typically attempted less than once per sweep, at a rate
proportional to the number N_ perim of edges at the perimeter of
the sheet. Precisely, in each sweep we propose a shell monomer
insertion move with probability N_ perim k_s^0 τ (where k_s^0 τ
is presumed to be much less than 1/N_ perim).
First, one of the N_ perim edges is
selected at random. Each of these edges corresponds to a particular
monomer i.  The type of growth move that is proposed depends on the
local connectivity and geometry of the sheet.



In a simple monomer insertion move (Fig. <ref> (A)), addition of a new shell monomer is
proposed by adding a single vertex that is connected to the two
existing vertices defining the selected edge. Because the vertices of
bound monomers are constrained to coincide exactly along the shared
edge, a new monomer defined by three vertices can be introduced by
adding just one new vertex.
The position of the new vertex is
generated by first identifying an “ideal” insertion site a_
  ideal, lying in the same plane as monomer i, at a distance
√(3)/3 l_0 from the midpoint of the selected edge, and equidistant
from the edge's two endpoints.
The position a of the new vertex is chosen randomly within a small
volume v_add centered on a_ ideal.  The new monomer's outward
normal vector is chosen to have a positive projection onto that of
monomer i.  The generation probability for such a simple monomer
insertion is

    gen(Γ→Γ') = N_ perim(Γ) k_s^0 τ1/v_add
        N_perim(Γ).


The reverse of a simple addition move deletes a shell monomer by
removing a single vertex and its two edges. We first select one of the
N_perim(Γ') edges at the perimeter of the sheet.
 A monomer added by simple insertion can be removed by selecting either
of its two exposed edges.
The generation probability for simple deletion is therefore

    gen(Γ' →Γ) = N_ perim(Γ') k_s^0 τ2/
        N_perim(Γ').


The grand canonical probability of Γ' relative to
the initial configuration is

    p̅(Γ')/p̅(Γ) = z_s K^2
      exp( -βΔ H),

The two factors of K account for constraining positions of two of
the vertices defining the new monomer.  The energy difference Δ H due to changes in
H_ shell+ H_ int
involves
only the elastic energy of the edges of the new monomer and the
interaction with nearby cargo lattice sites.

The Metropolis acceptance criterion for the simple monomer insertion
move is thus

    acc(Γ→Γ') = min[1, 2 z_s K^2
        v_add
        e^-βΔ H].

In the acceptance probability for simple monomer deletion,
acc(Γ'→Γ), the second argument of the min
function is inverted.

Some surface edges are not eligible for simple insertion in our
protocol. Specifically, in many configurations it is possible to add a
shell monomer by filling in a “wedge” at the sheet's perimeter, as
depicted in Fig. <ref> (B).  We do not attempt simple insertion moves at the
edges defining such a wedge. When such an edge is selected, we instead
propose the wedge insertion move described below.




 §.§ Wedge insertion/deletion


In some shell configurations a new monomer can be added without
introducing any new vertices at all. Geometrically, this corresponds
to filling in an empty wedge at the sheet's boundary, as sketched in
Fig <ref> (B).
We detect this possibility by checking
if two disconnected perimeter vertices share a neighbor and are within a distance l<l_max. Since the new
monomer is completely constrained in this case, the corresponding
generation probability is

    gen(Γ→Γ') = N_ perim(Γ) k_s^0 τ2/N_perim(Γ),

whose factor of 2 reflects that the same insertion can result from
selecting either of the two existing edges that define the wedge.
The generation probability for the reverse move is simply

    gen(Γ' →Γ) = N_ perim(Γ') k_s^0 τ1/N_perim(Γ'),


The grand canonical probability of Γ' relative to
the initial configuration is

    p̅(Γ')/p̅(Γ) = z_s K^3
      exp( -βΔ H).

As with the simple insertion move, factors of K account for
constraining positions of vertices defining the new monomer.  The
energy difference Δ H due to changes in H_
  shell+ H_ int again involves only volume exclusion, the
elastic energy of the edges of the new monomer, and the interaction
with nearby cargo lattice sites.

The resulting Metropolis acceptance criterion for wedge insertion
is

    acc(Γ→Γ') = min[1,z_s K^3/2 
        e^-βΔ H].

In the acceptance probability for wedge monomer deletion,
acc(Γ'→Γ), the second argument of the min
function is inverted.



 §.§ Vertex fusion/fission




The insertion moves described above represent events in which shell
monomers arrive from bulk solution and bind to the growing sheet.  A
separate move is required to allow for binding between monomers that
already belong to the sheet.  In the triangulated sheet
description, such an event fuses one or more existing vertices, so we
term them “fusion” moves (and “fission” in reverse).  At each
sweep we attempt a fusion move with probability N_ fuse k_
  fuse^0 τ, where N_ fuse is the number of vertices eligible for fusion
events.  The reverse moves, i.e., fission, are attempted with
probability N_ fiss k_ fuse^0 τ, where N_ fiss
is the number of eligible fission events. 
 In each case one of the eligible moves is selected at random and then proposed as detailed below.

Some fusion moves act to close gaps between monomers that already
share a common vertex, as depicted in Fig. <ref>.  In such a
“stitching” move, two vertices i and j, initially separated by a
distance less than l_ fuse, are merged into a single vertex. The
new vertex is placed
at the midpoint of i and j. The generation probability for a
stitching move is

    gen(Γ→Γ') = N_ fuse(Γ)
      k_ fuse^0 τ1/
        N_ fuse(Γ),

The corresponding change in grand canonical
weight is

    p̅(Γ')/p̅(Γ) = K
      exp( -βΔ H).


The reverse of a stitching move splits one surface vertex i into
two.  One of the new vertices is placed at a position randomly
distributed within a spherical volume v_ fuse= 4π l_
  fuse^3/3 centered on the existing vertex. The other new vertex is
placed opposite the original vertex, such that i coincides with the
midpoint of the new vertices.  The resulting generation probability

    gen(Γ'→Γ) = N_ fiss(Γ)
      k_ fuse^0 τ1/
        v_ fuse
        N_ fiss(Γ),

The acceptance probability for a stitching fusion move is thus

    acc(Γ→Γ') = min[1, K/v_ fuse
        e^-βΔ H].

In the acceptance probability for the corresponding fission move,
acc(Γ'→Γ), the second argument of the min
function is inverted.

For a shell that is sufficiently large and sufficiently curved, two
different parts of the perimeter can meet. When this happens, binding
of the two parts should be quite facile, since it is not necessary to
await arrival of a new monomer from the dilute solution. Physically,
this union event is just another instance of two shell monomers
binding. Algorithmically, it is more complicated, involving a highly
nonlocal change in vertex connectivity.

We permit such a nonlocal connection move only when two edges at the
sheet's perimeter nearly coincide. Denoting the vertices of one edge
as i_1 and i_2, and those of the other edge as j_1 and j_2, we
require that i_1 and j_1 are separated by a distance less than
l_ fuse and that i_2 and j_2 also reside within the
short distance l_ fuse. One new fused vertex is placed at the
midpoint of i_1 and j_1; the other is placed at the midpoint
of i_2 and j_2. The resulting generation probability is

    gen(Γ→Γ') = N_ fuse(Γ)
      k_ fuse^0 τ1/
        N_ fuse(Γ),

and the change in statistical weight is 

    p̅(Γ')/p̅(Γ) = K^2
      exp( -βΔ H).

In a nonlocal fission move, an edge is split in two. The new vertices
are placed as described for stitching fission, giving a generation
probability

    gen(Γ'→Γ) = N_ fiss(Γ)
      k_ fuse^0 τ1/
        v_ fuse^2
        N_ fiss(Γ).

The Metropolis acceptance probability for nonlocal fusion is finally

    acc(Γ→Γ') = min[1, K^2/v_ fuse^2
        e^-βΔ H].

In the acceptance probability for the corresponding nonlocal fission
move, acc(Γ'→Γ), the second argument of the
min function is inverted.

Counting fission moves is more involved than counting fusion moves:
For the local stitching scenario,
a vertex i on the perimeter is eligible for fission along an
edge ij if the vertex j is not on the perimeter.
In the case of nonlocal fission, a perimeter vertex i is eligible
to be split
if it has a
neighbor j which is also on the perimeter and if removal of the edge
does not lead to two topologically disconnected domains on the shell.
In practice, we test this condition efficiently using a breadth first
search.



 §.§ Cargo addition/removal






For a given configuration {σ} the cargo droplet possesses
N_ boundary sites adjacent to the droplet boundary. This set includes
unoccupied sites with one or more occupied neighbors, where a new
cargo monomer can be accommodated.  It also includes occupied sites
that are incompletely coordinated, where a monomer could be removed
without changing the droplet's connectivity.  During each sweep we
attempt with probability
N_ boundary k_c^0 τ to change the cargo configuration. We
select one of the boundary sites at random and attempt to change its
state, either inserting or deleting a cargo monomer. The Metropolis
acceptance probability for this procedure is

    acc(Γ→Γ') = min[1,
      e^-β (Δ H_ cargo + Δ H_ int)].





 §.§ Topological defects

In the absence of cargo, the ground state of a shell protein sheet
features only six-fold coordinated vertices except at the perimeter.
In order to close, however, the shell must develop either vacancies or
else undercoordinated vertices.  Five-fold coordinate sites can arise
as a consequence of insertion or fusion moves, and they play an
important role in encapsulation.

The Monte Carlo moves we have described could also produce
overcoordinated vertices. Because there is no evidence for seven-fold
(or higher) symmetric protein assemblies in the carboxysome
literature, we prohibit these defects from forming, with a single
exception.  A seven-fold defect is allowed to result from fusion that
causes two fivefold defects to share a neighbor.  This arrangement can
be viewed as accommodating a vacancy in the monomer sheet, rather than
creating a heptameric arrangement of shell monomers.  The large
favorable energy for such a fusion event makes the distinction between
vacancy and excess coordination thermodynamically unimportant.



 §.§ Constraints and parameterization

The constraints on connectivity we impose in these dynamics – there
is exactly one connected cargo droplet and exactly one connected sheet
– are artificial.  This is done as a matter of convenience, with
substantial computational savings.  The constraint should not
significantly affect assembly dynamics provided that nucleation of
cargo and shell are both slow compared to their growth.  The physical
conditions we examine fall in this regime.

The set of parameters governing our molecular simulations are listed
in Tables <ref> and <ref>.
One key feature is that elastic relaxation is very fast
relative to monomer addition, i.e., k_s^0 τ≪ 1 and k_c^0 τ≪ 1. Another is that shell growth is highly unfavorable in the
absence of cargo, as guaranteed by small values of z_s K^2 v_
  add and z_s K^3. The potent attraction between shell and cargo
makes simple shell monomer addition feasible but reversible when cargo
is present on the interior side of the new monomer. Wedge insertion
next to cargo is even more favorable and is typically accepted with
high probability.

In some cases, two disconnected edges can come together so as to
enclose a triangular domain, but one which was not explicitly occupied
by a monomer before the fusion event. In the triangulated sheet
representation, the newly enclosed domain is assumed to be occupied by
a shell monomer. Such a move therefore implicitly entails the arrival
and binding of a new monomer. Because this binding is strongly downhill in
free energy for the physical conditions of interest, we treat it
as inevitable and irreversible, and we do not represent it explicitly.



§ PHENOMENOLOGICAL FREE ENERGY




 §.§ Macroscopic equilibrium


In order to assess the equilibrium arrangements of cargo and shell
species, we view their assembly as a chemical transformation, in which N_s shell capsomers S combine with N_c cargo subunits C to form a complex A_N_s,N_c, which might or might not be a completely closed microcompartment.
This process can be represented as a simple kinetic scheme

    N_s S + N_c C ⇌ A_N_s,N_c.


At thermodynamic equilibrium, stability criteria require that the chemical potentials μ associated with the constituent materials are balanced.

    N_sμ_s + N_cμ_c  = μ_N_s,N_c.

Assuming these particles are present at very dilute concentration, the chemical potential of each species j can be written as a function of its density, the temperature T, and a density-independent standard state chemical potential μ^(0),

    μ_j(T) = μ_j^(0)(T) + T lnρ_j/ρ^(0).

Note that Boltzmann's constant has been absorbed into the unit of temperature, i.e., k_B = 1.

Combining the relationships above yields a law of mass action,

    ρ_N_s,N_c = ρ^(0)exp[ - βΔμ_N_s,N_c^(0) - N_s T lnρ_s/ρ^(0) -  N_c T lnρ_c/ρ^(0)].

The equilibrium constant relating densities of unassembled and assembled species is expressed here in terms of the reversible work required to assemble the components into a filled shell at standard state concentration,

    Δμ^(0)_N_s,N_c = μ_N_s,N_c^(0) - N_sμ_s^(0) - N_cμ_c^(0).

The simple relationship (<ref>) among particle densities allows us to investigate the equilibrium distribution of encapsulated structures. 
In particular, we ask at which value (N_s, N_c) the density ρ_N_s,N_c is maximal. 

Equivalently, we can find a minimum of the free energy

    G(N_s,N_c) = Δμ^(0)_N_s,N_c - N_s T lnρ_s/ρ^(0) -  N_c T lnρ_c/ρ^(0).

While an exact expression for G is not analytically tractable for
our molecular model, we can judge from the energy scales involved what
type of minima might exist.  In particular, it is clear that the
shell-shell and cargo-shell interactions have an energy scale
proportional to N_s.
while the cargo-cargo interactions have an energy
proportional to N_c. Using these intuitive scaling arguments, we can
postulate a simple expression for the free energy

    G(N_s,N_c) =  a N_s + b N_c + E_el(N_s,N_c),

 where the constants a and b are unknown functions of the densities, temperature, and binding free energies.
 
 The elastic energy, in the case of an idealized spherical geometry, does not scale with the size of the shell; it is a constant E_el(N_s, N_c) = 4πκ.
 In the more general setting of nonspherical geometries, we can approximate the free energy as
 
    G(N_s,N_c) =  a ζ N_c^2/3 + b N_c + const. + E_el(ζ),

       where ζ≡
     N_sN_c^-2/3 is a parameter that quantifies the asphericity of the
     structure.  Note that ζ is not the same quantity used to
     characterize the sphericity in Sec. S7.  A detailed analysis of
     the extrema of the free energy functional gives insight into the
     types of shapes that are thermodynamically preferred.  The free
     energy has a minimum value ζ_min=O(1) in which the
     structural is spherical.  Because the elastic energy
     E_el is an increasing function of
     ζ, the minimum
     of the free energy will have ζ = ζ_min for a
     fixed N_c, so long as a>0. As a result, global minima of the
     free energy can only occur at N_c=0 or N_c=∞, depending on
     the sign of b. Mathematically, negative values of a could
     produce other minima, but such a negative surface tension would result
     in highly distended shapes, which are not physically relevant for
     our model.  A similar scaling argument can be made for idealized
     icosahedral structures.  
       This accounting suggests that the metastable
     microcompartment structures seen in biological systems
     represent a departure from equilibrium statistics.



 §.§ Thermodynamic potential for compartment size and shell coverage

A more detailed phenomenological free energy functional can be
constructed by assuming that cargo-shell assemblies adopt the simple
geometry suggested by the analysis above.  Specifically, cargo
molecules are assumed to form a compact spherical globule, and the
shell proteins are assumed to coat the globule as a spherical cap.  At
the level resolved here, icosahedral structures would be governed by a
very similar potential.

It is convenient to write this functional in terms of the radius of the spherical cargo droplet, R, and the angle subtended by the spherical cap θ.
We let f_c denote the free energy per unit cargo and γ_c the surface tension associated with the condensed cargo droplet. 
Similarly, we define f_s and γ_s as the free energy per unit shell and the line tension around the growth edge of the shell. 
The variables ρ and ν denote the density and surface density of the cargo and the shell proteins, respectively.

The shell-cargo interaction energy is -ϵ per monomer. 
The elastic energy is determined by a bending rigidity κ.
Finally, the activity of the cargo monomers is z_c and the activity of the shell monomers is z_s.
Putting together all the terms, we have 

    F(R, θ) = 4/3π R^3 ρ (f_c - T ln z_c) + 2π R^2 (1-cosθ) ν (f_s - Tln z_s - ϵ) + 4π R^2 γ_c + 2π R sinθγ_s + 2πκ (1-cosθ).

The fully assembled microcompartment is stable when Δ F_close(R) ≡ F(R,π) - F(R,0) becomes negative.
If the rate of cargo growth is much slower than the rate of shell growth, the expected size of the assembly is

    R^* = √(κ/|ν(f_s-Tln z_s - ϵ)|)

The average size of an encapsulated structure thus depends strongly on the bending rigidity κ, suggesting one way to tune microcompartment
assembly; controlling monomers' detailed physical properties, however,
is likely a difficult engineering problem in practice.
The shell activity z_s is a much more pragmatic control parameter
for this purpose.


Both the shell and the cargo may need to overcome nucleation barriers before they begin growing spontaneously.
In practice, this means that the expected size of the completely assembled structure will depend on the height and location of these barriers.
The free energy F can be written in a form that highlights the connection to classical nucleation theory as

    F(R, θ) = a_c[ ( R/R_c)^2 - ( R/R_c)^3 ] + a_s[ {(R^*/R_s)^2 - (R/R_s)^2 }(1-cosθ) + R/R_ssinθ],

with 

    a_s   = 2πγ_s^2/ν (Tln z_s+ϵ -f_s)
    
    R_s   = γ_s/ν (Tln z_s + ϵ -f_s)
    
    a_c   = 36 πγ_c^3/ρ^2 (f_c - Tln z_c)^2
    
    R_c   = 3 γ_c/ρ (Tln z_c - f_c)
    
    R^*    = √(2πκ) R_s.

This expression highlights the characteristic sizes of the critical nuclei and their interconnected dependence on the microscopic parameters.
In particular, it should be noted that for the parameters that we are exploring the cargo has essentially no barrier to nucleation, whereas the shell has a significant one at small radius.
However, as the radius increases, the shell can grow with greater facility and encapsulation becomes kinetically accessible.
Fig. <ref> emphasizes that, even for large average radii, the distribution of sizes can be made tight by decreasing the rate of cargo addition.

In simulations of the minimalist model, we adopt a unit of length
ℓ = (2πν)^-1/2, choose a cargo density ρ =
3/(4πℓ^3), and set L_0=ℓ, simplifying many of the above
expressions. In particular, with these choices the monomer populations
can be written transparently in terms of geometric variables as
N_s = R^2 (1-cosθ) and N_c=R^3.

 

 





§ SIMULATIONS OF THE MOLECULAR MODEL


All simulations were performed using an implementation of the algorithm described in Section <ref> which we have open sourced. It is available on GitHub: [https://rotskoff.github.io/shell-assembly/]rotskoff.github.io/shell-assembly/.
The parameters used for all simulations are described in Tables S2 and S3.






All simulations were initialized with a hexamer of shell material which we did not allow to be removed.
This ensures that nucleation of the shell can occur. 
Cargo is initiated by populating all lattice sites within three lattice
spacings of the origin (corresponding to a total number 429 of cargo monomers).

The simulation completes when there are no remaining surface edges on the shell structure. 
At this point, growth moves are no longer proposed and the structure is equilibrated for 10^5 MC sweeps before exiting. 
An animation of a trajectory of the microscopic model is shown in Movie S1.


Simulations were visualized using the VMD software package <cit.>. 
In order to clearly show the shell structure, only occupied lattice sites that are above a distance cutoff of 1.5ł_0 are shown in our visualizations. 
The cargo is depicted either as density isosurfaces or spheres centered at the lattice sites. 
In the Movie S1, the shell is depicted as dynamic bonds with a distance cutoff of 1.5l_0, which leads to unconnected vertices appearing bonded in some frames.



§ FREE ENERGY CALCULATIONS FROM THE MOLECULAR MODEL


We performed a series of umbrella sampling calculations to compare the
thermodynamics of our molecular model with the minimalist model.
First, we considered a small cargo droplet N_ c=135,
corresponding to a sphere with diameter 4 l_0. As shown in
Fig. <ref>, the cost of deforming the shell makes
encapsulation uphill in free energy.  In the case of a droplet with
N_ c=429, a sphere of diameter 6 l_0, which is above the
nucleation barrier for cargo growth, we find that there is a modest
barrier to the initial growth which becomes downhill as N_ s
increases.  Finally, we calculated the free energy for N_
  c=1055, corresponding to sphere with diameter 10 l_0. In
addition to the free energy calculations performed for shell growth,
we also performed a free energy calculation to show that cargo growth
is downhill for the regime of interest.

All umbrella sampling calculations for N_ s were performed with a harmonic bias energies with force constant k=1.0, and centers separated by two monomers. 
The cargo umbrella sampling was performed with k=1.0 and windows separated by 5 cargo monomers. 
 






§ DYNAMICAL FLOWS FOR THE MINIMALIST MODEL


 

In very simple models, the routes of typical trajectories can be
inferred directly from the underlying free energy surface.  Kinetic
inference is not so straightforward for our minimalist model.
Dynamical flows are complicated by geometric factors,
e.g., the fact that cargo addition rates are proportional to the
number of available surface sites, which changes during the course of
assembly.  The free energy F(R,θ) in this case constrains
ratios of rates for each forward and reverse move, but does not give a
good sense for the directions along which the system typically
evolves.

In this section we develop a means to visualize dynamical tendencies for cargo droplet growth and shell envelopment in our minimalist model. 
Let ω(θ,R) be a vector field with components ω_R=Ṙ and ω_θ=θ̇ along R and θ that reflect the typical rate of change along these coordinates at each point in the θ,R-plane. 
First, we differentiate the relationships

    R(N_c)    = (3/4π(L_0 ρ)^-1 N_c)^1/3
    θ(N_s,N_c)    = cos^-1(1-N_sl_0/2 πν R(N_c)^2),

to express the time derivatives Ṙ and θ̇ in terms of R, θ, Ṅ_s and Ṅ_c, where dots represent differentiation with respect to time.
We then estimate average rates of change in N_s and N_c through the effective drift in these populations,

    Ṅ_s   = -k_s^+(N_s, N_c) GN_s
    Ṅ_c   = -k_c^+(N_s, N_c) GN_c,

which can in turn be written in terms of R and θ. 
Plotting ω as a function of R and θ yields a dynamical flow field that reveals the direction in which trajectories typically evolve from each point (θ,R).

The dynamical flow fields shown in Fig. <ref> highlight the importance of monomer arrival rates, particularly for the likelihood of shell closure. 
Panels (A) and (B) correspond to systems with the same free energy F(R,θ), but with different relative arrival rates, k_c^0, for cargo.  

For very fast cargo arrival, the nucleation of sufficient shell material to encapsulate the rapidly expanding surface of the cargo is kinetically impeded. 
The field shown in Fig. <ref> (A) corresponds to this case of “runaway” growth (see also, Fig. <ref>), with ω pointing primarily in the direction of increasing R throughout much of (θ,R) space.
In a typical trajectory, the radius of the cargo increases rapidly and the angle θ cannot grow to fully encapsulate the surface.


By contrast, when the shell arrival rate greatly exceeds that of cargo, the field points more strongly along increasing θ, promoting shell closure. 
As soon as the barrier to nucleating a shell disappears, the average dynamical flow leads trajectories rapidly and unequivocally towards encapsulation, as shown in Fig. <ref> (B). 
This representation of dynamical biases does not capture the possibility of barrier crossing, i.e., the spontaneous nucleation of a substantial shell even when it requires activation. 
In these systems, however, this nucleation barrier becomes comparable to k_BT, and therefore easily navigated, not far from the point where shell growth becomes thermodynamically downhill.

The location of the barrier θ^(R), for a given droplet
radius R, can be determined by finding the maximum of
F(R,θ)θ in the range [0,π]. Using the resulting
expression,

    θ^(R) = tan^-1( R/R_s/(R/R_s)^2 - ( R^*/R_s)^2),

we can estimate the radius at which the barrier to shell growth
becomes comparable to the thermal energy.  We judge the height of the
barrier,

    Δ F^(R) = F(R,θ^) - F(R,θ_0),

with reference to a typical value θ_0=0.2 of the spherical cap
angle prior to nucleation of a steadily growing shell. This activation
free energy becomes comparable to thermal energy at a radius R^
for which

    Δ F^(R^) = k_B T.

The critical radius R^ compares favorably with the point at
which shell nucleation occurs in kinetic Monte Carlo simulations, as shown in Fig. 2 (A). 



§ QUANTIFYING SIZE AND SHAPE DISTRIBUTIONS


 
 
 
  



In order to quantify the extent of shape fluctuations, we computed both a measure of sphericity and a measure of the ideality of the icosahedral shape. 
We measure the sphericity by quantifying deviations from the mean radius.
Let 𝐫_0 denote the center of mass of the vertices. 
We denote by d_i magnitude of the vector 𝐫_i which connects 𝐫_0 to a vertex indexed by i.
Then the sphericity is calculated as

    α(𝒱) = 1 - ∑_i∈𝒱( d_i^2 - d_i^2 ) / d_i.

 The number α is 1 on a perfect sphere.
 As the variance of the distance to shell grows, the value of the order parameter decreases.
 This parameter is not
 strongly influenced by the presence of sharp edges and vertices, but instead characterizes the extent to which a structure is isotropic rather than distended or elongated.

 
Bond orientational order parameters (cf. Ref.<cit.>) provide a  very sensitive measure of icosahedral symmetry.
We computed the distribution of two such parameters that distinguish icosahedral orientational order from face-centered cubic, body-centered cubic, hexagonal close packed, and simple cubic packings. 
Following Steinhardt et al., we define

    Q_6({r_1,…, r_N_defect}) = 1/N_defect∑_i=1^N_defect∑_l=-m^m Y_6m (r_i)

and 

    Ŵ_10 = ∑_m_1, m_2, m_3 = -l, …, l( l   l   l
    m_1   m_2   m_3) Q_l,m_1 Q_l,m_2 Q_l,m_3/( ∑_m=-l^l Q_lm)^3/2,

where the coefficient is the Wigner 3-j symbol.
For an ideal icosahedral structure

    Q_6 ≈ 0.66    Ŵ_10≈ -0.094

Previous studies have classified icosahedrality using a minimum of 0.6 for Q_6 or a maximum of -0.1 for Ŵ_10.
By these metrics, our typical structures deviate strongly from perfect icosahedra. 
The distribution of Ŵ_10 suggests that some structures
reasonably approximate
ideal icosahedra, but those structures are rare, in the tail of the distribution of Fig. <ref>.
In contrast, the Q_6 order parameter does not register icosahedrality in our structures, as seen in Fig. <ref>.

Though the structures are clearly not ideal icosahedra, they have the distinct faceted appearance of the structures typically seen in tomograms of carboxysomes and other bacterial microcompartments.
To quantify the extent to which the structure is faceted, it is useful to look at the extent of bending in the proximity of fivefold defects. 
We use the average bending angle of the five bonds around a defect as a faceting measure.
Let 𝒱_d denote the set fivefold defect vertices in a structure, of which there are N_d. 
We also define ℬ(v_i) to be the set of bonds connected to vertex v_i which has elements ij.
Then, the faceting parameter is

    θ_defect = 1/5 N_d∑_v_i∈𝒱_d∑_ij ∈ℬ(v_i) | θ_ij |.


We quantified the size of assembled microcompartments according to their “effective” radii.
We defined the effective radius as

    R_eff = (3 N_c/4π)^1/3,

which assumes an idealized spherical droplet.  While this measure of
the radius is only approximate, it readily provides insight into the
degree of size heterogeneity.  A distribution of assembled compartment
size is shown in Fig. <ref>, together with experimental data
reported in Ref. <cit.>.  It is
challenging to make a quantitative comparison in this case because the
number of experimental observations is small.





