Modeling the assembly of oppositely charged lock- and key-colloids

The interaction of oppositely charged lockand key-colloids is investigated using computer simulations. We show that indented spheres, i.e., lock-particles, can be specifically assembled with spherical key-particles using solely electrostatic interactions in addition to a hard overlap potential. An analytic description of the entropic and energetic contributions is derived and supported by simulations and explicit energy calculations, respectively. The analytic expression of the electrostatic contribution is further employed to build up a schematic model allowing for efficient large-scale Monte Carlo simulations. The influence of the charge/ionic strength, the degree of indentation, and the size/number ratio is discussed by analyzing the specific and unspecific associations from the simulations. Herein, both particle design and mixing conditions lead to the formation of stable specific clusters analogous to colloidal molecules whose valence is defined by the number of lock-particles associated with a key-particle. In addition, the approach is extended to the encapsulation of an excess of small key-particles in largely indented lock-particles. These two examples exemplify that highly specific pairwise interactions can be implemented by using solely oppositely charged particles with complementary geometries, which opens the road for a rational design of complex hierarchical self-assemblies of complementary building blocks.


Introduction
Self-assembly is a key construction principle that Nature uses so successfully to create often highly elaborate structures. 1 Proteins and peptides have, for instance, the intrinsic ability to spontaneously self-assemble into complex structures such as fibrils, shells or nanotubes, as observed in the formation of cytoskeletal filaments or viral assemblies. [2][3][4][5][6] While spherical colloids have been established as models to investigate various aspects of phase transitions and self-assembly processes, 7-10 colloidal particles with centrosymmetric interaction potentials are often unable to capture the intrinsic complexity of diverse biological materials. This represents a major limitation in any attempts to quantitatively understand their self-assembly and motivates the use of new synthetic colloids such as anisotropic colloids, patchy particles, or well-defined colloidal clusters [11][12][13][14][15][16][17] as novel model systems to study more complex associative processes and develop bioinspired applications and materials.
The lock and key assembly principle, introduced to describe the specific interactions of enzymes with substrates, [18][19][20][21] illustrates the exquisite ability of biological entities to fabricate and regulate their molecular machinery. 4 Herein, the molecular recognition mechanism stemming from shape complementary and pairwise interactions constitutes as well one of the foundations of modern supramolecular chemistry within the formation of specific host-guest and donor-acceptor complexes. 22 Besides, even at the molecular level, assemblies are often dynamic and the constituents may adapt to their local environment so as to maximize their interactions. This marks the evolution of the lock and key theory into the induced fit theory, 20,21 which can be seen in analogy to the more recent development of supramolecular chemistry towards adaptive chemistry. 22,23 A synthetic colloidal analog of a lock, or an acceptor, follows recent developments of indented spheres or bowl-shaped particles. These anisotropic colloids have been investigated in detail due to their rich phase diagram and with respect to their specific self-assembly with complementary spherical key-particles. 24 The prevailing strategies to trigger short range specific interactions driving the assembly of complementary shaped particles have so far relied on depletion interactions. The first experimental realizations could be rationalized by excluded volume theory, which has been extensively supported by computer simulations. [25][26][27][28] However, other interactions for colloidal lock and key selfassembly have rarely been explored so far, with the exception of some works on electric field induced lock-and key-interactions, 29 cell imprint interactions, 30 and point-charge models. 31 Recently, it has been demonstrated that thermoresponsive microgel-based lock-and key-particles may interact to form defined colloidal molecules simply by considering them as oppositely charged building blocks. 32 Other studies also point to electrostatics being an important factor where assembly between lock-and keyparticles is concerned. 33 By taking advantage of the fact that both the geometrical parameters and interactions could be tuned with temperature, the valence of the clusters can be controlled from a tetravalent ''methane''-to a divalent ''carbon dioxide''-like configuration. Oppositely charged microgels can, hence, be regarded as a colloidal analog of adaptive chemistry.
Those experimental findings have motivated this work to further model such complex lock and key assembly. Given the complexity of the experimental systems, we only focus on oppositely charged ''hard'' lock-and key-particles where we assume a uniformly distributed charge throughout the surface of the colloidal molecule, which in many respects makes this model a more realistic alternative to a simple point charge model. We present an analytic expression to describe the energetic contribution to the interaction potential and further build up a schematic model allowing for large scale Monte Carlo simulations of the assembly. The specificity of the association illustrated by the formation of stable specific colloidal molecules is then discussed in relation to the ionic strength, size and number ratio between the lock-and the key-particles, and the degree of lock-indentation.

