A general theory of DNA-mediated and other valence-limited interactions

We present a general theory for predicting the interaction potentials between DNA-coated colloids, and more broadly, any particles that interact via valence-limited ligand-receptor binding. Our theory correctly incorporates the configurational and combinatorial entropic factors that play a key role in valence-limited interactions. By rigorously enforcing self-consistency, it achieves near-quantitative accuracy with respect to detailed Monte Carlo calculations. With suitable approximations and in particular geometries, our theory reduces to previous successful treatments, which are now united in a common and extensible framework. We expect our tools to be useful to other researchers investigating ligand-mediated interactions. A complete and well-documented Python implementation is freely available at http://github.com/patvarilly/DNACC .


I. INTRODUCTION
The theoretical study of the statistical behavior of soft matter systems generally begins with a description of the interaction between their constituent particles. One exciting system that has emerged in recent years is that of DNA-coated colloids (DNACCs), where interactions between colloids are driven by partial hybridization of complementary ssDNA strands grafted on their surface [1][2][3][4] . To date, DNACCs have been used mainly for building highly sensitive biomedical probes 5 and to assemble binary nanoparticle crystals of numerous morphologies [6][7][8][9] .
The aim of this paper is to present a general theory of DNACC interactions that can be applied to any DNACC system, using only mild assumptions that are valid for common DNACCs. A key feature of our approach is that once the binding statistics of an arbitrary pair of DNA strands is known, all combinatorial and competitive effects can be computed separately. Disentangling these two aspects is what gives our theory its generality. Our efforts build on past theoretical attempts to understand DNACC experiments in terms of the statistical mechanics of bond formation 3,10-18 , which have all partially succeeded in modeling various experimental setups, differing for instance in coating densities, nucleotide sequences, strand lengths, stiffness and colloid sizes. To date, however, a rigorous and general solution to the problem has not been presented [17][18][19] .
A salient feature of DNACCs is that their interaction is mediated by surface-grafted ligands, each of which can form a single bond at a time. Owing to this restriction, we say that such interactions are "valence-limited" 20 . While we phrase our theory in the language of DNACCs, the treatment is equally applicable to any valence-limited interaction, including multi-valent particles for drug delivery, viruses and cells [21][22][23][24] .
The rest of the paper is structured as follows. In Section II, we present a statistical mechanical derivation of our theory in its most general form, followed by a specialization to DNACCs and a computationally simpler mean-field approximation. In Section III, we compare the results of our theory with those from more accurate and computationally expensive Monte Carlo simulations. We have chosen as examples specific DNACC systems that have been presented previously [15][16][17]25 , as well as two novel ones, to show that they can all be described using a single, unifying framework. In Section IV, we summarize our approach and reproduce the essential equations for reference purposes. Finally, in Section V, we discuss possible applications of our theory, both for DNACCs and for more general valence-limited interactions.
In Appendix A, we discuss an important technical aspect of our theory, namely, how to calculate the entropic cost of binding for two general polymeric ligands. We also summarize the results for short, stiff rods (which model short dsDNA tethers common in micron-sized DNACCs 3,4,14 ) and ideal chains (as used in Refs. 17 and 18). Appendix B explains how our theory is connected to the Local Chemical Equilibrium treatment recently introduced by Rogers and Crocker 19 (itself an improvement of the work of Biancaniello, Kim and Crocker in Ref. 3). We show how a physically motivated and simple change to the equations used by these authors corrects their theory's principal shortcoming, and that the corrected version reduces to a specialization of our theory. Finally, Appendix C presents an alternate derivation of the mean-field approximation to our theory using a saddle-point analysis, making explicit the connection between our theory and earlier successful but more specialized approaches 16,25 .
Although this paper explains our approach in full detail, we emphasize that the final results are simple and can be used widely by other researchers to calculate pair potentials and binding configurations of general DNACCs. To this end, we have made freely available a complete and well-documented Python implementation of the theory, which can be found at http: //github.com/patvarilly/DNACC. arXiv:1205.6921v3 [cond-mat.soft] 9 Aug 2012 II. THEORETICAL MODEL Under experimentally relevant conditions, DNAcoated colloids interact through entropic repulsion and hybridization-mediated attraction 3,12,15 . Whereas the repulsion is simple to estimate for common geometries, as shown in Appendix A, calculating the attraction is more complicated. The key subtlety is that the attraction is mediated by valence-limited binding, i.e., each DNA strand can bind at most one partner at a time 26 . In this section, we first present a general treatment of this attraction for generic valence-limited ligands. Then, we specialize our treatment to DNACCs. Finally, we construct a mean-field approximation that is suitable for geometries with translational symmetry.
A. Statistical Mechanics of valence-limited interactions Consider a system of N generalized ligands, depicted in Figure 1(a). At any given time, an individual ligand i can bind at most one other ligand j, with a resulting free energy change ∆G ij . If i and j cannot bind, either for geometric or chemical reasons, we set ∆G ij to infinity. Since i cannot bind to itself, we also set ∆G ii to infinity for all i. The set of all ∆G ij values completely characterizes the relevant chemical detail of the system. We explicitly assume that ∆G ij is independent of the presence of other ligands, which is accurate when the density of ligands is low enough that two ligands not bound to each other can be thought of as not interacting at all. This condition holds for the grafting densities and salt concentrations typically used in DNACC experiments with micron-sized colloids 3,12,14,27 .
Let φ be a binding configuration of the system, i.e. a list of bound ligand pairs (i, j). Taking the state where no bonds are formed as a reference state, the free energy of φ is the sum of the binding free energies of each bond. Summing over all φ, we obtain the exact partition function of the system, Z, relative to that of the reference state, Z 0 : For a particular ligand i, the terms in the sum above can be classified according to the binding partner of i, if any. Let Z −i be the partition function for a system where i is unbound, and let Z −i,−j likewise be the partition function for a system where both i and j are unbound. The sum of the terms in Eq. (1) where i is unbound is Z −i /Z 0 , whereas the terms where i is bound to a ligand j out of the remaining N − 1 ligands is exp(−β∆G ij ) Z −i,−j /Z 0 . With these definitions, it follows from Eq. (1) that where the sum is over the N − 1 ligands j distinct from i. Let p i be the probability that ligand i is unbound: In Eq. (2), the factor Z −i,−j /Z is the probability that both i and j are unbound. Noting that Z −i,−j /Z −j is the conditional probability for i to be unbound given that j is unbound, we can rewrite Z −i,−j /Z as Above, we have approximated the binding of i and j as uncorrelated events. This approximation is accurate if both i and j have many possible binding partners with comparable binding free energies. Using Eqs. (2)-(4), we obtain the following set of N self-consistent equations for the quantities {p i }: When rewritten in the following convenient form, these equations can be solved, for example, by self-consistent iteration: The probability p ij that ligands i and j are bound is then given by Eqs. (6) and (7) are a principal result of this paper. It can be shown that these results are a single-particle analog of the more familiar chemical equilibrium conditions between the concentrations of N species of particles, with binding constants proportional to exp(−β∆G ij ).
Eq. (6) can be recast as a minimization problem, which leads to a more efficient method of solution than fixedpoint iteration. A similar observation has recently been made when solving the equations that arise in umbrella sampling free energy calculations (WHAM/MBAR), which have a similar structure 28,29 . Additionally, this form allows us to prove that Eq. (6) has exactly one solution, which satisfies 0 < p i ≤ 1 for all i by construction. We define f ({p i }), the function to be minimized, as follows: where k ij = exp(−β∆G ij ). The stationary points of f are given by the solutions to the following system of equations, which are equivalent to Eq. (6). We now demonstrate that all stationary points of f are local minima, implying there is a single stationary point that is the global minimum. To this end, we prove that the Hessian of f at a stationary point, ∂ 2 f /∂p i ∂p j , is positive definite, implying that the stationary point is a local minimum, by showing that the quantity The non-negative quantities in square brackets are the elements of a strictly diagonally dominant, symmetric matrix with positive diagonal elements, which is thus positive definite 30 . Indeed, the i th element of its diagonal is equal to 1 + p 2 i k ii ≥ 1 (we allow nonzero k ii here, in anticipation of the mean-field self-consistent theory of Eq. (21) below), while the sum of the off-diagonal terms in the i th row satisfies where the equality and the strict inequality both follow from Eq. (6) at the stationary point. Hence, the righthand side of Eq. (11) is positive for every nontrivial choice of {v i }, so the Hessian of f is positive definite, as claimed.
Having calculated the probabilities of each bond forming (Eq. (7)), we now calculate the resulting attraction free energy, F att , given by This quantity is most easily estimated using thermodynamic integration. One possible integration path is to offset all interaction energies by an amount λ/β, i.e.
Once the free energy is known at one temperature, a more physical integration path can be used to obtain it at any other temperature. To proceed, we denote the dependence of bond energies on β explicitly with the notation ∆G ij (β), and similarly for p ij (β) and F att (β). We obtain (16) In calculating the integral in Eq. (16) from β 0 to β 1 , the free energy of binding at all intermediate temperatures is obtained as well.
Once F att is known, the free energy of entropic repulsion, F rep , must be added to obtain a full interaction free energy (see Appendix A for details). Additionally, it is usually convenient to offset all free energies to yield zero total energy when the colloids are an infinite distance apart, which we do in Section III.
In the limit of weak binding, when all the factors exp(−β∆G ij ) are very small, the above results simplify significantly, as all the complications of binding competition disappear. Discarding all terms of order exp(−2β∆G ij ) and higher, we obtain In particular, the last relation can be used to sidestep the thermodynamic integration in Eq. (15), and is related to a similar "Poisson" approximation used in the literature 3,10,13,17,19 . However, we have previously pointed out that under realistic experimental conditions, this weakbinding approximation is rather poor 18 , so its utility for experiments is limited.

