Predicting the size and morphology of nanoparticle clusters driven by biomolecular recognition

Abstract

Nanoparticle aggregation is a driving principle of innovative materials and biosensing methodologies, improving transduction capabilities displayed by optical, electrical or magnetic measurements. This aggregation can be driven by the biomolecular recognition between target biomolecules (analytes) and receptors bound onto nanoparticle surface. Despite theoretical advances on modelling the entropic interaction in similar systems, predictions of the fractal morphologies of the nanoclusters of bioconjugated nanoparticles are lacking. The morphology of resulting nanoclusters is sensitive to the location, size, flexibility, average number of receptors per particle f̄, and the analyte–particle concentration ratio. Here we considered bioconjugated iron oxide nanoparticles (IONPs) where bonds are mediated by a divalent protein that binds two receptors attached onto different IONPs. We developed a protocol combining analytical expressions for receptors and linker distributions, and Brownian dynamics simulations for bond formation, and validated it against experiments. As more bonds become available (e.g., by adding analytes), the aggregation deviates from the ideal Bethe's lattice scenario due to multivalence, loop formation, and steric hindrance. Generalizing Bethe's lattice theory with a (not-integer) effective functionality f_eff leads to analytical expressions for the cluster size distributions in excellent agreement with simulations. At high analyte concentration steric impediment imposes an accessible limit value f_acc to f_eff, which is bounded by f_acc < f_eff < f̄. A transition to gel phase, is correctly captured by the derived theory. Our findings offer new insights into quantifying analyte amounts by assessing nanocluster size, and predicting nanoassembly morphologies accurately is a first step towards understanding variations of physical properties in clusters formed after biomolecular recognition.

---

1 Introduction

While pioneer studies on colloids were applied to prevent their aggregates, the last two decades¹ has witnessed a growing interest in controlling colloid aggregation starting from the nanoscopic details up to the meso and macroscales. Self-assembly yields novel mesoscopic structures with applications in smart materials, self-healing structures,² drug delivery,^3,4 light-matter interaction^5,6 or biomolecular sensors. One route to programmable self-assembly consists on designing the directionality or spatial arrangement of surface patches for colloid–colloid interactions.⁷ Another route, considered hereby, consists on nanoparticle conjugation, where bonds are mediated by high-affinity molecular linkers (the analytes in biosensing) binding to receptors located at the nanoparticles surface. Clusters of bioconjugated nanoparticles driven by analyte–receptor interactions are usually based on antibody–antigen couples,⁸ protein pairs,⁹ or complementary ssDNA sequences.^10,11 Aggregates created by specific biomolecular recognition are opening new routes in biosensing, to detect from proteins to viruses.¹² As in nature's self-assembled biochemical networks (some of them recently reproduced in living cells using two-protein mixtures¹³), the concentration of free linker molecules in solution enables an extra thermodynamic parameter with subtle antagonistic effects, such as reversible re-entrant transition to fluid phase at high linker concentration¹¹ (owing to saturation of receptors) or a softening the otherwise strong dependence on gelation with temperature.¹⁴ Non-specific interactions are key to understand aggregation by analyte–receptor bonds¹⁰ as they lead to either effective repulsion between particles (e.g. excluded volume, steric effects),^11,14 or to attractive interaction (e.g. dispersion forces). While non-specific attraction is inherent to colloids, in these systems, the effective repulsion depends non-trivially on the functionality of the particle, the size and flexibility of the linker and receptor,¹⁰ and probably even their surface mobility.¹⁴ Interestingly, these non-trivial effects could in principle be tuned to modify the effective functionality (or valence) of the conjugated nanoparticles,¹⁰ an effect which we also study hereby. Recent theoretical and computational works on polymer-linked colloidal gels^15–17 report features which are qualitatively similar to those found here for protein-linked nanoparticles: high sensitivity to the receptor location,¹⁵ to the linker size and flexibility,¹⁶ to the average number of receptors per particle f̄, and to the analyte–particle concentration ratio. General features of the phase diagram, such as the re-entrant phase behaviour¹⁷ can be captured by thermodynamic approaches such as Wertheim's thermodynamic perturbation theory (ref. 18 and references therein). More specific analytical efforts based on statistical mechanics have focused on predicting the binding free energy (and mean field potentials) from the combinatorial and configurational entropy associated to analyte–receptor mixtures, also including effective repulsion.^11,14,19 Theoretical and experimental works show that the structures of these aggregates depend on how much nanoscopic details can be kept under control: from “addressable complexity” in weak (reversible) attractive interaction, leading to crystal phases,⁷ to aperiodic, amorphous or irreversibly gellified phases.¹⁴ Concerning the kinetics of aggregation, following the foundational works by Flory²⁰ and Stockmayer²¹ and subsequent developments for equilibrium polymerization,^22,23 theoretical and computational generalizations for different types of linker-mediated complexes have been more recently proposed by several groups.^{15,17,24–27} The aggregation strongly dependent on the relative size of the linkers and host particles, and usually treated as diffusion limited. Despite these advances, there is still a lack of studies on the fractal structure of linker mediated aggregates. For a certain parameter range we show that these aggregates also share some features with limited-valence colloids²⁸ with strong biding energy (E_bind/k_BT ∼ 10), which percolate at low packing fraction forming open structures, or “empty fluids”.²⁹

Biomolecular recognition-mediated nanoassembling is relevant in many applications related to medical care, based on nanoparticles (NP) of different compositions, sizes, and shapes.^30–32 Most of the available biosensing methodologies benefit from nanoparticle aggregation.³³ This is the case of photothermal³⁴ converter, colorimetric^34–37 optical,^5,38 electrical³⁹ or magnetic readouts.⁴⁰

Uniform nanoparticles can be now prepared with rather well defined surface characteristics, allowing for covalent stoichiometric attachment of oligonucleotides, nucleic acids, small molecules, peptides, proteins, antibodies, etc.⁴¹ The linkers might have a valence larger than one, capturing several receptors, while nanoparticles can be decorated with multiple receptors. In bioconjugation, the number of receptors per particle is a random variable f, whose statistics (binomial distribution) can be only fixed by experimentally controlling its average value f̄. The complexity of linker-mediated clusters thus lies some where between standard colloids with non-specific interactions and the addressable complexity of some patchy colloids (fixing f and the patches surface location). Their aggregates structures highly depend on the particle and analyte concentrations, the particle and receptor sizes and the average functionality f̄. We note however that f̄ is seldom experimentally controlled;⁴² suspensions are often saturated with receptors which physically adsorb to the particles forming a coating.

There is great practical interest in relating all these physico-chemical details with the resulting cluster size distribution and morphology, notably in the case of biosensors, either based on magnetic nanoparticles,^43,44 plasmonic resonance and colorimetric biosensors⁴⁵ and many other techniques where the transduced signal is enhanced by cluster formation.⁴⁶ We present a general analysis of this problem, which combines experiments, theory and simulations to study the morphology and abundance of clusters of bioconjugated nanoparticles, with experimentally controlled stoichiometric constraints. We show that the information about the average stoichiometry of the receptor–analyte bond is essential to control the cluster distribution, for instance by tuning the relative analyte–particle concentration.

