Surface-triggered cascade reactions between DNA linkers direct self-assembly of colloidal crystals of controllable thickness

TL;DR

This paper introduces a surface-triggered DNA cascade reaction system that enables the controlled self-assembly of colloidal crystals with tunable thickness, combining experiments, simulations, and theory to elucidate the underlying thermodynamics.

Contribution

It presents a novel mechanism where surfaces initiate DNA-mediated cascade reactions, allowing for the fabrication of finite, controllable-thickness colloidal structures.

Findings

Surface interactions activate colloids for assembly.

Thermodynamics of activation and binding are clarified.

Finite aggregates with controllable layers can be formed.

Abstract

Functionalizing colloids with reactive DNA linkers is a versatile way of programming self-assembly. DNA selectivity provides direct control over colloid-colloid interactions allowing the engineering of structures such as complex crystals or gels. However, self-assembly of localized and finite structures remains an open problem with many potential applications. In this work, we present a system in which functionalized surfaces initiate a cascade reaction between linkers leading to self-assembly of crystals with a controllable number of layers. Specifically, we consider colloidal particles functionalized by two families of complementary DNA linkers with mobile anchoring points, as found in experiments using emulsions or lipid bilayers. In bulk, intra-particle linkages formed by pairs of complementary linkers prevent the formation of inter-particle bridges and therefore colloid-colloid aggregation. However, colloids interact strongly with the surface given that the latter can destabilize intra-particle linkages. When in direct contact with the surface, colloids are activated, meaning that they feature more unpaired DNA linkers ready to react. Activated colloids can then capture and activate other colloids from the bulk through the formation of inter-particle linkages. Using simulations and theory, validated by existing experiments, we clarify the thermodynamics of the activation and binding process and explain how particle-particle interactions, within the adsorbed phase, weaken as a function of the distance from the surface. The latter observation underlies the possibility of self-assembling finite aggregates with controllable thickness and flat solid-gas interfaces. Our design suggests a new avenue to fabricate heterogeneous and finite structures.

Figures

a-c) are given by

a-c) are given by

Equations

$$Z$$

$$\frac{\partial }{\partial \{n\}}{\mathcal{F}}_{\mathrm{multi}}(\{n\}){|}_{\{n\}=\{\overline{n}\}}=0.$$

$${\overline{n}}_{ii}^{\mathrm{AB}}$$

$${\overline{n}}_{ij}^{\mathrm{AB}}$$

$${\overline{n}}_{ij}^{\mathrm{BA}}$$

$${\overline{n}}_i^{\mathrm{AC}}$$

$$\beta \Delta G_{ii}(\{\mathbf{r}\})$$

$$\beta \Delta G_{ij}(\{\mathbf{r}\})$$

$$\beta \Delta G_i^{\mathrm{s}}(\{\mathbf{r}\})$$

$$\beta F(\{\mathbf{r}\})$$

$$\beta {\mathbf{f}}_i$$

$$\Delta {\mathbf{r}}_i$$

$$k_{\mathrm{on}}^{ij}$$

$$k_{\mathrm{on}}^{ii}$$

$$k_{\mathrm{on}}^{i\mathrm{s}}$$

$$f_{\mathrm{multi}}(\Gamma )$$

$$P(\Gamma )$$

$${\mathcal{N}}_{\mathrm{select},\mathrm{A}}^i$$

$$n_i^{\mathrm{A}}$$

$${\mathcal{N}}_{\mathrm{select},\mathrm{B}}^i=\frac{N_{\mathrm{L}}!}{n_i^{\mathrm{B}}!n_{ii}^{\mathrm{AB}}!{\prod }_{j\ne i}{n}_{ij}^{\mathrm{BA}}!}$$

$${\mathcal{N}}_{\mathrm{select},\mathrm{C}}=\frac{N_{\mathrm{C}}!}{n_{\mathrm{s}}^{\mathrm{C}}!{\prod }_{i=1}^{N_{\mathrm{p}}}{n}_i^{\mathrm{AC}}!}$$

$${\mathcal{N}}_{\mathrm{select}}={\mathcal{N}}_{\mathrm{select},\mathrm{C}}\prod _{i=1}^{N_{\mathrm{p}}}\left[{\mathcal{N}}_{\mathrm{select},\mathrm{A}}^i{\mathcal{N}}_{\mathrm{select},\mathrm{B}}^i\right].$$

$${\mathcal{N}}_{\mathrm{pair}}=\prod _i\left[n_i^{\mathrm{AC}}!n_{ii}^{\mathrm{AB}}!\right]\prod _{i<j}\left[n_{ij}^{\mathrm{AB}}!n_{ij}^{\mathrm{BA}}!\right].$$

$${\mathcal{N}}_{\mathrm{comb}}$$

$$\mathcal{Z}$$

$$Z_{T=\infty }(\{\mathbf{r}\})$$

$${\Omega }_i(\{\mathbf{r}\})={\Omega }_0-e_i^{\mathrm{s}}(r_{i,z})-\sum _{j\in v(i)}{e}_{ij}(|{\mathbf{r}}_i-{\mathbf{r}}_j|)$$

$$\frac{\partial }{\partial \{n\}}{\mathcal{F}}_{\mathrm{multi}}(\{n\}){|}_{\{n\}=\{\overline{n}\}}=0.$$

$${\overline{n}}_i^{\mathrm{A}}$$

$${\overline{n}}_i^{\mathrm{B}}$$

$${\overline{n}}_{\mathrm{s}}^{\mathrm{C}}$$

$$V_{\mathrm{pp}}(r_{ij})$$

$$V_{\mathrm{ps}}(r_{i,\mathrm{z}})$$

$$V(\{\mathbf{r}\})=\sum _{i<j}{V}_{\mathrm{pp}}(r_{ij})+\sum _i{V}_{\mathrm{ps}}(r_{i,\mathrm{z}}).$$

$${\Omega }_{ij}(r_{ij},L)$$

Full text

I Introduction

Many recent contributions have unveiled the advantages of using complementary single-stranded (ss) DNA oligomers tethered to colloidal particles to program self-assembly Mirkin et al. (1996); Alivisatos et al. (1996); Jones et al. (2015). The selectivity of Watson-Crick base pairing underlies most of the functionalities and responsive behaviors achieved using DNA Jones et al. (2015). For instance, DNA has been used to self-assemble colloidal crystals lacking molecular analog J. et al. or featuring optical bandgaps Ducrot et al. (2017); Wang et al. (2017); Liu et al. (2016), or engineer bigels Var (2012) and re-entrant phase behaviors Rogers and Manoharan (2015); Angioletti-Uberti et al. (2012).

Recently, systems of particles functionalized by mobile linkers tipped by reactive sites received a lot of attention Beales et al. (2011); Pontani et al. (2012); van der Meulen and Leunissen (2013); Chakraborty et al. (2017a); Hadorn et al. (2012); Parolini et al. (2015); Hu et al. (2018). In these systems, DNA oligomers conjugated to hydrophobic tags are tethered, for instance, to lipid bilayers van der Meulen and Leunissen (2013); Chakraborty et al. (2017a). Lipid vesicles functionalized by mobile ligands are currently employed in nanotechnological platforms to mimic biological functionalities like cell adhesion and recognition Pontani et al. (2012); Parolini et al. (2015). In self-assembly, mobile ligands facilitate annealing of crystal defects van der Meulen and Leunissen (2013) and could allow remote control over the valency of the aggregates Feng et al. (2013); Angioletti-Uberti et al. (2014). In this respect, colloidal supported bilayers Troutier and Ladavière (2007); Rinaldin et al. (2018) functionalized with DNA linkers Rinaldin et al. (2018) are a new generation of a particularly versatile type of rigid and monodisperse building blocks with enhanced programmability given by the mobility of the binders Angioletti-Uberti et al. (2014).