Results and discussion
Lock-and key-particle definition Let us consider a lock-particle defined as a spherical particle with a radius R L indented by a spherical cavity with a radius R C as schematically depicted in Fig. 1 (left) for two different indentations. Herein, the two sphere centers, L and C, are separated by the distance c. The complementary spherical key-particle is defined using its center position K and its radius R K . Finally, we note that for c 4 R C + R L a lock-particle has no indentation and is equivalent to a sphere. In the following we will, from an analytic perspective, explore the entropic and energetic (electrostatic) contributions that will later be included in the simulations.

Entropic contribution
The free energy DG between a hard lock-and a hard key-particle was previously 32 shown to be where r = |r| = |K À L|, a = +PLK and b = +PLC, where P defines the contact point between the particles, k B is the Boltzmann constant, and T is the temperature. Here the angular dependence has been averaged, and thus the free energy describes the rotational entropic cost at a specific distance r around L. In Fig. 1 (right), DG(r) is presented for two distinct indentations and different R K -values using eqn (1). The free energy increases with an increase of R K and a local minimum can be observed in DG(r) for small R K . That is, for some r there is an attractive entropic force pushing the particles closer together. Worth noting is that a lock-particle can encapsulate a key-particle with no (rotational) entropic cost when c + R K o R C as shown in Fig. 1b (right) for R K o 0.5375R L . These analytic results are in agreement with simulations using lock-and key-particles with hard body conditions, 32 which is illustrated for two cases in Fig. 1.

