Ligand-receptor interactions are ubiquitous in biology and have become popular in materials in view of their applications to programmable self-assembly. Although complex functionalities often emerge from the simultaneous interaction of more than just two linker molecules, state of the art theoretical frameworks enable the calculation of the free energy only in systems featuring one-to-one ligand/receptor binding. In this Communication, we derive a general formula to calculate the free energy of systems featuring simultaneous direct interaction between an arbitrary number of linkers. To exemplify the potential and generality of our approach, we apply it to the systems recently introduced by Parolini *et al.* [ACS Nano **10**, 2392 (2016)] and Halverson and Tkachenko [J. Chem. Phys. **144**, 094903 (2016)], both featuring functionalized Brownian particles interacting via three-linker complexes.

## I. INTRODUCTION

The quantitative understanding of the ligand-receptor interactions is receiving much attention in view of the key role played in biology and their applications to the self-assembly of composite materials.

Biological cells respond to the presence of specific molecules via cell-surface receptors. Examples include *toll-like receptors*, triggering immune response to bacterial and viral activity,^{1} and *receptor tyrosine kinases*, involved in the regulation of several physiological processes.^{2} In order for the signals to be transduced across the cell membrane, the presence of the ligands typically triggers dimerization or oligomerization of the receptors, through interactions that involve multiple molecules.

Functionalizing Brownian units with specific linkers, often made of synthetic DNA molecules, are powerful tools to engineer the structure and response of self-assembled soft materials.^{3–11} Many functionalization schemes rely on one-to-one ligand-receptor interactions, but recently, designs featuring multi-linker complexation have been proposed to extend the accessible range of functionalities.^{9,12–15} In particular, Parolini *et al.*^{13} adopted three-linker complexes enabling toehold-mediated strand exchange reactions^{16} to control aggregation kinetics of lipid vesicles coated with DNA linkers. Halverson and Tkachenko^{12} also proposed the use of three-linker DNA complexes to program a cooperative behavior between functionalized particles, which could allow to control the sequence of binding events in the self-assembly.

Recently, Angioletti-Uberti *et al.*^{17} proposed an analytical expression for the free energy of systems featuring one-to-one ligand-receptor interactions that overcame some limitations of earlier approaches.^{18,19} In this Communication, we provide a more general framework to calculate the free energy of systems including multimeric complexes featuring an arbitrary number of ligand/receptors (see Fig. 1). We consider “particles,” e.g., biological cells or artificial colloidal units, functionalized by surface ligands/receptors (“linkers” or “molecules”). We assume that linkers can freely diffuse on the surface of the particles. An extension to immobile linkers can be derived following Ref. 20. Bonds can either involve linkers tethered to the same particle or to different particles. Excluded volume interactions between the molecules are neglected. Our results are exact in the limit of many linkers per particle.^{21,22} We envisage applications of our theory to the association of more complex molecules like DNA tiles^{23–25} or virial caspids.^{26}

In Sec. II, we derive our theory while in Sec. III we test it on the system introduced in Ref. 12, calculating the interaction free energy between particles and quantitatively justifying the postulated cooperative behavior. In Sec. IV, we examine the suspensions of DNA-functionalized vesicles of Ref. 13, discussing the thermodynamic ground state in relation to the kinetic behaviour characterized in the original publication.

## II. FREE ENERGY CALCULATION

We consider *c* families of different linkers, each with a number of units *N _{i}* (

*i*= 1, …,

*c*). The linkers can reversibly associate into complexes of

*m*units. For clarity we only consider complexes that never feature more than a single linker of each family (1 ≤

*m*≤

*c*, see Fig. 1 where

*c*= 3). In Sec. S1 of the supplementary material,

^{27}we show that relaxing this assumption does not change our main result (Eq. (7)). The state of the system is described by the number

*n*

_{i1,i2,…im}of all the possible complexes made by

*m*linkers of type

*i*

_{1},

*i*

_{2}, …,

*i*, with

_{m}*i*

_{1}<

*i*

_{2}< ⋯ <

*i*and 2 ≤

_{m}*m*≤

*c*.

We start by deriving an expression for the partition function *Z* of the system as the weighed sum over all the possible realizations of *n*_{i1,i2,…,im}. First we calculate the contribution of two-linker complexes {*n*_{i1,i2}} to *Z*, then we deplete the total number of linkers of each family *N _{i}* by the number of those involved in two-linker complexes and calculate the contribution from complexes with three molecules {

*n*

_{i1,i2,i3}}. This procedure is repeated recursively. When calculating the contribution of complexes with

*m*+ 1 linkers,

*N*has been reduced to $ N i ( m ) $ that is given by