B. Specializing to DNA-coated colloids
For a large class of ligands, binding occurs only in a small portion of the ligand. In DNACCs in particular, the DNA strands that are used consist of a long, inert tether capped by a reactive, short sticky end, shown in Figure 1(b). The thermodynamics of association of DNA oligomers is well understood 31 , so for a pair of strands, i and j, we can predict from their sequence the value of the hybridization free energy of their sticky ends in solution, ∆G 0 ij . Separating out this contribution to ∆G ij , we have where ∆G (cnf) ij is the free energy penalty that arises from the reduced configurational space available to two bound ligands with respect to two free ligands 12,15 . This quantity depends only on the polymer statistics of the tethers and the way in which they are attached to the colloids (e.g., tethered on the surface vs. mobile), but not on the sticky ends, and in many cases is a completely entropic effect. In Appendix A, we list explicit expressions for ∆G (cnf) ij applicable to short, rigid dsDNA tethers grafted on plates. We also show how to use common polymer statistics techniques to estimate ∆G (cnf) ij to arbitrary accuracy for more complicated tethers, such as long, flexible ssDNAs modeled as ideal chains 3,17 .

C. Mean field theory for DNA-coated plate interactions
In some regimes, notably micron-sized colloids with short DNA tethers, the range of the binding interactions is small compared to the radii of curvature of all the colloidal surfaces involved. In such cases, we can use the Derjaguin approximation 32 to calculate particle interaction energies, knowing only the interaction energy density of a pair of uniformly coated parallel plates as a function of plate separation. For this translationally invariant geometry, a substantial simplification is possible if the ligands are treated at a mean field level rather than individually.
The simplified system consists of two parallel plates separated by a distance h. We consider M different types of ligands, labeled by letters a and b below. Two ligand types can differ in their sticky end sequences, tether types or in the plate on which they are grafted. Corresponding to these M types, there are M (M + 1)/2 potentially different solution hybridization free energies, denoted by ∆G 0 ab . Let σ a be the grafting density of a-type tethers, which has units of tethers per unit area. Since we approximated the system as uniform, all a-type tethers have the same probability of being unbound, denoted by p a . By approximating the sum in Eq. (6) by a continuous integral, we obtain a mean-field, self-consistent set of equations where In the last equation, ∆G (cnf) ab is the configurational cost for a fixed a-type ligand to bind a b-type tether grafted at r b . The factor K ab has units of area, and can be interpreted as an average Boltzmann factor for tether binding times the area on plate B within reach of a fixed a-type tether. A relation analogous to Eq. (7) yields the density of bonds between a-type and b-type tethers, σ ab : Similarly, using a relation analogous to Eq. (15), we can use thermodynamic integration to estimate the free energy of attraction per unit area, f att (h), at fixed plate separation: where λ is the integration parameter defined in Eq. (14). Eqs. (21)-(24) are a second principal result of this paper. They constitute a self-consistent mean field theory for the DNA-mediated attraction of uniformly coated parallel plates. In Appendix C, we show that, when specialized to only one or two kinds of bonds, these equations reduce to the earlier self-consistent mean field theories of Refs. 16 and 25. In these particular cases, we also show how Eqs. (21)-(24) may be derived from a saddle-point approximation of a mean-field partition function.
For the case of short, rigid dsDNA tethers, the factors K ab have particularly simple expressions, derived in Appendix B. In particular, let tethers of type a and b have lengths L a and L b . Without loss of generality, let L a < L b . When a and b are on the same plate, whereas when a and b are on opposite plates, When the tethers are modeled as freely-jointed chains of segments of length , a choice that has been suggested for modeling ssDNA tethers 17 , the factors K ab can be approximated with the following expressions, derived in Appendix A 2. Let N a and N b be the number of segments in tethers of type a and b. When a and b are on the same plate, whereas when a and b are on opposite plates, In both of these equations, the error function, erf(x), is defined as As before, when all the factors K ab σ b and σ a K ab are very small, one can implement a weak-binding approximation to sidestep the thermodynamic integration in Eq. (15), which yields However, as pointed out previously, the utility of this approximation under realistic experimental conditions is limited. Frequently, DNACCs are covered uniformly at a fixed density σ with a single type of strand, and only the sticky ends vary from one colloid to the next. In this case, we can obtain an explicit mean-field expression for βf att between a plate grafted with a-type strands and one grafted with b-type strands. Defining J as the dimensionless factor 2K ab σ, we obtain so, This expression makes explicit the breakdown of the weak-binding approximation (in this case, βf att ≈ −σ ab ) when J 1. Unfortunately, even for the simple case of rod-like DNA strands, the above expression does not yield a simple analytic result for the pair potential between to DNA-coated spheres via the Derjaguin approximation. Since this difficulty only increases when we allow for uneven coverages and multiple potential interactions, we believe that our numerical thermodynamic integration procedure is justified.

