Dhananjay Ipparthi*^{a},
Andrew Winslow†^{b},
Metin Sitti^{cd},
Marco Dorigo^{a} and
Massimo Mastrangeli‡*^{cd}
^{a}IRIDIA, CoDE, Université libre de Bruxelles, Brussels, Belgium. E-mail: dhananjay.ipparthi@ulb.ac.be
^{b}Département d'Informatique, Université libre de Bruxelles, Brussels, Belgium
^{c}Physical Intelligence Department, Max Planck Institute for Intelligent Systems, Stuttgart, Germany
^{d}Max Planck ETH Center for Learning Systems, Stuttgart, Germany
First published on 19th September 2017
We investigate the parallel assembly of two-dimensional, geometrically-closed modular target structures out of homogeneous sets of macroscopic components of varying anisotropy. The yield predicted by a chemical reaction network (CRN)-based model is quantitatively shown to reproduce experimental results over a large set of conditions. Scaling laws for parallel assembling systems are then derived from the model. By extending the validity of the CRN-based modelling, this work prompts analysis and solutions to the incompatible substructure problem.
A set of mobile components capable of bonding with one another will tend to form assemblies in a bounded space. The formation of intercomponent bonds decreases the enthalpy of the system^{6} and the number of accessible component configurations.^{7} Under thermal conditions that make entropic contributions negligible, the assembly of components reduces the free energy of the system, as well as the number of further bonds that can be formed. However, modular assembling systems are characterised by a multitude of intermediate states, associated with local minima in the free energy landscape, besides a few degenerate target states that correspond to global minima.^{8} Assembling systems able to break out of local energy minima, thanks to, e.g., perturbing energy imparted to the system to compete with bond formation, are self-assembling systems.^{9,10} Assembling systems incapable of escaping local energy minima in finite time are aggregating systems.^{11} Irreversible intercomponent bonds arrest aggregating systems typically within one of their local energy minima. Once in a local minimum, aggregating systems are prevented from exploring further possible configurations to reach the target state(s). Self-assembly is often used in literature, albeit incorrectly, to describe aggregating systems as well.^{10} In the following, we use assembly to describe systems that aggregate.
An assembling system is in a depletion trap when no more components are available to advance the assembly, as all components are already employed in existing (sub)structures.^{12} In parallel assembling systems of homogeneous components, depletion traps are absorbing states of the dynamics and local energy minima corresponding to the formation of incompatible substructures, i.e., partial target structures which cannot complete assembly due to steric incompatibility.§ Fig. 1 illustrates an instance of the issue. By preventing the assembly of target structures, the formation of incompatible substructures significantly reduces the assembly yield in parallel homogeneous assembly, and therefore wastes resources.^{16} At the same time, the incompatible substructure problem is both a logical and a physical problem, and is thus amenable to both analytical and experimental study.
In this paper, we present a comprehensive study of model-based prediction of assembly yield for parallel assembling systems. We consider the parallel assembly of two-dimensional (2D), geometrically-closed target structures composed of homogeneous sectors of varying anisotropy. The results of an extensive set of parametric assembly experiments are quantitatively shown to closely match in most cases the assembly yield predictions obtained by a corresponding chemical reaction network (CRN)-based model. Consequently, this work significantly extends the validity of CRN-based modelling of assembling systems, and the analysis of the incompatible substructure problem provides a foundation for the study of dynamical aspects of parallel (self)-assembly.
The paper is organised as follows. Section 2 gives an account of prior art. The physical system and the CRN-based model used in this study are described respectively in Sections 3 and 4. The experiments conducted to compare the physical system and analytical model are presented in Section 5. The results are reported in Section 6. Scaling of system properties is discussed in Section 7. Finally, Section 8 provides conclusions and outlook.
In the context of parallel assembly, Hosokawa et al.^{27} studied the assembly of triangular components into hexagonal target structures. Hosokawa et al. described the negative impact of incompatible substructures on assembly yield in their experiments, and first proposed to study their parallel assembly system using the formalism of chemical reaction networks.
Klavins^{28} built a mechatronic version of the triangular components used by Hosokawa et al.^{27} Software embedded in their “programmable parts” defined rules, based on graph grammars, to guide the assembly of specified target structures. The programmability of the components was used to control the interactions, leading to coordinated behaviour of the system—an approach recently extended by Haghighat and Martinoli.^{29}
Miyashita et al.^{30} studied the influence of reversible reactions on the yield of an assembly system. Their latching self-propelled components floated on the surface of water and had similar shape to those of Hosokawa et al.^{27} Miyashita et al. conducted experiments to compare sequential aggregation, reversible assembly and random aggregation using these components. They used a CRN-based model to quantitatively describe their system.
Mermoud et al.^{31} and Haghighat et al.^{32} developed a magneto-fluidic system of passive, centimetre-sized components that self-assembled on water. The system was supported by a general, CRN-based, multi-level modelling framework for stochastic distributed systems of reactive agents. Using such framework, the authors were able to control the self-assembly of the water-floating components in real-time. Mastrangeli et al.^{33} developed a downscaled version of the system, tailored for the automated control of the acousto-fluidic self-assembly of microparticles.
Hacohen et al.^{34} proposed an algorithm to program the mechanical self-assembly of 3D macroscale objects. Taking inspiration from the fully specified assembly of DNA-based structures, Hacohen et al. showed that an arbitrary 3D object, properly dissected into components, can be re-assembled by unsupervised mechanical shaking when appropriate rules are uniquely encoded on the faces of the components. In their experiments, Hacohen et al. introduced enough components to assemble 2 target structures of 18 components in parallel, and were able to assemble only 1 target structure with no erroneous bonds.
Murugan et al.^{12} studied the formation of incompatible substructures in chemical reaction-based systems. Using a theoretical model relevant to DNA, proteins and colloids, the authors suggested that the incidence of incompatible substructures can be reduced by appropriately tuning the reaction rates by the stoichiometry of the system.
Borrowing the concepts of intercomponent bond and structure formation from chemical reactions, CRN-based formalisms and approaches^{35} are commonly used to model and simulate aggregating and self-assembling systems. However, experimental validation of CRN-based model predictions of assembly yield is, to date, mostly qualitative. For instance, Hosokawa et al.^{27} conducted a single experiment under 2 conditions (systems of 20 and 100 components) and 100 trials per condition, while Miyashita et al.^{30} conducted a single experiment under 2 conditions (systems of 6 and 7 components) with 10 trials per condition. While noteworthy demonstrations of the feasibility of CRN-based modeling, these contributions fall short of showing the details of the relationship between the CRN-based model and physical experiments. Such a lack of comprehensive and quantitative comparison motivates the work presented here.
Fig. 2 The 3D printed components with embedded magnets used in the experiments (dimensions overlaid). Positions and orientations of the magnets are marked. |
We designed target structures of equal area and three different shapes: circle, square and triangle. All component sets are homogeneous. Sectors of a circle with radius 25 mm compose the circular target structures. We use sets of sector-shaped components with 4 different sector angles: 45°, 60°, 72° and 90°. Magnets are embedded in the middle of the two straight faces. For square target structures, the components are 4 equal squares with side length of 22.2 mm, with magnets embedded in two adjacent lateral faces. For triangular target structures, the components are 3 identical isosceles triangles, the two equal sides and the base side measuring 38 mm and 67 mm, respectively. Magnets are embedded in the middle of the two equal sides.
(1) |
Given a large enough population, we study the probabilistic evolution of the system using a system of difference equations:
x(t + 1) = x(t) + F(x(t)) |
The geometry-specific values of P^{b} are calculated according to Hosokawa et al.^{27} We implement the calculation by placing the first substructure at the origin and the second substructure at infinity. The faces of the components form an angle of visibility. The number P^{b} is the probability that the angle of visibility of one substructure's face that includes the South polarity magnet (refer to Fig. 2) is within the angle of visibility of another substructure's face that includes the North polarity magnet. A more detailed explanation of the computation is available in Appendix A of Hosokawa et al.^{27} The values of P^{b} used in this work are presented in Table 1.
Chemical Reaction | X_{3} target | X_{4} target | X_{5} target | X_{6} target | X_{8} target | |
---|---|---|---|---|---|---|
1 | 2X_{1} → X_{2} | 0.389 | 0.438 | 0.460 | 0.472 | 0.484 |
2 | X_{1} + X_{2} → X_{3} | 0.222 | 0.375 | 0.420 | 0.444 | 0.469 |
3 | X_{1} + X_{3} → X_{4} | 0 | 0.188 | 0.320 | 0.417 | 0.453 |
4 | 2X_{2} → X_{4} | 0 | 0.250 | 0.340 | 0.389 | 0.438 |
5 | X_{1} + X_{4} → X_{5} | 0 | 0 | 0.160 | 0.278 | 0.438 |
6 | X_{2} + X_{3} → X_{5} | 0 | 0 | 0.240 | 0.334 | 0.406 |
7 | X_{1} + X_{5} → X_{6} | 0 | 0 | 0 | 0.139 | 0.328 |
8 | X_{2} + X_{4} → X_{6} | 0 | 0 | 0 | 0.222 | 0.375 |
9 | 2X_{3} → X_{6} | 0 | 0 | 0 | 0.250 | 0.359 |
10 | X_{1} + X_{6} → X_{7} | 0 | 0 | 0 | 0 | 0.219 |
11 | X_{2} + X_{5} → X_{7} | 0 | 0 | 0 | 0 | 0.281 |
12 | X_{3} + X_{4} → X_{7} | 0 | 0 | 0 | 0 | 0.313 |
13 | X_{1} + X_{7} → X_{8} | 0 | 0 | 0 | 0 | 0.109 |
14 | X_{2} + X_{6} → X_{8} | 0 | 0 | 0 | 0 | 0.188 |
15 | X_{3} + X_{5} → X_{8} | 0 | 0 | 0 | 0 | 0.234 |
16 | 2X_{4} → X_{8} | 0 | 0 | 0 | 0 | 0.250 |
The F_{i} are generally expressed as:
(2) |
Using only the mean values of x_{i} in our analysis would prevent us from accurately predicting the yield of a system. This is especially true for smaller values of x_{i}. Therefore, we use the master equation, which uses the probability distributions of each x_{i} instead of their mean values. The master equation is:
The master equation was solved numerically. The equation was iterated over interaction steps until only absorbing states and their corresponding probabilities were left. This inherently allowed listing the number of possible absorbing states.
We conducted 6 experiments. The experiment with 45°-sector shape components was conducted under 33 conditions parameterised by the number n of components, which varied from 8 (sufficient for 1 target structure) to 40 (sufficient for 5 target structures). 20 trials were conducted for each of the 33 conditions, and an additional 21 trials were conducted for 12 of the 33 conditions. Experiments for 60°,72° and 90° circular sectors and for square and triangular components were conducted under 13 conditions each. The number of components used were: N, N + 1, 2N − 1, 2N, 2N + 1, 3N − 1, 3N, 3N + 1, 4N − 1, 4N, 4N + 1, 5N − 1 and 5N, respectively. 11 trials were conducted per condition.
We conducted each trial as follows. We manually placed the individually separated components at various initial positions and orientations in the reactor, ensuring that the polarities on the right edge of every component were the same and that intercomponent bonds were not to instantaneously form upon inception of shaking. The shaking frequency was then set to 350 rpm for all circular sector-shaped and square components. Using the same shaking frequency and type of magnets, we observed lever-caused breakage in substructures of triangular components; thus in this case we ran the experiments at 300 rpm to eliminate the problem.
After inception of shaking, each experiment was run until the system reached an absorbing state, composed only of target structures and incompatible substructures (Fig. 3). The final population of structures was then manually recorded.||
Fig. 5 Number of possible and of observed absorbing states as a function of the number of components for the system of Fig. 4. |
Fig. 6 Waste, i.e., relative fraction of components not assembled into target structures, for the system of Fig. 4. |
The relation between absolute assembly yield and possible number of absorbing states is explicited in Fig. 16 of Appendix C. The curve is multi-valued since, as shown by Fig. 5, different values of n give rise to the same number of possible absorbing states. Nevertheless, a non-monotonic increase of Y with the number of possible states can generically be observed.
The absolute assembly yield data and error bounds for the 60°, 72°, and 90° sector conditions (Fig. 7–9 (left)) show non-monotonically increasing trends similar to the 45° sector case. The corresponding comparison between possible and observed absorbing states are shown in Fig. 7–9 (right). The relative yield data, presented in Fig. 17–20 of Appendix D, show the convergence of T to an N-dependent asymptotic value already observed for N = 8.
Yield data and number of absorbing states for systems of triangular components are shown in Fig. 11 (left), Fig. 22 and Fig. 11 (right), respectively. Also in this case, the absolute yield curves for experimental and model-based results have closely matching, non-monotonically increasing trends; not all possible absorbing states were observed; and T appears to saturate for large n.
Though time is not a parameter of our experiments, it is worth noting that the trials with the triangular components took 15 minutes on average to reach the absorbing state, compared to an average of 5 minutes for square- and sector-shaped components. This can be only partially justified by the use of a lower shaking frequency in the former case (see Section 5). Visual observations suggested that the triangular components tended to mutually align and stack themselves against the edge of the container. This ordering and packing effect induced by local component density^{36} sterically hindered the roto-translational motion of the components and significantly reduced the rate of interbond formation.
Beyond the third target structure, the formation of the (z − 1)th target structure does not seem to cause a significant drop in yield.
To improve the reliability of the experimental data, the absorbing state space of each assembly system should be completely explored, and a sufficient number of trials conducted to build a reliable empirical estimate of the distribution over the absorbing states.^{37} Even with the non-negligible number of experimental trials we conducted (see Section 5), ultimately limited by time constraints, we were able to explore only a subspace of the possible absorbing states for all experimental conditions. In the experiment with N = 8 we did increase the number of trials per condition from 20 to 41 for selected values of n, as evidenced in Fig. 15. The additional trials did marginally help explore more of the absorbing state space (Fig. 5).
The parallel and combinatorial nature of the assembly process and the existence of incompatible substructures let expect the number of absorbing states to be a power law function of the number of components. Fig. 12 plots the growth of the number of possible absorbing states with the number of components. It can be observed that different sets of components are associated with different scaling exponents; and that, in spite of the different geometry, pairs of sector-shaped components give rise to the same scaling exponent, respectively 1.00 in the case of 120° and 90° sector components and 1.92 for 72° and 60° ones. The origin of such subdivision remains to be investigated. A power law scaling of the absolute assembly yield with the number of components is derived from Fig. 13. In this case, all component types share a very similar scaling exponent (∼1.00). Scaling of the relative assembly yield with the number of components is shown in Fig. 14, which further evidences the convergence of T toward a fixed, N-dependent values for large values of n.
Fig. 12 Scaling of the number of possible absorbing states for large populations of homogeneous sector-shaped components of varying anisotropy. |
Fig. 13 Scaling of absolute assembly yield for large populations of homogeneous sector-shaped components of varying anisotropy. |
Fig. 14 Scaling of relative assembly yield for large populations of homogeneous sector-shaped components of varying anisotropy. |
This work preludes to future analyses of ways to avoid incompatible substructures. A trivial solution to the problem is the compartmentalisation of the assembly space in subregions, each bounding the exact number of components making up a target structure. An obvious approach is also to convert the assembling system into a self-assembling system, whereby intercomponent bonds are reversible and can be broken under specific conditions.^{24} In a system designed so that only the target structures do not break up upon collisions, the components population would be able to escape local energy minima. Another approach consists in templating the assembly of the target structures, for instance by using (non-homogeneous) sets of geometrically complementary components^{16} or addressable, uniquely-encoded intercomponent bonds.^{34,40,41} In yet another approach, originally proposed but never realised by Hosokawa et al.,^{27} mechanical conformational switching^{42} would direct the assembly kinetics of passive components with binary internal states and limit the number of possible coexisting nucleation sites.^{20}
Finally, the abstract model we used takes into account the combinatorial entropy of the system and a limited information on the geometry of the components. Additional entropic considerations might improve the model by integrating complementary spatial information pertaining to, e.g., component density in the assembly space.^{7,36} Developing a model that takes into account spatial entropy will be part of our forthcoming work.
In our study, we bootstrap the physical experiment outcomes in order to obtain a confidence interval for the mean yield (μ_{Y}) of the physical experiments. Using the obtained confidence interval, the error bound, i.e., maximum difference between the model prediction and experimental mean, is calculated according to Algorithm 2. The experimental results (experimentalResults) are in the form of a matrix with numOfConditions rows and numOfTrials columns. We iterate through all the experiment conditions using the index i. A vector (μ_{Y}) of null values is initialised; this vector stores the yield values of the resampled data. For each condition, represented as a row in experimentalResults, we randomly choose (with replacement, i.e., resampled) values numOfTrials times. The yield of the resampled data is computed and stored in μ_{Y}. The process of resampling is repeated a given number of times (numOfBootstrapSamples).
The vector μ_{Y} is sorted in ascending order. As numOfBootstrapSamples = 200 in our work, the lower bound (lowerBound) is the 6th element, and the upper bound (upperBound) is the 195th element respectively. The error bound is the maximum difference between the model yield and one of the bounds.
The code and data employed to estimate the confidence intervals are available to download.††
Fig. 15 Results for the systems with N = 8 (45°) sectors per circular target structure: comparison of model-predicted and experimental absolute assembly yield, including error bounds for 20 and 41 trials. Vertical dotted lines extend the markers on the x-axis for the fully formed target structures. See also Fig. 4. |
Fig. 16 Scatter plot of absolute assembly yield as a function of possible absorbing states for the system with N = 8 (45°) sectors per circular target structure. Experimental yield values are from 20 trials (refer to Fig. 4 and 5). |
Fig. 17 Relative assembly yield for the system of Fig. 4 (N = 8). |
Fig. 18 Relative assembly yield for the system of Fig. 7 (N = 6). |
Fig. 19 Relative assembly yield for the system of Fig. 8 (N = 5). |
Fig. 20 Relative assembly yield for the system of Fig. 9 (N = 4, sector-shaped components). |
Fig. 21 Relative assembly yield for the system of Fig. 10 (N = 4, square-shaped components). |
Fig. 22 Relative assembly yield for the system of Fig. 11 (N = 3, triangle-shaped components). |
Footnotes |
† Current address: College of Engineering and Computer Science, The University of Texas Rio Grande Valley, Edinburg, USA. |
‡ Current address: Department of Microelectronics, Delft University of Technology, Delft, The Netherlands. E-mail: m.mastrangeli@tudelft.nl |
§ The parallel assembly of target structures composed of pairs of components is not affected by incompatible substructures, since in this case the only substructures are single components, and there is a single absorbing state associated to each number of components. This limit case has important technological applications, such as, e.g., the sequential three-dimensional self-assembly and packaging of light emitting diodes^{13} and the self-assembly of polymeric microcapsules.^{14,15} |
¶ N = 8 represents a more complex case than the N = 6 case presented by Hosokawa et al.,^{27} Miyashita et al.^{30} and Klavins.^{28} In our experiments, N was limited to 8 by fabrication constraints. |
|| A sample video of an assembly trial is available at http://iridia.ulb.ac.be/~dipparth/site/video/. |
** When computing the yield considering also the additional 21 trials (i.e., on a total of 41 trials), the difference between the model prediction and the observed value becomes smaller; similarly, the corresponding bounds are reduced (see Fig. 15 in Appendix C). |
†† http://iridia.ulb.ac.be/~dipparth/site/. |
This journal is © The Royal Society of Chemistry 2017 |