_{i} where $ n i ( \u2113 ) $ is the total number of linkers of type *i* involved in complexes of size ℓ, and *τ* is the operator that orders *m* indices. $ N i ( c ) $ is the number of linkers of type *i* that are free, and will be also indicated by *n _{i}* below. The partition function is then expressed as

where the curly brackets indicate the ensemble of all the complexes formed by a given number of linkers. Note that in Eq. (2), *Z*^{(ℓ)} is a function of $ { N i ( \u2113 \u2212 1 ) } 1 \u2264 i \u2264 c $ and, as a consequence of Eq. (1), of the number of complexes with *m* ≤ ℓ.

Defining Δ*G*_{i1,…,im} as the free energy associated to the formation of a *i*_{1}⋯*i _{m}* complex,

^{28}we can define the contribution to the partition function from all the complexes of size

*m*as

where *β* = 1/(*k*_{B}*T*) and Ω^{(m)} accounts for the combinatorial factors. The latter can be written as

where the first product counts the number of ways one can choose the molecules belonging to the complexes {*n*_{i1,i2,…,im}} starting from $ { N i ( m \u2212 1 ) } $ free linkers, while the second term accounts for the number of independent ways to build such a set of complexes. Using Eqs. (4) and (3) in Eq. (2), we can calculate the partition function and the free energy *F* of the system

where the double curly brackets {{…}} indicate the ensemble of complexes of arbitrary size, and $A$ is a functional introduced for later convenience.

In the limit of large *N _{i}*, we can simplify Eq. (5) using a saddle point approximation. In particular the stationary point of $A$, given by $\delta A/\delta n i 1 , \u2026 , i m =0$, identifies the average number of complexes $ n \xaf i 1 i 2 \cdots i m =\u3008 n i 1 i 2 \cdots i m \u3009$. The stationary conditions for the functional $A$ as defined by Eq. (5) become

Note that Eq. (6) is the relation for chemical equilibrium expressed in terms of the total number of molecules. When considering tethered linkers (Fig. 1), the complexation free energy Δ*G*_{i1,…,im}^{28} also includes rotational and translational entropic costs controlled by the length of the spacers and by the size of the particles.^{7}

Using Eq. (6) in Eq. (5) to express Δ*G*_{i1,…,im} as a function of equilibrium number of complexes, we can evaluate the free energy of the system as $F=A ( { { n \xaf i 1 , i 2 , \u2026 , i \u2113 } } ) $. By considering only the dominant term in the second line of Eq. (5), and using Stirling’s approximation we find

where Eq. (1) has been used in the second equality to factorize the terms $log n \xaf i $, and in the third equality to express $ N i \u2212 n \xaf i $ in terms of higher order complexes. Finally we obtain the main result of this work

Being written in terms of the equilibrium number of complexes, Eq. (7) cannot be used to sketch free-energy landscapes,^{24,25} however it is applicable to calculate effective interactions between functionalized objects as demonstrated in Secs. III and IV. Note that in order to guarantee the extensivity of the functional $A$, Δ*G*_{i1,…,im} − (*m* − 1)log*N* (with *N _{i}* ∼

*N*) should be kept fixed when taking the

*N*→ ∞ limit (see also Eq. (6)). For the case of one-to-one interactions (

*m*= 2), Eq. (7) reduces to the result of Ref. 17.

## III. BINDING COOPERATIVITY IN DNA-FUNCTIONALIZED PARTICLES

As a first example, we examine the cooperative self-assembly scheme recently proposed by Halverson and Tkachenko,^{12} based on the possibility of forming three-linker complexes, dubbed spiders. As shown in Fig. 2(a), we consider particles of type *A* functionalized by 3*N* mobile DNA linkers equally distributed among three families, each carrying different single-stranded DNA sequences or *sticky-ends*, labelled as *α _{i}*,

*i*= 1, 2, 3. Such sticky-ends can hybridize to form three different families of intra-particle

*loops*(ℓ

_{i}), involving two out of three types of linkers, or spiders (

*s*), involving all three types (see Fig. 2(a)). We then consider three types of particles

*B*,

_{i}*i*= 1, 2, 3, each functionalized by

*N*identical linkers carrying a sticky-end sequence $ \alpha \xaf i $ complementary to

*α*. Linkers on particles

_{i}*B*can form inter-particle

_{i}*bridges*

*b*, with particles

_{i}*A*. In the following, we consider linkers constituted by double stranded DNA spacers of length

*L*= 10 nm and point-like sticky ends,

^{21}rigid particles of radius

*R*= 10

*L*,

^{21}and

*N*= 100. See Sec. S2 of the supplementary material