III. RESULTS AND DISCUSSION
In this section, we compare the results of detailed Monte Carlo simulations (MC) with those derived from our theory, described above. We refer to the results of the theory where each strand is modeled individually (Eqs. (6), (7) and (15)) as the SCT results, and those of the form where the strands are modeled at a mean field level (Eqs. (21)-(24)) as the SCTMF results 33 .
A. Colloids coated with short dsDNA strands The simplest possible DNACC is a flat plate that is uniformly covered by a single type of DNA strand. Interaction energies obtained from plate geometries can be used to approximate interaction energies in more complex geometries by use of the Derjaguin approximation. Although, our theory can be used in arbitrary geometries, in this section we focus on plates for clarity.
The interaction free energy per unit area of two such DNACCs has been calculated previously from Monte Carlo simulations 15,16 . In this case, the tethers were modeled as short, stiff rods of length L = 20 nm, grafted at a density of 1/(20 nm) 2 . Figure 2 shows the fraction of bonds formed between the two plates when separated by a distance h = L, as a function of the sticky end solution hybridization free energy, ∆G 0 . The results are calculated both from Monte Carlo simulations and from our theory, using an identical realization of the tether grafting positions on a periodic 200 nm × 200 nm plate. The figure inset shows the free energy per unit area of these plates, as a function of plate separation, for three different values of ∆G 0 .
For the grafting density used in Figure 2, each tether typically has around 9 possible binding partners, so the mean-field approximation is excellent. Small discrepancies arise only at low temperature and large plate separation, where the number of possible binding partners is low and Poisson fluctuations in local coverage density are relevant. Such fluctuations are naturally accounted for in the explicit self-consistent theory, whose predictions are thus nearly indistinguishable from MC simulations in all cases. Figure 3 shows the analogous results for a pair of DNACCs coated with two kinds of tethers which has been described previously 25 . We label the tethers on one plate by a and b, and those on the other plate by a and b . Each tether type is grafted at a density of 1/(20 nm) 2 . The sticky ends of a and a hybridize strongly, with free energy ∆G 0 a , forming bridges between the two plates. Conversely, the sticky ends of a and b, as well as those of a and b , hybridize weakly, with free energy ∆G 0 a + 5 k B T , forming loops within each plate. A competition thus arises between strong bridges and weak loops, which has been shown to result in interaction potentials with an unusual temperature dependence that should significantly enhance the ability of these DNACCs Density of strong and weak bonds formed between plates with competing-linkages, at plate separation 20 nm. Strong a-a inter-plate bonds (bridges) have sticky end hybridization free energy ∆G 0 a , and weak a-b and a -b intraplate bonds (loops) have sticky end hybridization free energy ∆G 0 a + 5 kBT . See text for details. Inset: potential of mean force per unit area as a function of plate separation, for ∆G 0 a = −5 kBT (upper), −15 kBT (bottom) and −25 kBT (middle). When the sticky end binding is strongest, the plates barely interact.
to form binary crystals 25 . All the features of these interactions are reproduced with quantitative accuracy by our self-consistent theory.
The key advantage of the present theory is its generality and comparative simplicity. With little additional effort, we can consider systems with many different kinds of grafted DNAs, each with different entropic costs for binding, and with different and potentially nonuniform coating densities. In the previous two cases, the resulting SCTMF theory is equivalent to earlier self-consistent mean-field theories that were less general and that were derived more heuristically (Appendix C). In our third and fourth examples, we go beyond what has been done previously to illustrate the generality of our approach. For our third example, we generalize the setup of Figure 3 to include a third type of tether on each plate, respectively c and c . The sticky ends of c and c hybridize only with each other and with the same strength as the a and a sticky ends. Furthermore, we set the lengths of the a, b, a and b tethers to 30 nm and those of the c and c tethers to 10 nm. Finally, we set the grafting density of the c and c tethers to 1/(30 nm) 2 . Thus, our third system features three different types of rods, of different lengths, with different coverages and with a nontrivial interaction matrix. Figure 4 summarizes the behavior of this system. Our choices result in the interaction energy per unit area having up to two different minima, depending on the specific temperature (binding strength). The example also illustrates one way in which the mean-field approach fails, and how the full theory behaves for such cases. The large discrepancies in the number of c-c bonds formed (black solid and dashed curves in Figure 4), and the resulting discrepancies in the potentials of mean force, result from the low coverage density of these tethers. The typical separation between grafting points that result from Poisson statistics is far larger than the length of the tethers, but these inhomogeneities are not captured by a mean field approach. Importantly, since the full theory retains detailed information on the grafting points of c and c tethers, the MC results are again recovered to quantitative accuracy.
Our fourth and final example in this section is a system of two parallel plates, each with a 100 nm-radius circular patch of dsDNA strands, distributed uniformly within the patch at a coverage density 1/(20 nm) 2 . As above, the strands are modeled as 20 nm rigid rods. The only hybridization allowed is between strands on opposite plates, and we consider solution hybridization energies of −8 k B T , −7 k B T and −6 k B T . Figure 5 shows the potentials of mean force of these plates as a function of plate separation and lateral displacement (zero lateral displacement means that the centers of the circular patches line up vertically). Such detailed potentials could be used to inform parameter choices for coarsegrained models of patchy particles 34 when the patches are made from patterned DNA coatings. For nanoscopic colloids, we suggest that such patches could be realized by building the particles through DNA origami 35 . However, to our knowledge, such precise experimental control of coating patterns for microscopic colloids remains elusive, although recent experiments with triblock Janus particles are very encouraging 36 .

