We computationally study how an icosahedral shell assembles around hundreds of molecules. Such a process occurs during the formation of the carboxysome, a bacterial microcompartment that assembles around many copies of the enzymes ribulose 1,5-bisphosphate carboxylase/ oxygenase and carbonic anhydrase to facilitate carbon fixation in cyanobacteria. Our simulations identify two classes of assembly pathways leading to encapsulation of many-molecule cargoes. In one, shell assembly proceeds concomitantly with cargo condensation. In the other, the cargo first forms a dense globule; then, shell proteins assemble around and bud from the condensed cargo complex. Although the model is simplified, the simulations predict intermediates and closure mechanisms not accessible in experiments, and show how assembly can be tuned between these two pathways by modulating protein interactions. In addition to elucidating assembly pathways and critical control parameters for microcompartment assembly, our results may guide the reengineering of viruses as nanoreactors that self-assemble around their reactants.


eLife digest

Bacterial microcompartments are protein shells that are found inside bacteria and enclose enzymes and other chemicals required for certain biological reactions. For example, the carboxysome is a type of microcompartment that enables the bacteria to convert the products of photosynthesis into sugars. During the formation of a microcompartment, the outer protein shell assembles around hundreds of enzymes and chemicals. This formation process is tightly controlled and involves multiple interactions between the shell proteins and the cargo – the enzymes and other reaction ingredients – they will enclose. Understanding how to control which enzymes are encapsulated within microcompartments could help researchers to re-engineer the microcompartments so that they contain drugs or other useful products.

Recent studies have used microscopy to visualize how microcompartments are assembled. However, most of the intermediate structures that form during assembly are too small and short-lived to be seen. It has therefore not been possible to explore in detail how shell proteins collect the necessary cargo and then assemble into an ordered shell with the cargo on the inside. Experiments alone are probably not enough to understand the process, especially since microcompartment assembly can currently only be studied within live cells or cellular extract. Within these complex environments it is difficult to determine the effect of any individual factor on the overall assembly process.

Perlmutter, Mohajerani and Hagan have now taken a different approach by developing computational and theoretical models to explore how microcompartments assemble. Computer simulations showed that microcompartments could assemble by two pathways. In one pathway, the protein shell and cargo coalesce at the same time. In the other pathway, the cargo molecules first assemble into a large disordered complex, with the shell proteins attached on the outside. The shell proteins then assemble, carving out a piece of the cargo complex. The simulations showed that many factors affect how the shell assembles, such as the strengths of the interactions between the shell proteins and the cargo. They also identified a factor that controls how much cargo ends up inside the assembled shell.

Perlmutter, Mohajerani and Hagan found that, in addition to revealing how microcompartments may assemble within their natural setting, the simulations provided guidance on how to re-engineer microcompartments to assemble around other components. This would enable researchers to create customizable compartments that self-assemble within bacteria or other host organisms, for example to carry out carbon fixation or make biofuels.

A future challenge will be to investigate other aspects of microcompartment assembly, such as the factors that control the size of these compartments.



Encapsulation is a hallmark of biology. A cell must co-localize high concentrations of enzymes and reactants to perform the reactions that sustain life, and it must safely store genetic material to ensure long-term viability. While lipid-based organelles primarily fulfill these functions in eukaryotes, self-assembling protein shells take the lead in simpler organisms. For example, viruses surround their genomes with a protein capsid, while bacteria use large icosahedral shells known as bacterial microcompartments (BMCs) to sequester the enzymes and reactions responsible for particular metabolic pathways (Kerfeld et al., 2010; Axen et al., 2014; Shively et al., 1998; Bobik et al., 1999; Erbilgin et al., 2014; Petit et al., 2013; Price and Badger, 1991; Shively et al., 1973; Shively et al., 1973; Kerfeld and Erbilgin, 2015). Within diverse bacteria, BMC functions have been linked to bacterial growth, carbon fixation, symbiosis, or pathogenesis (Kerfeld and Erbilgin, 2015). Other protein-based compartments are found in bacteria and archea (e.g. encapsulins (Sutter et al., 2008) and gas vesicles (Pfeifer, 2012; Sutter et al., 2008)) and even eukaryotes (e.g. vault particles (Kickhoefer et al., 1998)), while some viruses may assemble around lipidic globules (Lindenbach and Rice, 2013; Faustino et al., 2014). Thus, understanding the factors that control microcompartment assembly and encapsulation is a central question in modern cell biology. From the perspectives of synthetic biology and nanoscience, there is great interest in reengineering BMCs or viruses as nanoreactors that spontaneously encapsulate enzymes and reagents in vitro (e.g. Luque et al., 2014; Douglas and Young, 1998; Rurup et al., 2014; Patterson et al., 2014; Patterson et al., 2012; Zhu et al., 2014; Rhee et al., 2011; Rurup et al., 2014; Wörsdörfer et al., 2012; Comas-Garcia et al., 2014), or as customizable organelles that assemble around a programmable set of core enzymes in vivo, introducing capabilities such as carbon fixation or biofuel production into bacteria or other organisms (e.g. Kerfeld and Erbilgin, 2015; Bonacci et al., 2012; Parsons et al., 2010; Choudhary et al., 2012; Lassila et al., 2014). However, the principles controlling such co-assembly processes have yet to be established, and it is not clear how to design systems to maximize encapsulation.

In this article we take a step toward this goal, by developing theoretical and computational models that describe the dynamical encapsulation of hundreds of cargo molecules by self-assembling icosahedral shells. Although our models are general, we are motivated by recent experiments on a type of BMC known as the carboxysome (Kerfeld et al., 2010; Schmid et al., 2006; Iancu et al., 2007; Tanaka et al., 2008). Carboxysomes are large (40–400 nm), roughly icosahedral shells that encapsulate a dense complex of the enzyme ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO) and other proteins to facilitate the Calvin-Bensen-Bassham cycle in autotrophic bacteria (Price and Badger, 1991; Shively et al., 1973; Shively et al., 1973; Iancu et al., 2007; 2010; Kerfeld et al., 2010; Tanaka et al., 2008). Recently, striking microscopy experiments visualized β-carboxysome shells assembling on and budding from procarboxysomes (the condensed complex of RuBisCO and other proteins found in the interior of carboxysomes) (Cameron et al., 2013; Chen et al., 2013). Genomic analysis suggests that many BMCs with diverse functions assemble via similar pathways (Cameron et al., 2013; Kerfeld and Erbilgin, 2015). However, the mechanisms of budding and pinch-off to close the shell remain incompletely understood because of the small size and transient nature of assembly intermediates. Moreover, experiments suggest that α-carboxysomes (another form of carboxysome) assemble by a different mechanism, in which shell assembly encapsulates an initially diffuse pool of RuBisCO (Iancu et al., 2010; Cai et al., 2015). The factors determining which of these assembly pathways occurs are unknown.

BMC assembly is driven by a complex interplay of interactions among the proteins forming the external shell and the interior cargo. It is difficult, with experiments alone, to parse these interactions for those mechanisms and factors that critically influence assembly pathways, especially due to the lack of an in vitro assembly system. Models which can correlate individual factors to their effect on assembly are therefore an important complement to experiments.