So far DNA mediated interactions have been used almost exclusively to fabricate extended colloidal structures. However, the ability to constrain self-assembly spatially is central in many nanotechnological applications including encapsulation and development of point-of-care devices Chapman et al. (2015). Localized self-assembly is also important in biology where, for instance, many cell functionalities rely on dynamic compartmentalized environments (e.g. Refs. Conduit et al. (2015); Case and Waterman (2015)). Developing bottom-up methods for controlled surface coating is also a key problem in chemistry Cécile et al. .

In this contribution, we show how to use DNA to yield localized self–assembly of colloidal structures. In particular, we study a biomimetic system comprising colloidal supported bilayers functionalized by two types of complementary mobile linkers interacting with surfaces carrying a single type of DNA linker. In most systems, reactions between oligomers tethered to different particles are, to a good extent, independent events. Instead, in the present work, the surface initiates a cascade reaction between DNA linkers that propagate concomitantly with colloidal aggregation. In bulk, DNA linkers predominantly form intra–particle loops while the probability of forming inter–particle bridges remains negligible resulting in weak particle-particle interactions and stable gas phases, Fig. 1. Instead, the unpaired linkers tethered to the surface can easily stabilize colloidal particles through the formation of particle-surface bridges. Once bound to the surface, particles become activated and display a higher number of free DNA linkers appearing as a side product of the reaction leading to the formation of particle–surface bridges. Importantly, activated colloids can attract and activate others colloids from the bulk. This process triggers a domino effect leading to the self–assembly of colloidal crystals at the surface, Fig. 1. We use theoretical modeling, supported by Brownian dynamics simulations, to explain how enhanced particle–particle interactions at the surface is an entropic effect mainly controlled by the number of linkers per particle and the relative statistical weight of forming inter–particle and intra–particle linkages. Importantly, the domino effect does not propagate indefinitely. Instead, colloid-colloid interactions sharply decrease with the distance between newly activated colloids and the surface, Fig. 1. This observation explains the flatness of the solid-fluid interfaces, a missing result in wetting phenomena where roughness and thickness are correlated quantities Bonn et al. (2009). As compared to existing protocols leading to colloidal layer deposition (e.g. Refs. Trau et al. (1996); Reculusa and Ravaine (2003)), our method provides direct control over the number of deposited layers and does not require external intervention during aggregation. Such property arises from the possibility of controlling particle-particle interactions at different particle-surface distances (Fig. 1) using design parameters.

Our design is robust as proven by state-of-the-art simulations of self-assembly directed by ligand-receptor complexation. We prove the reliability of our methods using recent experiments that investigated the stability of suspensions of colloids featuring competition between bridges and loops in bulk Bachmann et al. (2016). Beyond addressing an important technological problem, this paper suggests new routes leading to the fabrication of localized structures and new responsive behaviors at functionalized surfaces.

II Theoretical and simulation methods

Figs. 2$a$ and 2$b$ present the design that we use to reproduce the peculiar collective behaviors anticipated by Fig. 1. We consider rigid colloidal particles of radius $R$ functionalized by two types of linkers (or ligands), A and B, freely moving on the surfaces of the particles. Each particle carries $N_{\mathrm{L}}$ ligands of each type, while the surface is decorated with $N_{\mathrm{C}}$ mobile linkers (or receptors). In our model, A can bind to B and C, while C and B do not pair. Therefore, the possible linkages featured by the system are inter-particle loops, particle-particle bridges, and particle-surface bridges (see Fig. 2$a$). Van der Meulen and Leunissen first studied silica particles coated with DNA linkers with a hydrophobized head immersed in the bilayers covering the colloids van der Meulen and Leunissen (2013). The group of Kraft is currently studying the potentialities of colloidal supported lipid bilayers functionalized by DNA for self-assembly Rinaldin et al. (2018); Chakraborty et al. (2017b). In our design, both ligands and receptors comprise short rods of double-stranded (ds) DNA of length $L$ tipped by reactive single-stranded (ss) DNA sequences (see legend of Fig. 2). In this study, we use $R=5\cdot L$ and $L$ is chosen as the unit length.

As it is often the case in soft-matter, interactions between complex building blocks include entropic contributions due to ensemble averages over all possible configurations of the system compatible with a given particle-particle (or particle-surface) distance. For multivalent interactions, entropic terms are due to the configurational constraints that surface anchored binders need to satisfy (coalescence of the endpoints and, for mobile constructs, colocalization of the tethering points), and combinatorial factors counting the possible ways of making a certain number of linkages Dreyfus et al. (2009); Martinez-Veracoechea and Frenkel (2011); Angioletti-Uberti et al. (2016). In Sec. II.1 we derive the multivalent free energy $F$ controlling particle–particle and particle–surface interactions Angioletti-Uberti et al. (2014); Jan Bachmann et al. (2016), while in Sec. S3: Simulation details we detail the simulation methods.

II.1 Derivation of the multivalent free energy

The partition function of $N_{\mathrm{p}}$ colloids carrying $N_{\mathrm{L}}$ A and B ligands facing a surface decorated with $N_{\mathrm{C}}$ receptors of type C (see Fig. 2) reads as

[TABLE]

where $\{\bf{r}\}$ denotes the cartesian coordinates of the particles ($\{\bf{r}\}=\{{\bf r}_{1}\cdots{\bf r}_{N_{p}}\}$) and $\{n\}$ ($\{n\}=\{n^{\mathrm{AB}}_{ii},\,n^{\mathrm{AB}}_{ij},\,n^{\mathrm{BA}}_{ij},\,n^{\mathrm{AC}}_{i}\}$ with $i,j=1,\cdots N_{\mathrm{p}}$) the ensemble of intra-particle ($n^{\mathrm{AB}}_{ii}$), inter-particle ($n^{\mathrm{AB}}_{ij}=n^{\mathrm{BA}}_{ji}$ and $i\neq j$), and particle-surface ($n^{\mathrm{AC}}_{i}$) linkages. $n^{\mathrm{AB}}_{ij}$/$n^{\mathrm{BA}}_{ij}$ is the number of bridges between particle $i$ and $j$ resulting from binding linkers of type $A$/$B$ tethered to particle $i$ with linkers of type $B$/$A$ tethered to $j$. ${\cal Z}$ and ${\cal F}_{\mathrm{multi}}$ are, respectively, the partition function of the system and the multivalent free energy at given $\{{\bf r}\}$ and $\{n\}$. ${\cal Z}$ comprises combinatorial terms, counting the number of ways of making $\{n\}$ linkages, and terms linked to the hybridization free energy of DNA reactive sequences SantaLucia (1998). In Sec. S1 of the SI we adapt the calculations of Refs. Angioletti-Uberti et al. (2014); Jan Bachmann et al. (2016); Martinez-Veracoechea and Frenkel (2011) to the system of Fig. 2 and report the explicit expressions of ${\cal Z}$ and ${\cal F}_{\mathrm{multi}}$.

At given colloid positions, $\{{\bf r}\}$, the most likely number of linkages featured by the system, $\{\overline{n}\}$, are calculated by maximizing the multivalent free energy

[TABLE]

Eq. 2, along with the expression of ${\cal F}_{\mathrm{multi}}$ (see SI Sec. 1), lead to the chemical equilibrium equations for the most likely number of linkages at given $\{{\bf r}\}$

[TABLE]