B. Colloids coated with long ssDNA tethers
To describe more flexible ssDNA tethers, we use the representation recently employed by Rogers and Crocker 17 to describe experimental measurements of DNACC interaction potentials. They model ssDNA tethers of contour length 40 nm as 8-segment freely jointed chains of Kuhn length 5 nm. The sticky ends that cap them are modeled as point-like, so that when two tethers bind, they form a 16-segment freely jointed chain. Appendix A 2 describes how to calculate ∆G for two such tethers i and j.
In Rogers and Crocker's experiment, there are two complementary types of sticky ends, A and B, whose sequences are predicted 37 by DINAmelt 38 to have binding enthalpy ∆H 0 = −47.8 kcal/mol and binding entropy ∆S 0 = −139.6 cal/(mol·K). At temperature T , the A-B binding free energy in solution is taken to be ∆G 0 = ∆H 0 − T ∆S 0 . Tethers with these sticky ends are coated on three types of spherical colloids, all of diameter 1.1 µm, A-type colloids are coated with 4800 ± 480 Atype tethers, B-type colloids are coated with 4200 ± 420 B-type tethers, and AB-type colloids are coated with 2400 ± 240 A-type and 3500 ± 350 B-type tethers. In all the following simulations, we have used the central value of all ranges.
We first address the interaction between A and B colloids. The coating densities used by Rogers and Crocker are sufficiently low that local grafting density fluctuations are relevant. We have thus prepared 14 independent configurations of uniformly coated A and B colloids. Figure 6 shows the number of bonds formed between two such colloids when separated by 5 nm to 30 nm versus sticky end binding strength ∆G 0 , as calculated by explicit Monte Carlo simulation and the explicit SCT theory, averaged over the 14 independent coating realizations. The agreement is nearly quantitative, both when averaged over grafting realizations and independently for each such realization (inset).  ure 6 to obtain the mean attraction between the A and B colloids. Further accounting for the entropic repulsion caused by confining the tethers between the two colloids results in a potential of mean force for this system, shown as dots in Figure 7. In this Figure, we also show the result of calculating the colloid interaction potential by calculating the interaction energy per unit area between two similarly coated plates, and then applying the Derjaguin approximation. The interaction free energy between plates is calculated in three different ways: by explicit MC simulations, by using SCT, and by using its mean-field approximation, SCTMF. Our results here help establish the Derjaguin approximation as valid for ssDNA tethers modeled as freely-jointed chains. For completeness, we show similar results for the AB-AB system explored by Rogers and Crocker in Figure 8.
It should be noted that the simple freely-jointed chain model of ssDNA used here and in Ref. 17 may not adequately capture the complications of the experimental system 18,19 . Our aim in this paper is limited to establishing the SCT as an accurate and fast alternative to detailed MC calculations, regardless of the underlying model used to calculate the configurational entropic factors associated with individual tethers binding.

IV. SUMMARY OF RESULTS
For ease of reference, in this section we summarize the main results of the paper. These results are implemented in the Python package DNACC, available at http://github.com/patvarilly/DNACC.