Previous experimental and theoretical studies of encapsulation by icosahedral shells, e.g. the assembly of viral capsids around their nucleic acid genomes (e.g. Hu and Shklovskii, 2007; Kivenson and Hagan, 2010; Elrad and Hagan, 2010; Perlmutter et al., 2013; 2014; Mahalik and Muthukumar, 2012; Zhang et al., 2013; Zhang and Linse, 2013; Hagan, 2008; Devkota et al., 2009; Dixit et al., 2006; Borodavka et al., 2012; Dykeman et al., 2013; 2014; Zlotnick et al., 2013; Johnson et al., 2004; Patel et al., 2015; Cadena-Nava et al., 2012; Comas-Garcia et al., 2012; 2014; Garmann et al., 2014a; 2014b; Malyutin and Dragnea, 2013), have demonstrated that the structure of the cargo can strongly influence assembly pathways and products. However, BMCs assemble around a cargo which is topologically different from a nucleic acid — a fluid complex comprising many, noncovalently linked molecules. We demonstrate here that changing the cargo topology leads to new assembly pathways and different critical control parameters.

We present phase diagrams and analysis of dynamical simulation trajectories showing how the thermodynamics, assembly pathways, and emergent structures depend on the interactions among shell proteins and cargo molecules. Within distinct parameter ranges, we observe two classes of assembly pathways, which resemble those suggested for respectively α- or β-carboxysomes. We find that tunability of cargo loading is a key functional difference between the two classes of pathways. Shells assembled around a diffuse cargo can be varied from empty (containing almost no cargo) to completely full, whereas assembly around a condensed, procarboxysome-like complex invariably produces full shells. While we find that the encapsulated cargo becomes ordered due to confinement, complete crystalline order in the globule before encapsulation inhibits budding. We discuss these results in the context of recent observations on carboxysome assembly, and their implications for engineering BMCs, viruses or drug delivery vehicles that assemble around a fluid cargo (e.g. Refs. [Kerfeld and Erbilgin, 2015; Parsons et al., 2010; Choudhary et al., 2012; Lassila et al., 2014; Luque et al., 2014; Douglas and Young, 1998; Rurup et al., 2014; Patterson et al., 2014; Patterson et al., 2012; Zhu et al., 2014; Rhee et al., 2011; Rurup et al., 2014; Wörsdörfer et al., 2012]).


Our model system is motivated by icosahedral viral capsids and BMCs (Tanaka et al., 2008; Kerfeld et al., 2010). Since icosahedral symmetry can accommodate at most 60 identical subunits, formation of large icosahedral structures requires subunits to assemble into different local environments. The subunits can be grouped into pentamers and hexamers, with 12 pentamers at the icosahedron vertices and the remaining subunits in hexamers. Viruses typically assemble from small oligomers of the capsid protein, which we refer to as the basic assembly unit (Hagan, 2014). Recent AFM experiments demonstrated that hexamers are the basic assembly unit during the assembly of BMC shell facets (Sutter et al., 2016), and the carboxysome major shell proteins crystallize as pentamers and hexamers (Tanaka et al., 2008). Motivated by these observations, our model considers two basic assembly units, one a pentamer and the other a hexamer, with interactions designed so that the lowest energy structure corresponds to a truncated icosahedron with 12 pentamers and 20 hexamers (Figure 1). While BMCs generally have more hexamers, our model is intended to explore the general principles of assembly around a fluid cargo rather than model a specific system. Further details of the model and a thermodynamic analysis are given in section 3 and the appendices.

Description of the model.

(A) Each shell subunit contains ‘Attractors’ (green circles) on the perimeter, a ‘Top’ (tan circle, ‘T’) in the center above the plane, and a ‘Bottom’ (purple circle, ‘B’ below the plane). (B) Interactions between complementary Attractors drive subunit dimerization, with the Top-Top repulsions (tan arrow) tuned to favor the subunit-subunit angle in a complete shell. Complementary pairs of attractors are indicated by green arrows in (A) for the pentamer-hexamer interface and in (B) for the hexamer-hexamer interface. (C) Bottom psuedoatoms bind cargo molecules (terra cotta circles, ‘C’), while excluder atoms (blue and brown pseudoatoms in (D)) placed in the plane of the pentagon experience excluded volume interactions with the cargo. (D) The positions of excluder atoms in the lowest energy shell geometry, a truncated icosahedron with 12 pentamers (blue) and 20 hexamers (brown).


To understand how assembly around multiple cargo molecules depends on the relative strengths of interactions between components, we performed dynamical simulations as a function of the parameters controlling shell subunit-subunit (εSS), shell subunit-cargo (εSC), and cargo-cargo (εCC) interaction strengths. All energy values are given in units of the thermal energy, kBT. We focus on parameters for which shell subunit-subunit interactions are too weak to drive assembly in the absence of cargo (εSS4.5). Except where mentioned otherwise, the cargo diameter is set equal to the circumradius of a shell subunit.

For the simulated density of cargo particles, the phase behavior (in the absence of shells) corresponds to a vapor at εCC=1.3, liquid-vapor phase coexistence for εCC[1.6,2.0] (the phase coexistence boundary is slightly below εCC=1.6), and a solid phase at εCC=3.0. We find that tuning εCC through phase coexistence dramatically alters the typical assembly process. Strong cargo interactions (εCC1.6) drive formation of a globule followed by assembly and budding of a shell, such as observed for β-carboxysomes (Figure 2A, Simulation Video 1), while under weak interactions (εCC<1.6) shell assembly usually proceeds in concert with cargo encapsulation (Figure 2B, Simulation Video 2), as suggested for assembly of α-carboxysomes. We now elaborate on these classes of assembly pathways, and how the resulting assembly products depend on parameter values.

Figure 2 with 2 supplements see all
Snapshots illustrating typical assembly trajectories.

(A) Multi-step assembly involving an amorphous globule of cargo and shell subunits. (B) Single-step assembly, in which shell assembly drives local cargo condensation. and (C) when shell-cargo interactions are too weak to condense the cargo. The values of the cargo-cargo (εCC), shell subunit-cargo (εSC), and subunit-subunit (εSS) interaction strengths are listed above each panel (all energies are in units of the thermal energy kBT), and the time (in units of 106 timesteps) is noted below each image. The color scheme here and throughout the manuscript is: Red=Cargo, Blue=Pentagon Excluder, Brown=Hexagon Excluder. Attractor and Bottom pseudoatoms are omitted to aid visibility. Videos of assembly trajectories are included below.

Video 1
Animation of a typical simulation showing assembly around a cargo globule.

Parameters are εCC=1.6εSC=7, and εSS=2.5.

Video 2
Animation of a typical simulation showing simultaneous assembly and cargo condensation.

Parameters are εCC=1.3εSC=9, and εSS=3.5.


Assembly and budding from a cargo globule

We begin by discussing assembly behavior when the cargo-cargo interactions are strong enough to drive equilibrium phase coexistence (εCC1.6). Near the phase boundary (εCC=1.6) a system of pure cargo particles is metastable on the timescales we simulate. However, for εSC>4, adding shell subunits drives nucleation of a cargo globule with shell subunits adsorbed on the surface. The subsequent fate of the globule depends on parameter values; typical simulation end-states are shown as a function of parameter values in Figure 3. For moderate interaction strengths (2.5εSS3.5) the globule grows to a large size, typically containing at least twice the cargo molecules that can be packaged within a complete shell. Adsorbed shell subunits then reversibly associate to form ordered clusters. Once a cluster acquires enough inter-subunit interactions to be a stable nucleus, it grows by coagulation of additional subunits or other adsorbed clusters. For the parameter set corresponding to Figure 2A, nucleation is fast in comparison to cluster growth, and thus two nuclei grow simultaneously. The last three images show the system immediately preceding and following detachment of the lower shell. Missing only one of its 32 subunits, the shell is connected to the remainder of the droplet only by a narrow neck of cargo. Insertion of the final subunit breaks the neck and completes shell detachment. The complete shell contains 120–130 cargo particles, which is slighty above random close packing (120 particles) but below fcc density (150 particles, see appendix 1.2).

Figure 3 with 3 supplements see all
Results of assembly around a cargo globule.

