Melting transition in lipid vesicles functionalised by mobile DNA linkers

Abstract

We study phase behaviour of lipid-bilayer vesicles functionalised by ligand–receptor complexes made of synthetic DNA by introducing a modelling framework and a dedicated experimental platform. In particular, we perform Monte Carlo simulations that combine a coarse grained description of the lipid bilayer with state of art analytical models for multivalent ligand–receptor interactions. Using density of state calculations, we derive the partition function in pairs of vesicles and compute the number of ligand–receptor bonds as a function of temperature. Numerical results are compared to microscopy and fluorimetry experiments on large unilamellar vesicles decorated by DNA linkers carrying complementary overhangs. We find that vesicle aggregation is suppressed when the total number of linkers falls below a threshold value. Within the model proposed here, this is due to the higher configurational costs required to form inter-vesicle bridges as compared to intra-vesicle loops, which are in turn related to membrane deformability. Our findings and our numerical/experimental methodologies are applicable to the rational design of liposomes used as functional materials and drug delivery applications, as well as to study inter-membrane interactions in living systems, such as cell adhesion.

---

1 Introduction

Ligand–receptor interactions play a crucial role in a large variety of biological processes, including cell adhesion and signalling.¹ In biology, the selective nature of ligand–receptor interactions enables functional behaviours and responsiveness to environmental stimuli, which can be replicated in biomimetic materials where the interactions between the unit components are mediated by supramolecular ligand–receptor complexes.² The eminent example is represented by systems of DNA-functionalised colloidal units,^3–6 where the artificial DNA linkers can be designed to control phase behaviour^7–13 and self-assembly kinetics,^14–19 as well as to engineer biological probes.^20–23

A complete understanding of the complex emerging phenomena observed in multivalent interactions is only possible through a combination of experiments and numerical/theoretical analysis. Modelling multivalent interactions²⁴ is however a challenging task, as it requires the calculation of ensemble averages over the many possible configurations of the supramolecular linker complexes.^20,25 Analytical methods capable of capturing the resulting entropic effects have been recently developed.^26–35 Such theories have been utilised to calculate effective potentials in DNA mediated interactions between solid^26–31 and deformable particles,^35–39 as well as to design superselective probes.^{20,25,32–35,40} Analytical models have however limited applicability to those biologically and technologically relevant scenarios where non-specific contributions significantly affect the interactions. Effects neglected by analytical models include steric repulsion between the linkers and the deformability of the interacting surfaces. The latter is a particularly crucial aspect when dealing with multivalent interactions between soft substrates such as biological membranes. Numerical approaches could provide a faithful description of specific and non-specific effects in multivalent interactions. However, in view of the large interval of relevant lengthscales, from the molecular scale of ligands/receptors to the micron scale of the interacting units, even coarse-grained models are unsuitable to simulate phase behaviour and material properties in ensembles of micron-sized multivalent objects.⁴¹ Overcoming the limitation of purely analytical and purely numerical approaches is a critical step towards the development of truly predictive theoretical methods to aid the design of synthetic multivalent materials and improve our understanding of emergent behaviours in multivalent biological systems. A suitable approach should be capable of describing simultaneously the behaviour of the individual ligand/receptor pairs and the global phase behaviour of the systems.

The validation of such multiscale theoretical framework requires dedicated experiments in which the global phase behaviour can be disentangled from that of the individual ligand/receptor pairs, and the effect of the latter on the former can be assessed as a function of the number and adhesive strength of the linkers.

In this work we investigate the self-assembly behaviour of a biomimetic system of DNA-functionalised lipid vesicles,^{18,35,36,42–44} where artificial DNA linkers play the role of ligands/receptors and membrane deformability can potentially affect the resulting multivalent interactions. DNA linkers can freely diffuse on the surface of the vesicles and form either intra-vesicle loops or inter-vesicle bridges, the latter being responsible for attractive interactions and driving vesicle aggregation. We study the response of the system to temperature changes, and clarify how the aggregation/melting transition of the liposomes is affected by the competition between loop and bridge formation and the non-selective free energy contributions related in turn to membrane deformability.

We propose a new “hybrid” framework to calculate the free energy of the interactions between such vesicles, that combines state-of-art analytical theories developed for solid particles^{29–31,35–39} with Monte Carlo simulations that account for configurational costs related to membrane deformability.

By exploiting a fully automated and programmable platform, we perform experiments that for the first time are capable of simultaneously monitoring the self-assembly state of the liposomes and the bound/unbound state of the DNA linkers through a combination of fluorescence microscopy and Förster resonance energy transfer (FRET) measurements.

In systems where the number of linkers per vesicle N is low, simulations are capable of quantitatively replicating the response to temperature changes observed in experiments, including the aggregation/melting temperature of large unilamellar vesicle (LUV) clusters and its correlation with the temperature-dependent fraction of formed DNA bonds. Such agreement is lost at high N, most likely due to the effects of linker–linker and linker–vesicle steric interactions, neglected by our current model. Our numerical results confirm the importance and the nature of the previously hypothesised entropic effects on the hybridisation free energy of surface-tethered linkers.^36,37,39,45 In cases where linkers can diffuse on the surface of the substrates, such as DNA-tethers on lipid vesicles or ligands/receptors on cells,⁴⁶ these entropic costs include the loss of translational freedom following the formation of a bond. Particularly severe are the effects experienced by linkers forming inter-vesicle bridges, which end up confined within the relatively small adhesion patch between the vesicles. With the present numerical method these contributions can be directly evaluated based on the observed size of the adhesion patch, and then compared to the overall configurational free energy costs evaluated via density of states calculations.^36,39

Experiments and simulations demonstrate that these repulsive free energy contributions, combined with the competition between loop and bridge formation, have a substantial effect on self-assembly behaviour of DNA-functionalised vesicles: a minimum number of linkers N_dim is required to stabilise adhesion.^11,12 If fewer linkers are present vesicles do not aggregate.^11,12 The fair agreement between experimental and predicted value of the threshold number of linkers validates our methods as an useful tool to design biomimetic self-assembling systems featuring complex functionalities.

The paper is structured as follows. Section 2 reports on experimental methods. In Section 3 we present our modelling framework. In particular in Sections 3.2 and 3.3 we present respectively the analytical and numerical methods. The latter section illustrates how the two are combined. In Section 4 we discuss the simulations' outcomes for a large set of different system parameters (listed in Section 3.6) and perform a detailed comparison with experimental results.

2 Experimental methods

In this section we present our sample preparation protocols and experimental setup designed to provide a complete characterisation of the self-assembly behaviour of DNA-functionalised LUVs. In particular we can simultaneously monitor temperature dependent vesicle aggregation and the binding/unbinding state of the membrane-anchored DNA linkers in multiple samples using fully programmable and automated microscopy/fluorimetry apparatus. Experimental information can then be directly compared to numerical results to validate the assumptions and highlight the limitations of the numerical/theoretical framework described in Section 3.

2.1 Experimental system

As sketched in Fig. 1(a), LUVs with diameter D_exp ≈ 400 nm are functionalised by hydrophobised DNA linkers featuring a cholesterol anchor partitioning within the lipid bilayer, a rigid double-stranded DNA (dsDNA) spacer of length [small script l] ≈ 10 nm and a single-stranded DNA (ssDNA) sticky end. Membranes are in a fluid state, which enables cholesterol anchors to freely diffuse onto the vesicles' surface. On average, half of the overall 2N linkers present on each vesicle carry a sticky end a, the second half carry the complementary sequence a′. One linker type can bind to the other forming either inter-vesicle bridges or intra-vesicle loops, the latter being possible because each vesicle carries both types of linkers. Bridges are responsible for inter-vesicle adhesive forces that ultimately lead to aggregation of the LUVs. Sticky ends a and a′ are labelled with fluorophores Cy3 and Cy5 respectively, so that Förster resonance energy transfer (FRET) can be used to assess bond formation.