A. Self-consistent theory (SCT)
Our central result describes the interaction details of a DNACC system where each strand i is modeled explicitly. The algorithm is: 1. Using a suitable microscopic model, calculate the binding free energy ∆G ij between every pair of strands i and j. When i equals j, or when i and j cannot bind, set ∆G ij to infinity. This binding free energy is the sum of two contributions: the solution hybridization free energy of the sticky ends, ∆G 0 ij , and the entropic cost of binding of the two tethers that hold the sticky ends, ∆G 2. Solve the self-consistent set of equations for each of the quantities p i , i.e., the probability p i that strand i is unbound: It is straightforward to solve these equations using fixed-point iteration, i.e., making an initial guess for {p i }, then repeatedly inserting the values of {p i } into the right-hand side of the above equations to obtain an improved estimate for {p i }. A solution can also be obtained through minimization of a function with no other local minima or saddle points, as detailed in Section II A.
3. Calculate p ij , the probability p ij that strands i and j are bound: The interaction free energy of a DNACC system can then be calculated through thermodynamic integration of the binding free energies ∆G ij , as follows: 1. Calculate the the free energy of repulsion, F rep , equal to the free energy of the system when none of the strands can bind, i.e., when all the ∆G ij are infinite. This is simply the sum of the free energies of repulsion of the individual strands. We assume this quantity can be calculated analytically or otherwise, as described in Appendix A.
2. Define a thermodynamic integration variable, λ, and denote by p ij (λ) the probability that strands i and j are bound when the all the binding free energies are respectively set to ∆G ij + λ k B T , calculated as above.
3. The free energy of the system is given by the following equation B. Self-consistent mean-field theory (SCMFT) When colloids are coated uniformly with a small number of different DNA species at high density, it is useful to model the colloid interaction in mean field. Assuming that the colloidal surfaces have radii of curvature significantly larger than the range of interaction of two strands, we may first calculate the interaction free energy of two parallel plates as a function of plate separation, h, then apply the Derjaguin approximation.
The details of the interaction between two parallel DNA-coated plates separated by a distance h are estimated as follows: 1. For each pair of DNA species, a and b, calculate a binding factor K ab as follows: where ∆G (cnf) ab is the configurational cost for a fixed a-type ligand to bind a b-type tether grafted at r b . If a DNA species is present in both parallel plates, treat the DNA on each plate as a separate species. Simple analytic expressions for K ab for rod-like and polymeric tethers are given in Eqs. (25) and (26). 2. Solve the self-consistent set of equations for each of the quantities p a , i.e., the probability that an a-type strand is unbound: where σ b is the number of b-type strands per unit area.
3. Calculate the number of a-and b-type bonds per unit are, σ ab , as follows: The free energy per unit area between two DNA-coated parallel plates separated by a distance h is estimated as follows: 1. Calculate the free energy density of repulsion, f rep (h), equal to the free energies of repulsion of the individual strand types, weighted by their coverage density (Eq. (A8)). Expressions for this quantity in the case of rod-like and polymeric ligands are given in Eqs. (A8), (A15) and (A33).
2. Define a thermodynamic integration variable, λ, and denote by σ ab (λ) the number of a-b bonds per unit area when the all the solution binding free energies are respectively set to ∆G 0 ab + λ k B T , calculated as above.
3. The free energy density of the system is given by the following equation Finally, to estimate the free energy of interaction between arbitrary colloids, use the Derjaguin approximation 32 . In the case of two spheres of radius R 1 and R 2 , respectively, with a distance of closest approach h, this approximation yields: In using the Derjaguin approximation, it is important to shift all free energies appropriately to ensure that f (h) is zero at infinite h.

V. CONCLUSIONS AND OUTLOOK
The theory for DNACC interactions that we have presented, both in its explicit form (Eqs. (6), (7) and (15)) and its mean-field approximation (Eqs. (21)-(24)) is capable of describing many different DNACC setups, varying in grafting density, sticky end interactions and tether details. The theory is robust and correctly incorporates the statistical mechanics of valence-limited binding, which plays a key role in the behavior of DNACCs, and it does so at much lower cost than full MC simulations. Extending the theory to treat general ligands within our approximation is straightforward, as exemplified in Appendix A for freely-jointed-chains.
When comparing our full theory to MC data in the context of DNACCs, we have found quantitative agreement. The mean-field approximation, which substantially generalizes earlier approaches, is in quantitative to semi-quantitative agreement with the exact results. Significant discrepancies arise only when the grafting density fluctuations that are neglected in the mean-field approach become relevant.
Our theory accounts for all the complications arising from valence limitation and combinatorial entropy. Nevertheless, the binding free energy ∆G ij of a pair of ligands i and j is taken as an external input, and the agreement of the theory's predictions with experimental results will depend significantly on the quality of the estimates for the values of ∆G ij . Simple models for generating these estimates, such as those used in this paper, can inform experiments, though apparently, they are not yet sufficiently accurate for quantitative agreement with experiments, as we have discussed in Ref. 18.
We expect our theory to be a significant aid in the search for viable approaches to design crystal structures a priori. Presently, existing free energy methods can establish the stable crystal structure of a system given its interaction potential 39 . Conversely, robust methods are being developed for solving the inverse problem, determining the interaction potentials that might yield a particular target structure 40 . Designing colloids that realize such interaction potentials experimentally is another problem. In this regard, DNACCs, with their highly tunable interactions and vast parameter space (e.g., particle size, ligand structure and flexibility, coating density, nucleotide sequence length and hybridization strength), are a promising tool to build colloids with programmable interactions. For this reason, our theory can complement existing tools in crystal structure design by quickly suggesting experimentally realizable constructs that result in particular colloid-colloid interaction potentials. Whereas for simple potentials one can envision a manual search of parameter space guided by intuition and experience, it might also be possible to recast the problem in terms of a global minimization algorithm that optimizes the DNACC structures to obtain interaction potentials that are close to some target. It is precisely to help others easily use our tools in these and other ways that we have written a comprehensive and well documented series of Python modules that implement our theory, available under the GNU license at http://github.com/patvarilly/DNACC.
While we have developed our theory and the accompanying programs with an eye towards applying them to DNA-coated colloids, other interesting applications are possible. In particular, our theory can help model any system where attraction is dominated by linker-receptortype bonding. These include patchy particles, where the patches are realized by a nonuniform ligand coat. More generally, many biological systems could be treated with our theory, for example, interactions between cells (or between cells and a substrate) that are mediated by membrane proteins. For these more general cases, care should be taken to include the effects of surface elasticity and how the entropic costs of ligand binding are, in turn, modified.

