How symmetry and complexity shape three-dimensional programmable self-assembly

Introduction

Assembling a complex structure from its components requires information that specifies how these individual components interact and organize1–5. In a top-down approach, assembly information is stored in a centralized blueprint that directs the assembly process, much like a robot assembling parts on an assembly line. But information can also be encoded locally in the interactions between components, allowing ordered structures to emerge spontaneously without external control. This bottom-up process, known as self-assembly, underlies many biological assemblies such as viral capsids6–8, lipid membranes9, microtubules10, and protein complexes4,5.

These biological systems have inspired the design of artificially programmed interaction rules between components so that target structures can be assembled, even at scales where individual components cannot be manipulated directly11. A particularly well-developed example is DNA nanotechnology12, where approaches such as DNA origami13, DNA bricks14, DNA wireframe structures15,16, and DNA cubes17,18 enable the construction of remarkably complex nanoscale structures19. Another promising approach includes DNA-coated colloids20, which enable the design of patchy-particle systems into complex structures21.

The information required to self-assemble a given finite structure can be locally encoded in two contrasting ways. In fully addressable designs, every component and interaction is uniquely specified4,22–24. Such designs provide precise control over assembly but are information-intensive, difficult to scale, costly to implement25, and need thermodynamic annealing protocols23,26 to assemble due to the high configurational entropy cost. At the other extreme are minimal designs, which generate complex structures by reusing the smallest set of components and interactions24,27,28. Such designs require substantially less assembly information and reduce synthesis costs. They may also have larger thermodynamic windows for high-yield assembly24. Nature frequently exploits this strategy. For example, in viral capsid assembly, large symmetric shells are constructed from repeated copies of a small number of protein subunits6,8. It is an open question whether these near-minimal solutions arise because they confer adaptive advantages6,29,30, or whether they occur preferentially simply because they are easier to find by evolutionary search28.

Some structures can be assembled from relatively simple minimal designs, reusing many copies of only a few component types, whereas others require designs that are close to fully addressable assembly. These differences allowed Cartwright and Mackay to introduce assembly complexity, defined as the minimum information required to specify the assembly of a structure31, as a general measure of structural complexity. Traditionally, the complexity of crystal structures has been characterized in terms of symmetry, with highly symmetric structures regarded as simpler than those with lower symmetry. However, many biological and other self-assembled structures are aperiodic, making symmetry an incomplete measure of structural complexity. Cartwright and Mackay argued that assembly complexity provides a more general description because it applies equally to periodic and aperiodic structures31,32. Related information-theoretic approaches proaches have also been used to quantify the structural complexity of inorganic crystals33. This work motivates the first question our paper will address: 1) How well does assembly complexity correlate with symmetry?

A second question follows from a simple intuition: structures whose assembly complexity is much lower than that of a fully addressable design should admit many more alternative assembly rules than structures whose minimal complexity is already close to this maximum. To quantify this intuition, we define the assembly rule probability of a structure as a measure of how frequently it is generated when assembly rules are sampled at random. This naturally raises the question: 2) How does assembly rule probability scale with assembly complexity and with symmetry?

For insight into this second question we first go down a more abstract route, and note that assembly complexity has parallels with Kolmogorov complexity $K(x)$ from algorithmic information theory (AIT)34, defined as the length of the shortest algorithm that outputs the string $x$ on a universal Turing machine (UTM). Unlike Kolmogorov complexity, which is uncomputable34, assembly complexity is grounded in the computable physical process of assembly itself.

A central concept in AIT is algorithmic probability. For a string $x$, the algorithmic probability $P_U(x)$ is defined as the probability that a (suitably chosen) UTM $U$ will produce $x$ when supplied with a randomly chosen input program. Since algorithmic probabilities for different UTMs differ only by machine-dependent multiplicative constants34, the UTM dependence can be omitted, and the probability written simply as $P(x)$. Levin’s coding theorem34,35 then predicts that $P(x) \approx2^{-K(x)}$. Intuitively, if an output string can be produced by a short input program, then it can be produced by exponentially many other longer programs as well. While this powerful theorem can only be proven for UTMs, Dingle et al.[36] showed that for many non-universal input–output maps, an analogue to algorithmic probability of AIT, the computable probability $P(x)$ of generating output $x$ upon random inputs, can be bounded from above as $P(x) \lesssim2^{-a\widetilde{K}(x)-b}$, where $\widetilde{K}(x)$ is an appropriate approximation to the true Kolmogorov complexity and $a > 0$, $b$ are map-dependent constants. This simplicity bias implies that high $P(x)$ outputs have low $\widetilde{K}(x)$, and high $\widetilde{K}(x)$ outputs have low $P(x)$, although low $\widetilde{K}(x)$ does not necessarily imply high $P(x)$. However, a statistical lower bound constrains the total probability weight of outputs that lie far below the upper bound[37], so that the majority of the probability weight should not be too far from the upper bound.

These AIT-based arguments are abstract but also broadly applicable. The predicted simplicity bias has been observed in diverse systems including two-dimensional polyomino self-assembly28, Richard Dawkins’ biomorphs[38], a model of tooth development[39], and RNA secondary structure[36] [28]. For functional non-coding RNA, the calculated $P(x)$ closely recapitulates the frequencies of folded structures observed in nature[40, 41, 42]. The consistency of this behavior across such different systems suggests that assembly rule probability should likewise decrease exponentially with assembly complexity in a broad class of self-assembling systems.

To address the two questions we posed above in a concrete setting, we study tile-based lattice models of self-assembly, which have the advantage of being tractable while still enabling a rich variety of structures to assemble. Tile-based lattice models trace back to Wang tiles, introduced by Hao Wang in 1961 as square tiles with labeled edges that can be placed only when adjacent edges match[43]. Originally developed to study mathematical tilings, Wang tiles later became the foundation of algorithmic self-assembly[44]. Building on Ned Seeman’s demonstrations of sequence-directed DNA self-assembly12, Winfree recognized that DNA tiles could function as programmable tiles and developed the abstract Tile Assembly Model (aTAM)[45, 46]. In the aTAM, structures grow from an initial seed by sequential tile addition governed by local glue-matching rules and binding-strength thresholds. Despite relying only on local interactions, the model can generate arbitrarily complex structures[47, 48] and even perform Turing-universal computation[45] [11], establishing a formal connection between self-assembly and Kolmogorov complexity[46, 47, 49].

A two-dimensional model related to the aTAM is the polyomino model of self-assembly in which rotatable square tiles with fixed-strength, non-cooperative interactions assemble into connected polyominoes[50, 51] [27]. Unlike the aTAM, the polyomino model is not Turing universal, but it provides a simple and physically motivated setting for studying programmable self-assembly. It has proved a useful model of protein quaternary structure assembly and provides a genotype–phenotype map in which genotypes encode assembly rules for square tiles (protein subunits), while phenotypes correspond to assembled polyominoes (protein complexes)[50, 51, 52, 53] [27]. Random sampling of these genotypes, as well as neutral evolution on the resulting genotype–phenotype map, produces, to first order at least, highly symmetric structures exponentially more often than asymmetric ones28. Remarkably, the resulting symmetry distributions closely match those of protein complexes in the Protein Data Bank (PDB), suggesting that the prevalence of symmetry may arise not only through adaptive selection but also through simplicity bias: highly symmetric structures generally have lower assembly complexity and are therefore discovered more readily by evolutionary search.

Much of the work on programmable tile assembly has focused on two dimensions. Yet many biological and engineered assemblies, including protein complexes, are inherently three-dimensional. Here we study a three-dimensional analogue of the polyomino model, in which cubic tiles with programmable interactions assemble into connected polycubes[54] [18, 24]. This allows us to investigate how assembly complexity, assembly rule probability, and symmetry are related in three-dimensional self-assembly.

Our main finding is that assembly rule probability decreases approximately exponentially with assembly complexity, consistent with the simplicity bias previously observed in two-dimensional polyomino systems28. Symmetry can reduce assembly complexity by allowing tile species and interface types to be reused, leading to a strong enhancement of certain highly symmetric structures under random sampling of assembly rules. However, symmetry alone does not determine assembly complexity or assembly rule probability. For example, a four-cube linear rod and a four-cube square have the same global symmetry, yet deterministic assembly of the rod requires substantially more information because its tile species and interface types cannot be reused as effectively. Assembly complexity therefore depends not only on global symmetry but also on how effectively local interaction rules can be reused during assembly.

The remainder of the paper is organized as follows. We first introduce the three-dimensional polycube model and define the assembly rule space, assembly rule probability, symmetry classification, and measures of assembly complexity. We then investigate how assembly rule probability depends on structure size, complexity, and symmetry. Finally, we examine why some symmetries reduce assembly complexity more effectively than others, and why structures with identical global symmetry can nevertheless differ substantially in their assembly rule probability.

Model and Methods

Polycube model overview

A polycube is a finite connected cluster of unit cubes joined face-to-face on a three-dimensional cubic lattice. The size of a polycube, denoted by $n$, is the number of constituent cubic tiles. In the polycube self-assembly model, structures emerge through local interactions between these tiles24. The complete specification of an assembly system is given by an assembly kit, consisting of $s$ tile species, with $c$ possible interface types assigned to their faces, and an interaction matrix defining which interface types bind.

Each cubic tile has six faces aligned with the Cartesian axes. Every face is assigned two attributes:

1. Interface type: Each face carries an integer label $i \in\{0,1,\ldots,c-1\}$ which defines its interface type. The label $i=0$ represents a neutral interface that does not bind to anything, while non-zero labels correspond to active interface types that form interacting pairs with their designated complementary interfaces as specified by the interaction matrix. If $c$ is even, the final label $c-1$ is also considered neutral. Following the convention of Wang tiles[43], these interface types are also referred to as colors.

2. Orientation: Each face is additionally assigned an in-plane orientation $r \in\{\mathrm{O},\mathrm{U},\mathrm{I},\mathrm{D},\mathrm{L},\mathrm{R}\}$, corresponding to Out ↙, Up ↑, In ↗, Down ↓, Left ←, and Right →. The orientation is defined with respect to a right-handed coordinate system in which the observer looks down the $+z$ axis towards the $xy$-plane (see Figure 1). This may be visualized as an arrow drawn on each face indicating its in-plane direction. Because each face lies in a distinct coordinate plane, only four of the six orientation states are geometrically admissible per face: North–South faces admit $\{\mathrm{L},\mathrm{R},\mathrm{I},\mathrm{O}\}$, East–West faces admit $\{\mathrm{U},\mathrm{D},\mathrm{I},\mathrm{O}\}$, and Front–Back faces admit $\{\mathrm{U},\mathrm{D},\mathrm{L},\mathrm{R}\}$. Orientation determines how complementary faces must align during attachment.

A tile species is completely specified by the interface types and orientations of its six faces listed in fixed NESWBF (North, East, South, West, Back, Front) order (see Figure 1):

where $s_i$ and $s_r$ denote the corresponding sequences of interface labels and orientations.

Tiles are free to rotate during the assembly process, and two tiles are regarded as the same species if they differ only by a rigid rotation. For example, consider the tile species

Self-assembly growth process

A rotation by $90^\circ$ clockwise about the $z$-axis yields