Fig. 1 (a) Experimental system. We consider LUV functionalised by two species of DNA linkers carrying complementary sticky ends a and a′ and a cholesterol anchor that partitions inside the lipid bilayer. Sticky ends a and a′ are labelled with fluorophores Cy3 and Cy5 as shown in the schematics. Two kind of bonds are possible: intra-vesicle loops and inter-vesicle bridges, the latter drive self-assembly of the suspension. (b) Simulation system. Lipid bilayers are modelled as triangulated meshes whose vertices are hard spheres (blue elements). Some of the vertices carry a reactive linker (green elements) that can eventually bind a complementary linker when closer than a distance L. Linkers forming a bond are labelled in red. In this snapshots loop formation is not allowed. The yellow elements highlight the vertices belonging to the adhesion patch between the vesicles.

2.2 Sample preparation

LUVs are prepared from dioleoyl-sn-glycero-3-phosphocholine (DOPC, Avanti Polar Lipids) doped with 0.8 molar percent Marina Blue 1,2-dihexadecanoyl-sn-glycero-3-phosphoethanolamine (DHPE, Life Technologies) for fluorescent imaging. Vesicles are prepared by extrusion in a 300 mM sucrose solution using a Mini Extruder (Avanti Polar Lipids) equipped with polycarbonate track-etched membranes with 400 nm pores (Whatman), as detailed in ref. 18.

Hydrophobised DNA linkers are pre-hybridised from the ssDNA

(i) 5′–CGT GCG CTG GCG TCT GAA AGT CGA TTG CGAAAA–3′ [Cholesterol TEG]

(ii) 5′–GC GAA TCG ACT TTC AGA CGC CAG CGC ACGA [Sticky End] A–3′ Cy3/Cy5,

where bases marked in bold belong to the dsDNA rigid spacer and unpaired A bases highlighted in italic are included to provide flexibility. Hybridisation is carried out in TE buffer (10 mM Tris, 1 mM EDTA, Sigma) with 100 mM NaCl as detailed in ref. 35.

Samples are prepared by mixing a 10 μl extruded vesicle solution with 90 μl iso-osmolar solution containing TE buffer, 87 mM glucose, 100 mM NaCl and a variable concentration of pre-assembled DNA linkers spanning from 12.7 μM to 0.25 μM, equally divided between a and a′ linkers. The nominal number of DNA linkers per vesicles is calculated assuming that all linkers partition into the bilayer,⁴⁷ that all the processed lipids form unilamellar vesicles with diameter of 400 nm, and that each lipid molecule contributes with 70 Ų to the bilayer area.⁴⁸ This estimate leads to a number of a and a′ linkers per vesicle between 71 ≤ N ≤ 3550. The assumption that all linkers partition onto the bilayer may breakup at high DNA density due to surface saturation, as discussed in Section 4.

For thermal processing and imaging, samples are injected into borosilicate glass capillaries (0.2 × 4 mm² rectangular inner section, CM Scientific), protected with a droplet of mineral oil at both ends, and permanently sealed with epoxy glue (Araldite).

2.3 Temperature cycling and imaging

Imaging and fluorescence spectroscopy are carried out on a fully automated Nikon Eclipse Ti-E inverted epifluorescence microscope equipped with a Nikon PLAN APO 20 × 0.75 NA dry objective and a Point Grey Research Grasshopper-3 GS3-U3-23S6M-C camera. Temperature is controlled with a tailor made Peltier stage hosting 7 capillaries with samples at different linker concentration. Temperature is repeatedly ramped up and down between 80 °C (or 90 °C) and 0 °C in steps of 1 °C and a rate of 0.5 °C min⁻¹, and measured with a precision <0.1 °C by a thermocouple located in close proximity of the capillaries.

At each temperature step, the motorised stage is moved to sequentially image all samples, recording snapshots with different fluorescence excitation wavelength and emission windows. For these different combinations of monochromatic LED sources (Philips Lumileds LUXEON UV, LUEXON Z), quad band fluorescence filter set (Semrock LED-DA/FI/TR/Cy5-A) and exchangeable motorized emission filters (Semrock FF01-600/37-25, FF01-731/137-25) are used. Specifically: Marina Blue (quad band filter set, excitation LED LHUV-0380-0200 at 380 nm); Cy3–Cy3 (quad band filter set, emission filter FF01-600/37-25, excitation LED LXZ1-PM01 at 530 nm); Cy3–Cy5 (quad band filter set, emission filter FF01-731/137-25, excitation LED LXZ1-PM01 at 530 nm); Cy5–Cy5 (quad band filter set, emission filter FF01-731/137-25, excitation LED LXZ1-PD01 at 627 nm). A single snapshot for each fluorescence channel is recorded on each sample, imaging an area of 563 × 352 μm² at a fixed height of a ≈8 μm from the bottom of the capillary. Dark frames where no excitation is used are acquired between each fluorescence image and used for background subtraction. All steps of the experiment are fully automated. A Perfect Focusing System (Nikon) enables correction for vertical drift.

2.4 Image analysis

Images collected on the Marina Blue channel are used for structural characterisation of vesicle aggregates, carried out as described in ref. 18. Briefly, Fourier analysis of the images enables the evaluation of a 2D projection of the samples' structure factor S(q), where q is the spatial frequency. The first moment of S(q) is then measured as an indicator of the aggregation state of the samples, going from high values for monomeric vesicles towards low values when aggregation takes place and S(q) develops a strong peak at low-q.

Images in the Cy3–Cy3, Cy3–Cy5 and Cy5–Cy5 channels are used to assess the efficiency of FRET between donor and acceptor fluorophores attached to a and a′ sticky ends respectively. The average intensities of each channel [scr I, script letter I]_Cy3–Cy3, [scr I, script letter I]_Cy3–Cy5, and [scr I, script letter I]_Cy5–Cy5 are extracted from the images, and used to evaluate the (ratio)_A

where , and the intensities and are measured in a reference sample that only contains Cy3 fluorophores. The (ratio)_A is linearly dependent on the FRET efficiency.⁴⁹ The Förster radius of the Cy3–Cy5 pair is ≈5–6 nm,⁵⁰ therefore high FRET efficiency is expected in bound a–a′ pairs, where donor and acceptor are kept at ≈2–3 nm from each other. The probability of FRET between unbound linkers is comparatively small, although not fully negligible in samples with the high DNA coverage, where the average distance between unbound a–a′ linkers goes down to ≈11.9 nm. As a function of temperature, (ratio)_A describes a sigmoidal curve from which the fraction of hybridised linkers ϕ(T) can be extracted aswhere B_H and B_L are linear fits of the high- and low-temperature plateaus of the sigmoidal curve. See discussion in Section 4.

3 Modelling framework

3.1 Interacting large unilamellar vesicles (LUV)

Vesicles are modelled as triangulated surfaces [see Fig. 1(b)] as described previously.^51,52 Following ref. 53–56 the vertices of the mesh are represented by hard spheres of diameter σ, taken as the unit length in our simulations. With the exception of the case described in Section 3.3.2, in this work we always simulate pairs of interacting vesicles. The Hamiltonian of two vesicles, each with N_v vertices is then given bywhere 〈α,β〉 and 〈i,j〉 denote neighbouring triangles and vertices respectively, with i, j ≤ 2N_v. In eqn (3)κ is the bending rigidity, n_α is the outward normal of the triangle labelled by α, and v_bnd(r_ij) is an infinite square-well potential that constraints two neighbouring vertices to be within a distance . v_excl(r_ij) is a hard core repulsion that constrains vertices i and j to stay at distance r_ij > σ. Note that v_excl acts also between vertices belonging to different vesicles.

We distribute 2N implicit linkers, N of type a and N of type a′, over the N_v vertices of each vesicle in a way that no more than one linker is allowed on the same vertex (2N < N_v). Two free complementary linkers can react if their distance is smaller than L. In the case of DNA linkers, L is equal to twice the length [small script l] of the spacers connecting the sticky ends to the cholesterol anchors. Accordingly, the Hamiltonian associated to linker–linker interactions is given by