where $\overline{n}_{i}^{\mathrm{A}}$, $\overline{n}_{i}^{\mathrm{B}}$, and $\overline{n}^{\mathrm{C}}_{\mathrm{s}}$ denote the number of free (unbound) linkers. In Eqs. 3 we defined the hybridization free energies of making inter–particle, intra–particle, and particle–surface linkages with $\Delta G_{ii}$, $\Delta G_{ij}$, and $\Delta G^{\mathrm{s}}_{i}$, respectively. The hybridization free energies comprise the binding free energy of the reactive oligomers free in solutions SantaLucia (1998) ($\Delta G_{0}$ and $\Delta G_{0}^{\mathrm{s}}$ for the dimerization of $A$ with $B$ and $A$ with $C$, respectively) augmented by configurational contributions as follows Jan Bachmann et al. (2016); Angioletti-Uberti et al. (2014)

[TABLE]

where $\rho_{0}$ is the standard concentration ($\rho_{0}=1\,$M$\cdot$liter*-1*) while $\Omega_{ij}(\{{\bf r}\})$ and $\Omega_{i}^{\mathrm{s}}(\{{\bf r}\})$ are the volume of the configurational space available to each linkage cross-linking particle $i$ with $j$, and particle $i$ with the surface, respectively. Note how configurational contributions depend on the position of the particles $\{{\bf r}\}$. In this work we consider reactive sequences tethered to particles’ surfaces through short, thin rods of double-stranded DNA of length $L$ (see Fig. 2). When $L$ is much smaller than the radius of the particles, the reactive sequences of unbound linkers are uniformly distributed within the layer of thickness $L$ surrounding the tethering surfaces Angioletti-Uberti et al. (2014). In this limit, the configurational space of bound sequences ($\Omega_{ij}$ and $\Omega^{\mathrm{s}}_{i}$) is the volume of the overlapping regions spanned by the reacting sequences before binding (see Fig. 2$b$). Similarly, the volumes available to unbound reactive sequences ($\Omega^{\mathrm{s}}$ and $\Omega_{i}$) are depleted by the volume excluded by the hard-core of the neighboring particles (respectively, $k^{\mathrm{s}}_{i}$ and $e_{ij}$ in Fig. 2$b$) and of the surface ($e^{\mathrm{s}}_{i}$ in Fig. 2$b$). We report the explicit expression of the terms appearing in Eqs. 4 in SI Sec. 1. When written in term of the stationary number of linkages, the multivalent free energy simplifies into a portable expression $F$ ($F={\cal F}_{\mathrm{multi}}(\{\overline{n}\})$) given by (see SI Sec. 1) Angioletti-Uberti et al. (2016, 2014); Jan Bachmann et al. (2016)

[TABLE]

where $F_{T=\infty}$ is the free energy without any linkage (as found at infinite temperature, $T$) and accounts for non-selective interactions and repulsive osmotic terms due to compression of the linkers in the contact region Dreyfus et al. (2009); Angioletti-Uberti et al. (2014, 2016). We use Eq. II.1 to sample colloidal configurations $\{{\bf r}\}$ in the Mean Field calculations (Sec. III.1) and in simulations (Sec. S3: Simulation details).

II.2 Simulation methods

Modelling self-assembly dynamics of particles forming reversible linkages requires an algorithm capable of evolving the number of linkages $\{n\}$ and colloids’ positions $\{{\bf r}\}$ in a concerted way Angioletti-Uberti et al. (2014); Jan Bachmann et al. (2016). At each step of our simulation scheme, we first upgrade the number of linkages between particles, $\{n\}\to\{n^{\prime}\}$, while keeping $\{\bf r\}$ fixed. We then calculate forces acting on each particle ${\bf f}_{i}$ using ${\cal F}_{\mathrm{multi}}$ (Eq. 1)

[TABLE]

where $\{\Delta G\}$ is the ensemble of possible hybridization free energies (see Eqs. 4). We provide the explicit expression of ${\bf f}_{i}$ and further details on the simulation method in SI Sec. 3. Using ${\bf f}_{i}$ we evolve particles’ position using a Brownian dynamics scheme, $\{{\bf r}^{\prime}\}=\{{\bf r}\}+\{\Delta{\bf r}\}$, with

[TABLE]

where $D$ is the particle diffusion constant, $\Delta t$ the integration step, and ${\bf{\cal N}}$ a normal distributed vector with covariance matrix equal to the unitary matrix.

We employ two schemes upgrading $\{n\}$ in different ways. In the implicit scheme (IMP), $\{n^{\prime}\}$ are taken equal to the solutions of Eqs. 2 (see SI Sec. 3 for the iterative procedure used to solve Eqs. 2 Angioletti-Uberti et al. (2014); Jan Bachmann et al. (2016)). Such scheme minimizes the multivalent free energy ${\cal F}_{\mathrm{multi}}$ at each step of the dynamics, implicitly assuming infinite reaction rates between DNA linkers ($k_{\mathrm{on}}$ and $k_{\mathrm{off}}$ in Fig. 2$c$). To probe the effect of finite reaction rates on the dynamics of self-assembly, we also developed a scheme in which we explicitly simulate linkages’ dynamics using the Gillespie algorithm (EXP) Gillespie (1977); Jan Bachmann et al. (2016). At a given $\{\bf{r}\}$, we start by calculating the rates at which inter–particle linkages, loops, and particle-surface bridges form

[TABLE]

where $k_{\mathrm{on}}^{0}$ is the dimerization rate of free oligomers in solutions. The previous equations follow from the definitions of the hybridization free–energies (Eqs. 4) while assuming that the $off$ rates of tethered DNA match the $off$ rates of free oligomers in solution Parolini et al. (2016); Jan Bachmann et al. (2016); Ho et al. (2009): $k_{\mathrm{off}}^{ij}=k_{\mathrm{off}}^{ii}=\rho_{0}\exp(\beta\Delta G_{0})k_{\mathrm{on}}^{0}$ and $k_{\mathrm{off}}^{i\mathrm{s}}=\rho_{0}\exp(\beta\Delta G_{0}^{\mathrm{s}})k_{\mathrm{on}}^{0}$. Once all $on$ and $off$ rates are known, we use the standard procedure of the Gillespie algorithm and sequentially sample one within all possible reactions along with the expected time ($\tau$) for it to happen. We increment a reaction clock ($\tau_{\mathrm{reac}}$) by $\tau$, $\tau_{\mathrm{reac}}=\tau_{\mathrm{reac}}+\tau$ and repeat the procedure until $\tau_{\mathrm{reac}}$ remains smaller than $\Delta t$. At that point, the Gillespie algorithm is arrested and a Brownian Dynamics step is performed (Eq. 7). Notice that in our scheme oligomers of the same type on different particles are treated as different reacting species given the fact that the $on$ rates (Eq. 8) are configuration dependent. We report further details on the Gillespie algorithm in SI Sec. 3.

In addition to Brownian dynamics simulations, we perform grand canonical moves allowing to use small simulation boxes without depleting the gas phase (see SI Sec. 3 for details).

Chosen unit of length and time are $L$ and ${L^{2}\cdot D^{-1}}$, where $D$ is the diffusion constant of diluted colloids. In particular, the $on$ rate of free reactive sequences in solution ($k_{\mathrm{on}}^{0}$) is expressed in unit of $LD$ while the hybridization free energies $\Delta G_{0}$ and $\Delta G^{\mathrm{s}}_{0}$ have been offset by a constant term equal to $k_{B}T\log(\rho_{0}L^{3})$ (see Eqs. 4).

III Results

III.1 Programming multibody interactions using mobile ligands