Electrostatic contribution
The Yukawa interaction-energy between uniformly surfacecharged lock-and key-particles was analytically derived when the key-particle is perfectly in front or behind the lock-particle cavity, i.e., for |(C À L)Á(K À L)|/(cÁr) = 1. An approximate formula for the same setup has been derived elsewhere 30 yet we present an exact expression in eqn (S1) in the ESI. † The energies in the specific LK configuration, U LK , and the unspecific LK configuration, U LK , were calculated for different ionic strengths, I, and size ratios, R K /R L , in Fig. 2. Note how the Debye-length k À1 ¼ e 0 e r k B T= 2e 2 N A I À Á À Á 1 2 is related to the ionic strength I, where e 0 is the vacuum permittivity, e r = 72 is the relative permittivity of the dispersing medium, T = 293.15 K, e is the elementary charge, and N A is the Avogadro constant. The calculations were performed considering R L = R C = c = 500 nm, oppositely charged lock-and key-particles with a constant absolute surface charge-density |s| = 3.183 Â 10 À4 e nm À2 . By defining the effective charge 34-36 of a lock-or key-particle as where X = {L, K} indexes the particle type, a simplified expression for U LK for such a setup can be retrieved as assuming L as the origin. Here R is a point on the rim of the lockparticle defined in the ESI, † and z = |K|. The analytic expression was found to be in agreement with explicit charge calculations performed at three different ionic strengths. Herein, 10 000 equally charged point charges were randomly distributed at the surface of the lock-particle, which yielded the same surface charge-density |s| as above, and the total electrostatic energy was retrieved by summing up all the point-charge interactions with the effective charge of the key-particle q K *. We note a maximum when the key-particle is in front of the lock-particle at R K = R L . So does the difference between the LK and LK configurations, DU. Hence, R K should be equal to R L in order to maximize the interaction-energy difference between specific LK and unspecific LK interactions. Such a configuration is quite costly from an entropic point-of-view, see Fig. 1(a) (right), and it is not straightforward to predict the optimal particle design for self-assembly. Simulations are thus needed to model the assembly process, which then also includes many-body effects.
However, whereas DG is solely defined by the geometries of the constituents, it is important to note that U LK , U LK , and DU are all strongly dependent on I and |s|, cf. DU p |s| 2 , and can thus be tuned over many units of k B T. A specific lock and key selfassembly is therefore expected, at least in the limit of low I and high |s|. In order to perform simulations of lock-and key-particles, one needs an angular-dependent potential accounting for the lock-particle anisotropy. This potential is constructed based on the exact expressions derived in the previous section. By denoting the exact potential in front of the lock-particle as V L (d LK (r,c)) and the analogue expression behind it as V L d LK ðr; cÞ À Á , where c = C À L, the angular-dependent potential V L (r,c) is approximated as Note that when V L is a function of a scalar distance, it is described by eqn (S1), ESI, † whereas when V L is a function of a vector distance (and c), it is described by eqn (4). Here l(r,c) A [0,1] is a function depending on the region where the potential is evaluated in, see Fig. 3 (top-right), and d LK (r,c) and d LK ðr; cÞ are distances to the surface of the lock-particle. A more detailed description of the potential can be found in the ESI. † In order to assess the validity of the approximate potential, we compared it to the potential of a lock-particle where the surface is decorated with 10 000 randomly positioned point-charges, see the schematic in Fig. 3 (top-left). The potentials stemming from these two respective models are displayed in the bottom panels of Fig. 3. Note that along the symmetry axis (x = 0), eqn (4) is exact. The same calculations were repeated for different ionic strengths and are provided in the ESI † ( Fig. S3-S5).
With the exception of some regions tangent to the cavity rim, the approximate potential, V L (r,c), was found to provide an accurate description of the lock-particle electrostatics with a reduced number Fig. 2 Interaction-energy between a lock-particle and a key-particle at contact as a function of the size ratio R K /R L calculated at different ionic strengths (solid lines). Energy calculations for the specific LK-configuration are shown in the top panel and for the unspecific LK-configuration in the middle panel, and the difference between the two energies is displayed in the bottom panel. The gray dotted lines refer to explicit charge calculations using a lock-particle decorated with 10 000 randomly distributed point-charges for ionic strengths of 10 À4 M, 10 À5 M and 10 À6 M. of elements as compared to the explicit charge potential, thus allowing for large-scale simulations. For the purpose of such simulations, we use the effective charge 34-36 of a lock-or key-particle, see eqn (2), to define their interaction energy as while between key-particles it is 4pe 0 e r e Àkjrj jrj : Finally, the interaction energy between two lock-particles was approximated using where the prime indicates the second lock-particle. Note that for this approximation the interaction-energy is slightly overestimated when lock-particles face each other at close distances. From this we gather that the model is at its worst for dense systems, and more appropriate for dilute ones, as it is the case in this work but also in other studies. 37

Methods
We performed Monte Carlo simulations using the Faunus software, 38 where we included particle translation-rotation moves and standard cluster moves centered around key-particles. 39 The probability of a single particle move was equal to the probability of a cluster move. The aim was to obtain simulations with typical acceptance in the range of 30-40%. This was achieved by using a maximum rotational displacement parameter of 0.5 radians for all moves, whereas the maximum translational displacement parameters varied between 0.01R L and 0.5R L for single particle moves and 0.5R L for cluster moves. Simulations were performed using hard lock-and key-particles and the electrostatic interactionenergy presented in eqn (5)- (7), and periodic boundary conditions. The hard overlap potential between lock-particles was modeled by checking for overlap between the surface-arcs of the particles. 40 This algorithm only checks for overlapping surfaces and thus lockparticles could in principle be enclosed within one another without any detectable overlap for polydisperse systems. In our simulations, however, the lock-particles were equally sized in each system. Encasement could, however, take place between lock-and keyparticles (using the previous algorithm), since the key-particles sometimes were small enough to entirely fit within a lock-particle.
To prevent this, we used the volume overlap code implemented in Faunus 38 between lock-and key-particles. The total number of particles in the system was N = 1320 and the number ratio between lock-and key-particles N L /N K varied from 0.5 to 10. The simulations were sampled in N Â 10 5 steps starting from energetically equilibrated configurations (over N Â 10 5 steps). The temperature, relative permittivity of the dispersing medium, and ionic strength were the same as mentioned earlier, and the particle number density was r = 0.1 mm À3 . The implemented model does not account for features arising from counterion-release, [41][42][43][44] and thus we chose the charge-densities s L = Às K = |s| so as to minimize its repulsive effect. Any other |s L | a |s K | would have the free energy minima shifted, which effectively corresponds to an inflated particle. A basic length-unit in the simulations was R L = 500 nm, which was kept constant throughout the simulations, in which other lengths such as R C , c, R K , and k À1 were scaled so as to get generic units. If not mentioned otherwise, then R C = R L . A bond is defined as specific (LK) when the distance between a lockand a key-particle is r r R K + R L ; that is, the key-particle is in the cavity of the lock-particle. Conversely, a bond is defined as unspecific (LK) when R K + R L o r o 1.1(R K + R L ). Different values than 1.1 ranging from 1.05 to 1.20 were tested and found to provide tantamount results within the simulated systems. The simulations in Fig. 1 used N L = N K = 1, no cluster moves, N Â 10 9 simulations steps, and only the hard overlap potential.