where m and n run over all linkers of type a and a′ respectively and ε_m,n = 1 if linkers m and n are bound, 0 otherwise. The hybridisation free energy ΔG_m,n is equal to ΔG_L if linkers m and n belong to the same vesicle and can form a loop, or to ΔG_B if the linkers are on different vesicles and can form a bridge. The form of ΔG_L and ΔG_B is discussed in Section 3.2. The overall Hamiltonian is then given by

[script letter H] = [script letter H]_ves + [script letter H]_link. (5)

If only [script letter H]_link is considered, analytical models used in our previous studies^35,36 enable the calculation of the partition function of two interacting membranes, neglecting non-specific membrane–membrane interactions and configurational contributions associated to membrane deformation. These calculations adapted to the present system are given in Section 3.2. The hybrid numerical/analytical framework discussed in Section 3.3 combines analytical calculations with Monte Carlo methods to sample the non-specific contributions described by [script letter H]_ves (see eqn (5)) to the overall partition function Z, which is then used in Section 3.4 (see eqn (17)) to investigate vesicle dimerisation.

3.2 Analytical modelling of multivalent interactions

Interactions between multivalent objects are strengthened by the combinatorial entropic contributions accounting for the many different ways of forming a given number of bonds starting from a set of linkers. The experimental system features identical vesicles functionalised by an equal number N of two complementary linkers, which can therefore form both inter-vesicle bridges and intra-vesicle loops [see Fig. 1(a)]. In this scenario, the contribution of [script letter H]_link (eqn (4)) to the partition function of two vesicles linked by n_B bridges, taking two non interacting vesicles as reference, becomes³⁶whereIn eqn (6) and (7) we sum over all the possible number of loops since, consistently with neglecting steric interaction between linkers, we assume that loops are not affected by the conformation of the vesicles. In eqn (6)Ω_L(N,N) is the “linker” partition function of an isolated vesicle featuring only loops. ΔG_L and ΔG_B are defined as the free energies for loop and bridge formation. For system of mobile linkers ΔG_B includes the dimerisation free energy of the reactive groups of the linkers when free in solution, indicated as ΔG⁰, and a term ΔG^rot accounting for the loss of rotational freedom following the binding of two linkers. In this study we focus on the case of rod-like dsDNA linkers tipped by reactive ssDNA sticky ends, therefore ΔG⁰ is simply the hybridisation free energy of the sticky ends as obtained using nearest-neighbour rules for the sequences reported in Fig. 1(a).^57–59 We estimate the rotational term as ΔG^rot = −k_BT log[1/(ρ_⊖L³)], where ρ_⊖ = 1 M is the standard concentration.^35,36 More accurate configurational terms can be estimated accounting for variations in the distance between the two binding linkers. We leave such refinement to future investigations in which linkers will be simulated as explicit objects. Differently from ΔG_B, ΔG_L should also include an entropic cost ΔG^confining_L = −TΔS^confining_L of confining the otherwise diffusive linkers to within a short distance from each other.^35–39 Including an analogous term ΔG^confining_B in ΔG_B is not necessary, as confinement entropy is already accounted for by our simulation procedure that explicitly enables surface mobility in bridge-forming linkers (see Section 3.3.1). We estimate ΔG^confining_L = −k_BT log(L²/A), where A is the area of a vesicle.^35,36 In summary, we obtain
Recently it has been reported that the molecular roughness of the bilayer can alter the affinity between complementary linkers.⁶⁰ When properly parametrised,⁶¹ such effects can be included into the definition of ΔG₀.

In this work we combine multivalent partition functions, like the one in eqn (6) or (8), with Monte Carlo estimates of the configurational free energy costs of pairs of vesicles linked by n_B inter-vesicle bridges. In our model such costs only depend on the number of bridges n_B. In particular, the interaction free energy does not depend on the number of formed loops (see Section 3.3.3). For computational efficiency we therefore use a simplified system where only bridges are possible, which could be realised experimentally using two families of vesicles each carrying only one of the two complementary linkers. In this case the contribution to the partition function due to [script letter H]_link for vesicles with n_B bridges and N linkers (per type) is simply given by Ω_B(N,n_B) (eqn (8)). In Section 3.3 we use Ω_† (e.g. † = B or † = B, L) to tag the contribution to the partition function due to the selective part of a generic multivalent system.

3.3 Numerical estimate of configurational effects

In this section we explain our strategy to combine the analytical partition functions of interacting multivalent vesicles Ω_† (with † = B or † = B, L, see Section 3.2) with the numerical calculation of the configurational costs due to non selective terms of the Hamiltonian ([script letter H]_ves in eqn (5)) as estimated via Monte Carlo. These terms account for the deformation of the membranes following the formation of flat adhesion patches in the interacting vesicles, for the steric repulsion between membranes, as well as for the entropic costs of confining linkers forming bridges within the patch area. The role played by elastic and entropic terms in non-selective adhesion of vesicles with different morphologies has been extensively studied in the past.^62–64 However a direct comparison with previous results is not straightforward in view of the many peculiarities of multivalent ligand/receptor interactions.^29,60

Fig. 2(a) sketches the thermodynamic path we employ for the simulations. First, as detailed in Section 3.3.2, we use a hit-or-miss algorithm to calculate the internal partition function z(Ω_†,n_B = 1) of two vesicles linked by a single bridge with respect to the reference state of two free vesicles. By definition, the internal partition function includes the contributions of all degrees of freedom but the centre of mass of the free/linked vesicles. In Section 3.3.3 we then use the results of Section 3.3.2 to calculate the internal partition function of vesicles featuring an arbitrary number (n_B) of bridges z(Ω_†,n_B). In Section 3.3.1 we briefly outline our simulation algorithm.

Fig. 2 Thermodynamic pathway for the determination of vesicle dimerisation partition function. (a) The configurational costs entering the partition function of dimers are calculated by first sampling the conformation of two vesicles that can feature at least one bridge (Section 3.3.2), and then by sequentially adding more bridges (Section 3.3.3). (b) The dimer partition function, also including the multivalent terms of Section 3.2, enables to calculate the probability to form dimers.

3.3.1 Simulation procedure. In a single simulation cycle we sample on average all degrees of freedom of the system by means of local Monte Carlo moves. We attempt to randomly displace the nodes of the triangulated membrane, including the linkers possibly present on them, using the Metropolis algorithm and the Hamiltonian given in eqn (3). The fluidity of the lipid membrane is simulated by bond-switch Monte Carlo moves described in ref. 55 and 65. Local MC moves result in slow relaxation with autocorrelation times of the slowest fluctuation mode of the order of t_autcorr = 10⁴ MC cycles for the systems studied in this work. To enable sufficient sampling all simulated trajectories are longer than 5 × 10⁵ MC cycles. Consistently with the fact of having a fluid lipid membrane, linkers are allowed to diffuse on the triangulated membrane. For computational efficiency, the diffusion of the linkers on the triangulated mesh has been implemented using random jumps in which a randomly chosen linker is first selected and then moved to a randomly selected free vertex. All displacements that take two bound linkers to a distance bigger than L are rejected. The sampling of the bonds is done using a heat-bath algorithm.⁴⁵ In particular we randomly choose a linker i and create a list [script L]_i of all possible complementary linkers that could potentially form a bond with i, eventually including the linker to which i is already bound. The selected linker has then a probability p_free of becoming (or remaining) free and a probability n([script L]_i)p_bound of getting (or staying) connected to a randomly chosen partner from the list [script L]_i, which counts n([script L]_i) elements. Consistently with eqn (4) and (5), the probabilities are defined as
Note that as explained in Section 3.2 for efficiency reasons we simulate explicitly only the formation of bridges.

When studying systems of adhering vesicles (Section 3.3.3) we also attempt moves in which a vesicle is randomly chosen and rigidly translated along a random vector. The move is rejected if it causes two vesicles to overlap or a formed bridge to stretch beyond its maximum allowed bond-length L.