Since this transformation preserves the relative arrangement of interface types and orientations, both representations correspond to the same tile species $T_1$.

Interactions between interface types are specified by a symmetric binary matrix $A_{ij}$ of dimension $c \times c$,

Each active (non-neutral) interface type is assigned a unique complementary partner such that odd integers bind exclusively to the subsequent even integers: $1 \leftrightarrow2, 3 \leftrightarrow4$, etc. Two faces bind only when their interface types are complementary and their orientations match.

Assembly proceeds through stochastic growth on the cubic lattice (Figure 1(b)). Each assembly run begins with the placement of an initial tile (seed) at the origin. In seeded assembly, the first tile species in the assembly kit is used as the seed, whereas in unseeded assembly the seed tile is chosen uniformly at random from the assembly kit.

Once the seed tile is placed, the growth of the polycube unfolds through the following steps:

1. All exposed faces of placed tiles carrying non-neutral interface types that can bind to at least one interface present in the assembly kit are identified as open interfaces.

2. One open interface is selected uniformly at random.

3. A face carrying the complementary interface type is selected from the assembly kit. When multiple such faces are available, on the same or different tile species, the selection is made uniformly at random.

4. The tile species containing the selected face is rotated into a compatible orientation and attached to the growing assembly.

Here, we assume infinite copies of these finitely many tile species throughout assembly. Steps 1–4 are repeated until no open interfaces remain, in which case growth has terminated, or until the structure exceeds the imposed bounds on its maximum linear extent or total tile count, in which case the simulation is stopped and the assembly is classified as unbounded (Figure 2(a)). If the assembly terminates in a bounded structure, we repeat the assembly process up to $k$ times for the same rule set, each time with different random choices of open interface (step 2) and of the face selected to bind to it (step 3). This tests whether the rule set uniquely determines the final structure despite variations in the order of attachment events. Rule sets for which all trajectories produce geometrically equivalent final structures, up to global cubic rotations and reflections, are classified as deterministic; if different trajectories yield non-equivalent outcomes, the rule set is classified as non-deterministic (see Sec. A of the Supplementary Material (SM) for details).

Together, boundedness and determinism define four possible classes of rule sets: bounded and deterministic (BD), bounded and non-deterministic, unbounded and deterministic, and unbounded and non-deterministic. A bounded and deterministic rule set reproducibly assembles into the same finite structure across independent assembly trajectories. A bounded and non-deterministic rule set can, for example, produce different finite structures along different assembly trajectories, as in multifarious assemblies such as the werewolf in Ref. 24. An unbounded and deterministic rule set will generate the same indefinitely growing structure, such as a periodic crystal. An unbounded and non-deterministic rule set can grow without bound and also yield different structures along different trajectories. We collectively classify these latter three classes as UND and restrict our subsequent analysis to BD rule sets. We refer to the finite structures they deterministically produce as BD assemblies.

Rule-space sampling

The polycube model defines an input–output mapping from assembly rule sets $I$ to assembled polycubes $O$. For BD rule sets, each rule set produces a unique output structure, yielding the deterministic mapping $f : I \to O$. Many different rule sets may produce the same output polycube, resulting in a many-to-one mapping.

For a fixed number of tile species $s$ and interface types $c$, the space of all rule-set encodings forms a discrete input space $S_{s,c}$. We use superscripts $2D$ and $3D$ to distinguish two-dimensional polyomino and three-dimensional polycube rule spaces, respectively. For example, $S_{s,c}^{3D}$ denotes the three-dimensional rule space with $s = 5$ tile species and $c = 31$ interface types, while $S_{s,c}^{2D}$ denotes the corresponding two-dimensional rule space.

In the two-dimensional polyomino model27, each of the four square tile edges is assigned one of $c$ interface types. Since edge interfaces have no rotational degree of freedom, no orientation is required, giving

In the three-dimensional polycube model, each of the six tile faces is assigned one of $c$ interface types and one of four geometrically allowed orientations, giving

In both models, interface interactions follow the fixed pairing convention encoded in $A_{ij}$. Even for modest parameters (e.g., $s = 2$, $c = 4$), the three-dimensional rule space exceeds $10^{14}$ rule sets, rendering exhaustive enumeration impractical. We therefore sample rule sets uniformly at random from $S_{s,c}$ and simulate the assembly process for each.

Characterizing the output structures

To characterize the output landscape generated by the polycube assembly model, each BD assembly is described by three main properties: its assembly rule probability, its geometric symmetry under cubic transformations, and its assembly complexity measured by the information required to specify its minimal assembly kit.

Assembly rule probability

Assembled polycubes are identified by their global geometry, independent of the specific arrangement of tile species. For example, the assemblies made with tile orders $T_1$–$T_2$–$T_1$ and $T_2$–$T_1$–$T_2$ are treated as the same structure because both correspond to the same three-tile linear polycube. Once tile identities are ignored, geometric structures may be identified under three standard equivalence conventions56,57:

• Fixed polycubes: two structures are considered distinct unless they differ only by translation.

• One-sided polycubes: two structures are considered equivalent under translations and rotations, but mirror figures are treated as distinct.

• Free polycubes: two structures are considered equivalent under translations, rotations, and reflections.

Throughout this work, outputs are classified as free polycubes, identifying structures up to all 48 symmetry operations of the cube (Table I). The total number of free polycubes of size $n$ is given by the On-Line Encyclopedia of Integer Sequences (OEIS) sequence A03811958. For two-dimensional assemblies, we similarly classify outputs as free polyominoes, with the total number of free polyominoes of size $n$ given by OEIS sequence A00010558 (see SM Table S1).

Class	Geometric Operation	#Conjugates
$I$	Identity	1
$A$	$\pi/2$ rotation about face axis	6
$B$	$\pi$ rotation about face axis	3
$C$	$\pi$ rotation about edge axis	6
$D$	$2\pi/3$ rotation about vertex axis	8
$E$	Reflection in face plane	3
$F$	Reflection in edge plane	6
$H$	$D \times K$	8
$J$	$A \times E$	6
$K$	Inversion through cube center	1

Given this equivalence criterion, multiple BD rule sets may produce the same geometric structure or satisfy the same structural property. We define the assembly rule probability, $P_{AR}(\cdot)$, as the fraction of BD rule sets that produce a specified outcome; UND rule sets are excluded. For a specific polyomino or polycube structure $p$,

whereas for assemblies of size $n$,

Thus, $P_{AR}(n)$ represents the fraction of all BD rule sets that produce an assembly of size $n$.

For comparisons within a fixed size class $n$, we write the size-normalized assembly rule probability as

where $X$ denotes a set of structures sharing a property, such as symmetry order $SO$, symmetry class $SC$, or assembly complexity $\tilde{K}$. Thus, $P_{AR}(X \mid n)$ is the fraction of BD rule sets producing size-$n$ assemblies that also satisfy the specified property $X$.

In contrast to $P_{AR}$, which is defined over the space of BD rule sets, the morphospace probability is defined over the space of geometrically distinct structures. For structures of size $n$, let $M_n$ denote the theoretical morphospace, i.e., the set of all geometrically distinct free polycubes or polyominoes of size $n$, with $|M_n|$ given by the corresponding OEIS enumeration. We define the size-normalized morphospace probability as

which is the fraction of geometrically distinct free structures of size $n$ that satisfy the specified property or belong to category $X$.

The theoretical morphospace need not coincide with the set of structures observed by finite rule sampling. Some structures in $M_n$ may be inaccessible within the chosen rule space, while others may be accessible but not encountered in a finite sample. We therefore define the sampled morphospace, $\hat{M}_n \subseteq M_n$, as the set of distinct size-$n$ structures observed during rule sampling. The corresponding sampled morphospace probability is

which gives the fraction of observed size-$n$ structures satisfying the specified property $X$. If all theoretically possible size-$n$ structures are discovered, then $\hat{M}_n = M_n$, and the sampled and theoretical morphospace probabilities coincide.

Symmetry classification

The geometric symmetry of each assembled polycube is determined by its invariance under the full cubic symmetry group $O_h$ (denoted $O_3$ in Lunnon’s notation[55]). This group contains 48 symmetry operations, including proper rotations, reflections, inversion, and improper rotations. These operations can be grouped into ten conjugacy classes. Operations belonging to the same conjugacy class represent the same type of geometric transformation but act about symmetry-equivalent axes or planes. Table I lists these ten conjugacy classes together with their geometric interpretation and the number of conjugates (symmetry-equivalent operations in each class). For example, class label $B$ refers to a $\pi$ rotation about a face axis and contains 3 conjugate operations, one about each Cartesian axis.

Each polycube is invariant under a subgroup $H \subseteq O_h$, where $H$ consists of all cubic transformations listed in Table I that leave the polycube unchanged. The number of operations in $H$ is the order of the subgroup, which we refer to as its symmetry order. Thus, a polycube with no non-trivial symmetry has symmetry order 1, corresponding only to the identity operation, whereas a polycube invariant under every operation in $O_h$ has symmetry order 48.