Fig. 2$a$ shows how the functionalized surface alters particle-particle interactions. We use Eq. II.1 to calculate the effective interaction $f_{12}$ between pairs of particles in bulk at different particle-particle distance $D$ (Case I) and in contact with the surface (Case II) or with a surface bound colloid (Case II). We offset $f_{12}$ by the value of $F$ at large $D$, $f_{12}=F(D)-F(D=\infty)$ (note that for isolated particles $F=F_{\mathrm{multi}}<0$ given that $\overline{n}^{\mathrm{AB}}_{ii}>0$, see Eq. II.1). In bulk, inter–particle loops dominate (Case I in Fig. 2$a$) resulting in weak particle-particle interactions $f_{12}$ (see full line in the graph of Fig. 2$a$) and a stable gas phase Bachmann et al. (2016). Because receptors (C) at the surface are not self–protected, particles will likely form particle–surface bridges at high enough receptor concentration ($\sigma_{\mathrm{R}}$) and strength (i.e., at low $\Delta G^{\mathrm{s}}_{0}$). When bound to the surface, a particle will free one linker of type B for each particle-surface bridge formed. Free B ligands will then stabilize another particle (see Case II in Fig. 2$a$). This second binding, as well as the attachment of the colloid to the surface, is thermodynamically more favorable than pairing two particles in bulk because it requires denaturating a single loop (instead of two) to form two inter–particle bridges. This is proven by the graph in Fig. 2$a$ showing how $f_{12}$ is magnified by a factor of ten when one colloid is in direct contact with the surface (see Case I and Case II). The graph of Fig. 2$a$ also shows how, when colloids are not in direct contact with the surface (Case III), effective interactions are much weaker. On the other hand, effective interactions between two surface–bound particles (not shown in Fig. 2a) are weaker than in bulk given the fact that both particles express free ligands of the same type. The change of the state of the ligands displayed by particles following the encounter with the surface underlies the possibility of controlling layer deposition.

III.2 Self-assembly of crystals with a finite thickness

We consider suspensions of colloids as in Fig. 2$a$ and report the formation of crystals at the surfaces. In particular, despite the weak in–plane interactions between particles, we never observe the formation of non–compact structures like colloidal chains. Here, we clarify the system parameters and the thermodynamic conditions controlling the morphology of the assemblies. To do so, we employ the results of a Mean-Field Theory (MFT) balancing the free energy gains of forming an aggregate due to multivalent interactions (Eq. II.1) with the entropic losses of caging particles into crystalline sites. Using simulations, below we prove that the proposed MFT is quantitative, therefore providing a predictive tool that will be useful to design future experiments. In particular, we used the MFT to fine tune the system parameters of all simulations presented in this work. The SI Sec. 2 reports details of the MFT calculations. The scripts implementing the MFT can be found at MFT under an MIT license.

We consider the thermodynamic equilibrium between particles in bulk at density $\rho_{\mathrm{id}}$ and crystalline structures made of $\Gamma$ layers with $N_{\mathrm{p,ads}}$ particles per layer ($N_{\mathrm{p,ads}}=N_{\mathrm{ads}}/\Gamma$, where $\Gamma$ is the number of layers and $N_{\mathrm{ads}}$ the number of particles in the crystal) for which we calculate the free energy $F(\Gamma)$ using Eq. II.1. In all cases studied using simulations, we report self-assembly of fcc (111) crystals comprised of hexagonal, stacked layers parallel to the surface. In the top panel of Fig. 3$a$, we then calculate the multivalent free energy gain per particle in direct contact with the substrate defined as

[TABLE]

where $F_{1}$ is the free energy of a single particle in bulk. We find that $f_{\mathrm{multi}}(\Gamma)$ is non-linear at small $\Gamma$, corresponding to magnified interactions between colloids closer to the functionalized surface (see the graph of $f_{12}$ in Fig. 2$a$). Such nonlinearity is more prominent at high values of $N_{\mathrm{L}}$. At larger values of $\Gamma$, $f_{\mathrm{multi}}(\Gamma)$ becomes linear and surface effects negligible.

At a given ligand/receptor strength and coating densities, the chemical potential of the colloid, $\mu$, controls the number of layers assembled. In diluted conditions, $\mu$ is proportional to the logarithm of the gas density, $\beta\mu\sim\log\rho_{\mathrm{id}}$. If we assume that the configurational space available to colloids in the crystal phase is $v_{0}$, the probability of forming crystals made of $\Gamma$ layers is then given by (see SI Sec. 2) Sear (1999); Charbonneau and Frenkel (2007)

[TABLE]

where $Z$ is a normalization factor. In the bottom panel of Fig. 3$a$, we report $P(\Gamma)$ for three different values of $N_{\mathrm{L}}$. We choose three different $\rho_{\mathrm{id}}$ resulting in a most likely number of layers equal to three. We verify the predictivity of $P$ using simulations (see below). Note that the definition of $P$ is conditional on having a finite normalization factor $Z$. The onset in parameter space at which $Z$ diverges corresponds to the fluid–solid phase boundary in bulk (see SI Sec. 2) Sear (1999); Charbonneau and Frenkel (2007). In Figs. 3$b$ and $c$ we then report the predicted bulk phase diagram in the $(\Delta G_{0},\,\rho_{\mathrm{id}})$ and $(\Delta G_{0},\,N_{\mathrm{L}})$ planes. The gas phase is stable at low $N_{\mathrm{L}}$ and $\rho_{\mathrm{id}}$ values. Note how at low $\Delta G_{0}$ (equivalently, at low temperature) the phase boundary does not depend on $\Delta G_{0}$. In this limit, the numbers of paired linkers in the gas and crystalline phase are equal, and combinatorial terms fully control the transition Zilman et al. (2003); Bozorgui and Frenkel (2008). Ref. Bachmann et al. (2016) studied experimentally the self-assembly of 200 nm diameter vesicles functionalized by two families of linkers as in Fig. 2. In Fig. 4, we adapt the MFT developed in the present work to the system of Ref. Bachmann et al. (2016) and report a phase diagram similar to the one of Fig. 3$c$. Our MFT allows predicting the entropic transition at low $T$ (low $\Delta G_{0}$), as well the phase boundary at high temperature, without any fitting parameter (see caption of Fig. 4 for details). The discrepancy between the theoretical and the experimental phase boundary at high $N_{L}$ (not relevant to this study) are due to steric interactions between linkers not included in our MFT Bachmann et al. (2016). Overall, Fig. 4 validates our modeling and suggests using (supported) vesicles in future experiments aiming at reproducing controllable colloidal layer deposition.

We now relax the approximations employed by the MFT and develop Brownian dynamics simulations. For three different values of $N_{\mathrm{L}}$, we run simulations at different gas densities $\rho_{\mathrm{id}}$ and verify the possibility of assembling crystals made of $\Gamma=1$, 2, and 3 layers. Fig. 5 reports snapshots obtained in steady conditions. The choice of $\rho_{\mathrm{id}}$ at different $N_{\mathrm{L}}$ has been guided by the MFT providing predictions of the averaged number of layers, $\langle\Gamma\rangle$, through Eq. 10. As anticipated before, most of our simulations reported the formation of fcc (111) crystals. Occasionally we observed fcc (100) crystals. Crystals of this type are arrested states, occasionally appearing when using small simulation boxes and high rates of particle insertion/deletion in the grand-canonical scheme. In particular, at high insertion rates, the second layer starts to form before the first layer could relax into a triangular lattice.

Given that our design is based on equilibrium considerations (see Eq. 10), different initial conditions lead to the same number of layers. Supplementary videos 1, 2, and 3 prove that the configurations reported in Fig. 5 are stationary. The insets of the panels in Fig. 5 report the histograms with the number of particles found at different distances $D_{\mathrm{s}}$ (expressed in units of $L$) from the plane averaged over steady configurations. Such histogram sharply transits from the maximum number of particles that can be fitted into a single plane to zero at the solid-fluid interface. For the system employing the smaller value of $N_{\mathrm{L}}$ the gas phase is denser because lower values of $N_{\mathrm{L}}$ result in higher $f_{\mathrm{multi}}$ (Eq. 9), and higher $\rho_{\mathrm{id}}$ are required to yield a given number of layer $\Gamma$ (see Eq. 10). Overall, Fig. 5 nicely demonstrates how the proposed system can be used to fabricate crystals with a prescribed number of layers without any direct intervention.