3.3.2 Configurational costs of forming the first bridge. In this section we calculate the internal partition function z(Ω_†,n_B = 1) of two vesicles bound by a single linkage by taking as a reference the internal partition function of two separate vesicles [first step in Fig. 2(a)]. Such partition function is given by the sum over all allowed conformations in systems of two vesicles, weighted by the number of possible inter-vesicle linkages in each conformation. We employ a Monte Carlo algorithm resembling what previously used to calculate the virial expansions of single chain observables of polymer suspensions^66,67 or pair interactions between particles functionalised by inert polymer brushes.⁶⁸ We divide the runs into M = 25 steps. At each step i (i = 1,…,M) we thermalise two non-interacting vesicles with N_therm = 10 000 MC cycles as defined in Section 3.3.1 and use the final configurations to estimate the internal partition function Q₁ of two vesicles each featuring a single linker. In particular we choose N_tries = 250 random displacements between the centres of mass of the two vesicles (r_j, j = 1,…,N_tries) uniformly sampled from the spherical shell of inner radius R_min,i and outer radius R_max,i. The radii R_min,i and R_max,i are chosen in a way to guarantee that all possible displacements between vesicle's centre of mass resulting in a linked configurations satisfy the condition R_min,i < |r| < R_max,i. Additionally a random rotation of the displaced vesicle around its centre of mass is performed. For each displacement and rotation, we loop over all the pairs of vertices belonging to different vesicles and count the number n^bonds_i of possible inter-particle bonds, with n^bonds_i = 0 if the two vesicles overlap. Q(1)_i is then estimated as
where 1/N_v² accounts for the probability of finding one linker on a given node of the membrane and the volume v_0,i encloses the displacement between the vesicles' centres of mass we sample, accounting for both changes in the absolute distance between the vesicles and rotations of one LUV around the other.

The final value of Q(1) is sampled using

The internal partition function of two vesicles featuring a given multivalent model Ω_† (see Section 3.2) linked by one bridge is then given by

z(Ω_†,1) = Q(1)Ω_†(N,1). (13)

The values of Q(1) for the systems considered in this work have been reported in Table 1.

Table 1 List of simulated systems featuring different vesicle size (number of vertices in the mesh) N_v, maximum linker length L and bending rigidity κ. In the fourth column we report the partition function of a pair of vesicles linked by a single bridge (eqn (12)). Relative statistical errors are ≈5–10%

Vertices Nv	Linker length L/σ	Bending rigidity κ/kBT	Q 1/σ^3
8004	3	20	0.15
4504	3	20	0.22
2004	3	20	0.40
4504	5	20	2.98
4504	10	20	67.14
4504	3	5	0.16
4504	3	30	0.23

---

3.3.3 Configurational costs of forming n_B bridges. The internal partition function of vesicles linked by n_B bridges z(Ω_†,n_B) is sampled using the successive umbrella sampling of Virnau and Müller⁶⁹ taking as bias parameter the number of bridges n_B [second step in Fig. 2(a)]. We successively constrain the algorithm to sample between states with n_B and n_B + 1 bridges. The ratio between the internal partition functions of the two states can then be directly evaluated as z(Ω_†,n_B + 1)/z(Ω_†,n_B) = [scr N, script letter N]_nB+1/[scr N, script letter N]_nB, where [scr N, script letter N]_nB+1 and [scr N, script letter N]_nB count the number of times the system visits states with n_B + 1 and n_B bridges respectively. Using z(Ω_†,1) as calculated in Section 3.3.2 we then obtain In Fig. 3(a) the ratio [scr N, script letter N]_nB+1/[scr N, script letter N]_nB is shown as a function of n_B for a bridge-only model († = B) with different amounts of linkers N. Similar to what done in eqn (13), we factorise the internal partition function using the selective term Ω_† and a configurational term indicated asfrom which we deriveNote that in view of eqn (15), ΔG^cnf is not a function of the hybridisation free energy of single linkers (ΔG_B and ΔG_L in eqn (9)). ΔG^cnf describes all the non-specific configurational effects that are not already included in the analytical multivalent partition function. Specifically, it accounts for membranes configurational entropy, and steric repulsion between the membranes, and confinement of the bridge-forming linkers to within the formed adhesion patch. For ideal linkers that do not interact sterically with each other or with the vesicles, ΔG^cnf is only a function of the number of bridges n_B, while it does not depend on the number of linkers N or the particular multivalent model used (as specified by † in Ω_†(n_B,τ)). This has been verified in Fig. 3(b) where we calculated ΔG^cnf(1 → n_B) by using a bridge-only model († = B, see Section 3.2) and three different values of N (N = 50, 100, 150). The results show an universal trend for ΔG^cnf. We observe that it increases almost linearly with n_B, which enables the extrapolation of to values of n_B above the simulated range.

3.4 Dimerisation probability

To study the melting transition of the DNA-functionalised vesicles we calculate the probability of forming liposome dimers starting from a dilute suspension of free vesicles [see Fig. 2(b)]. We start by evaluating the internal partition function of pairs of bound vesicles (Z) that is obtained by summing z(Ω_†,n_B) over all possible number of bridgesHere η is a symmetry factor equal to 2 if the vesicles are identical, to 1 otherwise. The dimerisation probability of two vesicles in a solution with a specific vesicle concentration c^exp can be written as

P_dim = (1 − P_dim)²Zc^expN_A, (18)

which can be solved for P_dim leading to the final expression Note that concentrations are here expressed in units of mol m⁻³, and N_A is Avogadro's number.

Fig. 3 Evaluating configurational free energy in pairs of interacting vesicles from MC simulations. (a) Ratio [scr N, script letter N]_nB+1/[scr N, script letter N]_nB between the number of times the system visits a configuration with n_B + 1 bridges and one with n_B bridges for different number of linkers (per type) N. Model parameters: ΔG₀ = 0 kJ mol⁻¹, N_v = 2004, and Ω_† = Ω_B. (b) Computed configurational free energy (1 → n_B) as a function of n_B for the same values of N as in panel (a). No significant dependence on N is observed (see inset). Note that accordingly to eqn (16) a different choice of ΔG₀ for the calculation of [scr N, script letter N]_nB+1/[scr N, script letter N]_nB would produce identical ΔG^cnf.

It is important to stress the portability of our method. In view of the universality of [eqn (16) and Fig. 3(b)], our procedure allows to calculate the internal partition function in pairs of vesicles at different temperatures or number of linkers N by re-weighting with the proper multivalent partition function (see eqn (15)), without the need of further MC simulations.

3.5 Other quantified observables

Besides the configurational contribution to the vesicle–vesicle interaction free energy ΔG^cnf, and the vesicle dimerisation probability P_dim, we quantify the area A_p of the adhesion patch, and the number of formed intra-vesicle bridges and inter-vesicle loops.

The adhesion patch between vesicles [see Fig. 1(b)] is defined as the region featuring vertices that could potentially bind to at least one vertex on the partner vesicle. Note that the exact position of the patch border depends on L. For numerical efficiency the patch area is once evaluated every 2500 MC steps. From the area A_p of the adhesion patch, we can evaluate the confining contribution to the configurational free energy of bridge formation as

where A is the total area of the vesicle. As discussed in Section 3.2, ΔG^confining_B describes the loss of translational freedom following the formation of a bridge, when two initially free linkers become confined to a distance L from each other and within the adhesion patch.

The average number of bridges is calculated as

If loops are present, their average number can be estimated using a saddle point approximation of the multivalent partition function in eqn (7), resulting in Using this expression along with eqn (21) the average amount of bound linkers (loops + bridges) in pairs of connected vesicles is given by Loops can also form in pairs of vesicles that are not connected. Accounting for this possibility we can estimate the total number of bound linkers as

〈n_B+L〉_tot = P_dim〈n_B+L〉_dim + (1 − P_dim)2n̄_L (24)

where P_dim is given in eqn (19).

Vesicle dimerisation temperature is estimated by evaluating the temperature dependence of P_dim. Starting from T = 0 °C and incrementing T in steps of 1 °C, the dimerisation temperature is defined as the first point where P_dim drops below 0.5.

Similarly, the melting temperature of the DNA linkers is defined as the temperature above which less than 50% of the linkers are hybridised.

3.6 Simulated systems