(A) The most frequently observed assembly outcome is overlaid on a color map of the theoretical free energy density difference Δfassem (Equation (3)) between assembled shells and the unassembled globule. Results are plotted against the shell-cargo adsorption strength εSC and the shell-shell interaction strength εSS for indicated values of the cargo-cargo interaction strength εCC. (B) Representative snapshots of the predominant assembly outcomes shown in (A).

Figure 3—source data 1

List of all simulation outcomes for Figures 3A,5A.

Figure 3—source data 2

Criteria used to categorize assembly outcomes.

The sizes of each cargo globule and shell assemblage, and associations between shell assemblages and cargo globules, were determined by clustering. The outcome was then categorized according to the criteria listed in this table.


Increasing the shell-shell interaction strength drives faster shell assembly and closure, thus limiting the size of the globule before budding. For the largest interaction strength we simulated (εSS=4.5) the globule typically does not exceed the size of a single shell, and multiple globules nucleate within the simulation box (Figure 2—figure supplement 1). This observation could place an upper bound on shell-shell interaction strengths, since multiple nucleation events were rare in the carboxysome assembly experiments (Cameron et al., 2013) (however, we discuss potential complicating factors within the cellular environment below). To quantify the relationship between assembly mechanism and parameter values, we calculate an assembly order parameter, defined as the maximum number of unassembled subunits adsorbed onto a globule during an assembly trajectory. The order parameter is shown as a function of the interaction strengths in Figure 4. For εCC1.6 and εSS3 we observe large values of the order parameter (e.g. >32, the red and yellow regions in Figure 4), which indicate formation of a large amorphous globule consisent with the procarboxysome precursor to carboxysome shell assembly (Cameron et al., 2013).

Figure 4 with 1 supplement see all
Dependence of assembly pathway on shell-cargo and shell-shell interaction strength.

The assembly order parameter, defined as the maximum number of unassembled shell subunits adsorbed on a globule at any point during a trajectory, is shown as a function of εSC and εSS for indicated values of the cargo-cargo interaction εCC. Large numbers of adsorbed unassembled subunits (>32) indicate the two step assembly mechanism (Figure 2A), whereas smaller values correspond to simultaneous assembly and cargo condensation (Figure 2B).


Other assembly products

Outside of the optimal parameter ranges, we observe several classes of alternative outcomes. Overly weak shell-shell interactions fail to drive assembly. For εCC=1.6 and εSC4 the cargo vapor phase is metastable, and the system remains ‘Unnucleated’ (with no cargo globule) on simulated timescales (we discuss alternative initial conditions below). Stronger cargo-cargo or shell-cargo interactions result in unassembled ‘Globules’, where a cargo globule forms but the shell subunits on its surface fail to nucleate. As εSS increases, we observe assembly on the globule, leading either to complete shells or two classes of incomplete assembly. In the first incomplete case, ‘Attached’, one or more shells almost reaches completion, but fails to detach from the droplet within simulated timescales. ‘Attached’ configurations occur for low εSC, when the subunit-cargo interaction does not provide a strong enough driving force for the last subunit(s) to penetrate the cargo and close the shell. Overly strong interactions drive the other class of incomplete assembly: ‘Over-nucleated/Malformed’, in which an excess of partially assembled shells deplete the system of free subunits before any shells are completed. In this regime it is also common to observe malformed structures, in which defects become trapped within growing shells.

As the cargo-cargo interaction increases (εCC1.8), multiple effects narrow the parameter range that leads to complete assembly and detachment. Firstly, cargo globules nucleate rapidly at multiple locations within the simulation box, increasing the likelihood of the ‘Over-nucleated’ outcome. Secondly, the threshold value of εSC required for cargo penetration increases, resulting in ‘Attached’ shells over a wider parameter range. We also observe a configuration we refer to as ‘Stalled’, in which shell assembly fails to penetrate the globule surface (and thus does not even proceed to the attached stage). The latter is especially prevalent for εCC=3.0, when the cargo crystallizes even in the absence of shell encapsulation. For both ‘Attached’ and ‘Stalled’ configurations, regardless of the initial number of nucleation events, we typically observe coarsening into a large globule.

Simultaneous shell assembly and cargo condensation

For εCC=1.3 the cargo forms an equilibrium vapor phase in the absence of shell subunits. However, above threshold values of εSS and εSC, the diffuse cargo molecules drive nucleation of shell assembly. The subsequent assembly pathway depends sensitively on the shell-cargo interaction strength. For low εSC (Figure 2C), assembly captures only a few cargo molecules, leading to complete, but nearly empty shells. For larger εSC (Figure 2B, and Simulation Video 2), the shell-cargo interactions drive local condensation of cargo molecules. Shell assembly and cargo complexation then proceed in concert, resembling the mechanism proposed for assembly of α-carboxysomes (Iancu et al., 2010). Thus, tuning the shell-cargo interaction dramatically affects cargo loading, with a sharp transition from empty to filled shells around εSC=2. This transition closely tracks the equilibrium filling fraction (Figure 5C), measured by simulating a complete shell made permeable to cargo molecules. This effect is comparable to the condensation of water vapor below its dew point inside of hydrophilic cavities. In contrast, assembly around a globule only generates full shells.

Figure 5 with 2 supplements see all
Results of assembly around a cargo with weak interactions (εCC=1.3kBT).

(A) The most frequently observed assembly outcome as a function of εSS and εSC. The distribution of outcomes for εSS=4 is shown in Figure 3—figure supplement 2, and a data file containing the outcome for each trial at each parameter set is included (Figure 3—source data 1). (B) Representative snapshots for the outcomes shown in (A). The complete shell outcomes are shown with the excluders rendered opaque (left) and transparent (right) to enable visualizing the encapsulated cargo. (C) The number of cargo molecules encapsulated by shells assembled in dynamics simulations (red symbols) is compared to the results of equilibrium simulations (black line). The dynamics results are averaged over all complete shells (for any εSS) assembled at each value of εSC, the error bars indicate 95% confidence intervals. Most simulations were performed for 3×108 timesteps; simulations with εSS=4.5, εSC4, and εCC=1.3 exhibited partially assembled shells at 3×108 timesteps, and were continued up to 7.2×109 timesteps.


Assembly of full shells (by either pathway, Figure 2A or Figure 2B) is typically about two orders of magnitude faster than assembly of empty shells (Figure 2C). This disparity demonstrates the key role that the cargo plays in promoting shell association, during all stages of assembly. Cargo molecules initially promote shell nucleation by stabilizing interactions among small, sub-nucleated clusters. Then, the presence of a condensed globule provides a large cross-section for adsorption of additional subunits, significantly enhancing the flux of subunits to the partial capsid, thus increasing its growth rate. The condensed cargo particularly facilitates insertion of the last few subunits, which are significantly hindered by steric interactions, as noted previously for simulations of empty virus capsids (Nguyen et al., 2007).

Figure 5A shows how the products of assembly around cargo with weak interactions depends on parameters. While moderate parameter values lead to complete assembly, overly weak εSC and εSS (lower left region of Figure 5A) prevent shell nucleation, leading to the ‘Unnucleated’ outcome. In the limit of large εSC but weak εSS the shell-cargo interaction stabilizes small disordered globules (50 cargo particles, lower right region of Figure 5A), while under strong subunit and weak cargo interactions (εSS=4.5, εSC<5) shells nucleate but cannot condense the cargo, leading to the complete but slow assembly just discussed. As for assembly around a globule, overly strong interactions lead to overnucleation and malformed shells. However, the predominant mode of malformation is now shell collapse. Because the cargo is below its dew point, the locally condensed globule leads to a negative pressure on the shell subunits, which can flatten the shell and thus prevent closure of a symmetric shell.