^{27}and Ref. 21 for details. Below we calculate the free energy

*F*(

*AB*) of clusters made by a single

_{n}*A*particle and a variable number

*n*of

*B*particles taken as in Fig. 2. We demonstrate a cooperative effect by which the free energy gain from binding the

*n*th

*B*particle Δ

*F*=

_{n}*F*(

*AB*) −

_{n}*F*(

*AB*

_{n−1}) is higher than the gain from binding the (

*n*− 1)th one, for

*n*= 2, 3. This is due to the necessity of breaking spider and loop complexes formed on the

*A*particle for the 1st or the 2nd

*B*particles to bind. Our theory allows to calculate the free energy gain for binding the 1st, the 2nd, and the 3rd

*B*particles, chosen as a model parameters in Ref. 12. The number of complexes at equilibrium are given by

^{21}

where $\Delta G \u2113 0 $, $\Delta G s 0 $, and $\Delta G b 0 $ are the hybridization free energies of the sticky-ends associated to loop, spider, and bridge formation respectively, *ρ*_{⊖} = 1M is the standard concentration, and *v*_{A/B/AB} are volume factors reported in Sec. S2 of the supplementary material^{27} that quantify the configurational entropic costs of binding mobile tethers (Refs. 7 and 21 and Sec. S1 of the supplementary material^{27}). Note that by using the generic notation of Sec. II, we would have $ n \xaf \u2113 k \u2261 n \xaf \alpha i \alpha j $, $ n \xaf s \u2261 n \xaf \alpha 1 \alpha 2 \alpha 3 $, and $ n \xaf b i \u2261 n \xaf \alpha i \alpha \u0304 i $. Different types of loops and bridges are assumed to have the same hybridization free energy. In Eq. (10) *ϵ _{i}* = 0, 1 specifies if

*B*is bound or not to

_{i}*A*. In particular $n= \u2211 i \u03f5 i $ indicates the valence of particle

*A*. Equations (8)–(10) are then closed by the conditions $N= n \xaf \alpha i + n \xaf \u2113 j + n \xaf \u2113 k + n \xaf s + n \xaf b i $ and $N= n \xaf \alpha \u0304 i + n \xaf b i $.

First we consider an isolated *A* particle and calculate the number of loop and spider complexes as a function of $\Delta G \u2113 0 $, choosing $\Delta G s 0 =3\Delta G \u2113 0 $. As shown in Fig. S1 of the supplementary material,^{27} when $\Delta G \u2113 0 =\u221210\u2009 k B T$ only spiders are present on *A*. We fix $\Delta G \u2113 0 $ to this value as a reasonable guess to maximize the cooperative behaviour. We then consider particle clusters *AB _{n}* (with

*n*= 0, 1, 2, 3), with distances between the centers of

*A*and

*B*particles equal to

*d*= 2

*R*+

*L*, and calculate the number of bridges

*n*

_{bi}as a function of $\Delta G b 0 $ (see Fig. 2(b)). We find that bridges form at higher values of $\Delta G b 0 $ when

*n*is higher. Finally we use Eq. (7) (contextualized to this system in Eq. (S14) of the supplementary material

^{27}) to calculate the free energy of the system including the repulsive part of the potential calculated accounting for the entropic compression of the DNA strands between the particles (see Eq. (S2) of the supplementary material

^{27}). We consider clusters in which all of the

*B*particles are at the same distance

*d*from the

*A*-particle, and for which

*B*-particles do not interact with each other (see Fig. 2(a)). Figure 2(c) shows the free-energy change associated to the binding of a single

*B*particle to a cluster as a function of

*d*. As expected, the free-energy gain obtained when adding the second

*B*particle is higher than that obtained by binding the first, and the gain achieved upon adding the third particle is significantly higher than both the former.

## IV. INTERACTION FREE ENERGY IN THE PRESENCE OF TOEHOLD-EXCHANGE-MECHANISM