Table 1 lists the parameters of the simulated systems. We explored the influence of vesicle size, the maximum distance between bound linkers L, and the bending rigidity of the membrane κ on different properties of interest. Simulation units are converted to physical units by comparing the diameter of the vesicles in experiments D^exp, with that of the simulated liposomes D^sim. In particular the relation D^simσ = D^exp allows to convert the simulation unit length σ and the linkers' binding range L^exp = Lσ in physical units. In the system with N_v = 8004 and L = 3σ, using D^exp ≈ 400 nm we find σ = 9.1 nm.

4 Results and discussion

4.1 Configurational free energy

In this section we discuss numerical estimates of the configurational contribution ΔG^cnf(1 → n_B) to the interaction free energy of pairs of liposomes linked by n_B bridges, and disentangle components deriving from membrane deformation and steric repulsion from those caused by the confinement of bridge-forming linkers within the inter-vesicle adhesion patch. Fig. 4 shows the configurational free energy ΔG^cnf(1 → n_B) and the area A_p of the inter-vesicle adhesion patch as a function of the number of formed bridges at different vesicle size [(a) and (b)], linker length L [(c) and (d)], and membrane bending modulus κ [(e) and (f)]. In all cases, we find that ΔG^cnf increases linearly with the number of bridges, and this is correlated with the formation of larger adhesion patches at higher n_B. However the patch area is non-linear in n_B with A_p that increases more rapidly for small number of bridges. Similar trends have been previously reported when studying adhesion between fluid vesicles and solid supports at different vesicle–substrate attraction.^64,70

Fig. 4 Simulated configurational free energy for the formation of vesicle dimers ΔG^cnf(1 → n_B) (left column) and area of the adhesion patch A_p (right column) as a function of the number of inter-vesicle bridges n_B. (a and b) Effect of changing vesicle size (number of vertices in the triangulated mesh) N_v. (c and d) Effect of changing maximum bond length L. (e and f) Effect of changing bending modulus κ. Circles in the left column mark the confinement contribution ΔG^confining_B to ΔG^cnf(1 → n_B) evaluated using eqn (20) (multiplied by the number of bridges n_B) and values of A_p shown in the right column.

As discussed in Section 3.2, ΔG^cnf implicitly includes the entropic cost ΔG^confining_B of confining each linker engaged in a bridge within the patch region. This term can be estimated using the measured patch area (right column of Fig. 4) and eqn (20), and it is shown in Fig. 4 (left column, circles). Interestingly, n_BΔG^confining_B accounts for the entire configurational free-energy costs of binding a pair of vesicles, regardless of the simulation conditions tested in Fig. 4. Deviations of n_BΔG^confining_B from ΔG^cnf are difficult to quantify due to uncertainties in the estimation of A_p. From this observation we deduce that translational entropic terms hindering bridge formation are largely dominant over other contributions to ΔG^cnf. In particular steric repulsion between the membranes, and membrane stretching following adhesion have negligible effect. Note that ΔG^cnf accounts for the overall configurational free energy of the system, while ΔG^confining_B as derived in eqn (20) accounts for the contribution of a single linker, hence we compare ΔG^cnf with n_BΔG^confining_B.

Nonetheless, changing physical parameters of the vesicles causes changes in the area of the adhesion patch, and thereby in ΔG^confining_B and ΔG^cnf. In Fig. 4(a) and (b) we explore the effect of changing vesicle size. The patch area increases with the number of vertices N_v, but the ratio L²A_p/A² decreases, causing the configurational penalty for bridge formation to become more severe, effectively weakening the attraction between membranes. Similar considerations can be used to understand the results of Fig. 4(c) and (d) where we study the effect of changing the bond length L. Increasing L results in more stable bridges (see eqn (20)) and therefore in larger adhesion patches. However, especially for large L, the increase in patch area is also due to the fact that bridges made by longer tethers can explore a wider portion of the curved membrane region at the periphery of the adhesion patch. In Fig. 4(e) and (f) we test the effect of changing membrane bending modulus κ. In the range of values we tested, centred around the experimental bending modulus of DOPC bilayers κ ≈ 19k_BT,^71–73 we observe little effect on the patch size and thereby on the configurational free energy. This evidence confirms that contributions arising from the elastic deformation and the suppression of membrane thermal fluctuations are overwhelmed by the entropic terms related to bridge formation.^62,63 Our findings are consistent with the observation that the effect of thermal fluctuations becomes negligible for liposomes in the size-range of LUVs.⁷⁴

4.2 Vesicle aggregation and DNA hybridisation

In this section we study temperature dependent self-assembly of vesicle suspensions by means of our model and experiments. We show that experiments and modelling agree on the fact that at low number of binders N self-assembly is suppressed. Numerically, self-assembly is characterised by means of the vesicle dimerisation probability P_dim (see Section 3.4). In experiments, we make use of the fully automated microscopy/fluorimetry platform described in Section 2.3 and characterise aggregation of vesicles through Fourier analysis of epifluorescence images and the binding/unbinding state of DNA linkers via FRET.

For a meaningful comparison with experimental data we choose model parameters that better match the experimental ratio between vesicle diameter and linker length, in particular N_v = 8004 L = 3σ and κ = 20 (see Section 3.6). As done in experiments (see Fig. 1) we consider identical vesicles functionalised by two types of complementary linkers featuring intra-vesicle loops and inter-vesicle bridges († = L, B in Sections 3.2 and 3.3.3). Fig. 5(a) shows the calculated dimerisation probability as a function of temperature and the number of linkers per vesicle and type N. The P_dim(T) curves describe a sigmoidal shape, shifting towards higher temperature and becoming more sharp as the number of available linkers increases, a characteristic behaviour of multivalent interactions already observed in DNA functionalised solid particles.²⁷ As N is decreased, P_dim(T) tends to converge to a low-temperature plateau smaller than 1, and eventually smaller than 0.5 for N = 125, effectively suppressing dimerisation. This is a unique characteristic of multivalent interactions featuring competition between loops and bridges,^11,12 also confirmed by our experiments. Fig. 5(b) shows the experimental curves of the first moment of the structure factor as measured upon heating up vesicle samples with different N from low to high temperature. The step-like features mark the sharp melting of the aggregates. Also in experiments we observe a threshold value of N below which aggregation is suppressed. In samples with N = 177 partial self-assembly is observed in 2 out of 5 nominally identical repetitions, while the other 3 samples showed no aggregation, indicating that N = 177 is close to the threshold value for vesicle clustering. All samples with N = 71 showed no sign of aggregation, while all samples with N = 355 aggregated. The presence of a threshold value N_dim in the number of linkers, below which dimerisation does not take place, can be rationalised by observing that at low temperature, where the overall number of bonds (loops or bridges) is maximised, the driving force for vesicle adhesion is purely entropic and due to the combinatorial advantage of having a fraction of the linkers forming bridges.^11,12 At low N such gain is not sufficient to overcome the repulsive configurational contributions to the free energy, thus vesicle adhesion becomes unfavourable. The agreement between the experimental and simulated value of N_dim is remarkable, particularly in view of the high sensitivity of this value to changes in the model parameters. To exemplify such sensitivity, in Fig. 5(c) we report the simulated N_dim as a function of a hypothetical unbalance δ(ΔG_B − ΔG_L) between the hybridisation free energy of forming loops and bridges (see eqn (9)), where δ(ΔG_B − ΔG_L) = 0 marks our original choice described in Section 3. A small bias of 3 kJ mol⁻¹ between the two free energies produces a theoretical threshold value of N_dim ∼ 4000, more than one order of magnitude larger than the experimental value. In the inset of Fig. 5(c) we highlight how the experimentally determined N_dim can only be captured by simulations for uncertainties in the estimated ΔG_B − ΔG_L smaller than ±1 kJ mol⁻¹. Note that an unbalance between ΔG_B and ΔG_L would also result from a wrong estimate of ΔG^cnf. In this respect Fig. 5(c) certifies the accuracy of our model in conditions of low N.