Thermodynamic model

The simple free energy model (Equations (1–2)) reproduces the threshold parameter values required for shell assembly with no adjustable parameters (color map in Figure 3). Since it is an equilibrium model and only considers the free energy difference between complete and unassembled configurations, it cannot distinguish between parameter values that lead to complete assembly or kinetic traps at the long but finite simulation times. However, the thermodynamic calculation does suggest that the simulations resulting in ‘Attached’ shells would eventually reach completion on a longer timescale. We do not show Δfassem in Figure 5A because the globule is always less favorable than assembled shells for εCC=1.3, but the yield of well-formed shells in our simulations roughly follows the prediction of the equilibrium theory (Figure 5—figure supplement 1).

Effects of varying other parameters or initial conditions

To investigate whether the results described above depend on assumptions within our model, we performed several sets of additional simulations. Firstly, we performed simulations in which the ratio between cargo diameter in shell subunit size was varied. As shown in Figure 5—figure supplement 2, assembly is most robust for our default cargo diameter (for which the model was parameterized), but productive assembly occurs for cargo diameters varied over a factor of four. Secondly, we performed assembly simulations with anisotropic cargo molecules with a shape motivated by the octomer structure of the RuBisCO holoenzyme (Figure 2—figure supplement 2).

Thirdly, we performed a set of simulations in which we pre-equilibrated the cargo globule before introducing shell subunits into the system (Figure 3—figure supplement 2, Simulation Video 3). Investigating this alternative initial condition was motivated by the fact that RuBisCO is present in the cell before induction of the carboxysome gene in the experiments of Ref. (Cameron et al., 2013), and by the observation that multiple carboxysomes bud sequentially in time from a single procarboxysome. For εCC=1.6 the results are very similar to those obtained without pre-equilibrating the cargo. However, for εCC>1.6, successful assembly and detachment is limited to more narrow ranges of shell-shell and shell-cargo interaction strengths than in Figure 3, due to an increased prevalence of ‘Attached’ and ‘Stalled’ configurations. The latter are particularly common for εCC=3, when the cargo forms a hexagonally close packed crystal which strongly resists deformation by shell protein assembly.

Video 3
Animation of a simulation with a pre-equilibrated cargo globule.

Parameters are εCC=1.6, εSC=6, and εSS=3.5.


Taken together, the results from both assembly protocols (Figure 3 and Figure 3—figure supplement 2) suggest that moderate effective cargo-cargo interactions are most consistent with the observations of shell assembly and budding in Refs. (Cameron et al., 2013; Chen et al., 2013). Such interactions are strong enough to drive cargo globule formation, but malleable enough to allow shell assembly to deform and eventually sever intra-globule interactions.

Organization of encapsulated cargo

Studies of assembled carboxysomes report varying degrees of order for the encapsulated cargo, ranging from none to paracrystalline order (Iancu et al., 2007; 2010; Kaneko et al., 2006; Schmid et al., 2006). We therefore studied the relationship between cargo order and interaction parameters using equilibrium simulations (see Figure 6 and Figure 6—figure supplement 1). Below εCC<3kBT, we do not observe true fcc order of the encapsulated cargo. However, for all parameters leading to significant filling, even those well below the cargo liquid-vapor transition, the cargo becomes organized in concentric layers (Figure 6). We observe similar cargo organizations within shells which have budded from cargo globules in dynamical simulations. These results demonstrate that ordering of the cargo does not require crystallinity of the initial globule. Moreover, the magnitude of ordering increases with cargo loading, but, for fixed loading, is essentially independent of the cargo-shell interaction strength εSC. We observe ordering within filled shells due to confinement, even if even if εSC is set to 0 (Figure 6—figure supplement 1), as previously noted by Iancu et al. (Iancu et al., 2007).

Figure 6 with 1 supplement see all
Order of the encapsulated cargo.

The spherically averaged density of cargo molecules inside a shell is shown as a function of radius for (A) εCC=1.6 and (B) εCC=1.3 for indicated values of the cargo-shell adhesion strength εSC, measured in equilibrium simulations. The density of the encapsulated cargo ranges from below random close packing to near hexagonal close packing density as εCC and εSC are increased (see Figure 3—figure supplement 3). A snapshot of cargo inside the shell is shown in Figure 5—figure supplement 2. The raw data for this figure is provided in Figure 6—source data 1.

Table 1

Description of the assembly outcomes presented in Figures 3,5.

Complete shell (full)Complete shell, full of cargo molecules
Complete shell (empty)rComplete shell, almost empty of cargo molecules
AttachedNearly complete shells attached to a globule by a neck of cargo
Over-nucleated/MalformedMultiple globules, with incomplete or malformed shells on their surfaces
×StalledLarge globule with multiple incomplete or malformed shells on its surface
GlobuleCargo globule with unassembled shell subunits on its surface
UnnucleatedDiffuse subunits and cargo molecules


We have described an equilibrium theory and a dynamical computational model for the assembly of shells around a fluid cargo. Our simulations show that assembly can proceed by two classes of pathways: (i) a multi-step process in which the cargo forms a dense globule, followed by adsorption, assembly, and budding of shell proteins, or (ii) single-step assembly, with simultaneous aggregation of cargo molecules and shell assembly. This result demonstrates that the minimal interactions included in our model are sufficient to drive both classes of assembly pathways, suggesting that they are a generic feature of assembly around a fluid cargo. Moreover, while we cannot rule out the existence of active mechanisms in biological examples such as carboxysomes, our model demonstrates that the same interactions which drive assembly of shells can also drive budding from and closure around an amorphous globule of cargo.

Our results suggest bounds on the relative strengths of interactions that drive BMC assembly in cells. The decisive control parameter determining the assembly pathway is the cohesive energy between cargo molecules, which could arise through direct cargo-cargo interactions or be mediated by auxiliary proteins (Cameron et al., 2013). Relatively weak cargo interactions lead to single-step assembly pathways, while stronger interactions favor formation of the cargo-shell globule. However, the strength of cargo-shell and shell-shell interactions also play a role. Strong shell-shell interactions cause assembly to proceed rapidly during globule formation, limiting the size of the globule. Moreover, if a large globule is already present (e.g. due to time-dependent protein concentrations within a cell), strong interactions tend to drive malformed assemblies. We find that an important functional difference between the two classes of assembly pathways is control over the amount of packaged cargo. While the multi-step assembly pathways always generate a shell filled with cargo molecules, shells assembling around a diffuse cargo can be tuned from nearly empty to completely full by controlling the strength of cargo-shell interactions.

These results have implications for reengineering BMCs to encapsulate new core enzymes. Recent works demonstrated that protein cargos can be targeted to BMCs via encapsulation peptides that mediate cargo-shell interactions. However, packaged amounts were much lower than for native core enzymes (Parsons et al., 2010; Choudhary et al., 2012; Lassila et al., 2014). Our simulations show that both cargo-shell and cargo-cargo interactions (direct or mediated) must be controlled to assemble full shells.

We also find that a general equilibrium theory describes the ranges of parameter values for which assembly occurs. However, the dynamical simulations demonstrate that, at finite timescales, there is a rich variety of assembly morphologies. Formation of ordered, full shells requires a delicate balance of cargo-cargo, cargo-shell, and shell-shell interactions, all of which must be on the order 5-10kBT. This constraint is consistent with previous studies on viruses and other assembly systems, which found that formation of ordered states requires multiple, cooperative weak interactions between subunits (Hagan, 2014; Whitelam and Jack, 2015). Outside of optimal parameter regimes, the simulations predict alternative outcomes, ranging from no assembly to various alternative trapped intermediates, with the morphology depending on which interaction is strongest. We find that assembly is least robust to parameter variations when the cargo crystallizes before shell assembly. The assembling shell is unable to deform or penetrate the cargo complex, leading to defect-riddled, non-budded complexes. Within the limits of our simplified model, this observation suggests that procarboxysome complexes are at least partially fluid prior to successful shell assembly. Moreover, we find that observations of ordered cargo within assembled shells may be explained by packing constraints.