III.3 Controlling thickness in colloidal layer deposition

In this section, we study in more detail the thermodynamic properties of the structures reported in Fig. 5. We consider the system with $N_{\mathrm{L}}=40$ at three different chemical potentials forming from one to three layers as obtained in the second column of Fig. 5. In Fig. 6 we use simulations to study the contribution ($f_{\mathrm{layer}}$) of particles found at a given layer (specified by the plane-particle distance $D_{\mathrm{s}}$) to the multivalent free energy $F$. Specifically, $f_{\mathrm{layer}}$ is calculated using the values of $\{\overline{n}\}$ sampled by the simulations along with Eq. II.1 in which the terms involving $\overline{n}^{\mathrm{AB}}_{jq}$ and $\overline{n}^{\mathrm{BA}}_{jq}$ are evenly split between particles $j$ and $q$. For perfect crystals, we have $f_{\mathrm{multi}}(\Gamma)=\sum_{i=1}^{\Gamma}f_{\mathrm{layer}}(i)$, with $f_{\mathrm{multi}}$ given by Eq. 9. Consistently with the behavior of $f_{12}$ in Fig. 2 and $f_{\mathrm{multi}}$ in Fig. 3, we find that $f_{\mathrm{layer}}$ sharply increases with the particle-surface distance, $D_{\mathrm{s}}$. As a result, when considering systems assembling multiple layers, colloids in the bottom layers are irreversibly adsorbed while at higher layers adsorption/desorption energies become thermally accessible. As predicted by Eq. 10, layers with $f_{\mathrm{layer}}$ greater than $k_{B}T\log(\rho_{\mathrm{id}}v_{0})$ are unstable (see dashed lines in Fig. 6$a$-$c$).

Although $f_{\mathrm{layer}}(i)$ is also a function of $\Gamma$ (given that the numbers of linkages $\{\overline{n}\}$ featured by the particles in the crystal are coupled variables), Fig. 6 shows that such dependency is weak. On the one hand, this observation suggests that $f_{\mathrm{layer}}(\Gamma)$ is well approximated by $\Delta f_{\mathrm{multi}}(\Gamma)\equiv f_{\mathrm{multi}}(\Gamma)-f_{\mathrm{multi}}(\Gamma-1)$ as confirmed by Fig. 6$a$-$c$. Moreover, Fig. 6 also suggests that adding an extra layer to the crystal has little effect on the number of linkages already formed. This motivated a simplification of the MFT in which, when calculating $\Delta f_{\mathrm{multi}}(\Gamma+1)$, we treat the number of free receptors featured by particles sitting in the layer $\Gamma$ as a constant parameter. In the low–temperature regime (low $\Delta G_{0}$ in Fig. 3$b$, $c$), this assumption allows validating an analytic expression, $\Delta f_{\mathrm{adi}}$, that matches well with the values of $\Delta f_{\mathrm{multi}}$ reported in Fig. 6. The explicit expression of $\Delta f_{\mathrm{adi}}$ and further checks of its validity are reported in SI Sec. 2. In general, we advise using the MFT to refine experimental designs MFT .

In Fig. 7 we study the number of inter–particle bridges featured by colloids found at different layers ($h=1,\,2$, and 3) of steady configurations of a system assembling a crystal with $N_{\mathrm{L}}=40$ and $<\Gamma>=3$ (corresponding to Fig. 6$c$). We decompose the total number of bridges emanating from particles belonging to a given layer (tagged with $h$ in Fig. 7) into in–plane, $n_{\mathrm{l}}$, and off–plane linkages, $n_{\mathrm{v}}$. In–plane bridges involve the six neighboring particles of a triangular lattice while off–plane bridges the three+three neighboring particles of an fcc (111) crystal belonging to layers $h\pm 1$ (when present) or the surface if $h=1$. We report simulation data points along with their averages and mean-field predictions. Simulations and theoretical predictions are in agreement. Close to the surface, all bridges are off–plane. This finding shows how two surface–bound particles exhibit an excess of free B ligands and few intra–particle or in–plane bridges. The number of off–plane and in–plane linkages converge when increasing $D_{\mathrm{s}}$, as expected in bulk conditions. However, for the system of Fig. 7, the total number of bridges goes to zero at high $D_{\mathrm{s}}$ in favors of loops consistent with the fact that, in bulk, the gas phase is stable, see Fig. 3$b$–$c$. Notice how particles (and occasionally dimers) can temporarily bind colloids in the third layer but rapidly detach as confirmed by the fact that they do not form any lateral bridge (see Supplementary video 3).

In Fig. 8 we summarize the number of self-assembled layers as a function of the chemical potential of the system, Fig. 8$a$, and the hybridization free energies of the reactive sequences, Fig. 8$b$. We compare numerical simulations with theoretical predictions and report a quantitative agreement. Compared to theoretical results, numerical predictions do not feature fractional values of $\Gamma$. In particular, in our simulations, we never observe half-filled top layers. SI Fig. 10 reporting the number of adsorbed particles at different chemical potentials versus time confirms this result. Line tension effects, not considered by the MFT, are likely to play a role in hampering the formation of terraced structures. Avoiding terracing is a further advantage of the present design that makes our approach potentially interesting in applications seeking at forming flat interfaces. Fig. 8 also investigates the effect of systematic errors of our evaluation of the configurational volume available to colloids in the crystal phase ($v_{0}$) on the MFT predictions.

We can use our quantitative MFT for design purposes. In SI Fig. 11 we extend the analysis of Fig. 8 and report an extensive study of the expected number of layers, $\langle\Gamma\rangle$, as a function of $\Delta G_{0}$, $N_{\mathrm{L}}$, and $\rho_{\mathrm{id}}$. As expected, $\Gamma$ diverges as we approach the bulk phase boundary (see Figs. 3$b$ and $c$) from the gas phase either by increasing the number of linkers, $N_{\mathrm{L}}$, or the chemical potential, $\rho_{\mathrm{id}}$. Moreover, we can also investigate the effect of changing the strength of particle–surface relative to particle–particle hybridization free energy, $K=\Delta G_{0}-\Delta G^{\mathrm{s}}_{0}$. SI Fig. 12 shows how, by increasing $K$, the system sharply transits from not forming any layer (the surface cannot denaturate loops) to a constant thickness (no loops left on the particles of the first layer). SI Figs. 11, 12 clarify how, generally, the averaged number of layers is not a function of $\Delta G_{0}$. This highlights the fact that layering is solely controlled by entropic factors (namely combinatorial terms and the ratio between the statistical weight of loops and bridges) as already observed for the bulk gas–solid transition (see Fig. 3$b$ and $c$).

IV Discussions and Conclusions

The results of the previous section show the robustness and flexibility of our design strategy to achieve colloidal layer deposition. Our system provides multiple tuning parameters allowing to control the thickness of crystals self-assembled at the functionalized surface. For instance, as shown in Fig. 5 and SI Fig. 11 and 12, using the number of linkers per particle, $N_{\mathrm{L}}$, as a design parameter could elude experimental limitations on the possible values of the chemical potential ($\rho_{\mathrm{id}}\sim\exp[\beta\mu]$). We notice how $N_{\mathrm{L}}$ has been already employed in recent experiments (Ref. Bachmann et al. (2016)) to control aggregation of DNA functionalized vesicles resulting in a phase diagram that we have successfully reproduced in the present study (see Fig. 4).