Simulations
Typical simulation snapshots generated with N L /N K = 1 at different ionic strengths including the uncharged conditions (Binfinite screening) are shown in Fig. 4 (top). The size of the key-particles was set to R K = 0.8R L , which is slightly smaller than the lock-particles and their cavities. Without effective electrostatics (far left), no assembly was observed (as expected for such hard lock-and key-particles). However, by decreasing I o 10 À4 M, i.e., by increasing the amplitude and range of the interactions, particles start to build up defined clusters. Fig. 4 (bottom) displays simulation snapshots captured at constant I = 10 À5 M for different N L /N K -values. Starting from the left, for N L /N K r 1, there are primarily 'dipolar' clusters consisting of a key-particle specifically associated in the cavity of one lock-particle. By increasing the number ratio (going right), higher order clusters start to form, primarily 'CO 2 ' molecules, i.e., two lock-particles specifically bonded to a single key-particle. Based on these initial results, we conclude that specific clusters can be formed and tuned using the number ratio and ionic screenings. The simulation results shown in Fig. 4 are quantified in Fig. 5 considering the average numbers of specific, unspecific, and total bonds per key-particle as functions of ionic strength/ Debye-length, for different indentations c. The dotted lines in all figures represent the maximas of specific bonds (assuming that each lock-particle is specifically bound to a single key-particle). Larger values thus refer to higher order clusters consisting of several key-particles. Though the total number of bonds seems to increase (fairly) monotonically with a decrease of I, the specific number of bonds does show a maximum around I E 10 À5 M, independent of the cavity center c. After this maximum, unspecific bonds are increasingly favored and the key-particles are effectively shifted from the cavities of the lock-particles toward their backs/ sides. Conversely, upon decreasing I the key-particles shift from unspecific to specific sites. This binding-feature has been shown to occur in experimental setups. 45 From geometrical arguments, we note that small c decreases the number of possible specific bonds per key-particle as compared to large c. However, for larger indentations, the area of the cavity is decreased and thus the interactions become energetically less specific. By going from left to right in Fig. 5, the maximum number of possible specific bonds increases. The top left graph shows a maximum around 2, which is in line with geometrical considerations, for which a key-particle can at most be specifically assembled with 2 lock-particles. By increasing c to 1.4R L , this number gets larger with a maximum for N L /N K = 10 at E2.33. The simulation results show that basically 2/3 of the clusters are linear ''CO 2 '' molecules and 1/3 have a trigonal planar geometry like SO 3 2À . For N L /N K = 2, however, most of the clusters present a ''CO 2 ''configuration. The ''valency'' of the key-particle centered clusters can therefore be tuned by adjusting the lock-indentation, the number ratio, and/or ionic strength. As has been noted before for spherical particles, 10 we confirm also for lock-and key-particle assembly of various degrees of indentation (by observing configurations) that the geometrical formation of the decorated particles most often diverges from the optimal positions retrieved from using k = 0 as is in the Thomson problem. 46,47 In Fig. 6, the number of bonds per key-particle is plotted as a function of the center of the cavity c. For N L /N K = 1 and c r 1.6R L , almost all the lock-particles are specifically bound to a key-particle. The total number of bonds, though, is larger than unity, indicating the formation of larger order structures, such as dipolar chains as shown in a simulation snapshot in Fig. S6 (left) (ESI †). Such dipolar chains are not observed using N L /N K = 2 in Fig. S6 (right) (ESI †), where rather different higher order structures appear. This points to the importance of the stoichiometry for the generation of defined hierarchical self-assemblies. The total number of bonds has a maximum at c = 1.6R L (the same as the maxima for the specific bonds) after which the contacts become mostly unspecific. These results indicate that the presence of an indentation may promote the assembly process and increase even the total number of bonds when compared to spherical particles. We rationalize this observation by noting that the entropic cost decreases seemingly more than does the attractive electrostatic interactions as the cavity is shifted from c = R L towards c = 1.6R L . After this point the slight energetic difference between a specific and an unspecific bond is too small to overcome even moderate entropic costs. Thus, for c E 1.6R L , the balance between the two contributions provides the maximum yield of specific and total bonds. In the ESI, † we present more simulation results analogous to Fig. 4-6 but using R K = R L and R K = 1.2R L .
In Fig. 7, the number of bonds per key-particle is shown as a function of the radius of the key-particle. The total number of bonds increases with an increase of R K . The surface-area of the key-particle increases quadratically with R K , which geometrically allows for more lock-particles to associate. Also, since the surface charge-density is constant, the total charge of the particle as well increases quadratically with R K . Yet, going from small to large R K , there is an initial increase of the number of specific bonds, which quickly decays toward zero after its peak value at R K E 1.0 À 1.1R L . This peak is in fact rather a quite broad plateau when N L /N K r 1, whereas for larger N L /N K -values, a single maximum is observed. Note that for N L /N K Z 1, a local plateau is as well visible for the transition around 2 specific bonds per key-particle, which is geometrically strongly constrained by the degree of indentation. Thus, depending on the number ratio, the impact of the size of the key-particles on the cluster conformation may be more or less pronounced.
In Fig. 8 we present configurations from simulations where the radius of the cavity perfectly matches the radius of the keyparticle, i.e. R C = R K (left), and when it does not, i.e. R C o R K (right). The total number of bonds per key-particle is roughly the same, B4; however, for a geometrically perfect match between the cavity and the key-particle, the specific number of bonds is B4, while the same number using R C o R K is B0. The R C /R K ratio defines the geometrical complementarity of the two components, and is thus of great importance for the assembly specificity, keeping in mind that the energetic contribution is maximum for matching curvatures, i.e. R C /R K = 1. Tuning this parameter therefore allows switching from specific    The left configuration is sampled by using c = R C = R K (i.e. perfect geometric match between the cavity and key-particle) and the right by using c = R C = R L . to unspecific bonds. Similar results were obtained for R K = 2R L and N L /N K = 10 as shown in Fig. S12 in the ESI. † Having established the influence of diverse parameters on the specific lock-and key-particle self-assembly, simulations can now be compared to recent experimental results. Herein, oppositely charged microgel-based particles were associated into specific colloidal molecules as shown in Fig. 9a. Due to their respective thermosensitivity, the size and shape of the particles could be adjusted with temperature as well as the valence and conformation of the colloidal molecules (Fig. 9a). We refer to a previous study for further information on the experimental systems. 32 Interestingly, the model presented in this work was able to capture most of the assemblies by fixing the surface charge-density and implementing a size ratio and degree of indention roughly corresponding to different deswellings, see Fig. 9b. In both experiments and simulations, trimeric, ''SO 3 2À ''-like and tetrameric ''CH 4 ''-like colloidal molecules are observed in the swollen regime corresponding to R L E R K and to small lock-indentations with a transition to ''CO 2 ''-like molecules for lower size ratios and larger indentations. Although the model captures most of the features, the experimental realizations differ in some aspects. The microgel based systems are soft and deformable. Consequently, the entropic cost is expected to be much lower than that in the simulations. In analogy to the induced fit theory, the microgels may as well adapt their configuration to maximize their contact area. Therefore, we speculate that the interactions for soft particles may become even more specific. In addition, such systems experience van der Waals and hydrophobic interactions at higher temperatures in their collapsed configurations, which thus promotes specific interactions even more. As an example, the keyparticle is fully enclosed within two lock-particles, forming an ''H 2 ''-configuration in the experiments performed at 48 1C, whereas simulation results exhibit a much looser ''CO 2 ''-configuration.
To conclude this study, we have tested the possibility of encapsulating small key-particles into highly indented lockparticles schematically depicted in Fig. 1b. We use R K = 0.3R L , which is small enough for the key-particle to be enclosed with no entropic cost as discussed earlier. Two simulation snapshots obtained for two different surface charge-densities at I = 10 À5 M with N K /N L = 4 are shown in Fig. 10. In the ESI, † we present similar results for N K /N L = 1 (Fig. S11, ESI †). At |s| = 3.183 Â 10 À4 e nm À2 , most of the key-particles were found to be freely dispersed. By charging up both types of particles by a factor of ffiffiffiffiffi 10 p , the entropic barrier for encapsulation was overcome, which leads to an almost complete encapsulation of all keyparticles by ''Pac-Man''-like lock-particles.
Finally, we comment on the effect of counterion-release, which can either promote the association 42 or add an additional repulsion in the case of particles presenting different absolute surface charge-densities and perfectly matching curvatures. 43 To avoid the latter, we used equally magnified and oppositely polarized charge-densities. Yet, an additional attraction stemming from the counterion-release is to be expected, contributing mostly at large distances and for high charge-densities. 44 Since at small distances this effect becomes negligible, such a mechanism is not the origin of the cluster stability. Fig. 9 (a) 2D CLSM micrographs of a dispersion using N L /V E 0.068 mm À3 and N L /N K E 6. To the right are schematic representations illustrating the specific self-assembly in colloidal molecules with a valency of: B4 (CH 4 ; top), B2 (CO 2 ; middle), and a collapsed state with B2 (H 2 ; bottom). 32 Scale bars in (a): 1 mm. (b) Results from simulations using N L /V E 0.086 mm À3 and N L /N K E 6 obtained for different size ratios and degrees of indentation.