An important limitation of the present study is that the model interactions are specific to the shell geometry shown in Figure 1 (containing 20 hexamers) because alternating edges on hexagonal subunits have attractive interactions only with pentagonal subunits. In reality BMCs contain many more hexamers (formed from multiple protein sequences) and thus must include a greater range of hexamer-hexamer interactions. Extension of the model to allow for this possibility would allow consideration of two important questions: (1) The mechanism controlling insertion of the 12 pentagons required for a closed shell topology. (2) The relationship between assembly pathway and BMC size polydispersity. In particular, experiments suggest that β-carboxysomes are more polydisperse than α-carboxysomes (Price and Badger, 1991; Shively et al., 1973; Shively et al., 1973; Iancu et al., 2007; 2010; Kerfeld et al., 2010; Tanaka et al., 2008). We speculate that in the case of assembly around vapor-phase cargo, the size of the assembling shell will be primarily dictated by the preferred shell protein curvature and thus relatively uniform. However, during assembly around a condensed globule, the shell protein interactions could be strained to accommodate a globule which is larger or smaller than the preferred curvature, causing the shell size to depend on a complex balance of intermolecular interaction strengths and variables such as the local RuBisCO concentration.

Our model is minimal, intended to elucidate general principles of assembly around a fluid cargo, and thus may apply to diverse systems including prokaryotic microcompartments, viruses, and engineered delivery vehicles. The predicted trends for how assembly mechanisms and morphologies vary with control parameters can be experimentally tested by microscopy experiments. Such testing will be most straightforward in vitro (e.g. Luque et al., 2014; Douglas and Young, 1998; Rurup et al., 2014; Patterson et al., 2014; Patterson et al., 2012; Zhu et al., 2014; Rhee et al., 2011; Rurup et al., 2014; Wörsdörfer et al., 2012), where subunit-subunit interactions can be tuned by varying solution conditions and the stoichiometries of shell and cargo species can be readily varied. While there is currently no BMC assembly system starting from purified components, our findings can be tested in vivo by mutations which alter known protein binding interfaces, or by altering expression levels of RuBisCO or carboxysome proteins.

We anticipate that our model can serve as a qualitative guide for understanding how such multicomponent complexes assemble in natural systems, or to reengineer them for new applications. More broadly, our results demonstrate that the properties of encapsulated cargo, such as its topology, geometry and interaction strengths, strongly influence assembly pathways and morphologies.

Materials and methods

Computational model

Shell subunits

Request a detailed protocol

We have adapted a model for virus assembly (Perlmutter et al., 2013; 2014; Perlmutter and Hagan, 2015a; Wales, 2005; Fejer et al., 2009; Johnston et al., 2010; Ruiz-Herrero and Hagan, 2015) to describe assembly of an icosahedral shell around a fluid cargo. Each subunit contains ‘Attractors’ on its perimeter that mediate subunit-subunit attractions (as in Ruiz-Herrero and Hagan, 2015). Attractor interactions are specific – complementary pairs of Attractors (see Figure 1A,B and appendix 1) have short-range interactions (modeled by a Morse potential), whereas non-complementary pairs have no interactions. A repulsive interaction between pairs of ‘Top’ (type ‘T’) pseudoatoms favors the correct subunit-subunit angle. The ‘Bottom’ (type ‘B’) pseudoatoms mediate short-ranged subunit-cargo attractions (e.g. due to interactions with shell ‘encapsulation peptides’ (Kinney et al., 2012; Cameron et al., 2013; Fan et al., 2010)), represented by a Morse potential. We also add a layer of ‘Excluders’ in the plane of the ‘Top’ pseudoatoms, which represent subunit-cargo excluded volume interactions. The strengths of subunit-subunit and subunit-cargo attractions are parameterized by potential well depths εSS and εSC respectively (appendix 1).


Request a detailed protocol

As a minimal representation of globular proteins, the cargo is modeled as spherical particles which interact via an attractive Lennard-Jones (LJ) potential, with well-depth εCC. The attractions implicitly model hydrophobic and screened electrostatics interactions between cargo molecules, as well as effective cargo-cargo interactions mediated by auxiliary proteins (e.g. the carboxysome protein CcmM (Cameron et al., 2013)).


Request a detailed protocol