The rigorous free-energy calculations presented in the previous sections certify the generality of the mechanism leading to controllable layer deposition. The only two indispensable ingredients underlying the sought effect are the mobility of the binders and the presence of competition between inter– and intra–particle linkages. Nevertheless, in systems of DNA mediated interactions, kinetic limitations may hamper yielding, for instance, of thermodynamically stable crystals in favor of disordered aggregates. In the last five years, several design strategies have been proposed leading to enhanced crystallization. In particular, Van der Meulen and Leunissen clarified the advantages of using mobile linkers when aiming at yielding regular structures as due to the possibility of bound particles to pivot around each other van der Meulen and Leunissen (2013).

Beyond defect annihilation, the finite rates at which DNA linkages form and break could affect self-assembly kinetics. We recently showed how, at low temperature, colloidal systems carrying mobile binders struggle to form compact aggregates as due to small DNA denaturation rates impeding re-distribution of interparticle bridges and effectively stabilizing low valency aggregates Jan Bachmann et al. (2016). To probe the effect of finite reaction rates, in Fig. 9 we report the outcomes of expensive simulations in which DNA reaction dynamics is explicitly simulated using a Gillespie algorithm fed by the rates of forming/breaking DNA linkages (see the EXP method Sec. S3: Simulation details). The reaction kinetics is parametrized by the $on$ rate of free reactive oligomers in solution, $k^{0}_{\mathrm{on}}$, expressed in units of $D\cdot\rho_{0}^{-1}\cdot L^{-2}$, where $D$ is the diffusion coefficient of diluted colloids and $\rho_{0}$ the standard concentration. Finite values of $k^{0}_{\mathrm{on}}$ sensibly slow down adsorption to the point that steady states cannot be reached using affordable simulations. For $k^{0}_{\mathrm{on}}\to\infty$ we recover the results found using Brownian dynamics simulations calculating the number of linkages using chemical equilibrium equations (IMP in Fig. 9). Importantly Fig. 9$b$ shows that, in the present system, finite reaction kinetics do not hamper the relaxation of the system toward the equilibrium state. We also notice that it is possible to design suitable DNA sequences capable of converting intra–particle loops into inter–particle bridges without the need of denaturating DNA. In particular, the scheme that has been recently validated by Parolini *et al * Parolini et al. (2016) or by Zhang et al Zhang et al. (2017), based on the toehold-exchange mechanism of Zhang and Winfree Zhang and Winfree (2009), could be readily implemented in the present system to speed up $k^{0}_{\mathrm{on}}$ by orders of magnitude Zhang and Winfree (2009).

In conclusions, in the present contribution, we presented a system capable of self-assembling fcc (111) crystals at functionalized surfaces with finite thickness fully controlled by design and thermodynamic parameters. Beyond addressing an important technological challenge, the proposed system is innovative in that it provides a general design principle to self-assemble finite and localized colloidal aggregates. Often self-assembly of finite structures relies on the use of directional interactions Patra and Tkachenko (2017); Halverson and Tkachenko (2013) involving many different particles Zeravcic et al. (2017, 2014). Such design is experimentally challenging, mainly due to the need of designing many orthogonal pair-interactions with comparable strength. Our strategy provides a valuable alternative to the use of multicomponent systems. Our design principle uses particles that are reprogrammed and interact differently when in contact with functionalized surfaces. We are therefore confident that our work will inspire new investigations and experiments leading to functional materials made of dynamic, trajectory-dependent building blocks.

Acknowledgements

This work was supported by the Fonds de la Recherche Scientifique de Belgique - FNRS under grant n∘ MIS F.4534.17. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CECI), funded by the Fonds de la Recherche Scientifique de Belgique - FNRS under grant n*∘* 2.5020.11.

Supplementary Figures

Supporting Information

S1: Calculation of the multivalent free energy of the system

We calculate the partition function of the system at a given colloids’ configuration $\{{\bf r}\}$ and number of linkages $\{n\}$ (see Main Eq. 1) leading to the multivalent free energy ${\cal F}_{\mathrm{multi}}(\{{\bf r}\},\{n\})=-k_{B}T\log{\cal Z}(\{{\bf r}\},\{n\})$. We start by computing the number of ways of making a certain set of linkages, $\{n\}=\{n^{\mathrm{AB}}_{ij},\,n^{\mathrm{AB}}_{ji},\,n^{\mathrm{AB}}_{ii},\,n^{\mathrm{AC}}_{i}\}$ with $i=1,\cdots N_{\mathrm{p}}$ and $i<j$. First, we count the number of ways, ${\cal N}_{\mathrm{select}}$, of selecting the binders (ligands and receptors) used to form a certain type of linkage (e.g., a bridge of type AB between particle $i$ and particle $j$, see left and central panel in Fig. 13). The number of ways of selecting $n^{\mathrm{AB}}_{ij}$, $n^{\mathrm{AB}}_{ii}$, and $n^{\mathrm{AC}}_{i}$ A ligands on particle $i$ to be used to form, respectively, the bridges with particle $j$ ($j=1,\cdots,N_{\mathrm{p}}$), the loops, and the surface–particle bridges is

[TABLE]

where $n^{\mathrm{A}}_{i}$ is the number of free A ligands

[TABLE]

Similarly, the number of ways of selecting the B linkers on particle $i$ and the receptors are

[TABLE]

where $n^{\mathrm{C}}_{\mathrm{s}}$ is the number of free receptors on the surface. Finally, ${\cal N}_{\mathrm{select}}$ reads as

[TABLE]

After selecting the ensemble of ligands/receptors to be used to form linkages $\{n\}$, we calculate the number of ways, ${\cal N}_{\mathrm{pair}}$, of reacting the pre–selected binders (see right panel of Fig. 13). Noticing that there are, for instance, $n^{\mathrm{AB}}_{ij}!$ ways of binding $n^{\mathrm{AB}}_{ij}$ A ligands on particle $i$ with $n^{\mathrm{AB}}_{ij}$ B ligands on particle $j$, we obtain

[TABLE]

Finally, the combinatorial term associated to configurations with a given set of linkage $\{n\}$ is

[TABLE]

To derive the final expression of ${\cal Z}$, we multiply ${\cal N}_{\mathrm{comb}}$ by the Boltzmann factor accounting for the hybridization free energy of binding a set $\{n\}$ of linkages (see Main text for the definitions)

[TABLE]

${\cal Z}_{T=\infty}$ denotes the partition function of the system when no linkages are formed (as found at high temperature). ${\cal Z}_{T=\infty}$ is related to entropic, repulsive forces due to the reduction of the configurational space of the DNA linkers compressed by two approaching surfaces. Neglecting excluded volume interactions between linkers, as often done when modeling DNA mediated interactions, and defining by $\Omega_{0}$ and $\Omega_{i}(\{{\bf r}\})$ the volume of the configurational space available, respectively, to a single ligand tethered to particle $i$ at infinite dilution and at finite density (similar definitions follow for the configurational space of receptors tethered to the surface, $\Omega^{\mathrm{s}}_{\mathrm{0}}$ and $\Omega^{\mathrm{s}}(\{{\bf r}\})$) we find

[TABLE]

In the previous expression, $V(\{{\bf r}\})$ accounts for additional (e.g. electrostatic) interactions between colloids. We use $V(\{{\bf r}\}))$ to regularize hard–core repulsions with negligible effects on the final results. Below (see Sec. S1.A), we report the expression of $V(\{{\bf r}\})$ that has been used in this work.

In this work we consider reactive sequences tethered to particles’ surfaces through short, thin rods of double-stranded DNA of length $L$. When $L$ is much smaller than the radius of the particles, the reactive sequences of unbound linkers are uniformly distributed within the layer of thickness $L$ surrounding the tethering surfaces. We then have $\Omega_{0}=4\pi R^{2}L$, $\Omega^{\mathrm{s}}_{0}={\cal A}L$ (where ${\cal A}$ is the area of the surface). Similarly, the configurational volumes defining the hybridization free energies (see Main Eqs. 4) read as