The system studied is illustrated in Fig. 1. Iron oxide nanoparticles (IONP) coated with dextran and carboxylic PEG are bioconjugated with the protein GST-MEEVF⁴⁷ which binds to the divalent protein dVFP-TPR2-MMY.⁴⁸ In experiments, we control the ensemble average number of MEEVF receptors per IONP (we consider f̄ = {4,8,12}), the nanoparticle concentration c_p and linker (analyte) concentration c_a. The system's modeling was conducted by integrating analytical techniques to solve the distribution of receptors over the IONPs (Fig. 1c) and the formation of primary analyte–receptor bonds (Fig. 1d), along with numerical simulations employing Brownian Dynamics (see Methods) to address the slow IONP–IONP polymerization 1e) process.⁴⁹ Simulations are successfully validated against experimental measurements of the average cluster size, polydispersity index (PDI) and the transition to a condensed colloidal phase (experimentally observed by pellet formation). Simulations indicate that this condensed phase correspond to a percolated network of protein-linked nanoparticles.‡ We analyze the effect of steric hindrance by neighboring particles in reducing their effective functionality and derive a simple expression for the accessible number of receptors (f_acc < f̄) at high bond concentration. We then investigate the morphology of the clusters (number of bonds per particle and fractal dimension) as a function of f̄. A generalization of Bethe's lattice analytic results, using a non-integer effective functionality f_eff is shown to successfully recover the cluster size distributions n(s). The effective valence varies with the analyte concentration: from f_eff ≈ f̄ (Bethe's limit at low bond density) to f_eff ≈ 〈f_acc〉 (high bond concentration). Using the high density limit, f_acc, into Bethes lattice theory we predict the locus of the gelation transition observed both in experiments and simulations. Concluding remarks are given after the discussion of the results.

Fig. 1 Sketch of the system illustrating (a) one iron oxide nanoparticle (IONP) of radius 15 nm with a magnetic core covered with a coating of dextran polymer of about 10 nm. (b) Atomic force microscopy (AFM) images of IONP assemblies compared with simulation configurations. (c–e) Illustration of the three steps of the experimental and theoretical protocol: (c) bioconjugation: a concentration c_p of IONPs are first incubated along with a concentration of c_r receptors (orange), fixing the average number of receptors per particle f̄. (d) Formation of primary receptor–analyte bonds (analytes in green). This process is fast and irreversible (high affinity, 13k_BT). (e) Formation of secondary (IONP–IONP) bonds and nanoclusters (slow kinetics solved by Brownian dynamics). At high analyte concentration non-specific (spurious) bonds (not involving receptors) might be formed due to adhesion of analytes to the dextran coating.

2 Methods

2.1 Experimental protocol

Experiments were carried out with magnetic iron oxide nanoparticles (IONPs) coated with dextran and carboxylic PEG synthesised by Micromode (product code is 104-56-701). The core diameter of the IONPs is 30 ± 4 nm size, the hydrodynamic diameter is 64 nm and the PDI is 0.08. Analyte–receptor recognition. IONPs were bioconjugated with a specific peptide sequence MEEVF,⁴⁷ that is recognized by the design TPR domain TPR2-MMY⁴⁸ that was used as target divalent analyte. Bioconjugation of magnetic nanoparticles. For the bioconjugation of IONPs coated with carboxylic PEG, three samples of 500 μL of IONPs at 3.5 mg of Fe per mL were incubated 4 hour at 37 °C with 150 mmol of EDC per g of Fe and 150 mmol of NHS per g of Fe. Samples were washed with amicon ultra centrifugal filters and redispersed in 10 mM sodium phosphate buffer pH 7.4 (PB buffer) at least 3 times. These pre-activated IONPs were incubated overnight at 37 °C with different concentrations of the receptor in order to estimate a ratio of 4, 8 and 12 receptors per particle. After that, the bioconjugated IONPs were purified by filtration through a sepharose 6 CLB column using PB buffer. Receptor-bioconjugated IONPs incubation. 50 μL of MEEVF-IONPs at 1 g L⁻¹ of Fe were incubated 1 hour at room temperature with different concentrations of divalent analytes (dVFP-TPR2-MMY) in PB buffer. Quantification of magnetic nanoparticles. We employed a Nanosight NS300 (Malvern Instruments, USA) in order to determine the number of IONPs per unit volume in magnetic suspensions with iron contents of 1 g Fe per L. The samples were diluted 1 : 5000 in DDW and injected into the instrument chamber using a 1 mL syringe. Camera settings were adjusted in order to focus the objective. The video data were collected for 60 seconds and repeated 3 times per sample. Colloidal Characterization of selected IONPs was performed through hydrodynamic size measurement by dynamic light scattering (DLS) in a Zetasizer Nano ZS equipment (Malvern Instruments, USA). In order to perform DLS measurements, the bioconjugated magnetic nanoparticles at 1 g Fe per L were incubated with increasing concentration of the analytes. After 60 minutes, each sample was diluted in buffer PB to a final concentration of 0.05 g Fe per L in a commercial cuvette. The energy source was a laser emitting at 633 nm, and the angle between sample and detector was 173°. DLS uses the Stokes–Einstein relation to determine the hydrodynamic size of solutes (either isolated particles or clusters), as D_H = k_BT/(3πηD_t), where D_t is the translational diffusion coefficient measured from DLS. Using dilute IONP suspensions we measured single-particle size D_H → σ = 64 nm. Moreover, DLS permits us to determine the nanoclusters average size and their polydispersity index (PDI). As explained below, above a critical value of the analyte concentration c^(g)_a, experiments show the formation and precipitation of a condensed (pellet) phase. In such case, we performed DLS measurements after having resuspended the solution (note however that such measurements are not representative).

2.2 Theoretical and computational modeling

We decompose the overall process of nanoassembly formation in three separate steps, illustrated in Fig. 1. First, we analytically derive the distribution of receptors per nanoparticle, for a given average number of receptors per particle f̄. This step is independent of the rest as in experiments the particles are first incubated with receptor proteins prior to linker molecules addition. Second, we solve (both analytically and by random sampling) the distribution of primary analyte–receptor links, which permit us to draw ensembles of particles with the expected distribution of loaded and void receptors for given values of f̄, particle concentration c_p and analyte concentration c_a. Third, for each particular realization of the ensemble of bioconjugated particles, we use Brownian dynamics (BD) simulations to model the aggregation via the formation of secondary (receptor–analyte–receptor) bonds (Fig. 1). Statistics of nanoassembly formation are gathered by repeating this three-step protocol. In the experiments, the second and third steps occur concurrently. However, as illustrated in Fig. 1, both stages can be independently solved due to the large time scale separation between the (fast) kinetics of formation of primary analyte–receptor bonds and the subsequent (slow) IONP–IONP polymerization.⁴⁹ Both processes are diffusion-limited and irreversible due to the high receptor–analyte affinity (binding free energies of about 13k_BT). This fact, and the large size discrepancy of the analyte proteins and the nanoparticles, ensures a solid separation between primary and secondary bond formation. In particular, the hydrodynamic diameter of the IONPs particles is 64 nm, while that of analyte proteins is around 4 nm, leading to diffusion times of τ_IONPs ≈ 1200 μs and τ_analytes ≈ 0.3 μs where τ = D_t/σ². We now provide details on each of the three steps of our modelling approach.

2.2.1 Distribution of receptors on particles. Experimentally we control the average number of receptors per particle f̄. However, the number of receptors on each IONP follows a binomial distribution. It can be assumed that the maximum number of receptors that each particle can have (that depends on the number of free carboxilic groups in the dextran coating) is much larger than the average number receptors per particle. Hence, the bonding of a receptor to a particle is independent of the number of receptors that the particle already had. Consequently, if there are N_p particles, each receptor can be bonded to a specific particle with a probability p_r = 1/N_p, and the probability of a particle having f receptors bonded follows a binomial distribution,where F = N_p·f̄ is the total number of receptors in the system. The number of particles in the simulations and especially in the experiments is very high, consequently p_r ≪ 1 and F ≫ 1 and the binomial distribution can be approximated by a Poisson distribution,

In practice, instead of using eqn (2) it is faster to use random sampling to distribute the receptors among the particles: we generate F = N_pf̄ random numbers uniformly distributed between 1 and N_p. Then, the number of receptors associated to a particle i is just the number of times that the number i appears in the list of generated numbers.

In Fig. S4 (ESI†) we have compared the predictions for P_r(f) using the Poisson distribution, the Binomial distribution and random sampling. An excellent agreement is found for the two limiting cases considered hereby (f̄ = 2 and f̄ = 12). Thus, for computational convenience, we deploy random sampling to populate the receptors on the IONPs.

2.2.2 Distribution of primary analyte–receptor links. The formation of primary analyte–receptor links can be treated as a stochastic process. In analogy with the random sampling for the distribution of receptors on the particles, we now assign one identity index “id” to each receptor and generate a list of N_a (number of analytes) random integers with a uniform distribution between 1 and F. This list of generated numbers represents the list of receptors that are initially occupied with analyte in the simulation.

The random sampling protocol samples the probability of a particle having f_occ occupied receptors, P(f_occ). To verify the numerical sampling we need to derive an analytical relation for P(f_occ). We derive in the ESI† (Section S.5) a compact form for P(f_occ) which consist on another Poisson distribution,

Somewhat surprisingly this distribution does not depends on f̄, provided that the number of receptors in the system is greater than or equal to the number of analytes c_a < f̄c_p. Note that when the system is saturated of analytes, c_a > f̄c_p, the distribution of occupied receptors is trivially equal to the distribution of receptors, eqn (2). In Fig. S5 (ESI†) we compare the distribution of particles having f_occ receptors initially occupied obtained by random sampling and by eqn (3). The comparison is basically perfect, so that to fasten up computations, random sampling is the preferred choice to set the initial conditions of the BD simulations.

2.2.3 Brownian dynamics for nanoassembly formation. Simulations of the aggregation process were carried out using Brownian Dynamics (BD) with translational and rotational displacements. IONPs are taken as spherical particles with diameter, σ = 2a, equal to the experimental value, 64 nm. The receptors are fixed to the surface and move rigidly with it. The position of the receptors on the particles surface is set randomly. The receptor “status” can be either “free” or “occupied” by an analyte. The initial receptor status is set according to the sampling protocol explained above. Recall that we consider divalent analytes and we consider that a secondary receptor–analyte bond is formed when the distance between the occupied and the free receptors is smaller than 2 times the longitudinal length of the receptor (approximately 4 nm). Newly formed bonds are simulated with harmonic springs, as explained shortly.

Simulations were performed in cubic periodic boxes of side L = 2000 nm ≈ 31σ. The translational diffusion coefficient D_t (given by the Stokes–Einstein relation D_t = k_bT/(6πηa) provides the physical time scale of the simulation. The time step was set to dt = 1 × 10⁻⁴ in units of τ = a²/D_t. Individual runs reached 1.5 × 10⁷ time steps, which is enough for a particle to diffuse through the full box τ_box = L_box²/D_t. As shown in Fig. S7 (ESI†), these simulation times are enough to guarantee that the number of bonds formed saturates to a stationary value. While two or three of these runs are enough to obtain a good estimation of the mean hydrodynamic size and PDI of the clusters, we need between 30 and 200 runs to extract the cluster-size distributions.

All simulations have been conducted using a particle concentration of c_p = 1.27 μM, which corresponds to 6111 particles in our simulation box, with the exception of those in which we aimed to study the effect of particle concentration on the average hydrodynamic size and the PDI of the nanoclusters. In these cases, we varied the concentration of MNP from 0.049 μM up to 2.9 μM (i.e., from 200 IONP to 12 000 IONP). Additionally, throughout all simulations, we used analyte concentrations ranging from 0.1 to 4.0 μM (i.e., 411 analytes to 19 270 analytes).

Brownian dynamics for the particle displacements were solved using the Euler–Maruyama method,⁵⁰ a generalization of the Euler method to stochastic differential equations. Not being interested in the precise aggregation kinetics, but rather on the statistics of the nanoclusters, we did not include hydrodynamic correlations between particles displacements, which might modify cluster formation dynamics but not their structure.⁵¹ The interaction forces between IONPs are just steric to prevent overlapping. The steric force was implemented using the Weeks–Chandler–Andersen (WCA) potential⁵² which is a purely repulsive, short range potential based on Lennard-Jones interaction, truncated at r = 2^1/6σ,

where r is the distance between two particles and ε = 0.11k_BT is the energy parameter. On the other hand, whenever a bond is formed between two IONPs, a harmonic spring is set between the two receptors linked by the same analyte. The harmonic force is:with d = r + a(u_j − u_i) (Fig. S8, ESI†). 2[small script l] = 4 nn is the equilibrium distance between the surface of two bonded particles and u_i, u_j are the unitary vectors pointing the position of the receptors on the IONP surfaces (see Fig. S8, ESI†). The spring constant k was set to k = 100k_BT/a², which yields bond fluctuations of about 2 nm, within the expected range for the sum of thermal variations in receptors, analyte–protein and dextran-coating (in any case, results were not observed to vary using larger values of k).

The translational displacement of each particle over time dt, is described by standard Brownian dynamics,

where M_t = 1/(6πηa) is the translational mobility, η is the viscosity of the solvent and dW a Wiener increment (three random components with zero mean, taken from a Wiener process, such that 〈dW_i²〉 = dt).

In analogy to the translational movement, the overdamped Langevin equation for the rotational motion includes a drag torque, a stochastic torque and all the external torques acting on the particles, which in our case they arise from the harmonic forces between the linked receptors,

τ_i = au_i × F_i (7)

We note that the steric forces do not exert any torque because they are central forces. The infinitesimal change in the orientation of the particles, dq is,^53,54where dq is the change in any set of generalised coordinates defining the orientation of the particles, M(q) a mobility matrix and dW_ϕ the components of the Wiener increment vector.

To integrate eqn (8), we choose a set of generalised coordinates q such that the angular velocity ω, is ω = dq/dt. Using this choice ∂_qU is the total torque that acts on the particles and the mobility matrix M(q) is the rotational mobility, M_r, that relates the torque with the angular velocity as t = M_rω. In addition if we consider quasi-spherical particles the stochastic drift term k_BT∂_qM(q) is 0 and the mobility tensor is a diagonal matrix. Using this assumptions eqn (8) can be rewritten as,

with M_r = 1/(8πηa³). From the angular velocity we can obtain a rotation vector dϕ = ωdt whose norm (θ) is the rotation angle and whose direction (u_ϕ) is the direction of the rotation axis.

In the literature, we can find multiple ways of parameterize the orientations and rotations of each particle⁵³ (i.e., using the Euler Angles, storing the directions of the three-coordinate axis of the particles…). We describe the orientation of each particle using quaternions (Q_i(t))⁵⁵ because they provide a singularity-free description of rotations and they are very efficient in the use of the memory (only 4 numbers are required to store the orientation of each particle).

The rotation vector that actualize the orientation of each particle can be easily transformed to a quaternion Q_ϕ,⁵⁶

Q_ϕ = (cos(θ/2), sin(θ/2)u_ϕ) (10)

The orientation of the particle in the time t + dt, Q_i(t + dt), can be obtained by multiplying their orientation in the time t, Q(t), by Q_ϕ,

Q_i(t + dt) = Q_ϕ·Q_i(t) (11)

The position of the receptors on the surface of the IONP was stored using the polar and the azimuthal angles in the reference system of the IONP. The code was written and compiled in CUDA-C/C++ using the library UAMMD⁵⁷ and simulations were performed in our own groups GPU cluster.

3 Results and discussion

A key difference between colloidal aggregation via nonspecific interactions and directed aggregation of bioconjugated nanoparticles mediated by free analytes in solution, is that the latter is limited by valence (stoichiometry) and bond competition. The total number of available receptors per unit volume is f̄c_p/2 where the factor 2 reflects two receptors per bond. The number of secondary (IONP-analyte-IONP) bonds depends non-trivially on the particle and analyte concentrations (c_p, c_a) and the IONP average functionality f̄. Yet, an upper bound for the density of bonds can be easily derived. Forming one bond requires one receptor–analyte complex (loaded receptor) encountering a free receptor. At low analyte concentration (c_a < c_pf̄/2) the number of bonds per unit volume increases linearly with the analytes c_bond = c_a. The number of secondary bonds reaches a maximum value at the stoichiometric or balance mixture . Within this range the probability of forming a bond is given by,so the maximum value of p corresponds to 2/f̄ and takes place in the balanced mixture (all the receptors may form a bond). In the range receptors start to saturate and the number of available bonds starts to decrease as c_bond = c_pf̄ − c_a, in this range p = 2[1 − c_a/(f̄c_p)]. Finally for c_a > f̄c_p all the receptors are filled with analyte and no bond can be formed. While this trend roughly determines the size of the nanoclusters, it is altered by several collective deviatoric effects: notably the presence of non-accessible receptors shadowed by neighboring particles (steric hindrance) and the formation of non-specific “fake” bonds between IONP not involving receptors (Fig. 1e). We first validate our theoretical approach against experimental results for the average size and polydispersity of the IONP nanoclusters. Then, we analyze in detail the structure of these aggregates and the (experimentally observed) gelation transitions. In doing so, we show that Bethe's theory can be satisfactorily generalized to explain the cluster distributions and percolation transitions for different f̄.

3.1 Nanoclusters size and dispersion: comparison with experiments

Experimental and simulation results for the average size (cluster hydrodynamic diameter D_H) and polydispersity index (PDI) of the IONP clusters as a function of analyte and IONP concentrations are illustrated in Fig. 2 for c_p = 1.27 μM (i.e. 7.64 × 10¹⁴ IONP mL⁻¹). In simulations, the hydrodynamic diameter of a cluster of mass s is obtained from ref. 58, (average over independent cluster configurations). Notably, for all the functionalities f̄ considered, simulation estimates for the average 〈D_H〉 over the cluster ensemble (we gather about a tens of thousand of individual clusters) and the corresponding polydispersity index PDI = Std[D_H]/〈D_H〉 are in quantitative agreement with the experimental measurements. These quantities are related to the first and second moment of the cluster distribution and it is convenient to analyze them concurrently for increasing analyte concentration c_a. As shown in Fig. 2, the case f̄ = 2 (two receptors per particle in average) perfectly follows the stoichiometric trend for the bond formation, explained above. In particular, the maximum average cluster size and the maximum polydispersity take place at with . For most receptors are loaded with analytes and secondary bond formation starts to saturate, leading to smaller clusters with smaller PDI.

Fig. 2 Mean hydrodynamic size (left column) and polydispersity index (PDI, right column) of the IONP nanoclusters as a function of the concentration of analyte c_a when the average number of receptors per nanoparticle is (a and b) f̄ = 2; (c and d) 4; (e and f) 8 and (g and h) f̄ = 12. The concentration of nanoparticles in all cases is c_p = 1.27 μM. Star symbols indicate precipitation (pellet formation) in the experiments and the formation of a percolated network (gel phase) in simulations. The shadowed regions indicate the formation of a percolating structure (gel) in simulations. For c_a > c^(g)_a ≈ 2 μM the suspension precipitates in the experiments so it makes no sense to compare with simulations.

At small analyte concentration, the case f̄ = 4 (where experimental results are available) shares features with f̄ = 2. Simulations perfectly capture the increase in average cluster size and PDI in experiments. However, for f̄ = 4 (in fact for any f̄ > 2) simulations predict the formation of a percolating cluster (corresponding to gel phase). For f̄ = 4 the onset of gelation takes place in the range 1.5 μM < c^g_a ≤ 1.75 μM and the gel phase persists up to c_a ∈ [3.0,3.5) μM. The observation of the gel phase is indicated in Fig. 2 with larger (black) star-shaped symbols and a shadowed region. The onset of the gel phase predicted by simulations coincides quite precisely with the formation of a pellet in the experiments (i.e. a solid-like phase of condensed IONPs which precipitates). This experimental event is indicated by red stars in Fig. 2. We note that in the experiments, pellet precipitation was observed for c_a ≥ 2.0 μM; hence in this range it is not possible to compare with simulations.

For f̄ = 4 the stoichiometric balance (for which the maximum number of available bonds is attained) takes place at . Unlike the f̄ = 2 case, the system gelates before the stoichiometric condition is attained as a consequence, a peak in average cluster size at is not to be observed. Above the transition, the fraction of particles composing the percolating cluster is 0.6. Thus, close to the gel transition, a high enough number of particles remain unbounded or forming finite clusters and the average cluster size does not significantly increase for c_a > c^g_a (Fig. 2c). By contrast, this effect creates a large cluster polydispersity and the PDI reaches a maximum value at (Fig. 2d). For the number of available bonds decreases, and so it does the PDI, as the infinite cluster grows up to 95% of the particles. Interestingly, for c_a ≥ 3.5 μM most receptors are saturated with analytes and simulations predict re-entrant melting (reverse gel–sol transition) as indicated in Fig. 2c (the gel phase is indicated with a shadowed area). Indeed, as the available secondary bonds scarce, the scenario becomes similar to . A similar re-entrant network formation (without phase separation) as a function of linker concentration and centered at the stoichiometric ratio, was recently observed in simulations of linker-mediated aggregation of colloids.¹⁷ This re-entrant transition is however not always easily observed in experiments. As shown in Fig. 2c, under experimental conditions, aggregates are still large enough to form precipitates at c_a ≥ 3.5 μM. A plausible hypothesis for such a deviation from simulations is the formation of non-specific bonds, whereby analytes physisorb to the dextran-coated IONP surfaces acting as bridges for another IONP and creating spurious bonds, as illustrated in Fig. 1e. These non-specific (spurious) bonds become more frequent as c_a surpasses the corresponding dissociation constant for non-specific protein-dextran adhesion, and when specific interacting partners are not available. To validate this hypothesis, experiments were repeated without receptors (f = 0) and fixed c_p: indeed, the formation of aggregates was observed for large enough values of c_a, confirming the presence of non-specific interactions in absence of receptors. We note that non-specific adhesion can be easily included in the theoretical model. However, we deliberately exclude this effect in the theory to validate, by comparison, the region where specific ligand–analyte bonds are dominant in experiments.

The trends for f̄ = 8 and 12 (Fig. 2e–h) share most of the features of f̄ = 4. Simulations and experiments observe the formation of a gel phase (percolating in simulations and precipitating in experiments) at a slightly larger analyte concentration c^g_a ∈ (1.75–2.0] μM. For c_a < c^g_a simulations quantitatively agree with experimental measurements of the average cluster size and PDI. In these cases, the stoichiometric balance would take place at high analyte concentration (respectively and 7.62 μM). But, the percolating cluster grows as c_a is increased above c^g_a and well before , the infinite cluster ends up by gathering the total number of particles in the simulation box, leading to a vanishing PDI at large enough c_a.

To further validate our theoretical and computational models for cluster formation, Fig. 3 compares two assays with fixed analyte concentration (c_a = 0.15 μM and 0.75 μM) and varying IONP concentration (average valence f̄ = 4). As expected, a rather sharp peak in the average cluster size is observed at the stoichiometric mixture . For c_a = 0.15 μM the peak takes place around 0.45 × 10¹⁴ IONP mL⁻¹ (c_p = 0.075 μM) while for c_a = 0.75 μM the maximum cluster size is found when 2.25 × 10¹⁴ IONP mL⁻¹ (or c_p = 0.37 μM). Experimental results are in good agreement with simulations, although small deviations are noticed for the c_a = 0.75 μM case. As an aside, it has to be stressed that details of these assemblies are rather sensitive to each parameter of the system. For instance, simulations showed deviations of about 15% in the cluster average size if the same functionality was used for all IONPs (f = f̄) instead of using eqn (2). More examples are provided below. As a relevant conclusion in Fig. 3, it is possible to maximize the average cluster size by working close to the optimal ratio c_p/c_a = 2/f̄. However, we shall later show that the system percolates when p = f̄c_p/c_a > p_c, where the threshold p_c is derived later (in eqn (20)). Thus working too close to p = 1 will lead to pellet formation, which is undesirable for biosensor applications. In particular, in Fig. 3 simulations with c_a = 0.75 μM and p ≈ 1, produced pre-percolated structures with large clusters (in some cases with 400 NPs). In the dilute case, c_a = 0.15 μM, and for p ≈ 1, we observed clusters with up to a 20% of the total number of NPs. While Fig. 3 should be a master curve (as even in dilute suspensions one should expect the formation of percolated structures around p ≃ 1) waiting times becomes prohibitively long as the concentration is reduced (this could actually be a benefit for biosensor experiments).

Fig. 3 Mean hydrodynamic size of the aggregates for increasing IONP concentration c_p and two values of analytes concentrations (c_a = 0.15 μM and 0.75 μM). The mean number of receptors per particle is f̄ = 4. The x-axis is normalized, clearly showing that the peak in cluster size takes place for the balanced mixture, c_pf̄/(2c_a) = 1.

3.2 Steric hindrance and accessible links

After validating our theoretical and computational approach against experiments we analyze the statistics of receptors and bonds obtained from simulations. To start with, we consider the effect of steric hindrance in limiting the number of available receptors. This effect is illustrated in Fig. 4, where it is seen that the area described by a solid angle θ_ex around one bond linking two particles (shadowed region in the bottom panel in Fig. 4a), cannot be reached by the receptor of a third. More precisely, for an excluded angle θ_ex, the area excluded around one IONP receptor is A_ex = A_sp(1 − cos(θ_ex))/2 ≈ θ_ex²A_sp/2 with A_sp = 4πa². Following the sketch in Fig. 4b is not difficult to derive with b = a + [small script l] and [small script l] the linker length (the bond length is 2[small script l]). In our case, the particle radius is a ≈ 32 nm and [small script l] = 2 nm, so one bond excludes a net area of A_ex = 0.11A_sp. As a simple estimation, a particle with f̄ receptors but having only f_acc of them accessible, will exclude a fraction f_accA_ex/A_sp of the particle area. Making the approximation that the ratio between accessible and total receptors equals the fraction of accessible area, f_acc/f̄ ≈ (1 − f_accA_ex/A_sp), one gets,Fig. 4a compares this estimation with results for the average 〈f_acc〉 obtained from test-particle random sampling. We generate 10⁵ configurations of a tagged bioconjugated particle having a random location and a number receptors on its surface according to the binomial distribution in eqn (1) (see Methods section), excluding f = 0 cases. To measure the maximum number of accessible receptors for each configuration f_acc, we consecutively bind neighbour particles to the tagged one until no more bonds can be formed due to receptor saturation or to steric hindrance (particle–particle overlaps). Thus, f_acc is the limit value of accessible bonds in a dense environment (e.g., the interior of clusters). If we take the limit of f̄ → ∞ in eqn (13), we find that, on average, the maximum number of bonds a particle can form is close to 9, which is larger than the close packing limit for the coordination number,⁵⁹ 6, probably due to the linkers excluded area and the finite separation between particles. Within the range of interest of f̄, the agreement with the approximation in eqn (13) is excellent, which suggests that, at least for small linker length [small script l], the excluded area can be understood as sum of two-body exclusions. The present analysis indicates the strong effect that the receptor length [small script l] might have on the number of accessible receptors f_acc. For instance, for a = 32 nm, and [small script l] → 0 nm, the fraction of excluded area per bond grows to 0.25, while for [small script l] > 13.25 nm, θ_ex = 0, so in that case no receptor will become unavailable due to steric effects (however, for long linkers -and depending on their flexibility-another way of “loosing” receptors is the formation of multiple bonds between two particles¹⁶). In any case, tuning the length of the linker might be an alternative¹⁰ way to control the effective valence of bioconjugated particles. For instance, linkers such as ssDNA strands will expand/shorten by modifying the pH or the electrolyte concentration, while proteins are sensitive to temperature or chaotropic agents (urea).

Fig. 4 (a) Average number of accessible receptors 〈f_acc〉: symbols obtained from random sampling (see text, error bars indicate the standard deviation) and the line is the theoretical estimation in eqn (13). The angle θ_ex excluded by two bonded particles to a third partner is illustrated in the inset. (b) Geometry to estimate θ_ex: using the triangles ABC and ADC we obtain sin(α) and cos(β) respectively, and then use α + β + θ_ex = 90° to derive θ_ex. The receptor length is [small script l] ≈ 2 nm and particle radius a ≈ 32 nm, providing θ_ex = 38.4°.

3.3 Bond statistics

Another important quantity is the number of bonds n_bonds, which is directly related to the average number of nearest neighbors 〈NN〉 in a cluster. Here we take the Bethe's lattice^60,61 as reference (see S.1. and Fig. S1, ESI†). Each node of the Bethe's lattice has a fixed number of neighbors (fixed valence f) and there is a fixed probability p of creating a link between two nodes. In a Betthe's lattice the number of bonds formed in a cluster with s particles is n_bonds = s − 1. This simple relation follows from the fact that cyclic structures (rings, etc.) are not allowed (Fig. S1b, ESI†). Consequently, in the Bethe lattice, on average, each particle will be connected to a number of neighbors given by,where the 2 arises from the fact that each bond connects two particles. Somewhat anti-intuitively, NN does not depends on f.

Compared with the Bethe's lattice (Fig. S1b, ESI†), nanoclusters of bioconjugated IONP will certainly present a larger number of bonds and neighbors because of the formation of cyclic structures, shown in the AFM images of Fig. 1b (see also Fig. S1c, ESI†). The deviations from the ideal Betthe's lattice trend increase with the number of bonds available and therefore with c_a (as long as ). This fact is clearly observed in Fig. 5a where the symbols represent simulation results for 〈NN〉 for increasing values of c_a without reaching the percolation of the system. Taking the analogy with the Bethe's lattice limit for s → ∞, we propose a helpful analytical relation for NN which only depends on 〈NN〉_∞; the mean number of neighbors that a particle would have in a cluster with infinite particles (note that 〈NN〉^Bethe_∞ = 2)

eqn (15) nicely recovers simulation results (see lines in Fig. 5a) and provides the fitting values of 〈NN〉_∞ (shown in Fig. 5b against f̄ for c_a = 0.5 μM, 1.0 μM and 1.5 μM). As an interesting outcome, at low bond concentration (c_a < 0.5 μM) the number of neighbours in the IONP assemblies approaches the ideal Bethe's lattice value NN → NN_Bethe. At high bond concentration and large clusters, 〈NN〉_∞ is maximum for f̄ = 3 or 4. For f̄ > 4, the number of options for the formation of a new bond increases with f̄ and this explains a decrease in NN. This gain in entropy leads to less compact clusters for f̄ > 4, as we reveal below by analyzing the Flory exponent.

Fig. 5 (a) Mean number of nearest neighbors that a particle has as a function of the number of particles composing the cluster for f̄ = 4 and c_p = 1.27 μM. (b) Average number of bonded particles per particle in large clusters a function of f̄.

3.4 Compactness of nanoclusters: Flory exponent

In fractal aggregates, the relation between the cluster size R_g(s) and the number of particles it contains s tends to a power law R_g ∼ s^ν for large enough s. The fractal dimension d_F = 1/ν is related to the Flory exponent ν which can be extracted from the scaling of the cluster gyration radius⁶² where r_ij is the distance between pairs of cluster particles. We gathered enough statistics from simulations to examine subtle details in the Flory exponent which provide information on the shape and the kinetics of the nanoclusters. The cluster size and shape have a strong influence on the cluster's transduction signal (e.g. its magnetic response³³) and for the applied sake, it is interesting to study how ν varies with the number of particles in the cluster, s. To that end we calculate a “local” exponent as,To perform the derivative in eqn (16) we first reduce the numerical noise by fitting R_g(s) with the function,

R_g(s) = −c₁s^−c₂ + c₃s^ν_∞ (17)

where c_i are non-negative fitting parameters and ν_∞ is the Flory exponent in the limit of an infinite cluster. Values of c_i and ν_∞ are given in the ESI.†Fig. 6 shows results for ν(s) for different values of f̄. The analysis was performed at c_a = 1.5 μM and c_p = 1.27 μM, which are values close to (slightly below) the gelation transition. In agreement with previous studies, we observe that the fractal dimension of clusters d_f = 1/ν increases with the number of particles.⁶¹ The limiting values of the Flory exponent for s → ∞ (which apply for about s > 10²) are plotted as a function of f̄ in Fig. 6c (see also Table S1 in ESI†).

Fig. 6 (a) An example of the scaling of the gyration radius of clusters against the number of particles in the cluster, illustrating the fit with eqn (17) (red line, details in ESI†). (b) The “local” Flory coefficient ν(s) defined in eqn (16). (c) The limiting value of the Flory exponent ν_∞ ≡ ν(s → ∞) against the average functionality f̄.

The variation of ν_∞ with f̄ is not large, yet, it reflects relevant aspects of cluster formation. Notably, values of ν_∞ lie between that of reaction limited cluster aggregation (RLCA) ν_RLCA = 1/2.1 ≈ 0.476 and a Gaussian freely jointed chain (FJC) ν_FJC = 0.5.

Clusters obtained for f̄ = 2 are Gaussian-like (ν_∞ ≃ 0.5). If the polymerization were linear, this would correspond to new bonds completely uncorrelated with the previous bond direction. However, according to the binomial distribution in eqn (1), for f̄ = 2, a percentage of the IONPs present f > 2, which leads to the formation of branches. Another subset presents f = 1, leading to dead-ends in some of the branches. To investigate the effect of branches and dead-ends, we performed simulations with a single functionality f = f̄ = 2 for all IONPs. In such a case we obtained an exponent corresponding to an expanded chain (self-avoiding walk) ν ≈ 0.6. Therefore, for a polydisperse functionality with f̄ = 2, the origin of more compact assemblies (ν ≈ 0.5) is the formation of branches f ≥ 3 and dead-ends f = 1. For f̄ = 3 and 4, we found even more compact structures with ν_∞ = 0.473, quite close to the RLCA value 0.476. RLCA corresponds to a bond formation which is hampered by a repulsive barrier, so that the probability of forming bond of two “colliding” particles is smaller than one. In the present scenario, such a barrier is purely entropic, corresponding to a smaller number of available spots to create a new bond in these clusters. As shown in Fig. 5, for f̄ = 3 (and 4), the clusters present the largest number of neighbors (e.g. 〈NN〉_∞ ≈ 2.4 for c_a = 1.5 μM). Remarkably, this number is quite close to the average accessible receptors per particle 〈f_acc〉 (Fig. 4). This implies that the chance to find a new available site in clusters with f̄ = 3 or 4, is smaller than for larger f̄. For instance, for f̄ = 12, a new particle joining a cluster may encounter up to 〈f_acc〉 ≈ 5 available receptors, while particles in a cluster only typically have up to 〈NN〉_∞ ≈ 2.25 neighbors (Fig. 5b). At very large f̄, the particle corona is filled with receptors and, for large bonding probability p ∼ 1, one expects to approach the diffusion limited (DLCA) scenario (colliding particles form a bond with probability one, leading to more open structures ν_DLCA = 1/1.8 ≈ 0.56). Consistently with this line, the value of ν_∞ gradually increases for f̄ > 4, reaching 0.493 for f̄ = 12 (and p = 0.2), as shown in Fig. 6.

3.5 Distribution of IONP clusters

A more detailed level of description concerns n(s), the distribution of clusters of mass s (S.3, eqn (S20), ESI†). To analyze n(s) we deploy the analytical predictions for a Bethe lattice⁶¹ (summarized in the S.1, ESI†). A Bethe's lattice is characterized by particles having an equal number of linkers f and a probability p of forming a link. Bethe's theory derives an analytical relation for the probability density of forming clusters with s particles,

n(s) = g_sp^s−1(1 − p)^2+s(f−2) (18)

where g_s is the number of combinations of placing s_f = 2 + s(f − 2) nodes in the cluster. For large s, it can be approximated using Stirling's relation. Above a certain critical binding probability p_c = 1/(f − 1) this theory predicts a transition to a gel phase (infinite cluster). Note that for f = 2 this corresponds to p_c → 1. Far from the critical percolation threshold, p < p_c, the distribution is exponential n(s) ∝ exp[−s/ξ(p)], where ξ(p) is the typical (average) number of particles in a cluster (see S.1, ESI†). Close to the percolation threshold p ≈ p_c, a Taylor expansion of ξ(p) around p_c indicates that for f > 2, the cluster typical size diverges as ξ(p) ∝ |p − p_c|⁻². For f = 2, the scaling is different ξ(p) ∝ |p − p_c|⁻¹. Close to gelation, Bethe also predicts the transition to a power-law dependence on the cluster size distribution, n(s) ∝ s^−τ (with τ = 5/2 = 2.5 for the Bethe lattice), which is a landmark of sol–gel phase transition. This critical scaling was first predicted by Fisher,⁶² and the exponent τ bears his name.

In order to use Bethe's theory as a template for our system, we fix the probability of forming a bond p, defined in eqn (12). Although steric hindrance limits the number of available receptors, these unusable linkers can also trap analytes, reducing the number of bond-forming analytes to a fraction, c^(acc)_a = c_a〈f_acc〉/f̄. Therefore, consideration of steric hindrance leads to the same binding probability as eqn (12), p = 2c^(acc)_a/(〈f_acc〉c_p).

In our system, deviations from Bethe's lattice are certainly expected. First, particles have a different functionalities f, determined by a binomial distribution with mean f̄. Secondly, structures with cyclic bonds are formed, and become significant as the available number of bonds increases. Finally, due to steric hindrance the average functionality is gradually reduced with the local particle and bond concentration, becoming 〈f_acc〉 < f̄ (Fig. 4). In order to compare with simulation results for n(s), we will abuse Bethe's theory and insert a real number as the “effective valence” f_eff in the analytical relation for n(s), eqn (18) (we note however that p is fixed, by eqn (12)). As p increases, cyclic-bonding structures become more frequent, as revealed by a significant deviation in the mean number of neighbors per particle, with respect Bethe's lattice (see 〈NN〉 in Fig. 5). At high bond concentration, steric hindrance gradually leads to the limit of accessible receptors. Therefore the effective functionality f_eff (fitting n(s) viaeqn (18) needs to decrease with c_a, within the range 〈f_acc〉 ≤ f_eff ≤ f̄. A way to estimate f_eff is to combine the information about the fractal dimension of the clusters R_g(s) ∼ s^ν and the cluster size distribution n(s;f_eff,p) in eqn (18): the effective functionality f_eff should recover the experimental trends for the average cluster size: , where P_s(s) = sn(s)/C is the probability of a particle being in a cluster of size s and . Numerical solution of this inverse problem leads to values of f_eff which are consistent with those derived from direct fitting of n(s) (see below). In particular, f_eff decreases with c_a, going from f̄ to 〈f_acc〉. Note that this procedure is also a fast and useful way to approximate n(s) from very limited information (e.g. from the average cluster size), an idea already used in ref. 63. However, we need to verify that the shape of n(s) is indeed recovered by eqn (18). To that end, we calculate n(s) from BD simulations providing a direct and more accurate estimation of f_eff.

3.5.1 Cluster distributions for f̄ = 2. The case f̄ = 2 is qualitatively different from f̄ > 2 and deserves a separate discussion. Fig. 7a–d shows simulation results for n(s) for different concentration of analytes: from low concentration c_a = 0.1 μM up to the balanced mixture c_a = c_p = 1.27 μM. Even at c_a = 0.1 μM (low binding probability) the IONPs cluster distribution for f̄ = 2 differ from the Bethe lattice with functionality f = 2 (blue lines), presenting a larger population of large clusters and a deviation from the exponential trend. For f̄ = 2, these deviations from the Bethe's scenario are due to the dispersion in the receptor number f. We confirmed this statement upon comparison with simulations of mono-disperse case (f = 2 for all IONPs) which showed n(s) ∼ exp[−s/ξ(p)] with ξ(p) in quite good agreement with Bethe's prediction up to c_a ≈ 1 μM (however, percolation was not observed due to loop formation). According to the binomial distribution of receptors P(f = 0) = e^−f̄ (see eqn (1), and inset of Fig. S4, ESI†) which for f̄ = 2, means that a significant fraction (13%) of particles have no linker f = 0, and will never polymerize. Excluding these impaired particles (f = 0) from the ensemble of f̄ = 2, leads to an average of f̄^(b) = 2.31 receptors per binding-viable particles (f > 0). Note that this effect is already negligible for f̄ = 4 as, f̄^(b)(4) = 4.07. One thus expects that f_eff (2) > 2, even at low c_a, which is consistent with the deviation of n(s) from the single exponential (Fig. 7a). Notably, the generalized Bethe relation with a real effective valence f_eff (red lines) nicely matches all the observed distributions. As expected, the values of f_eff obtained from the best fit to eqn (18) (Fig. 7a–d) decrease as c_a approaches ; being f_eff ≈ f̄^(b) at moderate c_a, and approaching the high cluster density limit 〈f_acc〉 ≃ 1.9 close to .

Fig. 7 Cluster distributions obtained for IONP concentration c_p = 1.27 μM and average receptors per IONP equal to f̄ = 2 (a–d) and f̄ = 4 (e–h). Symbols correspond to simulation results, blue lines to the Bethe's lattice (BL) with single functionality f = f̄ and red lines are the best fits using Bethe's analytical relation eqn (18) with an effective functionality f_eff. For f̄ = 4, results for c_a = 1.5 μM correspond to the onset of percolation.

For the balanced mixture the binding probability is p = p_c = 1 and Bethe's theory with f = 2, predicts the formation of a large unique chain in an infinite time.^22,23 However, for f̄ = 2 the bioconjugated IONPs reached an equilibrium state with no infinite clusters. Consistently with the absence of percolation transition, Fig. 7d shows that the distribution n(s) obtained for presents an exponential decay for s > 50. As verified by comparison with the mono-disperse f = 2, for f̄ = 2 we find formation of branches (particles with f > 2) and to a lesser extent loops. However, it should be noted that in the case of mixtures of patchy colloids (where saturation of available bonds is not a limiting factor, as it happens in linker-mediated aggregation), an infinitesimal fraction of f = 3 particles leads to percolation.²⁴

3.5.2 Cluster distributions for f̄ = 4, 8, 12. Fig. 7e–h shows simulation results for the distribution n(s) obtained from f̄ = 4 and compare them with the best fit to the generalized Bethe's relation eqn (18). For f̄ > 2 the effect of receptor poly-dispersion becomes less important (e.g. f̄ ≃ f̄^(b)). And in this line, for f̄ = 4 we find that when the analyte concentration is small (c_a = 0.1 μM, binding probability p ≈ 0.04), the cluster distribution is quite close to that produced by a Bethe lattice with f = f̄, with n(s) being a single exponential. For c_a = 0.5 μM (or p ≈ 0.2) steric hindrance and cyclic structures become relevant and n(s) starts to deviate from the Bethe lattice with f = 4. An excellent fit to n(s) is obtained using a smaller effective functionality f_eff = 3.55 into eqn (18). As c_a is further increased and larger clusters are formed, steric hindrance tends to further reduce f_eff, becoming f_eff ≃ 3.1 for c_a = 1 μM (p = 0.4), which is close to the high-concentration accessible limit 〈f_acc〉 ≈ 2.8. A similar outcome is found for the distributions for f̄ = {8, 12} (ESI†) which can be also mapped to the effective Bethe's distribution. As f̄ becomes larger, however, the formation of loops becomes more probable and the effective functionality f_eff starts to deviate from f̄ at smaller values of c_a.

Fig. 8 collects the results for the effective functionality against f̄, showing its upper f̄^(b) and lower 〈f_acc〉 boundaries. Values of f_eff for small analyte concentration c_a = 0.1 μM (purple diamonds) approach the Bethe's scenario, where loops are unlikely to be formed and most of the receptors are accessible. At c_a = 1.0 μM (green diamonds) loop formation and the steric hindrance within the larger clusters induce a decrease of f_eff to values much closer to the accessible limit f_acc.

Fig. 8 Variation of the effective functionality f_eff of the particles with the average number of receptors f̄ for two values of the analyte concentration c_a = 0.1 μM and 1.0 μM. Note that f_eff decreases with c_a: at low c_a it is close to the average number of receptors of binding-viable particles f̄^(b) (i.e. the ensemble with f > 0 receptors) while for large c_a, it approaches the high density limit of accessible receptors, f_eff ∼ f_acc (the shadowed region indicates the standard deviation of f_acc).

3.5.3 Percolation transition. While the observed cluster distribution n(s) are not generally pure exponential functions, for c_a < c^(g)_a (i.e. before percolation), they all present an exponential tail. Close to the sol–gel transition (observed for f̄ ≥ 4), the exponential tail becomes distorted and at the percolation threshold, n(s) fully exhibits a power law scaling. This transition is illustrated in Fig. 7h. Concerning the Fisher exponent τ governing the critical scaling n(s) ∼ s^−τ, we find τ = 1.75 (green line in Fig. 7h. This value is smaller than the Fisher exponent for the Bethes lattice, which (S.1, ESI†) turns out to be independent on the valence (τ_Bethe = 5/2). When compared with Bethes picture, our nanoclusters present a broader distribution close to the gel transition. This is expected considering that the possibility of forming multiple bonds between particles within the same cluster (rings, polyhedrons, etc.) enhances the diversity in cluster sizes, with respect the Bethe model.

It is possible to estimate the locus for the transition to gel phase by generalizing Bethe's prediction for the critical probability using the functionality at the accessible limit, p_c = (〈f_acc〉 − 1)⁻¹ = 2c^(g)_a/(f̄c_p). This relation leads to,

In the case of Fig. 3 (where f̄ = 4, 〈f_acc〉 ≈ 2.5 and c_p is varied at fixed c_a) we consistently observed formation of large structures only for p > p_c = 0.66. In terms of c_a (and fixed c_p = 1.27 μM) for f̄ = 4, 8 and 12 relation 20 provides the direct (sol–gel) transition at c^(g)_a = 1.41, 1.70 and 1.91 μM, respectively. These predictions are in quite good agreement with simulations and experiments. As shown in Fig. 2, the sol–gel transition takes place within the range c^g_a ∈ (1.5–1.75] μM for f̄ = 4 and 8, while for f̄ = 12 it lies within c^g_a ∈ (1.75–2.0] μM. Also, according to the prediction p_c ≈ 0.66 for f̄ = 4, the reverse gel–sol transition should take place for c_a ≈ 3.4 μM which is in the range of Fig. 2. In line with the prediction of eqn (20), as f̄ is increased, the sol–gel transition takes places at slightly larger values of the analyte concentration. One possible reason for the small discrepancies observed between the predictions of eqn (20) and the simulation results could be the formation of loops, as discussed by Lindquist et al. for linker-mediated aggregation.²⁷ These loops hinder the percolation of the system, thus requiring a higher concentration of analytes to reach the percolation transition.

To conclude, it is interesting to compare our aggregates with other linker-mediated aggregates. As an example, for f̄ = 4, we observe the gelation transition for IONP packing fraction of ϕ = c_p[m⁻³](4πa³/3) ≈ 0.085 with f_acc ≃ 2.8 and inverse binding energy k_BT/E_bind ≈ 0.08. In the case of the so called “empty fluids”, formed by open percolating structures formed by limited valence colloids²⁹ these values are quite consistent with the “empty fluid” phase of limited valence colloids.²⁹ In the case of linker-mediated colloids²⁷ where linkers and particles have the same size, for a volume fraction 0.09, the percolation thresholds were found around k_BT/E_bind ≈ 0.11. In the case of polymer-linked colloids,¹⁵ the spinodal line obtained using dumbbells as linkers and patchy colloids with f = 6 and orthogonally placed receptors, is indeed very close to our result k_BT/E_bind ≈ 0.06 for colloid volume fraction around 0.08.

4 Conclusions

We reported on an experimental and theoretical study on predicting the formation of nanoparticle clustering driven by biomolecular recognition. These clusters strongly depend on molecular details, such as the size and binding energy of the receptor–analyte interaction, as well as the average number of receptors per particle f̄, and the concentration of nanoparticles c_p and analytes c_a. Iron oxide nanoparticles displaying superparamagnetic behaviour were decorated with molecular receptors which specifically bind to the target molecule (analyte protein) with selectivity and affinity (a dissociation constant of about 1 μM). We proposed a theoretical and computational approach and successfully validate it against experimental results for the average size and polydispersity index (PDI) of the so-formed IONP nanoclusters. The theoretical route consists in three separate steps: (i) distribute the receptors amongst nanoparticles with a fixed (experimentally controlled) average valence f̄, (ii) distribute the analytes amongst the receptors according to the given set of parameters (c_p, c_a and f̄) and (iii) perform Brownian dynamics simulations for the formation of bonds between IONPs, in the limit of irreversible binding. The IONPs are considered spherical particles with excluded volume interactions, having translational and rotational degrees of freedom, while receptor–analyte–receptor bonds consists on harmonic springs.

Experimental results for the mean hydrodynamic size and PDI of clusters are reproduced with considerable accuracy. For f̄ = 2, a peak in the cluster size and PDI takes place at the stoichiometric mixture . For this ratio c_a/c_p = 2/f̄, the number of possible bonds is maximum. For f̄ > 2 (and using c_p = 1.27 μM) we found that the system gelates at a smaller analyte concentration . Notably, the numerical predicted sol–gel transition precisely coincides with the precipitation of a pellet phase in experiments. In the second part of this work we have used Bethe's lattice theory as a template for an analysis on the statistics of the IONP clusters. We consider deviatoric effects from Bethe's scenario, such as the effect of receptor poly-dispersity, steric hindrance (some receptors are not accessible due to the “shadow” of neighboring IONPs) and the formation of cyclic structures in the clusters. We analyzed the fractal structure of the clusters, providing a “local” Flory exponent ν(s). Interestingly, we found a close relation between the fractal dimension d_f = 1/ν, the number of accessible receptors, and the average number of neighbors per IONP. Finally, we have studied the distribution of aggregates n(s) as a function of c_a and average functionality f̄. Introducing a real effective functionality f_eff into Bethe's relation, allow us to excellently recover the measured cluster distributions. This effective valence f_eff decreases with the analyte concentration (i.e. with the binding probability) starting from f ≈ f̄ at low c_a towards f_eff ∼ 〈f_acc〉 as c_a is increased. This analysis enabled us to estimate the effect of polydispersion in receptor number (note that P(f) is a binomial distribution), and the creation of branches and cyclic structures. Notably, these two effects induce significant deviations between f̄ = 2 and the mono-valence f = 2 counterpart, even at low c_a. For f̄≥4 a sol–gel transition is observed, characterized by a power-law dependence for the cluster size distribution n_s(s) ∼ s^−τ. Generalizing Bethe's lattice prediction, we derived a relation (eqn (20)) which correctly forecast the locus of the sol–gel transition in the parameter space (c_a, c_p and f̄). Interestingly, for f̄ = 4, a re-entrant transition towards gas phase occurs as c_a is further increased, however this was not experimentally observed due to non-specific adsorption. Understanding the self-assembly of nanoparticles mediated by biomolecular recognition is essential to decipher their response under external fields, for instance the hysteresis cycles of bioconjugated magnetic nanoparticles under AC magnetometry. In a future work we will consider nanoclusters of bioconjugated ferromagnetic nanoparticles, where one expects that equilibrium structures will non-trivially depend on the interplay between the magnetic energy (dipolar interactions between particles) and the chemical energy of the bonds.

Author contributions

Pablo Palacios: conceptualization, methodology (theory and simulations), software, writing-original draft, review and editing. Elena Sand-de Diego: methodology (experiments), review and editing. Raúl Pérez Peláez: software. A. L. Cortajarena: resources (analytes), review. Francisco J. Terán: conceptualization, methodology (experiments), resources (nanoparticles), review and editing. Rafael Delgado-Buscalioni: conceptualization, methodology (theory and simulations), writing original draft, review and editing.

Conflicts of interest

The authors declare no conflict of interest.

Acknowledgements

This work has been partially funded by the Spanish Research Agencies (PID2020-117080RB-C51, PDC2021-121441-C21, PCI2019-103600,CEX2020-001039-S, PID2019-111649RB-I00 PID2020-117080RB-C53, MDM-2017-0720, RED2018-102626-T) and Comunidad de Madrid (NANOMAGCOST, S2018/NMT-4321). European COST Actions CA17115 (MyWave), and CA17140 (Nano2Clinic) are also acknowledged. ESD thanks Comunidad de Madrid for financial support (PEJ-2017-AI/BMD-7517). Authors thank Dr Patricia Pedraz for carefull acquisition of AFM images.

Notes and references

1. A. L. Hiddessen, S. D. Rodgers, D. A. Weitz and D. A. Hammer, Langmuir, 2000, 16, 9744–9753.
2. L. Zhai, A. Narkar and K. Ahn, Nano Today, 2020, 30, 100826.
3. A. R. Town, J. Taylor, K. Dawson, E. Niezabitowska, N. M. Elbaz, A. Corker, E. Garcia-Tuñón and T. O. McDonald, J. Mater. Chem. B, 2019, 7, 373–383.
4. Y. Gao, A. Ahiabu and M. Serpe, ACS Appl. Mater. Interfaces, 2014, 6, 13749–13756.
5. Y. Wang, Z. Gao, Z. Han, Y. Liu, H. Yang, T. Akkin, C. Hogan and J. Bischof, Sci. Rep., 2021, 11, 898.
6. I. Isnaeni and Y. Herbani, J. Phys.: Conf. Ser., 2017, 817, 012039.
7. N. A. Mahynski, E. Pretti, V. K. Shen and J. Mittal, Nat. Commun., 2019, 2028.
8. S. Neupane, Y. Pan, S. Takalkar, K. Bentz, J. Farmakes, Y. Xu, B. Chen, G. Liu, S. Qian and Z. Yang, J. Phys. Chem. C, 2016, 121, 1377–1386.
9. K. L. Gurunatha, A. C. Fournier, A. Urovas, M. Valerio-Lepiniec, V. Marchi, P. Minar and E. Dujardin, ACSNano, 2016, 10, 3176–3185.
10. S. Angioletti-Uberti, P. Varilly, B. M. Mognetti and D. Frenkel, Phys. Rev. Lett., 2014, 113, 128303.
11. W. Rogers, J. Chem. Phys., 2020, 153, 124901.
12. K. Wu, R. Saha, D. Su, V. D. Krishna, J. Liu, M. C.-J. Cheeran and J.-P. Wang, ACS Appl. Nano Mater., 2020, 3, 9560–9580.
13. M. Heidenreich, J. M. Georgeson, E. Locatelli and L. Rovigatti, et al. , Nat. Chem. Biol., 2020, 16, 939–945.
14. X. Xiuyang, H. Hu, M. P. Ciamara and R. Ni, Sci. Adv., 2020, 6, eaaz6921.
15. M. P. Howard, R. B. Jadrich, B. A. Lindquist, F. Khabaz, R. T. Bonnecaze, D. J. Milliron and T. M. Truskett, J. Chem. Phys., 2019, 151, 124901.
16. M. P. Howard, Z. M. Sherman, A. N. Sreenivasan, S. A. Valenzuela, E. V. Anslyn, D. J. Milliron and T. M. Truskett, J. Chem. Phys., 2021, 154, 074901.
17. T. Kwon, T. A. Wilcoxson, D. J. Milliron and T. M. Truskett, J. Chem. Phys., 2022, 157, 184902.
18. M. S. Wertheim, J. Stat. Phys., 1986, 42, 459.
19. P. Varilly, S. Angioletti-Uberti, M. M. Bortolo and D. Frenkel, J. Chem. Phys., 2012, 137, 094108.
20. P. Flory, J. Am. Chem. Soc., 1941, 63, 3083–3090.
21. W. H. Stockmayer, J. Chem. Phys., 1943, 11, 45–55.
22. R. M. Ziff and G. Stell, J. Chem. Phys., 1980, 73, 3492–3499.
23. P. G. van Dongen and M. Ernst, J. Stat. Phys., 1984, 37, 301.
24. D. de Las Heras, J. M. Tavares and M. Telo da Gama, Soft Matter, 2010, 7, 5615–5626.
25. J. M. Tavares, G. C. Antunes, C. S. Dias, M. M. Telo da Gama and N. A. M. Araújo, Soft Matter, 2020, 16, 7513–7523.
26. R. Teixeira, D. de Las Heras, J. Tavares and M. Telo da Gama, J. Chem. Phys., 2021, 155, 044903.
27. B. A. Lindquist, R. B. Jadrich, D. J. Milliron and T. M. Truskett, J. Chem. Phys., 2016, 145, 074906.
28. F. Sciortino and E. Zaccarelli, Curr. Opin. Colloid Interface Sci., 2017, 30, 90–96.
29. E. Bianchi, J. Largo, P. Tartaglia, E. Zaccarelli and F. Sciortino, Phys. Rev. Lett., 2006, 97, 168301.
30. D. Mendez-Gonzalez, M. Laurenti, A. Latorre, A. Somoza, A. Vázquez González, A. Negredo, E. López-Cabarcos, O. Calderón, S. Melle and J. Rubio Retama, ACS Appl. Mater. Interfaces, 2017, 9, 12272–12281.
31. L. Gloag, M. Mehdipour, D. Chen, R. Tilley and J. Gooding, Adv. Mater., 2019, 31, 1904385.
32. C. Pina-Coronado, Á. Martínez-Sobrino, L. Gutiérrez-Gálvez, R. Del Caño, E. Martínez-Periñán, D. García-Nieto, M. Rodríguez-Peña, M. Luna, P. Milán-Rois, M. Castellanos, M. Abreu, R. Cantón, J. C. Galán, T. Pineda, F. Pariente, Á. Somoza, T. García-Mendiola, R. Miranda and E. Lorenzo, Sens. Actuators, B, 2022, 369, 132217.
33. L. Gutiérrez, L. de la Cueva, M. R. del Árbol, E. Mazarío, S. de Bernardo, J. M. de la Fuente, M. P. Morales and G. Salas, Nanotechnology, 2019, 30, 112001.
34. W. Zhou, K. Hu, S. Kwee, L. Tang, Z. Wang, J. Xia and X. Li, Anal. Chem., 2020, 92, 2739–2747.
35. F. Ghasemi, M. Hormozi-Nezhad and M. Mahmoudi, Nanoscale, 2018, 10, 6361–6368.
36. P. Moitra, M. Alafeef, K. Dighe, M. Frieman and D. Pan, ACS Nano, 2020, 14, 7617–7627.
37. A. Latorre, C. Posch, Y. Garcimartín, S. Ortiz-Urda and A. Somoza, Chem. Commun., 2014, 50, 3018–3020.
38. M. Donolato, P. Antunes, R. Bejhed, T. Torre, F. Østerberg, M. Strömberg, M. Nilsson, M. Strømme, P. Svedlindh, M. Hansen and P. Vavassori, Anal. Chem., 2014, 87, 1622–1629.
39. A. Moyano, M. Salvador, J. Martínez-García, V. Socoliuc, L. Vekas, D. Peddis, M. Alvarez, M. Fernández, M. Rivas and M. C. Blanco-López, Anal. Bioanal. Chem., 2019, 411, 6615–6624.
40. J. Perez, F. J. Simeone, Y. Saeki, L. Josephson and R. Weissleder, J. Am. Chem. Soc., 2003, 125, 10192–10193.
41. J. M. Perez, L. Josephson, T. O'Loughlin, D. Högemann and R. Weissleder, Nat. Biotechnol., 2002, 20, 816–820.
42. K. Wu, J. Liu, D. Su, R. Saha and J.-P. Wang, ACS Appl. Mater. Interfaces, 2019, 11, 22979–22986.
43. I. Morales, R. Costo, N. Mille, J. Carrey, A. Hernando and P. de la Presa, Nanoscale Adv., 2021, 3, 5801–5812.
44. C. Iacovita, J. Hurst, G. Manfredi, P. A. Hervieux, B. Donnio, J. L. Gallani and M. V. Rastei, Nanoscale, 2020, 12, 1842–1851.
45. Y. You, S. Lim and S. Gunasekaran, ACS Appl. Nano Mater., 2020, 1900.
46. K. Wu, D. Su, R. Saha, J. Liu, V. K. Chugh and J. P. Wang, ACS Appl. Nano Mater., 2020, 3, 4972–4989.
47. M. E. Jackrel, R. Valverde and L. Regan, Protein Sci., 2009, 18, 762–774.
48. R. P. Ilagan, E. Rhoades, D. F. Gruber, H.-T. Kao, V. A. Pieribone and L. Regan, FEBS J., 2010, 277, 1967–1978.
49. C. S. Dias, C. A. Custódio, G. C. Antunes, M. M. Telo da Gama, J. F. Mano and N. A. M. Araújo, ACS Appl. Mater. Interfaces, 2020, 12, 48321–48328.
50. P. Kloeden and E. Platen, The Numerical Solution of Stochastic Differential Equations, Springer, Berlin, Heidelberg, 1992, vol. 23.
51. A. Vázquez-Quesada and R. Delgado-Buscalioni, Phys. Rev. Fluids, 2020, 5, 053301.
52. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
53. S. Delong, F. Balboa Usabiaga and A. Donev, J. Chem. Phys., 2015, 143, 144107.
54. T. Evensen, S. Naess and A. Elgsaeter, Macromol. Theory Simul., 2008, 17, 403–409.
55. W. R. Hamilton, London, Edinburgh Dublin Philos. Mag. J. Sci., 1844, 25, 10–13.
56. N. Rösch, Int. J. Quantum Chem., 1987, 32, 401.
57. R. P. Peláez, Universally Adaptable Multiscale Molecular Dynamics, https://github.com/RaulPPelaez/UAMMD.
58. J. G. Kirkwood, J. Polym. Sci., Part B: Polym. Phys., 1996, 34, 597–610.
59. G. T. Nolan and P. Kavanagh, Powder Technol., 1992, 72, 149.
60. M. Fixman, J. Chem. Phys., 1962, 36, 306–310.
61. D. Liu, W. Zhou, X. Song and Z. Qiu, Fract. Fractional, 2017, 1, 12.
62. M. E. Fisher, Phys. Physique Fizika, 1967, 3, 255–283.
63. J. Tavares, P. Teixeira, M. M. Telo da Gama and F. Sciortino, J. Chem. Phys., 2010, 132, 234502.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm00536d

‡ Here, we ascribe this percolation transition to a gel phase, although, viscoelastic characterization of the condensed phase were not performed.

---