We simulated assembly dynamics using the Langevin dynamics algorithm in HOOMD (a software package that uses GPUs to perform highly efficient dynamics simulations [Anderson et al., 2008]) and periodic boundary conditions to represent a bulk system. The subunits are modeled as rigid bodies (Nguyen et al., 2011). The simulations were performed using a set of fundamental units (URL. http://codeblue.umich.edu/hoomd-blue/doc/page_units.html), with 1du defined as the circumradius of the pentagonal subunit (the cargo diameter is also set to 1 du). Unless specified otherwise, each simulation contained enough subunits to form four complete shells (48 pentamers and 80 hexamers) and 611 cargo particles (a shell typically encapsulates 120–130 cargo particles) in a cubic box with side length 40du. The simulation time step was 0.001 in dimensionless time units, and dynamics was performed for 3×108 timesteps unless mentioned otherwise.

We performed two sets of simulations, using different initial conditions. In the first, simulations were initialized by introducing cargo particles and shell subunits simultaneously with random positions and orientations (except avoiding high-energy overlaps). The second set of initial conditions was motivated by the possibility that the cargo globule could form before shell subunits reach sufficient concentrations within the cell to undergo assembly. To model this situation, we pre-equilibrated the cargo by performing a long simulation with only cargo particles present. Shell subunits were then introduced with random positions and orientations (excluding high-energy overlaps). For εCC1.6, the assembly simulations thus began with a cargo globule already present. For εCC<1.6 the two protocols are equivalent, since no globule forms during cargo equilibration.

Sample sizes

Request a detailed protocol

To cover the largest range of parameter space possible given the computational expense associated with each simulation, we performed 5 independent simulations at most parameter sets. To assess statistical error and to estimate the distribution of different assembly outcomes, we performed 10 independent trials for one value of εSS at each value of εSC and εCC. We also performed additional simulations at parameter sets for which 5 trials did not result in a majority outcome, or when necessary to obtain better statistics on the number of encapsidated cargo particles. Based on these results, performing additional simulations at other parameter values would not qualitatively change our results. (It would increase the statistical accuracy of estimated boundaries between different outcomes; however, these boundaries correspond to crossovers rather than sharp transitions.)

Thermodynamics of assembly around a fluid cargo

Request a detailed protocol

To complement the finite-time simulations, we have developed a general thermodynamic description of assembly around a fluid cargo. We consider shells composed of species α=1,2,M, with nαshell subunits of species α in a complete shell, which encapsulates n0 cargo molecules (the index 0 refers to cargo molecules henceforth). Assembly occurs from a dilute solution of cargo molecules with density ρ0, shell subunits with density ρα for each species, and the density of assembled, full shells as ρshell. These are in equilibrium with a globule containing n0glob cargo molecules and nαglob subunits for each species α. We assume that, due to the asymmetric nature of the shell-cargo interaction, the shell subunits reside at the exterior of the globule (as we observe in our simulations). The globule containing unassembled shell subunits thus resembles a spherical microemulsion droplet (Safran, 1994). Minimizing the total free energy (see appendix 2) gives:

(1) v0ρshell=exp[(Gshellαnαshellμα)/kBT]

where Gshell is the interaction free energy of the assembled shell and μα are the chemical potentials of free cargo molecules and shell subunits, given by μα=kBTln(ραv0), with v0 a standard state volume and the globule composition given by

(2) Gglob({nαglob})nαglob=μαfor α=0M,

with Gglob(nsglob, n0glob) as the globule free energy.

(1) – (2) are the general equilibrium description for a system of assembling shells with a disordered-phase intermediate; application to a specific system requires specifying the forms of Gshell and Gglob. In appendix 2 we specify these equations for our computational model, allowing us to compare the equilibrium calculation with simulation results, using no free parameters.

To compare the relative stabilities of the globule and assembled shells, we also calculate the free energy difference

(3) Δfassem=ftot({nαglob=0})-ftot(ρshell=0),

where the first term on the right-hand side is the minimized free energy for a system containing shells and free subunits but no globule, while the second term corresponds to the minimized free energy for a system containing subunits and the globule, but no assembly.

Appendix 1: Model Details

1.1 Interaction potentials

Our subunit model is based on a model for viral capsid assembly, developed by Wales (Wales, 2005) and Johnston et al. (Johnston et al., 2010), which we have adapted to describe interactions with cargo molecules.

Each subunit contains ‘Attractors’ on its perimeter that mediate subunit-subunit attraction (as in [Ruiz-Herrero and Hagan, 2015]). Attractor interactions are specific – complementary pairs of Attractors have short-range interactions (modeled by a Morse potential), whereas non-complementary pairs have no interactions. For simplicity, complementarity is defined based only on the low-energy structure (Figure 1D); i.e., there is no attraction between pairs of pentagons. Complementary pairs of attractors are: for the hexagon-hexagon interaction, A4-A4, A5-A6, and for the hexagon-pentagon interaction A1-A4, A2-A8, A3-A7. The strength of attractive interactions is parameterized by the well-depth εSS. Because vertex attractors (A1, A4) have multiple partners in an assembled structure, whereas edge attractors have only one, the well-depth for A1-A4 and A4-A4 interactions is set to εSS/2, while all other attractor interactions use εSS. The ‘Top’ height, or distance out of the attractor plane, sets the Top-Top distance between interacting subunits, which determines the preferred subunit-subunit angle. We use a height of h=1/2rb, with rb=1 the distance between a vertex attractor and the center of the pentagon. The ‘Bottom’ (type ‘B’) pseudoatoms mediate subunit-cargo attractions, represented by a Morse potential with well-depth εSC. We also add a layer of ‘Excluders’ in the plane of the ‘Top’ pseudoatoms (positioned as in Figure 1), which represent subunit-cargo excluded volume interactions.

In our model, all potentials can be decomposed into pairwise interactions. Potentials involving container subunits further decompose into pairwise interactions between their constituent building blocks – the excluders, attractors, ‘Top’, and ‘Bottom’ pseudoatoms. It is convenient to state the total energy of the system as the sum of three terms, involving subunit-subunit (USS), cargo-cargo (ULJ), and subunit-cargo (UAds) interactions, each summed over all pairs of the appropriate type:

(A1) U=sub  isub j<iUSS+cargo  icargo j<iULJ+sub icargo jUAds

where sub isub j<i is the sum over all distinct pairs of subunits in the system, sub icargo j is the sum over all subunit-cargo particle pairs, etc.

Subunit-subunit interactions

The subunit-subunit potential USS is the sum of the attractive interactions between complementary attractors, and geometry guiding repulsive interactions between ‘Top’ - ‘Top’, ‘Bottom’ - ‘Bottom’, and ‘Top’ - ‘Bottom’ pairs. There are no interactions between members of the same rigid body. Thus, for notational clarity, we index rigid bodies and non-rigid pseudoatoms in Roman, while the pseudoatoms comprising a particular rigid body are indexed in Greek. For subunit i we denote its attractor positions as {𝐚iα} with the set comprising all attractors α, its ‘Top’ position {𝐭i}, and ‘Bottom’ position {𝐛i}. The subunit-subunit interaction potential between two subunits i and j is then defined as:

(A2) USS({aiα},ti,aj,tj)=εSSŁ(|titj|, σt,ij)+εSSŁ(|bibj|, σb)+εSSŁ(|bitj|, σtb)+α,βNai,NajεSS(|aiαajβ|, r0,ϱ,rcutatt)

where εSS is an adjustable parameter which both sets the strength of the subunit-subunit attraction at each attractor site and scales the repulsive interactions which enforce the geometry, Nai is the number of attractor pseudoatoms in subunit i, σtb=1.8rb is the diameter of the ‘Top’ - ‘Bottom’ interaction (this prevents subunits from binding in inverted configurations (Johnston et al., 2010), and σb=1.5rb is the diameter of the ‘Bottom’ - ‘Bottom’ interaction.

In contrast to the latter parameters, σt,ij the effective diameter of the ‘Top’ - ‘Top’ interaction, depends on the species of subunits i and j; denoting a pentagonal or hexagonal subunit as p or h respectively, σt,pp=2.1rb, σt,hh=2.436rb, and σt,ph=2.269rb. The parameter r0 is the minimum energy attractor distance, set to 0.2rb, ϱ is a parameter determining the width of the attractive interaction, set to 4rb, and rcutatt is the cutoff distance for the attractor potential set to 2.0rb. Since the interactions just described are sufficient to describe assembly of the shell subunits, we included no excluder-excluder interactions.

The function Ł is defined as the repulsive component of the Lennard-Jones potential shifted to zero at the interaction diameter:

(A3) Ł(x,σ)θ(σx)[(σx)121]

with θ(x) the Heaviside function. The function is a Morse potential:

(A4) (x,r0,ϱ,rcut)=θ(rcutx)×[(eϱ(1xr0)2)eϱ(1xr0)Vshift(rcut)]

with Vshift(rcut) the value of the potential at rcut.

Cargo-cargo interactions

The interaction between cargo particles is given by

(A5) ULJ({li},{lj})=i<jNlεCC(|litj|, σC,rcutc)

with the full Lennard-Jones interaction:

(A6) (x,σ,rcut)=θ(xrcut) ×{4[(xσ)12(xσ)6]Vshift(rcut)}

and εCC is an adjustable parameter which sets the strength of the cargo-cargo interaction, Nl is the number of LJ particles, σC is the cargo diameter set to 1.0rb except where mentioned otherwise, and rcutc is set to 3σC.

Subunit-cargo interactions

The subunit-cargo interaction is a short-range repulsion between cargo-excluder and cargo-‘Top’ pairs reresenting the excluded volume plus an attractive interaction between the cargo - ‘Bottom’ pairs. For subunit i with excluder positions {𝐱iα} and ‘Bottom’ psuedoatom {𝐛iα} and cargo particle j with position 𝐑j, the potential is:

ŁŁUAds({xiα},Rj)=αNxŁ(|xiαRj|,σex)+αNtŁ(|tiαRj|,σt)+αNbεSC(|ciαRj|, r0,ϱ,rcut)

where εSC parameterizes the shell-cargo interaction strength, Nx, Nt, and Nb are the numbers of excluders, ‘Top’, and ‘Bottom’ pseudoatoms on a shell subunit, σex=0.5rb and σt=0.5rb are the effective diameters of the Excluder - cargo and ‘Top’ - cargo repulsions, r0 is the minimum energy attractor distance, set to 0.5rb, ϱ is a parameter determining the width of the attractive interaction, set to 2.5rb, and rcut is the cutoff distance for the attractor potential set to 3.0rb.

Motivation for choice of interaction potentials

The choices we have made for potential functions (Morse or Lennard-Jones) between different classes of pseudoatoms are based on the need for tunability of the interaction length scale and the extent to which guidance on parameterization is available from the existing literature. In particular, the Morse potential enables controlling the interaction length scale independently from the particle excluded volume size, whereas the interaction length scale and excluded volume size are tuned by a single parameter in the Lennard-Jones potential. Our shell-shell interaction potential is based on previous models for viral capsid assembly (Wales, 2005; Johnston et al., 2010; Ruiz-Herrero and Hagan, 2015; Perlmutter et al., 2013; 2014; Perlmutter and Hagan, 2015b), and the choice of a Morse potential for attractor-attractor interactions and a Lennard-Jones potential for Top-Top interactions follows these previous works. The attractor interactions are modeled using a Morse potential because the length scale of their interaction strongly affects the subunit orientational specificity. We chose to model the cargo-cargo interaction using a Lennard-Jones potential because the phase behavior for this model has been extensively studied in the literature, thus limiting the need for model parameterization. However, we note that it could be of interest to study how the probability of shell detachment depends on the length scale of the cargo-cargo interaction; we speculate that a longer-range interaction would increase the probability of detachment by making it easier for shell subunits to penetrate into the globule. Finally, the shell-cargo interactions could have used either choice of potential; we elected to use a Morse potential due to its greater flexibility.

1.2 Maximum cargo loading

To give context to the densities of packaged cargo particles that we observe in simulations, we estimate the maximum possible cargo loading here. Our assembled shell has the geometry of a truncated icosahedron with an edge length of approximately 1.5du. Accounting for the volume occluded to cargo particles by the shell pseudoatoms, the interior volume is Vin109du3. The maximum number of cargo molecules that can be packaged (assuming hexagonal close packing) is thus NHCP154. However, this is an overestimate since the shell geometry is not commensurate with perfect hexagonal close packing. We thus estimate NHCP=148, the maximum number of packaged cargo particles seen in an equilibrium simulation. The maximum cargo loading for random close packing is then NRCP120.

Appendix 2 Thermodynamics of assembly around a fluid cargo

2.1 General theory

In this section we present a general thermodynamic description for assembly around a fluid cargo. The theory provides a description of phase behavior in terms of simple physical parameters, and enables evaluating the extent to which our finite-time dynamical simulations have approached equilibrium. We assume that the equilibrium distribution is dominated by three classes of system configurations: free cargo and shell subunits, a disordered globule of cargo molecules with unassembled shell subunits on its surface, and assembled shells filled with cargo molecules. Extension to consider partially assembled intermediates and partially filled shells is straightforward but would complicate the presentation; moreover, at conditions leading to productive assembly, concentrations of partially assembled intermediates are negligible at equilibrium (Hagan, 2009; 2014; Safran, 1994; Gelbart et al., 1994).

We consider shells composed of species α=1,2,M, with nαshell subunits of species α in a complete shell, which encapsulates n0 cargo molecules (the index 0 refers to cargo molecules henceforth). Assembly occurs from a dilute solution of cargo molecules with density ρ0, shell subunits with density ρα for each species, and the density of assembled, full shells as ρshell. These are in equilibrium with a globule containing n0glob cargo molecules and nαglob subunits for each species α. The total free energy density is then given by

(B1) ftot=α=0MkBTρα[ln(ραv0)1]+kBTρshell[ln(ρshellv0)1]+ρshellGshell+V1Gglob(n0glob,{nαglob})

where the sum runs over free cargo molecules and shell subunits, V is the system volume, v0 is a standard state volume, Gshell is the interaction free energy of the assembled shell, and Gglob(nsglob,n0glob) is the globule free energy. We then minimize ftot with respect to Nshell=Vρshell and {nαglob}, subject to the conservation of mass constraints:

(B2) ραT=ρα+nαglob/V+ρshellnαshellfor α=0M

where ραT denotes the total density of species α.

This results in Equations (1–2) of the main text.

2.2 Specification to our computational model

Equations (1–2) are the general equilibrium description for a system of assembling shells with a disordered-phase intermediate. To explore how assembly depends on the control parameters (εCC, εSC, εSS, ρsT, and ρsT) and to compare these equilibrium expressions against our simulation results, we now specify these relations to our computational model.

2.2.1 Globule and shell interaction free energies

We model the globule as a liquid droplet of Lennard-Jones (LJ) particles, with shell subunits adsorbed to its exterior surface. For simplicity, we treat shell subunit binding to the globule with the Langmuir adsorption model. To simplify the notation, we suppress dependencies on control parameters in the free energy expressions, but list them beneath. The free energy of the globule is then given by

(B3) Gglob(npglob,nhglob,n0glob)=γAglob(n0glob)+μliqn0glob+gAds(npglob+nh glob)+Gmix(npglob,nhglob,nmax(Aglob(n0glob)),

where γ(εCC) and μliq(εCC) are the bulk surface tension and chemical potential of a LJ liquid, gAds(εSC) is the shell subunit absorption free energy, npglob and nhglob are the numbers of adsorbed pentamers and hexamers respectively, Aglob=(34πρliq(εCC)n0glob)2/3 is the area of the globule, and ρliq(εCC) is the density of the LJ liquid. The final term is the mixing entropy of adsorbed subunits according to Langmuir adsorption, given by

(B4) Gmix(npglob,nh glob,nmax)/kBT=ln(nmaxnpglob,nhglob,nmax(npglob+nhglob)),

with nmax as the number of adsorbed subunits at saturation (calculated from simulations, see below).

For the free energy of shell assembly, we consider a shell comprised of npent=12 pentamers and nhex hexamers, which have nph pentamer-hexamer contacts with binding energy εph and nhh hexamer-hexamer contacts with energy εhh. For our T=3 model, nhex=20, nph=60, and nhh=30. The assembly free energy is then given by

(B5) Gshell=nphεph+nhhεhhT(npentspent+nhexshex +sconfig)+γAglob(n0glob)+μliqn0glob+gAds(npent+nhex),

with spent and shex the translational and rotational entropy penalty associated with binding of pentameric or hexameric subunits and sconfig accounting for the configurational entropy associated with subunit and shell symmetries. In our model the pentamers, hexamers, and capsid are 5-fold, 3-fold, and 60-fold symmetric, giving sconfig=kBln(5npent3nhex/60). Other parameters were calculated from simulations, as described next.

2.2.2 Determination of parameter values

Since our interactions are constructed from standard potential functions, some of the parameters discussed in the last section are known from the literature, and others can be calculated from simulations. Thus, it is possible to compare our equilibrium theory against simulation results with no fitting parameters. We present the parameter values and how they are obtained in this section.

Cargo parameters

The parameters characterizing the phase behavior of a Leonard-Jones fluid, γ, μliq, and ρliq can be obtained from the literature, but we performed fits specific to the parameter ranges of interest, 1.0εCC/kBT3.0. The surface tension γ was estimated using the approach of Mecke et al. (Mecke et al., 1997). We performed separate simulations containing only LJ particles, with numbers of particles and volume for each system set to achieve formation of a planar liquid vapor interface, and varying values of the LJ interaction strength εCC. We then calculated γ from the virial expression. For our LJ potential, truncated at rcut=3σ, we obtain (using the functional form of Ref. (Mecke et al., 1997)

(B6) γ(εCC)=2.936(1-εCC-11.3)1.688.

From the same simulations, we calculated the dependence of the bulk liquid density on εCC as

(B7) ρliq(εCC)=-1.439+2.165εCC0.115.

Although there are a number of empirical forms for the LJ equation of state available in the literature, they vary widely in complexity, number of fit parameters, and presumably accuracy over the parameter range we are interested in. We therefore estimated the liquid chemical potential μliq from the vapor-phase densities ρvap in LJ liquid-vapor coexistence simulations according to

(B8) μliq=kBTln(ρvapσ3)-Aγ/Nliq,

where A is the interfacial area, Nliq is the total number of particles in the liquid phase as a function of εCC, and γ is given by Equation B6. The results are fit well by the linear function

(B9) μliq(εCC)=3.13kBT-5.6εCC.

Shell subunit-subunit interactions

We estimated the subunit-subunit binding free energy values as functions of the well-depth parameter εSS by measuring the dimerization equilibrium constant in simulations of subunits only capable of forming dimers (Figure 1C). For both pentamer-hexamer and hexamer-hexamer dimers, we obtain binding free energies which are linear functions of the well-depth εSS. We interpret the y-intercept as the binding entropy, giving:


where the standard state volume is du3.

In Equation (B5) we then make the assumption that, because the interactions are orientationally specific, a subunit incurs its entire binding entropy penalty upon dimerization — because a bound subunit is already aligned to form additional interactions, these interactions do not lead to further entropy penalties. In reality, this is an under-prediction since some additional entropy losses occur on making additional bonds (Hagan and Chandler, 2006; Hagan et al., 2011), but these are not sufficiently large to qualitatively affect our results.

Appendix 2—figure 1
(A) Langmuir isotherms to estimate gAds(εSC).

(B) Estimate of the chemical potential for an equilibrated LJ system (before correcting for the finite size of liquid droplet). (C) Fit of the subunit dimerization free energies ghh(εSS) and gph(εSS) as a function of the well depth parameter εSS. (D) Fit of LJ droplet surface tension, including the tail correction.


Shell subunit adsorption onto globule

We estimated the shell subunit adsorption free energy by performing simulations of subunits which cannot assemble (εSS=0) in the presence of a cargo globule. We then measured the globule size and number of adsorbed subunits as functions of εSC. We found the results could be fit using the Langmuir adsorption model, with the adsorption free energy of a single subunit gAds as a fit parameter for each value of εSC. We assumed that the maximum number of adsorbed subunits (the number of lattice sites in the Langmuir model) does not directly depend on εSC, and hence fit this parameter globally, obtaining nmax=80 for a globule with n0glob=300 cargo molecules. In our calculations we assume that nmax is proportional to the globule surface area, consistent with observations from simulations. Our fit resulted in a linear relationship between the adsorption energy and free energy over the range of interest:

(B10) gAds=0.093kBT-1.17εSC.


    1. Bobik TA
    2. Havemann GD
    3. Busch RJ
    4. Williams DS
    5. Aldrich HC
    The propanediol utilization (pdu) operon of Salmonella enterica serovar Typhimurium LT2 includes genes necessary for formation of polyhedral organelles involved in coenzyme B(12)-dependent 1, 2-propanediol degradation
    Journal of Bacteriology 181:5967–5975.
    1. Borodavka A
    2. Tuma R
    3. Stockley PG
    (2012) Evidence that viral RNAs have evolved for efficient, two-stage packaging
    Proceedings of the National Academy of Sciences of the United States of America 109:15769–15774.
  1. Book
    1. Safran S
    Statistical Thermodynamics of Surfaces, Interfaces, and Membranes
    Addison-Wesley Pub.
    1. Shively JM
    2. Ball FL
    3. Kline BW
    Electron microscopy of the carboxysomes (polyhedral bodies) of Thiobacillus neapolitanus
    Journal of Bacteriology 116:1405–1411.
    1. Wales DJ
    (2005) The energy landscape as a unifying theme in molecular science
    Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 363:357–377.

Article and author information

Author details

  1. Jason D Perlmutter

    Martin Fisher School of Physics, Brandeis University, Waltham, United States
    JDP, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Farzaneh Mohajerani
    Competing interests
    The authors declare that no competing interests exist.
  2. Farzaneh Mohajerani

    Martin Fisher School of Physics, Brandeis University, Waltham, United States
    FM, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Jason D Perlmutter
    Competing interests
    The authors declare that no competing interests exist.
  3. Michael F Hagan

    Martin Fisher School of Physics, Brandeis University, Waltham, United States
    MFH, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9211-2434


National Institute of General Medical Sciences (R01GM108021)

  • Jason D Perlmutter
  • Farzaneh Mohajerani
  • Michael F Hagan

National Science Foundation (DMR-1420382)

  • Michael F Hagan

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


We are grateful to Maxim Prigozhin for illuminating discussions and for introducing us to the carboxysome assembly problem, and to Fei Cai, Cheryl Kerfeld and Charles Knobler for comments on the manuscript. This work was supported by Award Number R01GM108021 from the National Institute Of General Medical Sciences and the Brandeis Center for Bioinspired Soft Materials, an NSF MRSEC, DMR-1420382. Computational resources were provided by NSF XSEDE computing resources (Maverick and Keeneland) and the Brandeis HPCC which is partially supported by DMR-1420382. MFH performed part of this work while at the Aspen Center for Physics, which is supported by NSF grant PHY-1066293.

Version history

  1. Received: December 27, 2015
  2. Accepted: May 10, 2016
  3. Accepted Manuscript published: May 11, 2016 (version 1)
  4. Version of Record published: June 28, 2016 (version 2)
  5. Version of Record updated: July 12, 2016 (version 3)


© 2016, Perlmutter et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 2,763
    Page views
  • 516
  • 30

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Jason D Perlmutter
  2. Farzaneh Mohajerani
  3. Michael F Hagan
Many-molecule encapsulation by an icosahedral shell
eLife 5:e14078.

Share this article


Further reading

    1. Microbiology and Infectious Disease
    2. Physics of Living Systems
    Fabien Duveau, Céline Cordier ... Pascal Hersen
    Research Article

    Natural environments of living organisms are often dynamic and multifactorial, with multiple parameters fluctuating over time. To better understand how cells respond to dynamically interacting factors, we quantified the effects of dual fluctuations of osmotic stress and glucose deprivation on yeast cells using microfluidics and time-lapse microscopy. Strikingly, we observed that cell proliferation, survival, and signaling depend on the phasing of the two periodic stresses. Cells divided faster, survived longer, and showed decreased transcriptional response when fluctuations of hyperosmotic stress and glucose deprivation occurred in phase than when the two stresses occurred alternatively. Therefore, glucose availability regulates yeast responses to dynamic osmotic stress, showcasing the key role of metabolic fluctuations in cellular responses to dynamic stress. We also found that mutants with impaired osmotic stress response were better adapted to alternating stresses than wild-type cells, showing that genetic mechanisms of adaptation to a persistent stress factor can be detrimental under dynamically interacting conditions.

    1. Physics of Living Systems
    Josep-Maria Armengol-Collado, Livio Nicola Carenza, Luca Giomi
    Research Article Updated

    We formulate a hydrodynamic theory of confluent epithelia: i.e. monolayers of epithelial cells adhering to each other without gaps. Taking advantage of recent progresses toward establishing a general hydrodynamic theory of p-atic liquid crystals, we demonstrate that collectively migrating epithelia feature both nematic (i.e. p = 2) and hexatic (i.e. p = 6) orders, with the former being dominant at large and the latter at small length scales. Such a remarkable multiscale liquid crystal order leaves a distinct signature in the system’s structure factor, which exhibits two different power-law scaling regimes, reflecting both the hexagonal geometry of small cells clusters and the uniaxial structure of the global cellular flow. We support these analytical predictions with two different cell-resolved models of epithelia – i.e. the self-propelled Voronoi model and the multiphase field model – and highlight how momentum dissipation and noise influence the range of fluctuations at small length scales, thereby affecting the degree of cooperativity between cells. Our construction provides a theoretical framework to conceptualize the recent observation of multiscale order in layers of Madin–Darby canine kidney cells and pave the way for further theoretical developments.