ACKNOWLEDGMENTS
We thank M. Leunissen and F. Martinez-Veracoechea for useful discussions and a critical reading of the manuscript. We also thank G. Hummer for pointing out the analogy between Eqs. (6) and (7)  In this appendix we describe how to calculate the freeenergy change ∆G ij of binding two ligands, i and j, when the binding interaction is physically limited to a small, reactive sticky end. We begin with a general treatment, and later specialize to the two important cases of short, stiff polymers and of freely-jointed chains. These two models have been used to describe the DNA strands on DNA-coated colloids when these are short, rigid dsDNA 12,15,16,25 and when they are long, flexible ssDNA 17,18 . These models may also be useful to describe ligand-mediated binding in important biological systems, such as cells, where adhesion is, in many cases, mediated by different glycoproteins coated on their surface membrane 41 .
Let Q i , Q j and Q ij be the partition functions of the two unbound ligands and the bound complex that they form. In terms of these quantities, ∆G ij is given by As in the main text, we have assumed that the unbound ligands do not otherwise interact with each other. The Boltzmann factor exp(−β∆G ij ) is generally smaller than that of the sticky ends in solution at standard conditions, exp(−β∆G 0 ij ), owing to the entropic penalty of constraining the configurations available to ligands i and j when their reactive sticky ends are bound. This defines the configurational free energy cost of binding, ∆G (cnf) ij , as follows: A realistic model that reproduces ∆G ij for ligands such as DNA strands is expensive, as one would need to at least model each nucleotide of each strand explicitly 42 . It is computationally more convenient to separate the binding of the sticky ends from the configurational freedom of the rest of the ligand, and use a coarse-grained representation of the sticky ends where ∆G 0 ij is introduced as an explicit external parameter. Concretely, the internal partition functions of the free and bound sticky ends in solution, Z i , Z j and Z ij , are related to the equilibrium constant of binding 14 , where ρ 0 is the standard concentration. Using this relation in Eq. (A1) to factorize the unknown dependence of Q ij , Q i and Q j on the internal partition functions, we obtain are approximately independent of the chemical identity of the sticky end ("cnf" stands for "configurational"). It follows that We stress that this factorization is an approximation that holds only when the binding of the sticky ends is only weakly affected by the presence of the rest of the ligand, as occurs when the ligand is an ideal polymer or when the binding is very strong. However, it is a useful approximation because once ∆G (cnf) ij has been measured or calculated, one can predict the values of ∆G ij for families of ligands that differ only in their sticky ends, providing us with a powerful design tool.
We now specialize to the case of polymeric ligands with one end fixed and the other end sticky, such as polymers grafted at a surface. The main complication in treating polymers is dealing with excluded volumes. For instance, in DNACCs, the DNA strands cannot penetrate the colloids. In the discussion below, we denote by r i the grafting point of polymer i. Given two polymers i and j, we denote by r ij the vector r j − r i with length r ij . We denote by Ω the region of space that can be explored by the monomers of the polymers. Polymers restricted to Ω are deemed "confined", whereas those for which this restriction is absent are "unconfined".
We note that the free energy of repulsion due to a polymer i, denoted by F rep,i , is captured by the dependence of Q (A6) From this expression follows the total free energy of repulsion, Similarly, the mean-field free energy density of repulsion is given by where σ a is the coverage density of species a and F rep,a is the free energy of repulsion of any representative a-type strand.