Conclusion
We have derived an exact analytic expression for the electrostatic potential perfectly in front of and behind a lock-particle, and from this formula approximated the angular-dependent potential. Based on this potential, we have performed simulations of oppositely charged complementary shaped lock-and key-particles. The influence of the charge/ionic strength, the degree of indentation, and the size/number ratio was systematically investigated. We found that a combination of electrostatic interactions with a hard overlap potential is sufficient to specifically assemble and stabilize the two components and form colloidal molecules of different valences, as has been shown in experiments elsewhere. The energetic electrostatic contribution is maximized when the radius of the key-particle matches the size of the cavity; however, this also entails a high entropic cost. Nonetheless, we found that specific interactions are greatly promoted for particles presenting matching curvatures.
Our study provides an alternative to colloidal lock and key self-assembly, which usually relies on short range depletion (excluded volume) interactions. By solely altering the geometry of the particles, electrostatics can induce specific ensembles into well-defined clusters, a feature which further may be amplified or modified by the direct effects of charge-scaling and choice of screening-length. Therefore the model describes a versatile and multifaceted mechanism in understanding and directing selfassembly that could be easily applied to other pair-wise selective interactions and geometries. In addition to the model used for the simulation, the underlying analytic expressions derived in this work indeed allow extensions of the presented approach to a large variety of systems, among others multi-indented lockparticles or particles with an in-built surface roughness.

Conflicts of interest
There are no conflicts to declare. Fig. 10 Configurations from simulations using N K /N L = 4, I = 10 À5 M, c = 0.1R L , R C = 0.95R L , and R K = 0.3R L . The left configuration is sampled by using the absolute surface charge-density |s| as described earlier, whereas the right configuration used ffiffiffiffiffi 10 p jsj.