A simple analytical formula for the free energy of ligand–receptor-mediated interactions Stefano Angioletti-Uberti,1 Patrick Varilly,1 Bortolo M. Mognetti,1 Alexei V. Tkachenko,2 and Daan Frenkel1 1)Department of Chemistry, University of Cambridge, Lensfield Road, CB2 1EW Cambridge, UK 2)Center for Functional Nanomaterials, Brookhaven National Laboratory, Upton, New York 11973, USA Recently1, we presented a general theory for calculating the strength and properties of colloidal interactions mediatedbyligand-receptorbonds(suchasthosethatbindDNA-coatedcolloids). Inthiscommunication,we derive a surprisingly simple analytical form for the interaction free energy, which was previously obtainable only via a costly numerical thermodynamic integration. As a result, the computational effort to obtain potentials of interaction is significantly reduced. Moreover, we can gain insight from this analytic expression 3 for the free energy in limiting cases. In particular, the connection of our general theory to other previous 1 specialised approaches is now made transparent. This important simplification will significantly broaden the 0 scope of our theory. 2 n a We consider a general system of many linkers, such J as a solution of colloids coated with DNA strands that Nij=Nim=Njl=Nlm=0 Nij=Nim=Njl=Nlm=0 Nil=Njm=1 Nil=Njm=1 8 are capped with reactive sticky-ends. At any given time, eachlinkericanbindatmostoneotherdistinctlinkerj, N�i�j=1 N�i�j=0 ] l m l m withafreeenergychange∆G thatdependsonthepoly- t ij f mer statistics of the linkers (e.g., length, flexibility and o s grafting position). In many cases, including those of ex- t. perimental relevance, the probability that linker i is un- es a bound is approximately independent of whether or not pi i j i j m o any other linker is also unbound. Here, we show that in C l m l m - this limit, the free-energy of interaction of the system is N d given by n o (cid:88) (cid:88) c βF = lnp + p (1) att i ij i j i j [ i i<j 3 Ensemble realisations v where pi is the probability that linker i is unbound and 3 p is the probability that linkers i and j form a bond. ij 7 Previously1, we showed that these quantities are given FIG. 1. Two different realisations of an ensemble of two 8 by the unique physical solution to the following set of copies. Linkers are depicted as straight lines, and bonds are 1 self-consistent equations: shownasfilledcircles. Althoughthenumbersofbondsformed . in each ensemble, {N }, are equal, the number of copies 1 ij 1 pij =pipje−β∆Gij, (2) where both i and j are unbound, N−i−j, differs. 2 (cid:88) p =1 p . (3) 1 i ij − : j I. DERIVATION OF THE MAIN RESULT v i In what follows, we first motivate the free energy ex- X WeconsiderhereanensembleofN independentcopies pression in Eq. (1) through a calculation that closely r of the real system. In each copy, a different set of bonds resembles that for mixing entropy of solutions and a forms between the linkers (see Figure 1). Let N be the gases. This free energy is minimised for the values of i number of copies where linker i is unbound, and let N p and p that solve the self-consistent conditions ij i ij { } { } bethenumberofcopieswhereiandj areboundtoeach in Eqs. (2) and (3). We then show that the free en- other. Conversely, let N be the number of copies ergy in Eq. (1) is identical to that obtained through the i, j wherebothiandjareunb−ou−nd. Thesequantitiesarenot costly and numerical thermodynamic integration previ- independent: eachlinkeriiseitherunboundorunbound, ously proposed. In the discussion section, we compare so the performance of Eq. (1) with that of the thermody- namic integration, establish the explicit connection be- (cid:88) N + N =N. (4) i ij tween our Eq. (1) and a previous treatment of DNA- mediated colloid interactions2, and state the analogous j result to Eq. (1) for the mean-field system of plates dis- The fraction of copies where i is unbound (f ), where i cussed in Ref. 1. i and j are bound to each other (f ), or where they are ij 2 both unbound (f ) follow immediately: i, j − − Z( ,Nij, )N i, je−β∆Gij ··· ··· − − =Z( ,N +1, )(N +1). (8) fi =Ni/N, (5) ··· ij ··· ij fij =Nij/N, (6) The value of N i, j depends not just on the values of f =N /N. (7) N but on t−he−details of how those bonds are dis- −i,−j −i,−j {tribiju}ted between copies (see Figure 1). To remove this complication, we approximate the probability of j be- ing unbound as independent of whether or not i is also Let Z( N ) be the partition function of an ensemble { ij} bound. Hence, under the constraint that each pair of linkers i and j is boundtoeachotherinexactlyNij copies. Aclosed-form NiNj N =Nf Nf f = . (9) expression for Z( N ) can be constructed recursively, i, j i, j i j { ij} − − − − ≈ N by adding each bond one by one. For a given set of Nij , we need to work out how Z( Nij ) changes upon Neatly, this approximation allows us to treat N i, j as a{ddin}g one more i-j bond, which w{e do}as follows. We a function of only Nij . From the discussion ab−ov−e, we { } call the set of realisations of the ensemble with {Nij} obtain an expression for the increase in Z({Nij}) upon bondstheold ensemble,andthatofrealisationswithone adding one i-j bond to the system: more i-j bond, the new ensemble. In the old ensemble, there are N i, j copies where an i-j bond can be added. Z(··· ,Nij +1,···) e−β∆GijNiNj In doing th−is,−all the realisations of the new ensemble Z( ,Nij, ) ≈ N(Nij +1) ··· ··· areraeligseanteiornatse,do,nbeuwtitnhoatnuniiqjuebloyn.dFoinrceoxpaymXplbe,ugtinvoetnitnwYo =e−β∆Gij(N −(cid:80)kNik)(N −(cid:80)kNjk). orZ,andonewithani j−bondinYbutnotinXorZ,the N(Nij +1) − (10) same final realisation can be obtained by adding an i j − Thisrecursionrelation, andthefactthatZ =1whenno bond to Y in the former and X in the latter. Conversely, bonds form, allows us to write an approximate closed- inthenewensemble, wecangeneratearealisationofthe form expression for Z( N ), namely old ensemble in Nij +1 ways by removing one of the i-j { ij} bonds. For example, the old realisation with an i-j bond iincjopbyonXdbfruotmnoatnYeworrZeacliasnatbioenowbtitahinaendibydjebleotnindginthXe Z({Nij})≈(cid:89)i (N −(cid:80)N!kNik)!·N(cid:80)i1<jNij·i(cid:89)<j e−βNNiijj∆!Gij. − − andYbutnotZ,orfromonewithani j bondinXand (11) − ZbutnotY.Sincethenumberofwaysofgoingfromthe Using Stirling’s approximation and Eq. (11), the free en- old to the new ensemble is equal to the number of ways ergy per copy βf = (1/N)lnZ( N ) is then given a∗tt − { ij} ofgoingfromthenewtotheoldensemble,itfollowsthat by (cid:88) (cid:88)(cid:0) (cid:88) (cid:1) (cid:0) (cid:88) (cid:1) (cid:88) (cid:88) βf ( f )= f β∆G + 1 f ln 1 f + f lnf + f . (12) a∗tt { ij} ij ij − ij − ik ij ij ij i<j i j k i<j i<j Treating f as continuousin therange [0,1], theover- WhenN ,thevalues f arepreciselytheaverage ij ij { } →∞ { } all free energy per copy of the ensemble, Fa∗tt, follows values of fij . Eq. (14) implies that fij and pij { } { } { } from a saddle-point approximation: obey identical equations (Eq. (2)), and so are equal: βFa∗tt ≡−N1 ln(cid:90) (cid:32)(cid:89)Ndfij(cid:33)e−Nβfa∗tt({fij}) fij =pij. (15) i<j A. Connection to thermodynamic integration =βf ( f )+ (lnN/N) βf ( f ), (13) a∗tt { ij} O ≈ a∗tt { ij} where the(cid:80)integration is over all positive values of {fij} The free energy Fa∗tt, defined in Eq. (13), is equal to satisfying f 1 for all i and the values f are the free energy of the real system to the extent that the j ij ≤ { ij} obtained by minimising the free energy per copy, approximation in Eq. (9) is valid. Since this is the same approximation that we used previously1 to calculate the (cid:12) ∂βfa∗tt(cid:12)(cid:12) =0. (14) free energy in terms of a thermodynamic integral, Fatt, ∂fij (cid:12){fij} it is reasonable to suppose that Fa∗tt and Fatt are equal. 3 We now show this explicitly. Eqs. (12) and (13), we find that tivIenfroeuereonreigrginyaflopratpheer,rewael scyasltceumlatoefdlitnhkeeresxuascitnagttthraecr-- dβFa∗tt =(cid:88)∂βfa∗tt(cid:12)(cid:12) dpij + ∂βfa∗tt(cid:12)(cid:12) . (17) modynamic integration. Specifically, we replaced β∆Gij dλ i<j ∂fij (cid:12){pij} dλ ∂λ (cid:12){pij} byβ∆G +λ,whereupontheprobabilities p and p ij i ij becomefunctionsofλ. Wethenintegrated{the}appr{opri}- ThefirsttermvanishesowingtoEqs. (14)and(15), and atefreeenergyderivativeovertherange0 λ< ,and the second term follows immediately from Eq. (12). We obtained ≤ ∞ then have (cid:90) ∞ dβFatt (cid:90) ∞ (cid:88) dβdFλa∗tt =(cid:88)pij = dβdFλatt. (18) βFatt = dλ = dλ pij(λ). (16) i<j − dλ − 0 0 i<j Moreover, both F and F are zero when λ is infinite, a∗tt att so the two quantities are equal for all λ. The same replacement of β∆G with β∆G +λ can be The previous result, together with Eqs. (12)–(15), ij ij made in the ensemble of N copies. In that case, using yields the following closed-form expression for F : att (cid:88) (cid:88)(cid:0) (cid:88) (cid:1) (cid:0) (cid:88) (cid:1) (cid:88) (cid:88) βF = p β∆G + 1 p ln 1 p + p lnp + p . (19) att ij ij ij ik ij ij ij − − i<j i j k i<j i<j This expression is, in fact, equivalent to the much more compact Eq. (1). Concretely, (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) βFatt = pijβ∆Gij + (1 pij)lnpi+ pijln(pipje−β∆Gij)+ pij, (20a) − i<j i j i<j i<j (cid:88) (cid:88) 1(cid:88) 1(cid:88) (cid:88) = (1 p )lnp + p lnp + p lnp + p , (20b) ij i ij i ij j ij − 2 2 i j i,j i,j i<j (cid:88) (cid:88) (cid:88) (cid:88) = (1 p )lnp + p lnp + p , (20c) ij i ij i ij − i j i,j i<j (cid:88) (cid:88) = lnp + p . (20d) i ij i i<j II. DISCUSSION proximation” or the “weak binding regime”. However, in experiments with micron-sized DNA-coated colloids, Figure 2 reports the typical computational speedups this approximation can be significantly inaccurate6. At obtained by using Eq. (1) versus our original thermody- the nanoscale, where high bond strengths are commonly namic integral, Eq. (16), for systems of M linkers. The used, the Poisson approximation is expected to break speedup is higher for larger M and for stronger bonds down(seeFig.3forcomparison). Eq.(1)isinsteadquan- becauseeachevaluationofthethermodynamicintegrand titatively accurate for bonds of any strength1. involves solving a system of M equations, and the size Eqs. (1), (2)and(3)alsomakeexplicittheconnection of the integration domain scales linearly with the bond between our theory and previous treatments. For exam- strength. Typically,experimentallyrelevantregimesdeal ple, Dreyfus et al.2 model the attraction between two with hundreds to tens of thousands of strands, a regime DNA-coated spheres by first estimating the maximum that can now be treated exactly with Eq. (1). number Np of linkers on each sphere that could form a In the limit of weak-bonds, where p is close to 0, we bondwithalinkeronthesecondsphere,andthenassum- ij find that ing that each such linker can independently bind any of k linkers with an average free energy ∆F . For later (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) tether βFatt = lnpi+ pij = ln(1 pij)+ pij, convenience, we define a small expansion parameter x as − i i<j i j i<j (cid:88)pij +(cid:88)pij = (cid:88)pij. x≡ke−β∆Ftether, (21) ≈− − i,j i<j i<j andnotethattheexperimentsinRef.2takeplacemostly This approximate result has been widely used by pre- inwhattheauthorcallthe“weak-bindingregime”,where vious authors2–5 under the name of the “Poisson ap- x 1 and N x 1. In their model, the partition p (cid:28) (cid:29) 4 25 0 25linkers Presentmodel 44linkers Poissonapproximation 20 100linkers -0.5 200linkers up 15 Np -1 Speed 10 βF/att -1.5 5 -2 -2.5 -14 -12 -10 -8 -6 -4 -2 0 1 2 3 4 5 6 7 8 9 10 Bondstrength(kBT) x=kexp(−β∆Ftether) FIG. 2. Computational speedups obtained by using Eq. (1) FIG. 3. Free energy per linker in the “Poisson approx- vs. Eq.(16), forasystemofparallelplatescoatedwithcom- imation” and the present model. Higher values of x = plementary linkers. kexp(−β∆F ) lead to higher bonding probabilities, ei- tether therbecausebondsarestrongerorbecauselinkershavemore binding partners. The two models agree in the “weak bind- ing regime” (x(cid:28)1), but disagree when correlations between function of the system is neighbouring strands become significant. Z (1+x)Np, (22) ≈ from which the free energy F follows,7 III. ACKNOWLEDGMENTS att βF ln[(1+x)Np]= N ln(1+x) (23) ThisworkwassupportedbytheERCAdvancedGrant att p ≈− − 227758, the Wolfson Merit Award 2007/R3 of the Royal N (x x2/2+ ). (24) ≈− p − ··· Society of London and the EPSRC Programme Grant EP/I001352/1. P.V. has been supported by the Marie In the present framework, which does not treat the link- Curie International Incoming Fellowship of the Euro- ers as binding independently, every linker has the same pean Community’s FP7 PIIF-GA-2011-300045, and by probability p of being bound, given by the solution to a Tizard JRF from Churchill College. Research carried out in part at the Center for Functional Nanomaterials, 1 √1+4x 1 p= p= − . (25) Brookhaven National Laboratory, which is supported by 1+xp ⇒ 2x the U.S. Department of Energy, Office of Basic Energy Sciences, under Contract No. DE-AC02-98CH10886. The free energy then follows from Eqs. (1),(2): 1P.Varilly,S.Angioletti-Uberti,B.M.Mognetti, andD.Frenkel, βF =N p2x+2N lnp (26) “A general theory of DNA-mediated and other valence-limited att p p colloidalinteractions.”J.Chem.Phys.137,094108(2012). Np(x x2+ ) (27) 2R.Dreyfus, M.Leunissen, R.Sha, A.Tkachenko, N.C.Seeman, ≈− − ··· D.J.Pine, andP.M.Chaikin,“Aggregation-disaggregationtran- Thus, our theory recovers the results of Dreyfus et al. in sition of DNA-coated colloids: Experiments and theory,” Phys. Rev.E81,41404(2010). the weak binding regime. However, there is significant 3P.L.Biancaniello,A.J.Kim, andJ.C.Crocker,“Colloidalinter- disagreement already at second order in x, where linkers actions and self-assembly using DNA hybridization,” Phys. Rev. begin to compete for binding partners. Lett.94,58302(2005). Finally,usingthesameprocedureasinouroriginalpa- 4N. A. Licata and A. Tkachenko, “Statistical mechanics of DNA- mediatedcolloidalaggregation,”Phys.Rev.E74,41408(2006). per, we can directly write an attractive free energy den- 5W.B.RogersandJ.C.Crocker,“DirectmeasurementsofDNA- sity for a pair of plates, treated at a more approximate, mediatedcolloidalinteractionsandtheirquantitativemodeling,” spatial mean-field level. In the notation of Ref. 1, P.Natl.Acad.Sci.U.S.A.108,15687–15692(2011). 6B. M. Mognetti, P. Varilly, S. Angioletti-Uberti, F. J. Martinez- 1(cid:88) (cid:88) Veracoechea, J. Dobnikar, M. Leunissen, and D. Frenkel, “Pre- βf = σ p K p σ + σ lnp . (28) att α α αβ β β α α dictingDNA-mediatedcolloidalpairinteractions,”P.Natl.Acad. 2 α,β α Sci.U.S.A.109,E378–E379(2012). 7InRef.2,thefreeenergyofattractionincludesonlystateswhere Thisresultalsofollowsfromalarge-arealimitofEq.(19) atleastonebondisformed,thusβFattthereis−ln[(1+x)Np−1]. with random grafting points 1,8. We expect the simplifi- Our definition of Fatt is the free energy of the system of linkers, andincludestheconfigurationwherenobondsareformed. cationprovidedinthiscommunicationwillboosttheuse 8B. M. Mognetti, M. E. Leunissen, and D. Frenkel, “Controlling of our model for calculating interactions free-energy for the temperature sensitivity of DNA-mediated colloidal interac- general ligand-receptor mediated systems. tionsthroughcompetinglinkages,”SoftMatter8,2213(2012).