Short, stiff polymers as ligands
A short, stiff polymer, such as a short piece of dsDNA, can be modeled as a rigid rod, as has been done previously in Refs. 15, 16, and 25. We gather the results here for completeness and to clarify the origin of a Jacobian prefactor absent in Ref. 15 but addressed in Ref. 16. For rigid rods, the configurational partition functions are: where L i and L j are the lengths of rods i and j. The first integral is most easily evaluated in cylindrical coordinates, with the z axis parallel to r ij . Special care must be taken to isolate the Jacobian factor that arises when dealing with the product of the two δ-functions. After evaluating the integrals, we obtain, Here, Ω ij is the portion within Ω of the circumference of the circle explored by the common endpoint of the hybridized rods, Ω i and Ω j are the areas within Ω of the sphere that the endpoints of the free rods can explore, and γ i and γ j are the angles that the hybridized rods make with r ij . The factor of 1/ sin(γ i + γ j ) was obtained in Ref. 16 using two different representations of the δfunction, but the procedure here highlights the generality of the result. When the colloids are plates, one at z = 0 and one at z = h, the region Ω is simply that where 0 ≤ z ≤ h, and one can write an explicit expression for Ω ij . For the case where the rods are grafted on opposite plates, such an expression is given in Appendix A of Ref. 16 (see also Figure 8 therein). When the rods are grafted to the same plate, Ω ij is zero unless |L i − L j | ≤ r ij ≤ L i + L j , in which case The quantity Ω i also has a compact expression, The free energy of repulsion of a strand i follows immediately from Eq. (A6), 2. Long, flexible polymers as ligands Following Refs. 17 and 18, we can model a polymer i as a freely-jointed chain (FJC) with N i segments of length . Let r i,s be the endpoint of segment s, with 1 ≤ s ≤ N . For convenience, let r i,0 = r i , and let u i,s = r i,s − r i,s−1 . The chain segments are confined to the region Ω, but do not otherwise interact. It follows that where, abusing notation somewhat, the integration domain confines the polymer to the region Ω. A completely analogous treatment yields Q (cnf) j for a FJC with N j segments of length , grafted at r j .
One way of modeling the binding of i and j is to coalesce the endpoints of the two chains, so that r i,Ni = r j,Nj , resulting in a FJC with endpoints fixed at r i and r j and with N hyb = N i + N j segments {u s }, where 1 ≤ s ≤ N hyb . However, if the persistence length of the polymers is short, then the binding site may have a size comparable to the segment length . This is the case for ssDNA polymers with sticky ends that are 6 to 12 nucleotides long. Hence, we also consider a more general model for binding of i and j in which the final segment of the two polymers coalesces, so that r i,Ni = r j,Nj −1 and r i,Ni−1 = r j,Nj , resulting in a FJC of N hyb = N i + N j − 1 segments 43 . The partition function of the bound system is whereQ (cnf) ij is the partition function of a confined FJC with fixed end points, r i and r j . The additional factor in Q (cnf) ij results from fusing the end segments of the individual polymers.
Although the integrals in Eqs. (A16) and (A18) cannot be evaluated analytically even in simple geometries, as was done in the previous section, one can calculate their value numerically using standard polymer techniques. Let Q 1 (N ) = (4π ) N be the partition function of an unconfined N -segment FJC fixed at one endpoint, and let Q 2 (N, R) be the partition function of the same FJC with fixed end-to-end distance R. Their ratio yields the probability density for the end-to-end distance of the FJC to be R, which is known in closed form 44,45 . Using these partition functions, we can writeQ (cnf) ij in the following convenient form:Q The ratio on the right-hand side of Eq. (A20) is the probability that a random walk from r i to r j lies wholly within Ω (i.e., does not penetrate the colloid). To compute it, we use Rosenbluth sampling 39 , a common technique to generate a properly weighted ensemble of such FJCs. Specifically, to grow each segment s of one chain that starts at r i , we first generate K possible segments u (k) s+1 (1 ≤ k ≤ K), sampled with probability density p(N hyb − (s + 1), |r j − r s+1 |) using rejection sampling 46 . We then pick at random one of the m s segments in {u (k) s+1 } wholly within Ω. The Rosenbluth weight of that chain is defined as The fraction of FJCs that lie wholly within Ω is then the average Rosenbluth factor, where the subscript implies that the average is computed as described above. We can estimate Q (cnf) i and Q (cnf) j analogously, except that at each step in growing the chain, now of length N i or N j , the possible segments {u (k) s+1 } are chosen uniformly on the surface of a sphere of radius , since the end of the chain is free. We obtain Substituting Eqs. (A17)-(A24) into Eq. (A5), we obtain the following expression for the configurational entropy cost of binding i and j: In Figure 9 we show the configurational entropy cost of hybridizing two FJCs (made of 8 segments with = 5 nm) tethered to parallel plates. Interestingly, if ∆G (cnf) ij is expressed in terms of the distance between tethering points, r ij , and the plate separation, h, the dependence on the latter is marginal. This suggests that, for ideal chains, the cost of stretching a chain (p(N hyb , r ij ) in Eq. (A25)) is more important than the excluded volume terms (W ). Hence, we may approximate ∆G (cnf) ij analytically by neglecting all chain degrees of freedom but the sticky end. In particular, using a Gaussian approximation for p(N, r) 45 , we can model the ligand as a reactive point-like sticky end that is confined between the two walls, and is attached to the tethering point via a harmonic spring with spring constant k = 3 k B T /N 2 . The partition function of an unhybridized polymer is then given by where the integration domain, Ω, is [−∞ : ∞] 2 × [0 : h], and the error function, erf(x), is defined as in Eq. (27). Similarly, the partition function of the hybridized polymer is given When i and j are tethered to opposite plates, so that the bound pair forms a bridge, we obtain whereas when i and j are tethered to the same plate, so that the bound pair forms a loop, we obtain (A30) From these formulas, using Eq. (A5), we obtain analytical estimates of ∆G (cnf) ij . In Figure 9, these estimates are shown in the solid curve. We see an overall agreement with the exact result up to r ij = 50 nm, where the finite extensibility of the chain becomes important.
The above formulas can be used to derive an approximate expression for the coefficients K ab used in the meanfield treatment of Section II C (Eq. (26)). Indeed, for bridges, and for loops, In both of these cases, we have assumed that Q (cnf) ab goes to zero sufficiently quickly for the integrals to be insensitive to their upper bounds.
Finally, the above formulas also yield simple approximate expressions for the free energy of repulsion, F rep,i , of a strand i. From Eq. (A6), we obtain ligands (DNA strands) bind to one another. We can take a complementary "sticky-end-centric" approach to obtain the same results, which sheds light on the recent Local Chemical Equilibrium treatment of Rogers and Crocker 17 and corrects its main deficiency. As in the main text, we first treat all strands separately, then make continuum approximations for simplified geometries. Following the notation of Ref. 17, let C 0 i (r) be the normalized probability density of finding the sticky end of tether i at r, before any hybridization occurs. For this treatment only, we assume that the strands interact "ideally", i.e., with an interaction energy ∆G 0 ij when their sticky ends coincide and none otherwise, so that the configurations of the tether of i with sticky end at r have the same statistical weight whether the tethers are bound or unbound. In terms of partition functions, strands i and j interact ideally if their partition function when bound, Q ij , is given by where Q i (r) is the partition function of unbound strand i with sticky end at r, and Q j (r) is defined similarly. The restriction of ideal binding, absent in the general treatment of the main text, is satisfied when the tethers are ideal chains, for example. With this additional assumption, we rewrite Eq. (1) as follows: where ρ 0 is the standard concentration, 1 M. The analog of Eq. (2) is: As in the main text, let p i = Z −i /Z, and approximate Z −i,−j /Z by p i p j . As in Ref. 17, let C i (r) be the concentration of free sticky end i at r, equal to p i C 0 i (r) by comparison to the equation above. There being only one such sticky end, C i (r) is also the probability density that the sticky end of tether i is free and is at r. We obtain (B4) More concisely, let C ij (r) be the concentration of bound sticky ends i and j, given by We can then finally write the self-consistent condition for the quantities {p i } as follows: Comparing the above condition to the equations in Ref. 17, we conclude that the "local chemical equilibrium" approximation introduced by the authors is equivalent to our approximation Z −i,−j /Z ≈ p i p j , but that the self-consistent condition is significantly different from the condition that they proposed, Physically, in the treatment of Rogers and Crocker, depletion of free sticky ends i at a point r establishes a density gradient in C i (r) beyond that implied by the polymer statistics of tether i, but this density gradient is not allowed to relax, effectively imposing a nonequilibrium condition. Our treatment, summarized in Eq. (B6), corrects this shortcoming. As in the main text, for those cases where there is translational invariance, we may approximately treat together tethers with similar binding behavior. In the concrete case of two plates perpendicular to the z axis, and in the absence of binding, the concentration of a-type sticky ends at a height z is: where P a (z) = dx dy C 0 i (x, y, z), and i is any representative tether of type a. The corresponding mean field self-consistent equations are Rearranging these equations to obtain explicit equations for {p a } yields Comparing this last equation to Eq. (21), we find that This is a specialization of Eq. (22) under the additional assumption of ideal binding interactions. For the concrete case of short, rigid-rod-like tethers of length L, the probabilities P a (z) are particularly simple, and depend only on rod length, plate separation, and the plate on which a-type tethers are grafted. In particular, if they are grafted on the plate at z = 0 and the other plate is at z = h, then given in Appendix A of this paper and Appendix A of Ref. 16. We have verified that both routes give identical results when evaluated numerically for L b = L a and L b = 1.7L a .
In the absence of translational invariance or some similar symmetry, it is nontrivial, in general, to write down self-consistency conditions similar to Eq. (B6) where liketype tethers are grouped together. Unlike the incorrect condition of Eq. (B7), the correct self-consistency conditions are inherently non-local from the point of view of sticky ends: the change upon binding in sticky end concentration of a-type sticky ends at r depends on the precise identity of the a-type tether to which that sticky end is attached. Indeed, one needs to retain the level of detail inherent in Eq. (B6). As such, there is no advantage in adopting a "sticky-end-centric" point of view as opposed to a "tether-centric" point of view when dealing with arbitrarily shaped colloids, polymer statistics or non-uniform coverage. Moreover, the tether-centric equations are simpler, since they result from integrating over the irrelevant degrees of freedom describing the position of the sticky ends, and are more general, since they do not assume that binding is ideal.
To conclude this section, we present a numerical comparison between Monte Carlo results, our self-consistent theory and the LCE theory that illustrates the extent to which these approximate theories deviate from an exact treatment. We consider two periodic parallel plates, each with 500 tethers, at a density comparable to that used in Ref. 17 (4500/(4π × 550 2 ) = 0.001183 tethers/nm 2 ). The tethers are modeled as ideal 8-segment chains of Kuhn length 5 nm that, when bound, form a 16-segment ideal chain, as described in Appendix A. Figure 10 shows the number of bonds formed in this system as a function of the solution hybridization free energy, ∆G 0 , of the sticky  ends of tethers on opposite plates, when the plates are separated by 15 nm and by 25 nm. As is observed, our self-consistent theory closely tracks the Monte Carlo results, whereas the estimates obtained using the LCE approach of Rogers and Crocker differ appreciably from the exact ones.
Appendix C: Self-consistent equations from a Saddle Point Approximation In this appendix, we present an alternate derivation for the self-consistent mean field theory in Eqs. (21) and (22) in two specific setups, based on a saddle point approximation of a mean-field partition function. The two geometries we explore are the first two examples in Section III, which we refer to as the single linkage system and the competing linkages system.