[TABLE]

where $e_{ij}$ and $e^{\mathrm{s}}_{i}$ are, respectively, the volume excluded to the reactive sequence of a linker tethered to particle $i$ by the presence of particle $j$ and the surface (see Main Fig. 2). $k^{\mathrm{s}}_{i}$ is the volume excluded to a reactive sequence tethered to the surface by the presence of particle $i$ (see Main Fig. 2). $v(i)$ and $v_{\mathrm{s}}$ are the lists of particles in direct contact with ligands on particle $i$ and receptors on the surface. Similarly, the configurational space of bound sequences ($\Omega_{ij}$ and $\Omega^{\mathrm{s}}_{i}$) is the volume of the overlapping regions spanned by the reacting sequences before binding (see Main Fig. 2). We report the explicit expression of the terms appearing in Main Eqs. 4 and Eqs. 17 in Sec. S1.B.

At given colloid positions, $\{{\bf r}\}$, the most likely number of linkages featured by the system, $\{\overline{n}\}$, are calculated by maximizing the multivalent free energy

[TABLE]

The previous set of equations, along with the definition of ${\cal F}_{\mathrm{multi}}$, Eq. 15, lead to the chemical equilibrium equations reported in Main Eqs. 3. Main Eqs. 3 become equivalent to the following set of equations for the number of unbound linkers

[TABLE]

that are used to implement self-consistent calculations in our simulations (see Sec. S3: Simulation details). When written in term of the stationary number of linkages, the multivalent free energy simplifies into the expression reported in Main Eq. 5.

S1.A: Modeling hard–core repulsion

We use smooth pair potentials to regularize particle–particle and particle–surface hard–core repulsions ($V(\{{\bf r}\})$ in the Methods section of the main text). Following previous investigations, the repulsion between colloids is modeled using

[TABLE]

where $r_{ij}$ is the distance between particle $i$ and $j$, and $e_{ij}$ is defined below (Eq. 29). Similarly, if $r_{i,\mathrm{z}}$ is the distance of particle $i$ from the surface, the surface-particle repulsion is modeled using

[TABLE]

where $e^{\mathrm{s}}_{i}$ is defined in Eq. 30. Finally

[TABLE]

Such smooth regularizations allow using larger integration steps $\Delta t$ (see the Methods section of the main text) with negligible effects on the results of the manuscript. The latter claim follows from the fact that the typical surface–to–surface distance is comparable with $L$ while $V_{\mathrm{pp}}(r_{ij})>0$ and $V_{\mathrm{ps}}(r_{i\mathrm{z}})>0$ only when $r<0.75\cdot L$. $V_{\mathrm{pp}}$ and $V_{\mathrm{ps}}$ are effective potentials resulting from coating particles with inert strands (not carrying sticky ends) of length $0.75\cdot L$. This observation justifies the particular choice of $V$ given the fact that, in experiments, inert binders are often used to screen non-selective attractions (e.g., van der Waals forces).

S1.B: Configurational terms

Below we report the analytic expressions of the overlapping volumes defined in Main Fig. 2$b$ and used to calculate the hybridization free energies and the $on$ rates (see the Methods section in the main text)

[TABLE]

where $v_{\mathrm{par,par}}(r,R_{1},R_{2})$ and $v_{\mathrm{par,surf}}(r,R)$ are, respectively, the overlapping volume between two spheres of radius $R_{1}$ and $R_{2}$ placed at distance $r$ and the volume of a spherical cap of radius $R$ with base placed at distance $r$ from the center of the sphere:

[TABLE]

S2: Mean Field Theory

We detail the calculation of $f_{\mathrm{multi}}$ (see Main Eqs. 9) used to predict the probability of self-assembling crystals with $\Gamma$ layers (see $P(\Gamma)$ defined in Main Eq. 10). $f_{\mathrm{multi}}(\Gamma)$ is the free energy of an fcc (111) crystallite comprising $\Gamma$ layers as compared to a reference state in which particles are isolated and only feature loops, normalized by the number of particles in direct contact with the functionalized surface. We first calculate the number of free binders and free receptors on the particles belonging to layer $i$ along with the number of free receptors ($\overline{n}^{\mathrm{A}}_{i}$, $\overline{n}^{\mathrm{B}}_{i}$, and $\overline{n}^{\mathrm{C}}_{\mathrm{s}}$ in Fig. 14). We consider particles distributed on an fcc (111) crystal with fixed particle–particle and particle–surface distance. Therefore, the free energy of making inter–particle bridges ($\Delta G^{\mathrm{b}}$) is constant and is calculated using Main Eq. 4. Similarly, $\Delta G^{\mathrm{l}}$ and $\Delta G^{\mathrm{b,s}}$ are the hybridization free energies of forming loops and particle–surface bridges. We calculate $\Delta G^{\mathrm{b,s}}$ using a surface area $A^{\mathrm{(MF)}}$ corresponding to the averaged surface per surface–bound particle as sampled in a representative simulation (notice that $A^{\mathrm{(MF)}}$ affects $\Delta G^{\mathrm{b,s}}$ through $\Omega^{\mathrm{s}}$, see Main Eq. 4). Accordingly, we fix the number of receptors $N^{\mathrm{(MF)}}_{\mathrm{C}}$ to $N^{\mathrm{(MF)}}_{\mathrm{C}}=\sigma_{R}\cdot A^{\mathrm{(MF)}}$, where $\sigma_{R}$ is the receptor density.

Using Eqs. 19 we write the total number of free binders on the surface and on the particles belonging to the first layer as

[TABLE]

where in the expression of $\overline{n}^{\mathrm{A}}_{1}$ and $\overline{n}^{\mathrm{B}}_{1}$ we set $\overline{n}^{\mathrm{A}}_{2}=\overline{n}^{\mathrm{B}}_{2}=0$ when calculating the free energy of single–layer crystals ($\Gamma=1$ in Fig. 14). For particles in the intermediate layers $1<i<\Gamma$ we have

[TABLE]

while for the last layer (if $\Gamma>1$)

[TABLE]

We define by $\overline{n}^{\mathrm{A}}_{0}$ and $\overline{n}^{\mathrm{B}}_{0}$ ($\overline{n}^{\mathrm{A}}_{0}=\overline{n}^{\mathrm{B}}_{0}$) the number of free binders of particles isolated in bulk (therefore featuring $N_{\mathrm{L}}-\overline{n}^{\mathrm{A}}_{0}$ loops). $\overline{n}^{\mathrm{A}}_{0}$ or $\overline{n}^{\mathrm{B}}_{0}$ is calculated by setting $\Delta G^{\mathrm{b}}=+\infty$ and $\Delta G^{\mathrm{b,s}}=+\infty$ in one of the previous equations.

The free energy of $\Gamma$ isolated particles in bulk reads as (see Main Eq. 5)

[TABLE]

where $N_{\mathrm{L}}-\overline{n}^{\mathrm{A}}_{0}$ is the number of loops featured by each particle. On the other hand, the free energy per surface–bound particle of crystals with $\Gamma$ layers (see Fig. 14) is given by (see Main Eq. 1)

[TABLE]

Notice that in the configuration of Fig. 14 particles are sufficiently distanced that $F_{T=\infty}=0$. Finally, $f_{\mathrm{multi}}(\Gamma)$ (Main Eq. 9) used in the definition of $P$ (Main Eq. 10) and $\Delta f_{\mathrm{multi}}(\Gamma)$ (see Main Fig. 5a-c) are given by

[TABLE]