Fig. 5 Temperature-dependent dimerisation of vesicles. (a) Simulated dimerisation probability as a function of temperature and number of linkers (per type) N. Model parameters: N_v = 8004, L = 3σ. (b) Experimental first moment of the structure factor in aggregating LUVs as a function of temperature and N. Data are recorded while heating up aggregated samples. Straight lines (solid) are used to fit the curves just below and just above the melting transition, and identify the melting temperature (star symbol). Samples with N = 71 never show aggregation, while samples with N = 355 always aggregate. Samples with N = 177 aggregate in 2 out of 5 independent repetitions of the experiment (not shown for clarity), indicating that N = 177 is close to the threshold value for vesicle dimerisation N_dim. First moment curves collected on cooling show a smoother trend due to the slow aggregation kinetics, and are unsuitable to accurately locate the phase transition. If multiple cooling/heating ramps are performed, we observe a systematic shift of the melting temperature towards lower T. We ascribe this effect to vesicle degradation resulting for the repeated exposure to high temperatures. (c) Predicted N_dim as a function of an hypothetical unbalance δ(ΔG_B − ΔG_L) between bridge and loop-formation free energies, where δ(ΔG_B − ΔG_L) = 0 marks our original choice. In the inset the shaded region highlights the experimentally determined range for the threshold value 71 < N_dim < 355.

In Fig. 6 we study the fraction of DNA strands engaged either in loops or bridges as a function of temperature, experimentally measured on the same systems of Fig. 5. Panel (a) reports simulation results obtained using eqn (24). Panel (b) shows the experimental FRET (ratio)_A, linearly dependent on the FRET efficiency between the fluorophores on complementary linkers. The sigmoidal decrease of FRET efficiency marks the melting of the DNA bonds, while the low and high-temperature plateaux represent the regimes where all the linkers are bound or free, respectively. The fraction of hybridised linkers can therefore be easily estimated by fitting such plateaux with linear baselines as explained in Section 2.4 and shown in Fig. 6(b).³⁵ Note however that for high DNA coverage (N = 3550, 1775), a second small drop in FRET efficiency is observed at high T (∼80 °C) as highlighted in the inset of in Fig. 6(b). In these samples DNA linkers are densely packed, and the average distance between them is comparable to the Förster radius of the Cy3–Cy5 pair, causing a non-zero chance of energy transfer also between unbound linkers. When the temperature is increased to ∼80 °C, the dsDNA spacers melt and the single-strands carrying the sticky ends and the fluorescent labels, no longer bound to the cholesterol anchors, are released in solution. The detachment of fluorophore-carrying DNA from the membranes causes the suppression of the small FRET signal ascribed to high-density coating. In samples where the double-transition is present the high-temperature baseline is chosen to fit the plateau observed before the final FRET-efficiency drop, which particularly for N = 3550 spans a relatively small temperature range, possibly leading to uncertainty in baseline determination. Similarly uncertainties are found in the samples with the lowest DNA coverage (N = 71) and thereby the lowest melting temperature, where the low-temperature plateau spans a small temperature range [see Fig. 6(b)]. The fraction of hybridised DNA linkers as extracted according to this procedure is shown in Fig. 6(c) as a function of temperature and N. The DNA melting temperature is extracted as the point where the fraction of hybridised DNA is equal to 0.5, determined via linear interpolation. Both simulation and experimental results show the broad melting transition typical of short oligomers.

Fig. 6 Melting of DNA links. (a) Simulated fraction of formed DNA links (loops + bridges) as a function of temperature for systems with variable number of linkers per vesicle. Model parameters: N_v = 8004, L = 3σ. (b) Experimental FRET (ratio)_A as a function of temperature. Solid lines indicate high- and low-temperature baselines fitted by straight lines (see Section 2.4). In the inset: a detail of the high temperature plateau in samples with high N showing a second drop in FRET efficiency at ∼80 °C that marks the disassembly of the dsDNA spacers connecting the linkers to the membranes (see Section 4.2). The curves are obtained averaging over 2 consecutive temperature ramps collected on heating and cooling. No hysteresis is observed. At high temperature [scr I, script letter I]_Cy3–Cy5 (see eqn (1)) becomes very small in samples with low DNA concentration and its detection is affected by background noise, e.g. sample autofluorescence. This causes the high temperature value of (ratio)_A to increase for samples with low N. (c) Experimental fraction of formed DNA bonds extracted from the ramps of panel (b) as detailed in Section 2.4. Legend applies also to panel (b).

Fig. 7 summarises and compares the dimerisation/aggregation temperature of the vesicles and the melting temperature of DNA linkers, as computed by our simulations [panel (a)] and measured experimentally [panel (b)]. Simulations predict a linear increase as a function of N for both the DNA melting temperature and the vesicle dimerisation temperature. However, the slope of the two curves is different, with the vesicle dimerisation temperature increasing more sharply than the DNA melting temperature. As a result, a crossover temperature exists, above which the vesicles can dimerise even when less than 50% of the DNA bonds are formed, and below which more than 50% of the available linkers are needed for dimerisation. In Fig. 7(a) we also investigate the effect of uncertainties in the hybridisation free energy of the sticky ends ΔG⁰. For all results presented in this work, ΔG⁰ has been calculated using conventional nearest neighbours rules^57–59 applied on the sticky-end sequences shown in Fig. 1(a), including also the attractive effect of the dangling bases (of type A) neighbouring the hybridising duplexes. However this modelling choice is far from being univocal. On the one hand, it has been demonstrated that the presence of inert DNA tails emanating from the duplexes, like the dsDNA spacers in the present system, can have substantial destabilising effect.^18,75–77 On the other hand, Cy3 and Cy5 fluorophores have a stabilising effect, decreasing the hybridisation free energy by ≈2 kJ mol⁻¹.⁷⁸ The dashed lines in Fig. 7(a) exemplify the effect of including this attractive free energy contribution to the vesicle dimerisation and DNA melting temperatures. In both cases we observe a shift of about 10 °C, which demonstrates how sensitive the present results are to uncertainties in the estimation of ΔG⁰.

Fig. 7 Melting temperatures for DNA bonds and vesicle aggregation. (a) Simulated melting temperatures as a function of the number of linkers per vesicle (per linker type) N. Data shown and solid lines are obtained using hybridisation free energy between DNA sticky ends ΔG⁰ evaluated using nearest-neighbours thermodynamic rules on the sequences shown in Fig. 1(a).^57–59 To produce dashed curves ΔG⁰ has been corrected with an attractive term (2 kJ mol⁻¹) to account for the stabilising effect of Cy3 and Cy5 fluorophores.⁷⁸ (b) Experimental melting temperatures obtained by averaging over 4 independent experiments. Errorbars indicate standard errors.

Experimental data in Fig. 7(b) show that the DNA melting temperature and the vesicle aggregation temperature approach linear dependence on N only at low DNA coverage, and tend to plateau at high N. We ascribe this trend to the possible saturation of the lipid membranes when very high concentration of DNA linkers is added in solution. Saturation would result in an actual number of DNA linkers per vesicle smaller than the nominal value calculated as described in Section 2.2. High concentration of hydrophobised linkers may also promote the formation of stable DNA-cholesterol micelles, which would decrease the available number of linkers. As they approach the linear regime at low N, also the experimental curves for DNA melting temperature and vesicle aggregation temperature have a different slope, with the latter being more steep. The difference in slope is however less pronounced as compared to simulation results, which causes the vesicles aggregation temperature to reach lower values at high N. This discrepancy is possibly due to steric repulsion between linkers, neglected in simulations. In particular, the adhesion patch between two vesicles features an even higher linker concentration than the surrounding free-standing membrane due to the recruitment of bridge-forming linkers.³⁶ Steric repulsion within the adhesion patch would therefore limit the number of possible bridges and weaken vesicle–vesicle adhesion.

A more insightful view in this effect is given by Fig. 8 where we compare the number of formed DNA bonds, including both loops and bridges, evaluated at the vesicle melting temperature. Simulation results [panel (a)] predict that the number of DNA bonds at vesicle melting does not depend on N for a broad range of model parameters. If steric interactions between the linkers are neglected, multivalent theories predict that the ratio between loops and bridges is independent on temperature or the total number of linkers, and determined only by vesicle geometry.³⁶ Thus the constant trends in Fig. 8(a) imply that a fixed number of bridges is required to overcome entropic repulsion and bind vesicles to each other.