Single linkage system
In the single linkage case, two parallel plates of area A are respectively decorated with N a and N b strands of two complementary kinds of DNA, a and b, whose surface concentration is thus σ a = N a /A and σ b = N b /A. Without loss of generality, we assume that N a ≤ N b . Let n b ← a be the average number of b-type strands that a fixed a-type strand can bind to. When no bonds have yet formed n b ← a = σ b A b ← a , where A b ← a is the area that encloses all possible grafting points of b-type binding partners of a fixed a-type strand. However, after having formed n ab bonds between a-and b-type strands the average number of available linkers is reduced to 16 The average Boltzmann factor for binding a-and b-type strands is K ab /A b ← a , where K ab is defined in Eq. (22). With this results and the previous one, we can write a mean-field partition function of the system as follows Z = Na n ab =0 N a n ab n ab −1 j=0 K ab A b ← a n b ← a (j) , = Na n ab =0 N a n ab We now isolate the dependence of Z on A in order to take the limit A → ∞, which we then evaluate through a saddle-point approximation. Let s ab = n ab /A, and approximate the sum in Eq. (C2) by an integral with respect to s ab . Assuming the maximum of the integrand is far from the integration limits, we obtain where, using Stirling's approximation, S(s) = σ a ln σ a + σ b ln σ b − s ln s − s − (σ a − s) ln(σ a − s) − (σ b − s) ln(σ b − s) + s ln K ab . (C4) In the limit A → ∞, the partition function will be dominated by the values of s that maximize S(s). Setting dS/ds to 0 at s = σ ab , we obtain the following equation for the average density of hybridized strands, σ ab : which restates Eq. (23) after identifying σ a p a with σ a − σ ab and likewise for σ b p b . Using Eq. (21) one can also rewrite Eq. (C5) as which is the expression used in Ref. 16.

Competing linkages system
In this section we generalize the saddle point approach to the competing linkages system. Like in Section III, we have two parallel plates of area A, one grafted with aand b-type strands, the other with a -and b -type strands. The only allowed interactions are a-a , a-b and a -b . An equivalent treatment applies when, as in Ref. 16, allowed interactions are a-a , a-b and a -b, so that loops (intraparticle bonds) are forbidden. The only differences between the two systems are the labels in the equations and the values of the K ab factors below.
Using the same notation as in the previous section, we find that the average number of binding partners available to a fixed strand of each type is n a ← a (n aa , n a b ) = N a − n aa − n a b N a σ a A a ← a , Assume for simplicity that N a ≤ N a , N a ≤ N b and N a ≤ N b . A nearly identical argument applies for a different relative populations of strands. A mean-field partition function for the system is N a n aa n aa −1 k=0 K aa A a ← a n a ← a (k, 0) · N a − n aa n ab n ab −1 As in the previous section, we isolate the dependence of Z on A to take the A → ∞ limit. Let s xy = n xy /A, and approximate the sums in Eq. (C8) by integrals with respect to s aa , s ab and s a b . Assuming the maximum of the integrand is far from the integration limits, we obtain where the first, constant term does not depend on s xy .
The point (σ aa , σ ab , σ a b ) at which S is maximized, found by setting all its partial derivatives to zero, yields the average bond densities for each pair of species. They satisfy σ aa = (σ a − σ aa − σ ab )K aa (σ a − σ aa − σ a b ) = σ a p a K aa p a σ a , (C11a) These equations correspond to the mean-field selfconsistent theory of Eqs. (21) and (23)) for the competing linkages model.
We now show that the saddle points equations (C11) are also equivalent to the expressions derived in Ref. 16 (once σ ab and σ a b are replaced with σ ab and σ a b respectively). Using Eq. (C11c) two times, we have .
Substituting this result into Eq. (C11a), together with the definition of p a in Eq. (21), we find that = σ a K aa (σ a − σ aa ) · 1 + K a b (σ b − σ a b ) + K aa (σ a − σ aa ) + K ab (σ b − σ ab ) This is exactly the expression used in Ref. 16 to compute the density of a-a linkages, derived there by applying an heuristic self-consistent prescription to an approximate estimate of σ aa . A similar argument recovers the expressions for σ ab and σ a b used in Ref. 16.