As a second example, we examine the system studied experimentally in Ref. 13. Let us consider a suspension of identical micron-size lipid vesicles, functionalized by three families of mobile DNA linkers with sticky ends here labelled as *γ*, *δ*_{1}, and *δ*_{2}. As shown in Fig. 3(a), sticky end *γ* is made of three domains of equal length, *x*, *y*, and *z*. Sticky ends *δ*_{1} have two domains $ x \u0304 $ and $ y \u0304 $, complementary to *x* and *y*, whereas *δ*_{2} features domains $ y \u0304 $ and $ z \u0304 $. Linker *γ* can bind to *δ*_{1} and *δ*_{2}, with comparable hybridization free energy. A three-linker *γδ*_{1}*δ*_{2} complex is also possible, where *δ*_{1} and *δ*_{2} bind to domains *x* and *z*, respectively, and compete to occupy domain *y*. *δ*_{1} does not bind to *δ*_{2}. Two- and three-linker complexes can form either among linkers tethered to the same vesicles (loop-like) or between different vesicles (bridge-like). At a sufficiently high temperature, all of the linkers are unbound. If the suspension is quenched to low temperature, the formation of intra-vesicle loop-like complexes is kinetically favored over the formation of bridges, effectively sequestrating all of the available *γ* linkers. The aggregation of the liposomes, mediated by the formation of inter-vesicle bridges, is therefore limited by the opening of the intra-vesicle loops, which seldom occurs at low temperatures (Fig. 3(b), top). Through a *Toehold-Exchange Mechanism* (TEM),^{16} the formation of three-strand complexes mediates the swap between stable loops and stable bridges without the need for thermal denaturation. In particular, the toehold domain *z* (*x*) causes a free *δ*_{2} (*δ*_{1}) linker to transiently bind to an existing *γδ*_{1} (*γδ*_{2}) bond, facilitating the reaction *γδ*_{1} + *δ*_{2}⇋*δ*_{1} + *γδ*_{2} (Fig. 3(b), bottom). We indicate with 3*N* the total number of linkers per vesicle, *N* of which are of type *γ*, *Nχ* of type *δ*_{1}, and *N*(2 − *χ*) of type *δ*_{2}. The parameter *χ* ∈ [0, 2] controls the stoichiometric ratio between *δ*_{1} and *δ*_{2} and thereby the effectiveness of the TEM process. For *χ* = 0 or 2 three-strand complexes are not possible and the bridge formation and aggregation kinetics are dominated by the slow opening of formed loops. For *χ* = 1, TEM is most effective and aggregation kinetics is found to speed up by more than one order of magnitude at *T* = 15 °C.^{13}

We use our framework to calculate the free energy of the system, and demonstrate that, despite the large effect on aggregation kinetics, changing *χ* has little consequences on the thermodynamic ground state of the system. The DNA tethers are again modelled as freely pivoting rigid rods of length *L* = 10 nm, with freely diffusing tethering points and point-like sticky ends. For simplicity, we model two interacting vesicles as flat planes of area *A* = 0.5 *μ*m^{2} kept at a distance of *h* = 1.4*L* from each other. We chose *N* = 360. Hybridization free-energies between the sticky ends are taken from Ref. 13.

Explicit expression for the equilibrium distributions of all the possible complexes are shown in Eqs. (S15)-(S20) of the supplementary material.^{27} In the supplementary material (Eqs. (S21) and (S22)^{27}), we provide the expression for the interaction free energy between two vesicles (*per* *γ* strand), shown as a function of *χ* and *T* in Fig. 3(c). We observe that regardless of temperature, the free energy decreases by less than 10% when going from *χ* = 0, 2 to *χ* = 1, supporting the claim that with the architecture proposed in Ref. 13 aggregation kinetics can be substantially changed with little consequences on the thermodynamic ground state. The weak dependence of the overall free energy on *χ* is a direct consequence of the small number of three-strand complexes, always involving less than 10% of all *γ* linkers, as demonstrated in Fig. 3(d).

## V. CONCLUSIONS

We provide an analytical expression for the free energy of systems of ligand/receptors that can form complexes featuring an arbitrary number of molecules. Our framework can be applied to biologically relevant situations, where cell-surface receptors form trimers or oligomers, or to suspensions of colloidal particles functionalized by synthetic DNA ligands: an increasingly popular strategy to achieve controlled self-assembly of complex soft materials. To exemplify the versatility of our approach, we re-examine the artificial systems recently proposed by Halverson and Tkachenko^{12} and Parolini *et al.*,^{13} both featuring DNA-functionalized Brownian particles interacting through the formation of three-linker complexes. For the former, we are able to quantify the cooperative effects in the interaction free energy between the particles, taken as model parameters in the original publication. For the system of Parolini *et al.* we study the interaction free energy between vesicles with different linker stoichiometry. Our theory demonstrates that despite the substantial effect on aggregation kinetics observed experimentally, coating stoichiometry has a comparatively small effect of the thermodynamic ground state of the suspension.

## Acknowledgments

L.D.M. and L.P. acknowledge support from the EPSRC Programme Grant CAPITALS No. EP/J017566/1. L.D.M. acknowledges support from the Oppenheimer Fund and Emmanuel College Cambridge. S.B. and B.M.M. are supported by the Université Libre de Bruxelles (ULB). The python programs used to calculate the results presented here are available at http://dx.doi.org/10.5281/zenodo.47204.