The 98 subgroups of $O_h$ fall into 33 conjugacy classes[55]. Two subgroups are conjugate within $O_h$ if they contain the same geometric arrangement of symmetry operations but differ in the spatial directions along which those operations act. These 33 subgroup conjugacy classes enumerate all possible symmetry types of free polycubes[55], and we refer to them as symmetry classes. Table 2 lists each symmetry class together with its subgroup order, number of conjugate subgroups (#Conj.), corresponding Schoenflies symbol (Sch.), and class structure in terms of the operation classes of $O_h$ from Table I.

Class	Order	#Conj.	Sch.	Class structure
$I$	1	1	$C_1$	$I$
$A$	4	3	$C_4$	$I+B+2A$
$B$	2	3	$C_2$	$I+B$
$C$	2	6	$C_2$	$I+C$
$D$	3	4	$C_3$	$I+2D$
$E$	2	3	$C_s$	$I+E$
$F$	2	6	$C_s$	$I+F$
$H$	6	4	$S_6$	$I+K+2D+2H$
$J$	4	3	$S_4$	$I+B+2J$
$K$	2	1	$C_i$	$I+K$
$BB$	4	1	$D_2$	$I+3B$
$BC$	4	3	$D_2$	$I+B+2C$
$BE$	4	3	$C_{2h}$	$I+B+E+K$
$BF$	4	3	$C_{2v}$	$I+B+2F$
$CE$	4	6	$C_{2v}$	$I+C+E+F$
$CK$	4	6	$C_{2h}$	$I+C+F+K$
$EE$	4	3	$C_{2v}$	$I+B+2E$
$CD$	6	4	$D_3$	$I+2D+3C$
$FF$	6	4	$C_{3v}$	$I+2D+3F$
$AB$	8	3	$D_4$	$I+2A+2C+3B$
$AE$	8	3	$C_{4h}$	$I+B+E+K+2A+2J$
$BFF$	8	3	$D_{2d}$	$I+2F+2J+3B$
$CJ$	8	3	$D_{2d}$	$I+B+2C+2E+2J$
$EEE$	8	1	$D_{2h}$	$I+K+3B+3E$
$EF$	8	3	$C_{4v}$	$I+B+2A+2E+2F$
$EFF$	8	3	$D_{2h}$	$I+B+E+K+2C+2F$
$BD$	12	1	$T$	$I+3B+8D$
$CF$	12	4	$D_{3d}$	$I+K+2D+2H+3C+3F$
$BBC$	16	3	$D_{4h}$	$I+K+2A+2C+2F+2J+3B+3E$
$CCC$	24	1	$T_d$	$I+3B+6F+6J+8D$
$DEE$	24	1	$T_h$	$I+K+3B+3E+8D+8H$
$R$	24	1	$O$	$I+3B+6A+6C+8D$
$G$	48	1	$O_h$	$I+K+3B+3E+6A+6C+6F+6J+8D+8H$

It is important to distinguish between conjugacy classes of individual symmetry operations, listed in Table I, and conjugacy classes of subgroups describing the overall symmetry of a polycube, listed in Table 2. Following Lunnon’s notation[55], we use the same labels ($A$, $B$, $C$, …) for both levels of classification, but denote operation-based conjugacy classes in bold in Table I (e.g., A, B, C) to distinguish them from the subgroup-based conjugacy classes in Table 2. For example, in Table 2, symmetry class $A$ has order 4, three conjugate subgroups, and Schoenflies symbol $C_4$. Its class structure is $I+B+2A$. These operations correspond to rotations by 0, $\pi/2$, $\pi$, and $3\pi/2$ about a common fourfold face axis, i.e., a $C_4$ axis in Schoenflies notation. The three conjugate subgroups arise by placing this fourfold axis along any one of the three mutually perpendicular face-axis directions of the cube, and therefore represent the same symmetry class of a free polycube.

Each assembled polycube is assigned to the symmetry class corresponding to its invariant subgroup $H$. In the analysis, we further represent a symmetry class $SC$ with symmetry order $SO$ as $SC_{SO}$. For example, $A_4$ denotes symmetry class $A$ with symmetry order 4, and $G_{48}$ denotes symmetry class $G$ with symmetry order 48.

Similarly, the symmetries of two-dimensional polyominoes are described by the dihedral group of the square, $D_4$, whose subgroups fall into eight conjugacy classes. Grouping axis-parallel and diagonal mirror lines together, we classify polyominoes into six symmetry types: $C_1$ (no non-trivial symmetry), $C_2$ (twofold rotation), $C_4$ (fourfold rotation), $D_1$ (one mirror line), $D_2$ (two perpendicular mirror lines), and $D_4$ (full square symmetry).

Complexity measures

We quantify assembly complexity using three complementary measures, each inspired by Kolmogorov complexity $K$. The minimal number of distinct tile species is given by $\widetilde{K}_s$. The minimal number of interface types (including the neutral type) is given by $\widetilde{K}_i$. The third measure, $\widetilde{K}_{LZ}$, is obtained by applying a symmetrized version of the Lempel–Ziv 1976 compression algorithm to the binary encoding of the input rule set[59], as done in Ref. 36. Prior to complexity calculation, assembly rules are simplified to remove redundant information. Details of the simplification procedure and complexity calculations are provided in SM Sec. A.

For structures that can be generated by multiple BD rule sets, each complexity measure is defined as the minimum value over all corresponding BD inputs. In practice, we cannot always exhaustively enumerate all possible rules, and so the minimum is taken over all rules that are found by sampling. That implies that our measures may be upper bounds. However, simpler assembly descriptions typically occupy larger regions of rule space because there are more ways to add redundant rules. Therefore they are more likely to appear upon random sampling, and so we expect that for the most frequently assembled structures, observed minima provide a good approximation to the true minima. As a check, we compared the sampled minimum complexities of all tetracubes with those obtained using the SAT-based inverse-design solver of Ref. 24 and found that the two agree.

Assembly simulation parameters

Polycube assembly simulations were performed in the three-dimensional rule space $S^{3D}_{5s,31c}$, consisting of $s=5$ tile species and $c=31$ interface types. Interface type 0 was neutral, while all non-zero interfaces formed unique complementary pairs ($1 \leftrightarrow2$, $3 \leftrightarrow4,\ldots,29 \leftrightarrow30$) according to the interaction matrix $A_{ij}$ defined earlier. The rule space $S^{3D}_{5s,31c}$ contains more than $10^{62}$ rule set encodings, from which $10^9$ were sampled uniformly at random.

Similarly, polyomino assembly simulations were performed in the corresponding two-dimensional rule space $S^{2D}_{5s,31c}$, using the same numbers of tile species and interface types and the same interface-pairing convention. This rule space contains more than $10^{29}$ rule set encodings, from which $10^9$ were likewise sampled uniformly at random.

Unless otherwise specified, simulations in both dimensions were performed in seeded mode, with the first tile species initiating the growth process. To classify assemblies as bounded or unbounded, the maximum linear extent was set to 16 tiles and the maximum total assembly size to 100 tiles. Rule sets producing bounded assemblies were then tested for determinism using up to $k=25$ independent stochastic trajectories. Further details and justification of these simulation parameters are provided in SM Sec. F.

Results and Discussion

Sampling $10^9$ rule sets from $S^{3D}_{5s,31c}$ under seeded assembly yielded $1.7 \times10^7$ (1.7%) BD rule sets, while the remaining 98.3% were unbounded, non-deterministic, or both (UND). Of the BD rule sets, 56% produced only a single cube, and the remaining $7.3 \times10^6$ (0.73% of all sampled rule sets) assembled structures of two or more cubes. Together, the BD rule sets generated $2 \times10^5$ distinct free polycubes, reflecting the many-to-one mapping from rule sets to assembled structures.

For the corresponding two-dimensional system, sampling $10^9$ rule sets from $S^{2D}_{5s,31c}$ under seeded assembly yielded $3.79 \times10^8$ (37.9%) BD rule sets, while 62.1% were UND. Of the BD rule sets, 31% produced only a single tile, and the remaining $2.6 \times10^8$ (26% of all sampled rule sets) assembled polyominoes of two or more tiles. These BD rule sets generated $2.4 \times10^4$ distinct free polyominoes. As shown in Figure 2(c), the fraction of BD rule sets is substantially lower in three dimensions than in two.

Having identified the BD rule sets and their corresponding assemblies, we next examine how assembly rule probability is distributed with respect to assembly size, assembly complexity, and symmetry.

Size distribution of assembled structures

The number of mathematically possible free polycubes grows rapidly with polycube size $n$ (OEIS58 sequence A038119, see SM Table S1). As shown in the inset of Figure 3(a), the number of distinct free polycubes increases from just 2 at size 3 to 112 at size 6, reaching $6.4 \times10^{15}$ at size 22, the largest size for which exact enumeration is available. A conservative extrapolation suggests that by size 100, the cumulative number of free polycubes exceeds $10^{80}$, comparable to the number of atoms in the observable universe. This illustrates that these spaces are hyper-astronomically large[60].

Figure 3(a) shows how much of this enormous morphospace is recovered by random sampling of BD assembly rules. Assemblies range from polycube size 1 to 98, approaching the imposed size cutoff of 100. All 146 possible free polycubes of sizes up to 6 are recovered. Beyond this size, however, the fraction of the theoretical morphospace recovered decreases rapidly with increasing polycube size.

Figure 3(b) shows the assembly rule probability $P_{AR}(n)$, defined as the fraction of BD rule sets that generate polycubes of size $n$. Assembly rule probability generally decreases with increasing polycube size, indicating a strong bias towards smaller structures despite the much greater number of geometrically possible polycubes at larger sizes.

Figure 3 also exhibits pronounced local peaks at polycube sizes that are predominantly multiples of four, both in the number of distinct polycubes recovered and in their assembly rule probability. These peaks are absent from the corresponding distribution of polycubes in the theoretical morphospace shown in the inset of Figure 3(a). We return to the origin of these peaks later in our discussion of assembly complexity and symmetry.

To examine the effect of sampling depth, SM Fig. S1 compares size distributions obtained from $10^6$ to $10^9$ sampled rule sets. Increasing the number of samples reveals additional distinct polycubes, particularly at larger sizes, while the overall distribution of assembly rule probability across sizes changes little. The overall size distributions also show similar patterns under seeded and unseeded assembly conditions and for two-dimensional polyomino assemblies (SM Fig. S2).

Beyond finite sampling, the number of tile species and interface types available in our rule space $S^{3D}_{5s,31c}$ also limits which structures can be assembled, regardless of their existence in the theoretical morphospace. In SM Sec. E, we derive the minimum assembly complexity for deterministic assembly of two simple structures. A linear rod of size $n$ requires $\lceil n/2 \rceil+ 1$ tile species, whereas an L-shape of the same size requires $n - 1$ species, close to the fully addressable limit of $n$. Consequently, a rule space with $s$ tile species cannot deterministically assemble L-shaped structures larger than $s+1$, and therefore cannot recover the complete theoretical morphospace at larger sizes, even with exhaustive sampling.

Assembly rule probability decreases exponentially with assembly complexity

Figure 4 plots assembly rule probability against the three complexity measures introduced in Sec. IID. Single-tile structures are excluded from these plots since they do not involve assembly. Across two and three dimensions, under both seeded and unseeded assembly, the observed patterns are consistent with the simplicity bias predicted by Dingle et al.[36]: structures with high assembly rule probability have low assembly complexity, while structures with high assembly complexity have low assembly rule probability. This bias towards simplicity can be quantified by an upper bound on assembly rule probability that decreases exponentially with assembly complexity, with most of the probability weight concentrated relatively close to the bound on a logarithmic scale (see SM Fig. S8). These trends are observed for all three complexity measures, indicating that each of these measures captures something salient about the amount of information needed to assemble a structure.

The color gradient in Figure 4 highlights that larger polycubes generally occupy the high-complexity, low-probability region of the plots. In other words, they typically require more assembly information to specify their deterministic growth and therefore admit fewer compatible rule sets than smaller polycubes. This difference in assembly complexity helps explain the pattern observed in Figure 3, where larger structures typically arise less frequently under random sampling of assembly rules, even

though there are exponentially more of them in the theoretical morphospace.

However, a closer examination of the color gradient reveals substantial variation in both assembly rule probability and assembly complexity among polycubes and polyominoes of similar size. We next examine how assembly rule probability depends on complexity when assembly size is held fixed.

Assembly rule probability at fixed assembly size

Figure 5(a) and (b) show the size-normalized assembly rule probability against assembly complexity for polycubes of sizes 6 and 20, respectively. At both sizes, assembly rule probability decreases with increasing complexity, with an approximately exponential upper bound, demonstrating that simplicity bias persists at fixed size. The corresponding empirical upper-bound fit for size 20 is shown in SM Figure S10.

The color gradient in Figure 5(a) and (Figure 5b), indicating symmetry order, shows that the low-complexity, high-probability region is predominantly occupied by highly symmetric polycubes. Figure 5(c) shows the 30 most frequently assembled polycubes of size 20. All 30 have symmetry order 4 or higher, indicating that highly symmetric polycubes are strongly favored under random sampling of assembly rules. A particularly striking example is the polycube with full cubic symmetry (order 48), indicated by the red arrows in Figure 5(b) and circled in (c). Among the approximately $10^{14}$ possible free polycubes of size 20, this is the only structure possessing all 48 symmetries of the cube. Yet it is generated by approximately 1 in 30 BD rule sets producing size-20 polycubes. It is therefore more than twelve orders of magnitude more likely to occur under random sampling of assembly rules than under uniform sampling of free polycubes of size 20.

This strong preference for symmetric structures is not limited to a single polycube. At size 20, the 29 most probable polycubes together account for 50% of the assembly rule probability, as indicated by the black dashed lines in Figure 5(b). All 29 possess non-trivial symmetry and belong to just five symmetry classes: $A_{4}$, $AE_{8}$, $BBC_{16}$, $EF_{8}$, and $G_{48}$. Together, these classes represent only about $6 \times10^{-12}$ of the theoretical morphospace (see SM Table S2). Thus, half of the assembly rule probability is concentrated on structures drawn from an extraordinarily small fraction of all possible size-20 polycubes. This demonstrates the strong enrichment of symmetric structures under random rule sampling.

size	$I_{1}$	$A_{4}$	$B_{2}$	$C_{2}$	$D_{3}$	$E_{2}$	$F_{2}$	$H_{6}$	$J_{4}$	$K_{2}$	$BB_{4}$	$BC_{4}$	$BE_{4}$	$BF_{4}$	$CE_{4}$
1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1
4	0	0	0	1	0	1	0	0	0	0	0	0	1	0	0
5	4	0	0	2	0	6	3	0	0	0	0	0	1	0	3
6	46	0	3	5	0	28	6	0	0	1	0	0	5	1	4
7	394	0	4	17	1	126	26	0	0	1	0	0	5	0	12
8	3025	0	37	49	0	520	72	0	0	20	0	0	21	2	13
9	22707	0	52	138	1	2099	241	0	0	17	0	0	26	2	48
10	167732	0	342	374	9	8429	623	0	1	192	1	1	93	12	61
11	1241417	0	502	1062	1	33676	2028	0	1	175	1	2	116	9	186
12	9221624	2	2836	2851	7	135201	5374	1	6	1632	5	7	392	37	231
13	68936674	0	4343	8010	68	543248	16781	0	3	1522	3	5	497	32	727
14	518574100	2	22622	21432	9	2195182	45094	1	10	13088	13	17	1616	134	941
15	3925132946	2	35405	60142	61	8893547	138971	0	12	12339	16	23	2087	114	2855
16	29878869619	22	176176	161386	473	36196788	378196	0	63	102846	58	78	6690	426	3707
17	228629549175	14	281141	453253	81	147739046	1156028	0	41	97374	47	62	8656	375	11293
18	1757697391087	31	1363112	1222110	440	605138811	3183390	10	97	803363	136	181	27493	1393	14909
19	13570818452472	41	2205171	3435654	3316	2480501070	9679303	4	124	760594	181	218	35797	1236	44955
20	105182527335313	213	10527712	9316409	614	10233284681	26929415	6	553	6274638	546	717	113728	4398	59760
21	818093680980786	182	17221126	26231463	3315	42234214910	81575661	1	401	5928747	521	591	147838	3934	180092
22	6383353461322488	321	81462884	71496106	23537	174693725151	228958790	0	821	49133524	1247	1606	469824	13883	241827

A complementary way to illustrate the uneven distribution of assembly rule probability among structures of the same size is through rank–probability plots, in which structures are ordered by decreasing probability. Figure 6 shows $P_{AR}(p \mid n)$ against rank for polycubes and polyominoes of sizes 6, 11, and 20. The corresponding distributions for all sizes up to 20 are shown in SM Figure S11 and Figure S12. At each size, assembly rule probability decreases rapidly with rank, spanning several orders of magnitude. When structures of all sizes are combined, the rank–probability distribution for polycubes approximately follows a power law (SM Figure S11), whereas that for polyominoes appears not to (SM Figure S12). For individual sizes, the distributions exhibit more complex behavior in both systems. These results further demonstrate that even among structures of the same size, the assembly map explores morphospace highly unevenly, with a strong bias towards certain structures over others.

The gray dashed lines in Figure 6 mark the rank at which the cumulative assembly rule probability reaches 50%; further examples are shown in SM Figure S11 and Figure S12. In general, it is remarkable how tiny a fraction of all structures accounts for half of the BD rule sets. It is also interesting to note, as can be seen in SM Figure S11 and Figure S12, that sizes that are multiples of four need the fewest structures to account for 50% of the BD rule sets. This effect is more pronounced for polycubes than for polyominoes. It is likely caused by the fact that symmetries are more efficient at reducing the assembly complexity for structures whose sizes are multiples of four, which are most commensurate with the square or cubic geometry of polyominoes and polycubes. Note that this effect also explains why size 11 has more structures in the top 50% set than size 20 does in Figure 6. We return to assembly rule probabilities for sizes that are multiples of four in Sec. III F.

For classes whose structures are all recovered by our sampling, we mark the rank at which each symmetry class first appears (triangles) and the rank at which all of its possible structures have been found (squares) in Figure 6. At size 6, all symmetry classes are fully recovered in both dimensions. At larger sizes, complete recovery is observed mainly for highly symmetric classes, despite their rarity in the theoretical morphospace.

Together, Figure 5 and Figure 6 show that polycubes with high assembly probability typically have a relatively low complexity, and a relatively high symmetry. Symmetry provides a mechanism for reducing assembly complexity by allowing the same assembly instructions to be reused across symmetry-related regions of a structure. Symmetric structures that require less assembly information may therefore be generated by more rule sets and occur more frequently under random rule sampling.

However, the broad variation in assembly complexity and probability among structures with the same symmetry order in Figure 5, as well as differences between symmetry classes of the same order in Figure 6, suggests a more complex relationship between symmetry, complexity, and assembly rule probability. Before examining this relationship in more detail, we first consider how assembly rule probability is distributed across either complexity or symmetry.

Distribution of assembly rule probability over complexity

So far, we have focused on the assembly rule probabilities of individual structures, finding that highly probable structures often have low assembly complexity and high symmetry. However, the number of possible structures grows exponentially with size (SM Table S1), and larger structures typically have higher assembly complexity, so we expect there to be many more complex structures than simple ones. Even at fixed size, low-symmetry structures greatly outnumber highly symmetric ones (SM Table S2). Since lower symmetry is generally associated with higher assembly complexity, this suggests that complex structures may also be more numerous at a given size. Therefore, although individual complex structures tend to have lower assembly rule probabilities, their much greater number may compensate for this. This raises the following question: Which effect dominates the total assembly rule probability $P_{\mathrm{AR}}(\widetilde{K})$ at a given complexity: the higher individual probabilities of simpler structures or the greater abundance of complex structures?

Arguments from AIT34 suggest that the number of binary strings with Kolmogorov complexity $K$ grows as $2^K$, which would approximately counteract the $2^{-K(p)}$ individual probability for a string $p$ expected from the AIT coding theorem35, resulting in a nearly flat probability distribution over complexity, $P(K) \approx O(1)$. An open question is how well this argument works for computable input-output maps[61] with suitable approximations for complexity $\widetilde{K}$. Interestingly, for deep neural networks on Boolean functions, which can be represented as strings, and where high expressivity means that all functions can be found, it has been shown that the distribution $P(\widetilde{K})$ obtained upon random sampling of the weights (the analogue of sampling over assembly rules here) is nearly flat, at least on a log scale[62].

By contrast, for RNA secondary structures, which can be represented as ternary strings, many low- and high-complexity strings do not correspond to valid folded secondary structures. Consequently, the distribution of complexities in the RNA morphospace does not follow the $2^{\widetilde{K}}$ expectation for simple strings. This is one reason why, for RNA of length $L$, the probability distribution $P(\widetilde{K}\mid L)$ has an approximately Gaussian form, peaking at an intermediate complexity rather than remaining flat[42].

What $P_{AR}(\widetilde{K})$ to expect for polycubes and polyominoes is difficult to predict a priori, because the complexity distribution $P_M(\widetilde{K})$ of the full morphospace is difficult to estimate analytically. However, we expect both very simple and very complex rule sets to be less likely to produce BD outcomes. Thus, $P_{AR}(\widetilde{K})$ may peak at intermediate complexity, as observed for RNA secondary structures, although more complex behavior is also possible.

The left panel of each pair in Figure 7 (hatched bars) shows the complexity distributions of the sampled morphospace $P_{\hat{M}}(\widetilde{K}\mid n)$ for polycubes of sizes 6 and 20, as well as $P_{\hat{M}}(\widetilde{K})$ over all sizes. Corresponding distributions for additional polycube sizes are shown in SM Figure S13. At size 6, all 112 possible free polycubes are recovered, so the sampled morphospace coincides with the theoretical morphospace. At size 20, however, only 4678 distinct polycubes are recovered, compared with approximately $10^{14}$ in the theoretical morphospace. Although not all these structures are accessible in the rule space $S^{3D}_{5s,31c}$, many more should be accessible than are found in our sampling. Since high-complexity structures generally have lower assembly rule probabilities, this means that the sampled morphospace distributions $P_{\hat{M}}(\widetilde{K}\mid20)$ and $P_{\hat{M}}(\widetilde{K})$ will most severely underestimate the higher complexity parts of the full distribution. Nevertheless, even the sampled distributions show that there are typically many more complex structures than simpler ones.

The right panel of each pair in Figure 7 (solid bars) shows the corresponding assembly rule probability distributions $P_{AR}(\tilde{K} \mid n)$ and $P_{AR}(\tilde{K})$. Unlike the morphospace distributions, these can be estimated without exhaustively sampling all possible structures. The required sampling depth is largely determined by the least probable complexity values, which generally receive sufficient samples in our simulations for reasonably accurate estimates. Compared with the sampled morphospace, the assembly rule probability distributions are shifted towards lower complexity, reflecting the strong simplicity bias towards structures requiring less assembly information. The overall structure of $P_{AR}(\tilde{K} \mid n)$ illustrates the general principle that while simplicity bias implies that the most probable individual structures tend to be simple, this does not necessarily mean that random sampling of rules will predominantly produce simple structures, since there are many more low probability complex structures.

Distribution of assembly rule probability over symmetry orders and classes

Just as for assembly complexity, we can examine how assembly rule probability is distributed across different symmetries. Figure 8(a) compares the morphospace probability $P_M(X \mid n)$ with the assembly rule probability

$P_{AR}(X \mid n)$ across symmetry orders ($X = SO$) and symmetry classes ($X = SC$) for polycubes of size $n = 20$. Corresponding distributions for other sizes are shown in SM Figure S14. The theoretical morphospace is overwhelmingly dominated by low-symmetry structures (see SM Table S2). For example, among the approximately $10^{14}$ free polycubes of size 20, nearly all (99.99%) belong to the asymmetric class $I_1$; only $\approx10^{10}$ have higher symmetry. There is only one with full cubic symmetry (class $G_{48}$). But, rather strikingly, under random sampling of BD rule sets its probability is within a factor of about four of that of all $I_1$ structures. This enormous enhancement of class $G_{48}$ is linked to the high symmetry and low assembly complexity of this structure.

We also note that the probabilities of obtaining structures from all symmetry classes are quite similar (on a log scale), see also other polycube sizes in SM Figure S14. The reason for this pattern is that the rarer higher symmetry structures are more enhanced than the more common lower symmetry structures. This composite effect of two separate trends is analogous to what was observed for the $P_{AR}(\widetilde{K})$ distributions.

Figure 8(b) shows the probability enhancement factor $P_{AR}(X \mid n)/P_M(X \mid n)$ (dark green), which measures how much more likely a randomly sampled BD rule set is to produce a structure belonging to a given symmetry order or class than would be expected from uniformly sampling the theoretical morphospace of that size. Only the lowest symmetry $I_1$ structures are underrepresented, whereas all other symmetries found at this size are enhanced, with the enhancement factor generally increasing with symmetry order, reaching over twelve orders of magnitude for the fully symmetric class $G_{48}$. However, symmetry classes of the same order can differ in their enhancement by several orders of magnitude. For example, among order-4 classes, the enhancement ranges from approximately $10^{6}$ for $EE_{4}$ to $10^{11}$ for $A_{4}$, showing that the total number of symmetry operations, or symmetry order, alone cannot explain the large differences in assembly rule probability between symmetry classes.

Figure 8(b) also shows the fraction of theoretically possible structures in each symmetry order or class recovered by our sampling (light green). This coverage is extremely low for asymmetric structures, with only $3 \times10^{-11}$ of class $I_{1}$ recovered, but much higher for strongly enhanced symmetry classes. For example, we recover 71% of class $A_{4}$, 56% of $BBC_{16}$, and all structures belonging to $AB_{8}$, $AE_{8}$ and $G_{48}$. This shows that random sampling of assembly rules strongly favors symmetric structures, recovering a substantial fraction of some highly symmetric classes despite these classes representing only a minuscule fraction of the theoretical morphospace.

The strong enhancement of certain symmetry classes can be understood through their lower assembly complexity. Figure 8(c) shows the mean LZ complexity for each symmetry order and class, averaged either over BD rule sets $\langle\widetilde{K}_{LZ} \rangle_{AR}$ (dark blue) or over distinct sampled structures $\langle\widetilde{K}_{LZ} \rangle_{\hat{M}}$ (light blue). The most strongly enhanced classes, like $A_{4}$, $AB_{8}$, $EF_{8}$, $AE_{8}$, $BBC_{16}$ and $G_{48}$ have lower mean complexity. SM Figure S16 extends this analysis to show the range of LZ complexities and the complexity of the most frequently assembled structure in each symmetry order and class. Corresponding results for unseeded assembly are shown in SM Figure S17.

However, the relationship between symmetry and assembly complexity is not a strong or necessarily straight-forward one when averaged over all structures with a particular symmetry. As discussed above, symmetry classes with the same order can differ substantially in both mean complexity and probability enhancement. Moreover, as we shall see below, even within a single symmetry class, structures can vary significantly in complexity and assembly rule probability. For nearly every symmetry class, $\langle\widetilde{K}_{LZ} \rangle_{AR}$ is lower than $\langle\widetilde{K}_{LZ} \rangle_{\hat{M}}$, indicating that simpler structures within the class account for a greater share of its assembly rule probability. Understanding these patterns requires examining in more detail how global symmetry relates to assembly complexity across structures.

Why some symmetries reduce assembly complexity more effectively

The extent to which symmetry can reduce the amount of information needed to specify the assembly rules depends on whether they permit the same local interactions to be reused during deterministic assembly. That may not always directly correlate with the number of possible symmetry operations.

For example, rotations and reflections can differ in how effectively they allow local assembly rules to be reused. In our model, tile species can be rotated but not reflected during assembly. Proper rotations can therefore allow repeated motifs to be generated using rotated copies of the same tile species and interactions, whereas improper operations, such as reflections and rotoreflections, generally do not permit the same reuse without additional constraints on the interface patterns. Consider symmetry classes $A_{4}$ and $J_{4}$ (Table II). Both contain four symmetry operations, but $A_{4}$ is generated by a fourfold proper rotation, whereas $J_{4}$ is generated by a fourfold rotoreflection. This difference helps explain why structures belonging to $A_{4}$ generally have lower assembly complexity and much higher assembly rule probability than those belonging to $J_{4}$ (Figure 8). Among the order-4 symmetry classes, $A_{4}$ has the fewest structures in the theoretical morphospace of size-20 polycubes, yet the highest assembly rule probability. Of the 30 most frequently assembled size-20 polycubes shown in Figure 5(c), 18 belong to this class.

However, proper rotational symmetry does not guarantee efficient reuse of assembly rules. If a symmetry operation exchanges two cubes across their shared bond, assigning the same tile species to both positions would place identical interface types on the two contacting faces. Since our model allows binding only between distinct complementary interface types, such a bond cannot form. These constraints can limit the reduction in assembly complexity even when the structure possesses several proper rotational symmetries. This is relevant to the symmetry classes $BB_{4}$ and $BC_{4}$, which contain only proper rotations but show much smaller probability enhancements than $A_{4}$. Similarly, for planar polycubes, reflection through the plane containing all cube centers leaves every cube in place. Although this reflection contributes to the global symmetry order, it provides no additional opportunities for tile-species reuse.

We now return to the peaks in assembly rule probability at polycube sizes that are multiples of four, observed in Figure 3. They arise in particular from the strong preference for structures possessing fourfold rotational symmetry, particularly class $A_{4}$. More generally, the relatively enhanced preference for higher symmetry BD structures is also reflected in an elevated mean symmetry order (SM Figure S3) for these sizes. More specifically, when a fourfold rotation axis passes between cubes rather than through them, the cubes occur in groups of four related by rotation. The same tile species and interface types can be reused across these groups, reducing their assembly complexity. Structures that can take advantage of this have sizes that are multiples of four, and this explains the peaks in Figure 3.

Why structures with the same symmetry class can vary significantly in their assembly complexity and probability

The preceding discussion illustrates how different symmetry operations can vary in how efficiently they reduce assembly complexity. However, even structures belonging to the same symmetry class can require substantially different amounts of assembly information.

A striking example is provided by the seven free tetracubes shown in Figure 9(a). The square and the linear rod (outlined in orange) both belong to symmetry class $BBC_{16}$ and therefore possess the same global symmetry. Nevertheless, deterministic assembly of the rod requires at least three tile species and five interface types, whereas the square requires only one species and three interface types. Consistent with this difference, the square has the highest assembly rule probability among tetracubes, while the rod has the lowest.

The difference arises from variations in how their symmetry operations act on the individual cubes. In the rod, the fourfold rotation axis passes through every cube, leaving their positions unchanged. By contrast, the fourfold axis of the square passes through its empty center, mapping each cube to the other three positions. The square can therefore reuse the same tile species at all four positions, whereas the rod cannot exploit its fourfold symmetry in this way. Furthermore, the rod’s central bond prevents its two central cubes from sharing a species under the complementary binding rules discussed previously (see SM Sec. E for an analytic calculation of the assembly complexity of rods).

Another example of variation within the same symmetry class can be seen among size-20 polycubes belonging to class $A_{4}$. Of the 213 theoretically possible structures in this class, our sampling recovers 151. Although all possess the same fourfold rotational symmetry, their assembly rule probabilities span almost three orders of magnitude (Figure 9(b)). This variation is accompanied by quite substantial differences in their assembly complexity (Figure 9(c)). This illustrates how structures with identical global symmetry operations can require different amounts of assembly information.

Conversely, local rule reuse does not always require global symmetry. Figure 9(d) illustrates this using polycubes consisting of a central cube with single-cube arms. The T-tetracube, tripod, flat cross, and three-dimensional cross have symmetry orders 4, 6, 16, and 48, respectively, yet all can be assembled using just two tile species: one for the central cube and one for the arms. For the T-tetracube, global symmetry does not relate the stem to the side arms, but all three can share a species because they form the same local connection to the central cube. Similarly, the non-planar four-armed and five-armed hubs require only two species, even though their global symmetries would suggest three (the three colors in Figure 9(d)). These examples show that local rule reuse can reduce assembly complexity beyond what global symmetry would predict.

We have given a few examples showing that local rule reuse can depend on other structural properties beyond symmetry, such as compactness and the constraints of deterministic growth. Predicting more generally which features besides symmetry influence assembly complexity and assembly rule probability, and how their relative importance varies with structure size, remains an important direction for future work.

Conclusion

We studied how assembly complexity and symmetry shape the distribution of structures generated by a three-dimensional polycube model of programmable self-assembly. By randomly sampling assembly rules and analyzing the resulting bounded-deterministic structures, we find that polycubes with high assembly probability have low assembly complexity, consistent with the exponential enhancement predicted by simplicity bias for input-output maps[36, 37]. This extends earlier findings for two-dimensional polyominoes28 to a substantially larger and geometrically richer three-dimensional assembly space.

A striking finding is that the set of structures with the highest assembly rule probabilities is hugely enriched in high symmetry structures. For example, the size-20 polycube with full cubic symmetry is generated more than twelve orders of magnitude more frequently than would be expected from its abundance in morphospace. For these structures, symmetry provides an important mechanism for reducing assembly complexity, as symmetric structures can reuse tile species and interfaces, allowing them to be specified by shorter assembly rules. As a result, they can occupy disproportionately large regions of rule space.

Nevertheless, assembly complexity is underdetermined by symmetry. Structures with the same symmetry order, or even belonging to the same symmetry class, can require very different amounts of assembly information. This is because not all symmetry operations allow tile species and interfaces to be reused equally effectively under the local binding rules. Moreover, repeated local patterns can reduce assembly complexity even when the corresponding positions are not related by global symmetry.

Our results suggest that the application of Cartwright and Mackay’s31 idea that assembly complexity provides a more general measure of structural complexity than symmetry alone is not so straightforward for finite-sized assemblies. Their main application was to extended structures, which would be unbounded deterministic structures in our language of assembly, and it would be interesting to study these in more detail for polycubes and polyominoes. One possible application would be to the protein crystals that create surface layers coating bacteria and archaea[63].

Future work should also investigate which structural features beyond global symmetry allow tile species and interfaces to be reused. These may include local structural motifs, compactness, and the number of possible assembly pathways, which could be studied using graph-theoretic approaches[64]. Another important direction is to extend our analysis beyond the cubic lattice to other particle geometries and continuous space, where different interaction constraints may influence the relationship between symmetry and assembly complexity.

Understanding the links between geometry, assembly complexity, and assembly rule probability may help guide the design of programmable nanostructures, colloidal assemblies, and other engineered self-assembling systems where minimizing design complexity is both scientifically and technologically desirable. Recent studies of programmable shells have similarly shown that the arrangement of symmetry elements can influence the trade-off between design complexity and assembly fidelity[65]. More broadly, relaxing strict deterministic assembly constraints may offer further opportunities to reduce the number of distinct building blocks, as explored in semiaddressable self-assembly25. These findings may also have implications for evolution, where biases towards structures requiring simpler assembly rules could contribute to non-adaptive explanations for the prevalence of symmetry and other structural features in protein complexes28.

Acknowledgments

We thank John Mason for helpful communications about and contributions to the OEIS enumeration of free polycubes by symmetry class, and Abdurrahman Hadi Erturk and Andrew Turberfield for helpful discussions.

Data Availability

The code and data supporting this study are openly available through Zenodo at https://doi.org/10.5281/zenodo.22045750. Figs. 4 and 5 can be interactively visualized at https://akodiat.codeberg.page/polycubeViewer/.

Determinism testing and complexity calculations

Determinism testing

A rule set is classified as deterministic if all tested trajectories terminate in the same free polycube, i.e., if the resulting structures are equivalent under the operations of $O_h$. Final structures are compared as free polycubes, i.e. including reflections, although tiles only rotate. This is justified because mirror-image structures require the same assembly information: any rule set assembling a structure has a mirror-image counterpart, with the same numbers of tile species and interface types, that assembles its mirror image.

To test determinism, we assemble the rule set along $k = 25$ growth trajectories, generated by different random choices of open and complementary interfaces. The possible growth trajectories of a rule set can be represented by a choice tree, in which branching corresponds to alternative choices of open and complementary interfaces during assembly. For seeded assembly, we construct this tree and explore all possible trajectories when their number does not exceed $k$. Otherwise, the test is limited to $k$ sampled trajectories. For unseeded assembly, the choice of seed species introduces an additional branching level. We therefore do not construct the choice tree and instead perform $k$ independent assembly runs, selecting the seed species randomly in each run.

When only a subset of trajectories is sampled, non-determinism can be detected but cannot be excluded. Consequently, some rule sets classified as bounded-deterministic (BD) may admit unsampled trajectories leading to different structures. The sensitivity of the classification to the simulation parameters is examined in Sec. F.

Simplification of assembly rules and complexity calculations

Randomly sampled rule sets can contain information that does not affect the assembly process, such as non-neutral (active) interface types whose complementary partner is absent from the assembly kit. Retaining such redundant information would inflate the measured complexity without changing the assembled structure. We therefore calculate complexities from a reduced and simplified rule set, constructed separately for each of the trajectories used in the determinism test.

Different trajectories of the same BD rule set can use different subsets of the tile species and interface types present in the assembly kit, even when they produce the same final structure. For example, an interface type may have its complementary partner present in the kit but never participate in a bond because the corresponding attachment site becomes sterically blocked before that partner can attach. For each trajectory, we therefore retain only the tile species incorporated into the assembly, listed in order of first incorporation starting with the seed, and the interface types that form at least one bond; all other faces of the retained tiles are set to the neutral type. The reduced rule set is then simplified as follows:

1. Label 0 denotes the neutral interface. When the number of interface types $c$ is even, the final interface type $c - 1$ is also neutral under our interaction convention and is therefore relabeled as 0. For the odd values of $c$ used in the simulations reported here, this step has no effect.

2. Every active interface type whose complementary partner is absent from the assembly kit is set to the neutral type 0, since it cannot participate in a bond. The orientation of every neutral face, including faces that were neutral before simplification, is set to a fixed reference orientation, 0 for the N, E, S, and W faces and U for the B and F faces, so that neutral faces carry no orientational information.

3. The remaining active interface types are relabeled consecutively, reading the tiles in order of first incorporation and the faces of each tile in NESWBF order. When an interface type is first encountered, it is assigned the next available odd label and its complementary partner the following even label, so that the complementary-pair structure of $A_{ij}$ is preserved.

For example, with orientations omitted for clarity, suppose the incorporated tiles of a rule set in $S_{2s,5c}^{3D}$ carry the interface sequences $\{343002, 002434\}$, where $1 \leftrightarrow2$ and $3 \leftrightarrow4$. Types 3 and 4 form bonds, whereas type 2 cannot, because its partner 1 is absent. Setting type 2 to 0 gives $\{343000, 000434\}$, which becomes $\{121000, 000212\}$ after consecutive labeling.

From the simplified rule set of each trajectory, $\widetilde{K}_{s}$ is the number of tile species and $\widetilde{K}_{i}$ the number of interface types, including the neutral type. Since active interfaces are retained only in bonded complementary pairs, $\widetilde{K}_{i} = 1 + 2m$, where $m$ is the number of complementary pairs used.

To calculate $\widetilde{K}_{LZ}$, tile species whose faces are all neutral are first removed, since they carry no information relevant to assembly. The simplified rule set is then converted to a binary string by listing the tile species in order of first incorporation and the six faces of each tile in NESWBF order. Each interface label is represented using $\lceil\log_{2} c \rceil$ bits, where $c$ is the number of interface types in the rule space, and each orientation is represented using three bits: $O \to000$, $U \to001$, $I \to010$, $D \to011$, $L \to100$, and $R \to101$. For example, in a polycube rule space with $c = 5$ interface types, each face requires $3 + 3 = 6$ bits and each tile 36 bits. For 2D polyominoes, orientations are omitted and only interface types are converted to binary. Following Dingle et al.1, for a binary string $x$ of length $n$,

where $N_{w}$ is the number of words in the Lempel–Ziv 1976 parsing2. Evaluating the complexity for both the forward and reversed strings removes the dependence on reading direction. Unlike $\widetilde{K}_{s}$ and $\widetilde{K}_{i}$, $\widetilde{K}_{LZ}$ depends on the order in which tile species are incorporated and on the rotational representative of each tile, so it can differ between trajectories that use the same species.

For a given BD rule set, each complexity measure is taken to be its minimum over the $k$ tested trajectories. For a given structure, $\widetilde{K}_{s}$, $\widetilde{K}_{i}$, and $\widetilde{K}_{LZ}$ are then each defined as the minimum over all sampled BD rule sets that generate it. Each reported complexity thus involves two successive minimizations: first over the tested trajectories of each BD rule set, and then over all sampled BD rule sets generating the structure.

Polycube size distributions

Dependence on sampling depth

Figure 3 in the main text shows the size distribution of BD polycubes obtained from $10^{9}$ rule sets sampled from $S^{3D}_{5s,31c}$. To assess how this distribution depends on sampling depth, Fig. S1 compares results obtained from $10^{6}$ to $10^{9}$ sampled rule sets. Increasing the sampling depth recovers progressively more distinct polycubes (Figure S1). The gain is most prominent at large polycube sizes, where structures typically have low assembly rule probability and therefore require more samples to be found. Complete recovery of the theoretical morphospace is achieved only at small sizes: all 112 free polycubes of size 6 are recovered at sampling depths of $10^{8}$ and $10^{9}$. Figure S1(b) shows the distribution of assembly rule probability across polycube sizes at different sampling depths. Across all panels, small polycubes account for the largest fraction of BD rule sets, followed by a rapid decrease in assembly rule probability with increasing size, with characteristic local peaks at particular sizes (typically multiples of four). Thus, although deeper sampling continues to uncover additional rare structures, the qualitative form of the assembly rule probability distribution is robust to sampling depth.

3D polycubes
Size	Number of polycubes	Size	Number of polycubes
1	1	12	9,371,094
2	1	13	69,513,546
3	2	14	520,878,101
4	7	15	3,934,285,874
5	23	16	29,915,913,663
6	112	17	228,779,330,204
7	607	18	1,758,309,223,457
8	3,811	19	13,573,319,825,615
9	25,413	20	105,192,814,197,984
10	178,083	21	818,136,047,201,932
11	1,279,537	22	6,383,528,588,447,574

2D polyominoes
Size	Number of polyominoes	Size	Number of polyominoes
1	1	12	63,600
2	1	13	238,591
3	2	14	901,971
4	5	15	3,426,576
5	12	16	13,079,255
6	35	17	50,107,909
7	108	18	192,622,052
8	369	19	742,624,232
9	1,285	20	2,870,671,950
10	4,655	21	11,123,060,678
11	17,073	22	43,191,857,688

Dependence on dimension and seeding

Figure S2 compares the number of distinct structures obtained in two and three dimensions under seeded and unseeded assembly. In all four cases, the number of distinct assembled structures initially increases with size, reaches a maximum at intermediate sizes, and then decreases, despite the rapid growth of the theoretical morphospace with size (see Table S1). Consequently, an increasingly small fraction of the available morphospace is recovered at larger sizes.

Seeded assembly produces more distinct structures than unseeded assembly in both dimensions, consistent with the larger fraction of BD rule sets under seeded growth (main-text Figure 2(c)). The three-dimensional system generates more distinct structures at larger sizes than the two-dimensional system, but recovers a substantially smaller fraction of the corresponding theoretical morphospace. The distributions also exhibit local peaks at particular structure sizes, several of which occur at multiples of four.

Size and symmetry

Although assembly rule probability $P_{AR}(n)$ decreases with size $n$ overall, this decrease is punctuated by pronounced local peaks (main-text Figure 3). As discussed in the main text, a natural explanation is that highly symmetric structures, which can have low assembly complexity, are available at these sizes. Figure S3 compares $P_{AR}(n)$ with the mean symmetry order of the assembled polycubes, weighted by the number of BD rule sets generating each structure.

Several local peaks in $P_{AR}(n)$, particularly at larger sizes, coincide with elevated mean symmetry order. At large sizes, however, few structures are sampled, so the mean symmetry order at a given size can be dominated by a small number of structures.

We also resolve the size distributions by symmetry classes and orders. We write $P_M(n \mid SC)$ and $P_{AR}(n \mid SC)$ for the morphospace and assembly rule probabilities of size $n$ within symmetry class $SC$, and $P_M(n \mid SO)$ and $P_{AR}(n \mid SO)$ for the corresponding probabilities within symmetry order $SO$. The morphospace probabilities are calculated from the OEIS enumerations in Table S2. Since these are available for all 33 symmetry classes only up to size 22, both distributions are restricted to $n \le22$ and normalized over this range.

Figure S4 shows these distributions for each symmetry class. The morphospace distributions vary substantially between classes. Low-symmetry classes are generally increasingly dominated by larger polycubes, reflecting the rapid growth in the number of geometrically possible structures with size. In contrast, several higher-symmetry classes, such as $H_6$, $AB_8$, and $CCC_{24}$, are only possible at particular polycube sizes. Three classes, $BD_{12}$, $R_{24}$, and $DEE_{24}$, have no possible structures up to size 22. These features reflect geometric constraints on the sizes at which particular symmetries can be realized. The distributions obtained by random sampling of assembly rules (red) differ markedly from those of the underlying morphospace (blue). BD assemblies are generally biased towards smaller polycubes, while several symmetry classes, including $A_4$, $J_4$, $FF_6$, and $BBC_{16}$, show enhanced probabilities at particular sizes.

Figure S5 shows the corresponding distributions grouped by symmetry order. At low orders, the morphospace is increasingly dominated by large sizes, while BD assemblies are biased towards small sizes. At high orders, only certain sizes support the required symmetries, so both distributions are concentrated at a few sizes. In addition, the BD distributions for orders 4, 8 and 16 show peaks at sizes that are multiples of four, and those for order 3 show peaks at multiples of three, whereas the morphospace distributions show no such peaks. At these sizes, a fourfold or threefold rotation axis can pass between cubes rather than through them, so every cube belongs to a group of four or three symmetry-related positions that can share a tile species. At other sizes, some cubes must lie on the axis and gain no such reuse, so the corresponding structures require more assembly information and are generated less often.

Morphospace of size-6 polycubes and polyominoes

Size 6 is the largest size at which the full morphospace is recovered by sampling $10^{9}$ rule sets in both two- and three-dimensional rule spaces, $S^{2D}_{5s,31c}$ and $S^{3D}_{5s,31c}$, respectively. Figure S6 and Figure S7 show all 112 free polycubes and all 35 free polyominoes of this size, respectively, ordered by decreasing assembly rule probability and labeled by their complexities and symmetries.

Assembly rule probability versus assembly complexity

Dingle et al.¹ showed that, for a broad class of computable input–output maps, the probability $P(p)$ of obtaining an output $p$ upon random sampling of inputs is bounded from above by a function that decreases exponentially with the complexity of $p$, $P(p) \lesssim2^{-a\widetilde{K}(p)-b}$, where $\widetilde{K}(p)$ is a computable approximation to the Kolmogorov complexity and $a > 0$ and $b$ are map-dependent constants. Moreover, most of the probability weight is expected to lie relatively close to this bound⁴. Here we examine both properties for the polycube self-assembly map.

Empirical upper bounds

Figure S8 shows $\log P_{AR}(p)$ against species complexity $\widetilde{K}_s$, interface complexity $\widetilde{K}_i$, and Lempel–Ziv complexity $\widetilde{K}_{LZ}$ for three-dimensional (top two rows) and two-dimensional (bottom two rows) assembly, each under seeded and unseeded growth. In all panels, the highest probability at each complexity decreases approximately linearly with complexity on the logarithmic scale, consistent with an exponential upper bound on $P_{AR}$. The dashed lines are fitted upper bounds that lie on or above all plotted points and minimize the mean vertical distance to them. Structures generated by fewer than three BD rule sets are omitted since their probabilities cannot be estimated reliably from so few occurrences.

Distribution of probability weight below the bound

To see how probability is distributed below the bound, we sort the structures in Figure S8 at each complexity value by decreasing $P_{AR}(p)$ and color in blue the most probable structures that together account for at least 50% of the probability weight at that value; the remaining structures are shown in yellow. At most complexities, the blue points lie close to the bound (on a log scale), while the more numerous yellow points extend over several orders of magnitude below it. Half of the probability weight is therefore carried by a few structures near the bound.

Dependence on sampling depth

Figure S9 shows $\log P_{AR}(p)$ against complexity for $10^{6}$, $10^{7}$, $10^{8}$, and $10^{9}$ sampled rule sets in $S^{3D}_{5s,31c}$ under seeded assembly. The main effect of shallower sampling is a shortening of the low-probability, high-complexity tail of the distribution, for two reasons. First, $\widetilde{K}_s$, $\widetilde{K}_i$, and $\widetilde{K}_{LZ}$ are each estimated as the minimum complexity among the sampled BD rule sets generating a given structure, so fewer samples make it less likely that the true minimal rule set for a given structure has been found, inflating its apparent complexity. Second, and more importantly, a structure observed only a few times at low sampling depth is intrinsically hard to assign a reliable $P_{AR}(p)$, and such structures are systematically under-represented until sampling is deep enough to resolve them (structures generated by fewer than three BD rule sets are excluded from each panel for this reason).

Simplicity bias at fixed polycube size

Larger structures generally have higher assembly complexity and lower assembly rule probability (main-text Figure 4), so the observed relationship between complexity and probability may partly reflect differences in structure size. To separate these effects, Figure S10 plots the logarithm of the size-normalized assembly rule probability, $\log P_{AR}(p \mid n)$, against assembly complexity for polycubes of fixed size $n = 20$ in $S^{3D}_{5s,31c}$ under seeded assembly. The upper envelope decreases approximately linearly with complexity on the logarithmic scale, while the blue points, which together account for half of the assembly rule probability at each complexity, generally lie close to this envelope. Thus, the approximately exponential upper bound characteristic of simplicity bias persists even when structure size is held fixed.

Rank, complexity and symmetry distributions for assemblies of fixed size

Rank distributions

To examine how assembly rule probability varies among structures of the same size $n$, we rank the distinct structures by their size-normalized assembly rule probability $P_{AR}(p \mid n)$. Figure S11 shows the resulting rank distributions for polycubes of sizes 2–20 in $\mathcal{S}_{5s,31c}^{3D}$, and Figure S12 shows the corresponding distributions for polyominoes in $\mathcal{S}_{5s,31c}^{2D}$. The last panel of each figure combines all sizes, with $P_{AR}(p)$ normalized over all BD rule sets except those producing a single-tile assembly.

At sizes $n \le6$ the distributions contain every possible structure from the theoretical morphospace; at larger sizes only a decreasing fraction of the morphospace is recovered. For $n \ge6$, $P_{AR}(p \mid n)$ decreases rapidly with rank and spans several orders of magnitude, so some structures are generated by many more BD rule sets than others and the assembly map does not sample the morphospace uniformly.

To characterize the decay, we compare it with power laws $P(r) = r^{-\alpha}/H_{N,\alpha}$, where $r$ is the rank, $N$ is the number of distinct structures in the panel, and $H_{N,\alpha} = \sum_{r=1}^{N} r^{-\alpha}$ is a normalization factor ensuring that the reference probabilities sum to one. For $\alpha= 1$, this gives Zipf’s law, $P(r) = 1/(rH_N)$, shown as the gray dashed line. For the last panel in Figure S11, we also show a steeper power law with $\alpha= 1.3$ as a black solid line.

At fixed size, the distributions are not described by a single power law. The flat plateaus at the lowest probabilities correspond to structures generated by only one or a few BD rule sets and reflect the resolution of finite sampling rather than the true form of the distribution. When all polycube sizes are combined, the rank distribution follows an approximate power law with $\alpha\approx1.3$ over about five decades, decaying faster than Zipf’s law.

The red vertical lines mark the rank at which the cumulative assembly rule probability reaches 50%, that is, the number of most probable structures that together account for half of the BD rule sets producing structures of that size. This number is strikingly small at sizes that are multiples of four. For polycubes, half of the BD rule sets at sizes $n = 8$, 12, 16 and 20 are carried by fewer than 30 structures, compared with hundreds to thousands at sizes 9–19 that are not multiples of four. As discussed in main-text Sec. III F, this reflects the efficient reuse of tile species in structures with fourfold rotational symmetry, which is possible when the rotation axis passes between cubes, and hence when the size is a multiple of four. The same pattern is present but weaker for polyominoes, where the number of structures needed is small at most sizes.

Assembly complexity distributions

Figure S13 shows assembly complexity distributions for polycubes of sizes 6, 8, 10, 11, 12, and 20; sizes 6 and 20 are also shown in the top and middle rows of main-text Figure 7. Results are obtained from $10^9$ randomly sampled rule sets in $\mathcal{S}_{5s,31c}^{3D}$ under seeded assembly. For each size, the sampled morphospace probability $P_{\hat{M}}(\tilde{K} \mid n)$ is the fraction of distinct recovered polycubes in each complexity bin, so every structure counts once, whereas the assembly rule probability $P_{AR}(\tilde{K} \mid n)$ is the fraction of BD rule sets producing polycubes in each bin, so every structure counts as often as it is generated. The three column groups show species complexity $\tilde{K}_s$, interface complexity $\tilde{K}_i$, and Lempel–Ziv complexity $\tilde{K}_{LZ}$, the last in bins of width 5. At size 6 all 112 polycubes are recovered, so $P_{\hat{M}}$ equals the theoretical morphospace distribution; at larger sizes it describes only the structures found.

For every size and every complexity measure, $P_{AR}(\tilde{K} \mid n)$ is shifted towards lower complexity than $P_{\hat{M}}(\tilde{K} \mid n)$. At size 6, for example, $P_{\hat{M}}(\tilde{K}_s \mid n)$ peaks at $\tilde{K}_s = 4$ whereas $P_{AR}(\tilde{K}_s \mid n)$ peaks at $\tilde{K}_s = 3$; moreover, $\tilde{K}_s = 5$ accounts for about 10% of the structures but less than 1% of the BD rule sets. At size 20, most assembled structures have the maximum species complexity, $\tilde{K}_s = 5$. Complex structures therefore greatly outnumber simple ones, but simple structures are generated by many more BD rule sets than complex ones. Since sampling also misses rare high-complexity structures, the true morphospace is expected to be at least as skewed towards high complexity as the sampled one, so the shift of $P_{AR}$ relative to the full morphospace is expected to be at least as large as shown.

Symmetry distributions

Main-text Figure 8(a) compares morphospace and assembly rule probability distributions across symmetry orders and symmetry classes for polycubes of size 20. Figure S14 extends this comparison to polycube sizes 15–20. Here $P_M$ (light bars) is the theoretical morphospace probability, calculated from the OEIS enumerations (Table S2), and $P_{AR}$ (dark bars) is obtained from $10^9$ randomly sampled rule sets under seeded assembly.

Across these sizes, $P_M(SO \mid n)$ falls by many orders of magnitude from order 1 to the highest orders, reflecting the scarcity of highly symmetric structures in the morphospace. In contrast, higher symmetry orders account for a disproportionately large fraction of the assembly rule probability $P_{AR}(SO \mid n)$. This enrichment is particularly pronounced for symmetry order 48: at sizes where such structures occur (i.e., 18–20 here), they constitute only a small fraction of the morphospace but account for a substantially larger fraction of the assembly rule probability. Highly symmetric structures are therefore strongly enriched under random sampling of assembly rules.

Resolving the distributions by symmetry class reveals variation that symmetry order alone does not capture. Classes of the same order can differ in $P_{AR}$ by several orders of magnitude: at size 20, for example, $A_4$ has by far the highest probability among the order-4 classes, whereas $J_4$ lies several orders of magnitude lower. The relative probabilities of classes also change with size. Class $A_4$ receives a large share of the BD rule sets at sizes 16 and 20 but only a very small share at size 19, and none of its possible structures of size 15 was found in our sampling. At sizes that are not multiples of four, a fourfold rotation axis must pass through some cubes, which gain no reuse from the symmetry. How much a symmetry reduces assembly complexity therefore depends not only on the symmetry itself, but also on whether the size and topology of the structure allows every cube to share its tile species with symmetry-related copies.

Figure S15 shows the dependence of symmetry order distribution on sampling depth for sizes 6, 8, 10, 11, 12, and 20, from $10^5$ to $10^9$ sampled rule sets. Symmetry order is fixed by the geometry of a polycube once it is identified, so, unlike minimal assembly complexity, it does not change when further rule sets generating the same structure are found. The sampled morphospace distribution $\hat{P}_M(SO \mid n)$ (left) nevertheless changes strongly with depth at larger sizes, because deeper sampling keeps recovering additional, mostly low-symmetry, structures. By contrast, $P_{AR}(SO \mid n)$ (right) is essentially unchanged at the deepest depths, so $10^9$ sampled rule sets are sufficient to determine it.

Figure S16 extends the symmetry analysis of main-text Figure 8 by showing the minimum and maximum Lempel–Ziv complexities and the complexity of the most frequently assembled structure within each symmetry order and class. The corresponding results under unseeded assembly are shown in Figure S17.

Calculating assembly complexity of simple structures

In this section we derive the minimal species complexity $\widetilde{K}_{s}(N)$ and interface complexity $\widetilde{K}_{i}(N)$ (including neutral interface type) for deterministic assembly of two simple structures of size $N$: a straight rod, and an L-shape assembly obtained by moving one end cube of the rod perpendicularly. In the main text these quantities are estimated as minima over sampled rule sets, so they are upper bounds on the true minima; here we obtain the exact values.

Setting. As in the main text, each non-neutral (active) interface type binds only to its complementary partner, which is distinct from itself. Tile species are identified up to the 24 proper rotations of the cube (classes $I$–$D$ in main-text Table I). Orientations play no role here: orientation matching only selects the rotation of an attaching cubic tile about the bond axis and never prevents two complementary faces from binding. Both structures considered here are planar, so the results hold in both the two- and three-dimensional models, and the arguments apply to both seeded and unseeded assembly.

Both structures are paths of $N$ tiles, in which each tile is bonded to its neighbors along the path and to nothing else. For growth to stop at the target structure, every face that does not form a bond must be neutral. Each end tile therefore has one active face and each interior tile has two, so an end tile can never be the same species as an interior tile.

Linear assembly

Proof. Fully addressable starting point. Giving each of the $N-1$ bonds its own complementary pair, and setting all other faces to neutral interface type 0, makes every tile a distinct species. This construction is deterministic, but it uses $2(N-1)+1$ interface types and $N$ species, and ignores the symmetry of the rod.

Using the symmetry. A rod is unchanged by a half-turn about an axis through its midpoint, perpendicular to the rod. This half-turn swaps the tile at position $p$ with the tile at position $N+1-p$, and swaps the corresponding bonds. Since tiles can rotate, the two tiles of each symmetry-related pair can be the same species, and the two corresponding bonds can use the same complementary pair.

Odd $N$. The central tile is mapped onto itself, and the $N-1$ bonds form $(N-1)/2$ symmetry-related pairs. The rod can therefore be assembled using $(N-1)/2$ complementary pairs, that is, $N$ interface types including the neutral type, and $(N+1)/2$ species. For example, for $N=5$ (Fig. S18(b)), with $a\leftrightarrow b$ and $c\leftrightarrow d$, and listing only the west and east faces,

Here $T_{5}$ is a rotation of $T_{1}$ and $T_{4}$ a rotation of $T_{2}$. The central tile carries $d$ on both faces, which works because both its neighbors bind to it through $c$.

Even $N$. Now the central bond, between the two middle tiles, is mapped onto itself. It needs its own complementary pair, with one member on each of the two middle tiles, so these two tiles are different species. The rod can therefore be assembled using $N/2$ complementary pairs, that is, $N+1$ interface types, and $N/2+1$ species. For example, for $N=4$ (Fig. S18(a)),

where $T_{4}$ is a rotation of $T_{1}$, while $T_{2}$ and $T_{3}$ are distinct.

In both cases, growth is deterministic. The interface types are placed so that whenever an exposed interface can be bound by more than one tile or face, the alternatives are related by the half-turn symmetry. In the $N=4$ example, for instance, an exposed $a$ can bind either $T_{2}$ or $T_{3}$, and growth then proceeds as $T_{1}T_{2}T_{3}T_{1}$ or as $T_{1}T_{3}T_{2}T_{1}$, which is the same rod viewed from the other end. Every growth pathway therefore produces the target rod, from any seed.

Why no further reuse is possible. Consider what happens if we try to reuse more in the $N=4$ example (Fig. S18(c)). Setting $c=a$ (and hence $d=b$) turns $T_{2}$ into $(b,a)$, which binds to copies of itself, so the rod can grow forever. Setting $c=b$ (and hence $d=a$) turns $T_{2}$ into $(b,b)$, so the trimer $T_{1}T_{2}T_{1}$ can also form and stop.

These two failures illustrate the general problem. If the same complementary pair is used at two bonds that are not related by the symmetry, tiles intended for one bond can also attach at the other, creating alternative growth pathways. One of three things can then happen:

• the part of the rod between the two bonds is skipped, giving a rod that is too short (as for $c = b$);

• that part is repeated, giving a rod that is too long or grows without bound (as for $c = a$);

• if the two bonds carry the pair in opposite order, growth can turn around and run back along the rod, giving a rod whose length is correct only if the two bonds are symmetry related (Figure S18).

The seed may lie to the left of both bonds, between them, or to the right of both, and in each case at least one of these alternatives is available. The argument therefore holds for any seed position, and so applies to both seeded and unseeded assembly. Hence a complementary pair can be shared only by symmetry-related bonds, and the $N - 1$ bonds need at least $\lfloor N/2 \rfloor$ distinct pairs. With the neutral type, $\widetilde{K}_i \ge2\lfloor N/2 \rfloor+ 1$.

The same reasoning applies to tiles: two interior position can bear the same species only if they are symmetry related, and, as noted above, an end tile can never share a species with an interior tile. For even $N$ there is one further restriction. The two middle tiles are symmetry related, but they cannot be the same species. If they were, the central bond would join two copies of the same face, and hence two faces of the same interface type, which cannot bind. Hence $\widetilde{K}_s \ge\lfloor N/2 \rfloor+ 1$, and the constructions above are minimal.

This central constraint is the origin of the parity dependence in Eq. (8): extending an odd rod by one tile requires an additional species and complementary pair, whereas extending an even rod by one tile requires neither.

L-shaped assembly

Comparison. The rod and the L-shape differ only in the position of a single cube, yet their minimal complexities differ substantially. For large $N$, the rod needs about $N/2$ species and $N$ interface types, whereas the L-shape needs about $N$ species and $2N$ interface types. The difference comes from symmetry: the half-turn symmetry of the rod lets tiles and interfaces be reused in symmetry-related pairs, whereas the L-shape has no such symmetry and lies close to the fully addressable limit of $N$ species and $2N-1$ interface types.

Rule-set diversity of the square and rod tetrominoes

Although each structure is assigned a single assembly complexity, based on the simplest BD rule set that generates it, many BD rule sets of different complexities may produce the same structure. Structures that admit compact assembly descriptions are expected to be compatible with a larger and more diverse set of rule sets than structures that require more specific assembly instructions, and therefore to have higher assembly rule probability. To examine this, Figure S19 shows the distributions of complexities of all individual BD rule sets generating the two-dimensional square and rod of size 4, sampled from $S^{2D}_{4s,17c}$ under seeded assembly. The square is generated by substantially more rule sets, spanning a wider range of complexities, than the rod, consistent with its much higher assembly rule probability. The minimal complexities found by sampling, three species and five interface types for the rod, agree with Proposition 1.

Robustness of our results to simulation and polycube model parameters

Assembly size and spatial cutoffs

To identify unbounded growth, we impose a maximum linear extent of $D_{\max}=16$ tiles and a maximum assembly size of $N_{\max}=100$ tiles. This spatial cutoff was also used in previous polyomino assembly studies5. We expect the linear cutoff to be sufficient for identifying most unbounded assemblies. We found empirically that, for the $10^9$ samples we use, it is extremely rare to find structures close to these cutoffs in size. We therefore expect the effect of this cutoff to be small for the properties we are studying here. However, for larger numbers of allowed species and interfaces, or for longer runs, these cutoffs may need to be adjusted.

Sensitivity to determinism testing

The determinism test uses up to $k=25$ stochastic assembly trajectories per rule set. This value was chosen as a compromise between sensitivity to non-deterministic growth and computational cost, given that up to $10^9$ rule sets are sampled per rule space. A comparable number was used by Johnston et al.5, who tested determinism using 20 repeated assembly trajectories.

To assess the sensitivity of our classification to the number of trajectories, we re-tested all rule sets initially classified as BD in a sample of $10^6$ rule sets from $S^{3D}_{5s,31c}$, using an independent set of up to $k=100$ trajectories. Only 0.25% of these rule sets were reclassified as non-deterministic. These reclassifications occurred predominantly among rule sets generating rare structures, while the assembly rule probabilities of frequently observed structures remained essentially unchanged. Thus, increasing $k$ has little effect on the results, although non-determinism along unsampled trajectories cannot be ruled out.

Dependence on rule-space parameters

The main results are obtained from the rule space $S^{3D}_{5s,31c}$, which is a compromise between expressivity (how many kinds of polycubes can be made) and computational efficiency; larger rule-spaces are more expensive to sample. To test the effect of varying the number of tile species and interface types, we also sampled rule sets from $S^{3D}_{4s,17c}$. As shown in Figure S20, the relationships between assembly rule probability and both assembly complexity and symmetry are quite similar in the two rule spaces. The qualitative trends should not be that sensitive to our choice of $s=5$ and $c=31$. It would be interesting to do these analyses with larger numbers of species and interface types, but that is also computationally more expensive.

Dependence on sampling depth

The full morphospace of polycube structures is incredibly large. Sampling inevitably only captures a small fraction of the full morphospace of possibilities. To assess the effect of sampling depth, we compare results from $10^6$ to $10^9$ randomly sampled rule sets in $S^{3D}_{5s,31c}$ under seeded assembly. Increasing the number of samples recovers additional distinct polycubes, particularly larger and less frequently assembled structures, while the overall assembly rule probability distribution across sizes changes little (Figure S1). Similarly, the assembly rule probability distributions across symmetry orders at fixed size become comparatively stable at the largest sampling depths (Figure S15).

Minimal rule-set calculations

We estimate assembly complexity by taking the minimum complexity among sampled BD rule sets that produce a particular structure. This approach implies that longer sampling can improve estimates of the minimum assembly complexity, as additional rule sets may provide simpler descriptions of structures already observed (see, e.g., Figure S9). Empirically, we find that increasing the sampling depth has little effect on the estimated minimum complexities of the most frequently assembled structures, partly because simpler rule sets tend to be sampled more frequently than more complex ones. We expect sampling errors to be more significant for less frequently assembled structures. Formal methods, such as the SAT-based inverse-design solver of Ref. 6, could be used to determine the true minimum complexities, but are computationally expensive. Since the most frequently assembled structures account for much of the probability weight in our analysis, we do not expect more exhaustive minimization to substantially affect our main findings.

Dependence on seeded and unseeded assembly

The main results are obtained under seeded assembly, where the first tile species initiates growth. We also perform unseeded assembly, where the initial tile species is chosen randomly from the assembly kit. Unseeded assembly produces a significantly smaller fraction of BD rule sets (main-text Figure 2(c)), but the overall relationships between assembly rule probability, complexity, and symmetry remain qualitatively similar. The probability–complexity distributions retain their approximately exponential upper bounds (Figure S8), while the size distributions show similar patterns under both growth conditions (Figure S2). The corresponding symmetry distributions for size-20 polycubes are compared in Figure S16 and Figure S17.

Effects of cubic lattice geometry

The cubic lattice restricts the possible symmetries of assembled structures to those compatible with its geometry. The relative abundance of different symmetry classes in morphospace and their assembly rule probabilities will depend on the underlying lattice. Investigating the relationships between symmetry, assembly complexity, and assembly rule probability for other lattices or off-lattice systems remains an open direction for future work.

Limitations of the assembly dynamics

Our model considers sequential, irreversible attachment of tiles, whereas more general models could allow reversible attachment and detachment, variable bond strengths, and kinetic effects, as in the kinetic Tile Assembly Model (kTAM)7. Similarly, our model allows binding only between distinct complementary interface types, whereas more general interaction rules could permit self-binding interfaces through the interaction matrix $A_{ij}$, as well as cooperative attachment, as in the abstract Tile Assembly Model (aTAM)7. It would be interesting to investigate the relationships between assembly rule probability and complexity under these more general interaction rules and assembly dynamics.