Fig. 8 Fraction of formed DNA bonds evaluated at the melting temperature of vesicle aggregates as a function of the number of DNA linkers per vesicle N. (a) Simulation results. Curves display a constant trend in a broad range of model parameters (κ = 20k_BT). (b) Experimental results. The number of formed DNA bonds at vesicle melting temperature is extracted by linear interpolation of curves of the type shown in Fig. 6(c). Data are averaged over 4 independent experiments. Errorbars indicate standard errors. The curve shows an increasing trend as a function of N, approaching a plateau at low N and a linear asymptote (N/5) at high N.

Experiments predict an increase of the fraction of hybridised DNA at vesicle melting as a function of N, seemingly approaching a plateau at low N (compatibly with simulations) and a linear asymptote at high N. Since it is reasonable to argue that a fixed number of bridges is required to drive vesicle aggregation, we speculate that the observed trend may be caused by a dependence of the bridge/loop ratio on N. This could be again ascribed to steric repulsion between linkers within the crowded adhesion patch: although at high linker concentration the number of formed DNA bonds at vesicle melting is higher, excluded volume effects between linkers in the patch would result in a smaller fraction of bridges and a higher fraction of loops with respect to the ideal scenario in which linker–linker steric interactions are negligible.

5 Conclusions

Multivalent interactions are at the forefront of many applications in nanotechnology and underpin several functional behaviours in biology. A satisfactory understanding of all aspects of the problem is lacking, particularly in biologically relevant scenarios where substrate deformability couples with configurational free energy contributions and steric interactions. A numerical/theoretical framework validated by dedicated experiments and capable of assessing the consequences of such interplay is therefore of pivotal importance for biological and technological applications. Modelling multivalent interactions requires a significant coarse graining effort to bridge atomistic scales at which non-covalent bonds form with mesoscopic scales at which functionalised objects (vesicles, surfaces, polymers or nanoparticles) interact as a result of many possible reversible supramolecular linkages. Analytical theories have been developed to calculate the free energy of multivalent interactions.^{26,27,29–35} Considering a set of interacting multivalent objects, these methods attempt to enumerate all possible configurations of the supramolecular network of linkers mediating the interactions. The success of these theories enabled the development of efficient simulation schemes for characterising the phase behaviour of multivalent objects, coarse-grained as simple particles with effective pair interactions. Analytical multivalent theories, however, neglect non-selective interactions arising from particle deformability or excluded-volume effects between linkers.

In this work we addressed some of the limitations of current modelling approaches to multivalent interactions, particularly the impossibility of accounting for deformable particles. Working on a system of soft liposomes functionalised by linkers made of synthetic DNA, we propose a method that combines state-of-art multivalent theories with Monte Carlo methods. In particular we used a triangulated model of the lipid bilayer together with free energy calculations to estimate the configurational free energy cost of having a given number of bridges formed between interacting vesicles. We clarify how such contributions are mainly due to the entropic penalty of confining bridge-forming linkers within the flat contact area formed between adhering deformable vesicles. We then characterise the response of pair of interacting vesicles to changes in temperature, and in particular the temperature-dependent dimerisation probability and the melting curves of DNA linkers. Simulation results are compared to experiments on DNA functionalised large unilamellar vesicles, where the temperature-dependent vesicle aggregation state and DNA hybridisation can be independently monitored by fluorescence microscopy and Förster resonant energy transfer thanks to a fully automated and programmable setup. Both simulations and experiments confirm that a minimum number of linkers per vesicle is required to overcome configurational entropic costs for membrane deformation and produce stable aggregation. We observe deviations between simulated and experimental trends at high density of the DNA linkers. We argue that this disagreement is caused by excluded volume effects between pairs of linkers and between the linkers and the membranes, neglected in the present contribution to maximise the portability of the model. These deviations deserve future investigation.

Our experimental and numerical results highlight the importance of configurational free energy costs arising from the deformability of objects interacting via multivalent interactions. This is an ubiquitous scenario in biological contexts, where deformable cells adhere to each other or to the extra-cellular matrix thorough membrane ligand/receptors, but also in bio-nanotechnology and nano-medicine, where multivalent synthetic probes are designed to selectively target cells. The novel numerical approach we developed to reach these conclusions combining state-of-art analytical modelling with Monte Carlo simulations provides a valuable tool for further investigations of specific biological and nano-technological problems including tissue dynamics, cell sorting, cell–cell and cell–substrate adhesion, tissue scaffolding, and designing multivalent probes for drug and gene delivery.

Acknowledgements