To calculate the equilibrium layer distribution $P(\Gamma)$, we also need to estimate the entropic loss of caging colloids from the gas phase into a site of the crystal. We employ a cell model in which we assign a configurational space volume equal to $v_{0}$ to each particle in the solid phase. Following Main Refs. [42,43] we use $v_{0}=(L/2)^{3}$, where we identify $L$ with the interaction range of a square well potential. We study the sensitivity of our results to $v_{0}$ in Main Fig. 8. The probability of assembling $\Gamma$ layers from a diluted colloidal suspension at density $\rho_{\mathrm{id}}$ reads as $P(\Gamma)={1/Z_{\mathrm{MFT}}}\cdot\exp[-\beta f_{\mathrm{multi}}(\Gamma)](\rho_{\mathrm{id}}v_{0})^{\Gamma}$, where $\rho_{\mathrm{id}}\sim\exp[\beta\mu]$, $\mu$ is the chemical potential of the particles, and $\rho_{\mathrm{id}}$ is small enough to justify an ideal representation of the gas phase.

$Z_{\mathrm{MFT}}$ is a normalization factor that is well defined (i.e. $Z_{\mathrm{MFT}}=\sum_{\Gamma=0}^{\infty}\exp[-\beta f_{\mathrm{multi}}(\Gamma)](\rho_{\mathrm{id}}v_{0})^{\Gamma}<\infty$) conditional on the gas phase being stable in bulk. To extract the phase boundary in bulk, we notice that for $\Gamma\to\infty$ surface effects are negligible and $f_{\mathrm{multi}}(\Gamma)$ reads as $\Delta f(\infty)\cdot\Gamma$, where $\Delta f_{\mathrm{multi}}(\infty)$, $\Delta f_{\mathrm{multi}}(\infty)=\lim_{\Gamma\to\infty}f_{\mathrm{multi}}(\Gamma)-f_{\mathrm{multi}}(\Gamma-1)$, is the multivalent free energy per particle in an fcc crystal as compared to the gas phase. The gas–solid boundary in bulk (see Main Figs. 3$b$, 3$c$, and 4) is given by the relation $\exp[-\beta\Delta f_{\mathrm{multi}}(\infty)]\rho_{\mathrm{id}}v_{0}=1$ or $\beta\Delta f_{\mathrm{multi}}(\infty)=\log(\rho_{\mathrm{id}}v_{0})$ that, in the diluted limit $\rho_{\mathrm{id}}\to 0$, matches existing cell models Main Ref. [43].

S2.A: Analytic predictions

In this section we extract compact analytic expressions allowing to estimate $\Delta f_{\mathrm{multi}}$ (see Eq. 36). We consider the low temperature regime in which $\Delta G_{0}\to-\infty$. Using Main Eq. 5, we calculate the statistical weight of, respectively, particle–surface and particle–particle bridges divided by the statistical weight of intra–particle loops as

[TABLE]

In the previous expressions, we considered that for the colloidal arrangement used in the MFT (see Fig. 14) the configurational space of free binders is not excluded by any colloid or the surface ($\Omega_{i}=\Omega_{0}$ and $\Omega^{\mathrm{s}}=\Omega^{\mathrm{s}}_{0}$). Note that $\Omega_{ij}$ and $\Omega^{\mathrm{s}}_{i}$ are not a function of the specific particle $i$ given that particle–particle and particle–surface distances are kept constant (see Eq. 27). Once the first layer of particles is formed, each particle will present a number of free ligands equal to

[TABLE]

We now assume that when particles are added to the second layer, the number of particle–surface bridges as well as the number of lateral bridges between particles in the first layer do not change. Such approximation allows calculating $\Delta f_{\mathrm{multi}}(2)$ only using $M_{1}$. Similarly, by calculating the number of free linkers featured by particles in the second layer ($M_{2}$, which is only function of $M_{1}$), we can re-iterate the calculation of $\Delta f_{\mathrm{multi}}(3)$ for particles in the third layer. In general, the free-energy gain of adding layer $i$ when particles in layer $i-1$ express $M_{i-1}$ free linkers is

[TABLE]

where

[TABLE]

The number of free ligands per particle expressed by layer $i$ before attaching layer $i+1$ reads as

[TABLE]

Fig. 15 shows that the results of the analytic theory match the MFT predictions.

S3: Simulation details

We first calculate the force acting on particle $i$, ${\bf f}_{i}$ in Main Eq. 6 and 7. When using the implicit (IMP) scheme (see Main Sec. 2.2), the numbers of possible linkages are fixed to their most likely values, $\{n\}=\{\overline{n}\}$ (see Main Eqs. 12, 13). The force acting on particle $i$ then reads as

[TABLE]

where the second equality follows from the saddle point equations (Main Eq. 2) and the fact that the only direct dependency of ${\cal F}_{\mathrm{multi}}$ on $\{{\bf r}\}$ is due to $F_{T=\infty}$ (see Eqs. 15 and 16). In particular, we find

[TABLE]

By using the definitions of $\Delta G$ we find

[TABLE]

so that Eq. 43 becomes

[TABLE]

Using Main Eqs. 11, along with the fact that $\Omega^{\mathrm{s}}_{j}$ and $\Omega_{jq}$ are only function, respectively, of ${\bf r}_{j}$ and $\{{\bf r}_{j},\,{\bf r}_{q}\}$ we find

[TABLE]

When using the EXP algorithm (see Main Sec. 2.2), the linkages are evolved using the Gillespie algorithm (see next section). In this case, ${\bf f}_{i}$ is calculated as in the last equality of Eq.  42 in view of the fact that $\{n\}$ are not a function of $\{{\bf r}\}$. Therefore, the expression of ${\bf f}_{i}$ used in the EXP case is identical to Eq. 49 when replacing $\{\overline{n}\}$ with the actual values of the linkages $\{n\}$.

Concerning the grand-canonical Monte Carlo algorithm, insertion/removal acceptances are given by

[TABLE]

where $\Delta{\cal F}_{\mathrm{rem}}$/$\Delta{\cal F}_{\mathrm{ins}}$ is the change of the system free energy after removing/inserting a colloid in the simulation box. To optimize the acceptances, we set $\Delta{\cal F}_{\mathrm{rem}}$ and $\Delta{\cal F}_{\mathrm{ins}}$ to zero when the inserted/removed particle does not interact with any other particle or the surface.

S3.A: Gillespie algorithm

In this section, we detail the implementation of the Gillespie algorithm that has been used to simulate sticky-ends reactions in the EXP method.

At a given colloid configuration, $\{{\bf r}\}$, we start calculating all $on$/$off$ rates of making/breaking linkages

[TABLE]

as derived in Main Eq. 8. Notice that, for instance, $k_{\mathrm{on}}^{ij}$ is the $on$ rate of making a bridge between $i$ and $j$ either using an A sticky end tethered to particle $i$ or $j$. While $off$ rates are only function of $\Delta G_{0}$ or $\Delta G^{\mathrm{s}}_{0}$, $on$ rates also include configurational terms that are function of $\{{\bf r}\}$ (see Main Eq. 8). Accordingly, the list of all possible reactions is specified by the following affinities

[TABLE]

where, for instance, $a_{\mathrm{on,AB}}^{ij}$ refers to the possibility of forming a linkage between $i$ and $j$ using an A sticky end tethered to particle $i$. We then fire one within all possible reactions with probability

[TABLE]

where $a_{\mathrm{tot}}$ is the total affinity. Along with the type of reaction we sample the time for it to happen ($\tau$), distributed as $P(\tau)=\exp[-\tau/a_{\mathrm{tot}}]/a_{\mathrm{tot}}$, and increment a reaction clock $\tau_{\mathrm{reac}}$ by $\tau$. If $\tau_{\mathrm{reac}}<\Delta t$ (where $\Delta t$ is the simulation step), we update $\{n\}$, recalculate the affinities (Eq. 51), and fire a new reaction until reaching $\Delta t$.