SJB and BMM are supported by the Université Libre de Bruxelles (ULB). JK, LP, LDM and PC acknowledge support from the EPSRC Programme Grant CAPITALS number EP/J017566/1. LDM acknowledges support from the Oppenheimer Fund, Emmanuel College Cambridge, the Leverhulme Trust and the Isaac Newton Trust through an Early Career Fellowship. AS thanks the Human Frontier Science Program and Emmanuel College Cambridge. BMM thanks Ignacio Rondini for contributing to the early stages of this work. Computational resources have been provided by the Consortium des équipements de Calcul Intensif (CECI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11. In compliance with the requirements of EPSRC a representative dataset underlying this publication is available for download at http://dx.doi.org/10.17863/CAM.1301.

References

1. N. J. Gay, M. F. Symmons, M. Gangloff and C. E. Bryant, Nat. Rev. Immunol., 2014, 14, 546–558.
2. E. Sackmann and A.-S. Smith, Soft Matter, 2014, 10, 1644–1659.
3. C. A. Mirkin, R. C. Letsinger, R. C. Mucic and J. J. Storhoff, Nature, 1996, 382, 607.
4. A. P. Alivisatos, K. P. Johnsson, X. Peng, T. E. Wilson, C. J. Loweth, M. P. Bruchez and P. G. Schultz, Nature, 1996, 382, 609.
5. L. Di Michele and E. Eiser, Phys. Chem. Chem. Phys., 2013, 15, 3115–3129.
6. M. R. Jones, N. C. Seeman and C. A. Mirkin, Science, 2015, 347, 1260901.
7. S. Y. Park, A. K. R. Lytton-Jean, B. Lee, S. Weigand, G. C. Schatz and C. A. Mirkin, Nat. Mater., 2008, 451, 553–556.
8. D. Nykypanchuk, M. M. Maye, D. van der Lelie and O. Gang, Nat. Mater., 2008, 451, 549–552.
9. Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck and D. J. Pine, Nature, 2012, 491, 51–55.
10. J. D. Halverson and A. V. Tkachenko, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 062310.
11. S. Angioletti-Uberti, B. Mognetti and D. Frenkel, Nat. Mater., 2012, 11, 518–522.
12. W. B. Rogers and V. N. Manoharan, Science, 2015, 347, 639–642.
13. J. T. McGinley, I. Jenkins, T. Sinno and J. C. Crocker, Soft Matter, 2013, 9, 9119–9128.
14. M. E. Leunissen, R. Dreyfus, F. C. Cheong, D. G. Grier, R. Sha, N. C. Seeman and P. M. Chaikin, Nat. Mater., 2009, 8, 590–595.
15. Q. Xu, L. Feng, R. Sha, N. C. Seeman and P. M. Chaikin, Phys. Rev. Lett., 2011, 106, 228102.
16. W. B. Rogers, T. Sinno and J. C. Crocker, Soft Matter, 2013, 9, 6412–6417.
17. J. D. Halverson and A. V. Tkachenko, J. Chem. Phys., 2016, 144, 094903.
18. L. Parolini, J. Kotar, L. Di Michele and B. M. Mognetti, ACS Nano, 2016, 10, 2392–2398.
19. Y. Wang, Y. Wang, X. Zheng, é. Ducrot, J. S. Yodh, M. Weck and D. J. Pine, Nat. Commun., 2015, 6, 7253.
20. F. J. Martinez-Veracoechea and D. Frenkel, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10963–10968.
21. N. A. Licata and A. V. Tkachenko, Phys. Rev. Lett., 2008, 100, 158102.
22. T. A. Taton, C. A. Mirkin and R. L. Letsinger, Science, 2000, 289, 1757–1760.
23. J.-M. Nam, C. S. Thaxton and C. A. Mirkin, Science, 2003, 301, 1884–1886.
24. M. Mammen, S.-K. Choi and G. M. Whitesides, Angew. Chem., Int. Ed., 1998, 37, 2754–2794.
25. P. I. Kitov and D. R. Bundle, J. Am. Chem. Soc., 2003, 125, 16271–16284.
26. R. Dreyfus, M. E. Leunissen, R. Sha, A. V. Tkachenko, N. C. Seeman, D. J. Pine and P. M. Chaikin, Phys. Rev. Lett., 2009, 102, 048301.
27. R. Dreyfus, M. E. Leunissen, R. Sha, A. Tkachenko, N. C. Seeman, D. J. Pine and P. M. Chaikin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 041404.
28. W. B. Rogers and J. C. Crocker, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 15687–15692.
29. P. Varilly, S. Angioletti-Uberti, B. Mognetti and D. Frenkel, J. Chem. Phys., 2012, 137, 094108.
30. S. Angioletti-Uberti, P. Varilly, B. Mognetti, A. Tkachenko and D. Frenkel, J. Chem. Phys., 2013, 138, 021102.
31. L. Di Michele, S. J. Bachmann, L. Parolini and B. M. Mognetti, J. Chem. Phys., 2016, 144, 161104.
32. N. B. Tito, S. Angioletti-Uberti and D. Frenkel, J. Chem. Phys., 2016, 144, 161101.
33. H. Xu and D. E. Shaw, Biophys. J., 2016, 110, 218–233.
34. R. De Gernier, T. Curk, G. V. Dubacheva, R. P. Richter and B. M. Mognetti, J. Chem. Phys., 2014, 141, 244909.
35. S. Shimobayashi, B. M. Mognetti, L. Parolini, D. Orsi, P. Cicuta and L. Di Michele, Phys. Chem. Chem. Phys., 2015, 17, 15615–15628.
36. L. Parolini, B. M. Mognetti, J. Kotar, E. Eiser, P. Cicuta and L. Di Michele, Nat. Commun., 2015, 6, 5948.
37. L. Feng, L.-L. Pontani, R. Dreyfus, P. Chaikin and J. Brujic, Soft Matter, 2013, 9, 9816–9823.
38. L.-L. Pontani, I. Jorjadze, V. Viasnoff and J. Brujic, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 9839–9844.
39. S. Angioletti-Uberti, P. Varilly, B. M. Mognetti and D. Frenkel, Phys. Rev. Lett., 2014, 113, 128303.
40. G. V. Dubacheva, T. Curk, R. Auzély-Velty, D. Frenkel and R. P. Richter, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 5579–5584.
41. S. Angioletti-Uberti, B. M. Mognetti and D. Frenkel, Phys. Chem. Chem. Phys., 2016, 18, 6373–6393.
42. M. Hadorn, E. Boenzli, K. T. Sørensen, D. De Lucrezia, M. M. Hanczyc and T. Yomo, Langmuir, 2013, 29, 15309–15319.
43. R. J. Banga, N. Chernyak, S. P. Narayan, S. T. Nguyen and C. A. Mirkin, J. Am. Chem. Soc., 2014, 136, 9866–9869.
44. P. A. Beales and T. K. Vanderlick, J. Phys. Chem. A, 2007, 111, 12372–12380.
45. M. E. Leunissen and D. Frenkel, J. Chem. Phys., 2011, 134, 084702.
46. A.-S. Smith, K. Sengupta, S. Goennenwein, U. Seifert and E. Sackmann, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 6906–6911.
47. I. Pfeiffer and F. Höök, J. Am. Chem. Soc., 2004, 126, 10224–10225.
48. N. Kučerka, S. Tristram-Nagle and J. F. Nagle, J. Membr. Biol., 2006, 208, 193–202.
49. R. M. Clegg, A. I. Murchie, A. Zechel and D. M. Lilley, Proc. Natl. Acad. Sci. U. S. A., 1993, 90, 2994–2998.
50. F. Yuan, L. Griffin, L. Phelps, V. Buschmann, K. Weston and N. L. Greenbaum, Nucleic Acids Res., 2007, 35, 2833–2845.
51. A. Šarić and A. Cacciuto, Phys. Rev. Lett., 2012, 108, 118101.
52. A. Šarić and A. Cacciuto, Soft Matter, 2013, 9, 6677–6695.
53. Y. Kantor, M. Kardar and D. R. Nelson, Phys. Rev. Lett., 1986, 57, 791–794.
54. J.-S. Ho and A. Baumgärtner, Europhys. Lett., 1990, 12, 295.
55. A. Baumgärtner and J.-S. Ho, Phys. Rev. A: At., Mol., Opt. Phys., 1990, 41, 5747–5750.
56. G. Gompper and D. M. Kroll, J. Phys.: Condens. Matter, 1997, 9, 8795.
57. J. SantaLucia, Proc. Natl. Acad. Sci. U. S. A., 1998, 95, 1460–1465.
58. J. SantaLucia Jr and D. Hicks, Annu. Rev. Biophys. Biomol. Struct., 2004, 33, 415–440.
59. R. M. Dirks, J. S. Bois, J. M. Schaeffer, E. Winfree and N. A. Pierce, SIAM Rev., 2007, 49, 65–88.
60. J. Hu, R. Lipowsky and T. R. Weikl, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 15283–15288.
61. G.-K. Xu, J. Hu, R. Lipowsky and T. R. Weikl, J. Chem. Phys., 2015, 143, 243136.
62. R. Lipowsky and U. Seifert, Mol. Cryst. Liq. Cryst., 1991, 202, 17–25.
63. R. Lipowsky and S. Leibler, Phys. Rev. Lett., 1986, 56, 2541.
64. T. Gruhn and R. Lipowsky, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 011903.
65. H. Noguchi, J. Phys. Soc. Jpn., 2009, 78, 041007.
66. B. Li, N. Madras and A. D. Sokal, J. Stat. Phys., 1995, 80, 661–754.
67. S. Caracciolo, B. M. Mognetti and A. Pelissetto, J. Chem. Phys., 2008, 128, 065104.
68. B. M. Mladek and D. Frenkel, Soft Matter, 2011, 7, 1450–1455.
69. P. Virnau and M. Müller, J. Chem. Phys., 2004, 120, 10925–10930.
70. G. Linke, R. Lipowsky and T. Gruhn, Eur. Phys. J. E: Soft Matter Biol. Phys., 2007, 24, 217–227.
71. W. Rawicz, K. C. Olbrich, T. McIntosh, D. Needham and E. Evans, Biophys. J., 2000, 79, 328–339.
72. B. Sorre, A. Callan-Jones, J.-B. Manneville, P. Nassoy, J.-F. Joanny, J. Prost, B. Goud and P. Bassereau, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 5622–5626.
73. S. A. Rautu, D. Orsi, L. Di Michele, G. Rowlands, P. Cicuta and M. S. Turner, 2015, arXiv preprint, arXiv:1511.05064v1.
74. U. Seifert and R. Lipowsky, Phys. Rev. A: At., Mol., Opt. Phys., 1990, 42, 4768.
75. L. Di Michele, B. M. Mognetti, T. Yanagishima, P. Varilly, Z. Ruff, D. Frenkel and E. Eiser, J. Am. Chem. Soc., 2014, 136, 6538–6541.
76. N. Srinivas, T. E. Ouldridge, P. Šulc, J. M. Schaeffer, B. Yurke, A. A. Louis, J. P. Doye and E. Winfree, Nucleic Acids Res., 2013, 41, 10641–10658.
77. C. Wang, J. H. Bae and D. Y. Zhang, Nat. Commun., 2016, 7, 10319.
78. B. G. Moreira, Y. You and R. Owczarzy, Biophys. Chem., 2015, 198, 36–44.

---

Footnote

† These authors contributed equally to this work.

---