Analyte‐Driven Clustering of Bio‐Conjugated Magnetic Nanoparticles

The dynamics of bio‐conjugated magnetic nanoparticles suspended in buffer‐saline solutions containing target proteins (i.e., analytes) is investigated numerically on a mesoscopic level. To simulate the dispersion of magnetic nanoparticles the dissipative particle dynamics is employed, which allows to study rather large systems, while still retaining important microscopic nanoparticle properties. In addition, the method is coupled to the Landau–Lifshitz–Gilbert equation, describing the dynamics of the magnetic nanocrystals within the macrospin approximation. The binding of multivalent analytes to the recognition ligands of the nanoparticles leads to the formation of clusters of magnetic nanoparticles, which in turn drastically changes the macroscopic magnetic response of the solution. Such colloidal changes are experimentally observable, allowing to explore new approaches to quantify the analyte amount. The ratio of the concentrations between the analytes (biomarkers) and the recognition ligands on the nanoparticles is found to play an important role in the formation and hydrodynamic size of the clusters. The proposed computational framework has great potential to be integrated with experimental efforts to advance the development of nanoparticle‐based biosensors for disease diagnostics and other biomedical applications.


Introduction
Early diagnosis and monitoring the progression of a disease are crucial for its successful medical treatment. In many cases the method of choice for the detection are immunoassays, which magnetorelaxometry it was shown [5,6] that relaxation amplitudes for a dispersion of biotinylated magnetic nanoparticle increase as one increases the concentration of avidin (analyte). Furthermore, saturation magnetization was observed to increase as one increases the concentration of c-reactive proteins. [8] Remanence magnetization measurements [9] have also been performed to quantify the biomarker concentration.
In case of immunomagnetic reduction assays (IMR) the AC susceptibility measurements [3,[10][11][12][13] are performed. Due to the clustering and the consequently increased rotational viscosity one generally observes a decrease in the magnetic susceptibility as one increases the concentration of the analyte. Another related quantity that similarly decreases as one increases the concentration of the analyte is the area of the dynamic hysteresis.
While there have been many experimental studies on IMR, numerical simulations, which could provide additional insight into the mechanisms involved, are so far still lacking. In this paper, we couple a mesoscopic particle-based simulation method (DPD) with the spin dynamics (SD) of the nanoparticle core. SD methods coupled to molecular dynamics (MD), have been used in studies of coupling of spin and lattice dynamics. [14,15] See ref. [16] for a symplectic algorithm based on Trotter-Suzuki decomposition. A related MD-SD hybrid method has been developed in studies of giant magnetoresistance effect of nanoparticles suspended in a liquid gel, [17][18][19] where in contrast to our work, pointparticle approximation was used. In addition we distinguish between the dynamics of the orientation of the nanoparticle and its magnetic moment, which is crucial in a certain range of amplitudes and frequencies of the external field. [20,21] The presented hybrid DPD-SD method serves as a first step toward realistic IMR simulations and other biomedical applications based on magnetic nanoparticles. [22] We lay the experimental groundwork by characterizing the iron oxide nanoparticles (IONPs) to be used in subsequent studies, that is, we measure their size distribution and magnetic response, see Methods. In this work we establish the computational model of bioconjugated nanoparticles suspended in a solvent containing target proteins. We first discuss the behavior of a single nanoparticle under nonequilibrium conditions induced by a pressure difference across the system to gain insight into its rheological properties. In the next step we study the cluster formation as a function of the ratio of the analyte and recognition ligand concentrations. Using the developed simulation technique, we aim to show the connection between the amount of analyte and the change in the structure of the nanoparticle suspension. We characterize this change with the hydrodynamic radius, on which the magnetic response strongly depends.

Results and Discussion
We first discuss the rheological properties of the solvent. Then we consider a single coated nanoparticle with superparamagnetic behavior (i.e., zero magnetization in absence of external magnetic field) suspended in the solvent. We calculate its hydrodynamic radius and study the dynamics of the nanoparticle under nonequilibrium conditions. Finally, we consider a suspension of multiple nanoparticles with complete magnetic interaction and study the process of nanoparticle clustering.

Solvent
To put things into a physical context we first examine the physical properties of the DPD solvent, that is, water. We determine the self-diffusion constant of the water beads as well as its viscosity. The viscosity is an important parameter as it influences the response of the dispersion of magnetic nanoparticles under oscillatory magnetic fields. [20,23] To calculate the self-diffusion constant we determined the mean squared displacement r i (t) (MSD) of the water beads in a simulation box with periodic boundary conditions and vol- where N w is the number of DPD water beads and r c is the DPD length scale, Equation (20). We obtain D w = (0.0507 ± 0.0002)r 2 c ∕ , with the DPD time scale Equation (24). This is equal to D w = (5.96 ± 0.02) × 10 −9 ,m 2 s −1 in physical units, which is about 2.6 times larger than experimental self-diffusion of water at 25 • . [24] The dynamic viscosity is calculated using the open boundary molecular dynamics (OBMD) algorithm, described in Methods. Here we use a simulation box with dimensions 30r c × 20r c × 20r c and buffer sizes 9r c , Figure 1. A pair of oppositely equal forces along the y axis is introduced into the left and right buffers. In this way, a shear flow is induced across the sample, Figure 2. As can be seen, the velocity profile is linear in the region of interest and flat in the buffers with relatively small error bars, except near x = 0 and x = 30r c , which is due to the low density of the particles.
We then measure the shear rate in the middle of the simulation box at different values of the imposed shear stress, Figure 3, and fit the numerical values corresponding to low shear stress xy < 1.0 k B T∕r 3 c to the equation wherė= v y x is the shear rate. We obtain = (5.23 ± 0.04) k B T∕r 3 c = (0.202 ± 0.002) mPa s, which is 4.4 times smaller than real viscosity of water at 25 • . Such low viscosity is consistent with the standard DPD model, which is based on a rather soft potential. [25] Figure 3 shows a noticeable  shear-thinning abovė≈ 0.2 −1 . Since the clustering of the nanoparticles is limited by diffusion, the typical velocities of the nanoparticles are rather low. One can therefore neglect this phenomenon.

Nanoparticle Hydrodynamic Radius
The magnetization dynamics of superparamagnetic nanoparticles dispersed in solution is determined by two relaxation processes. The first is the Néel rotation of the magnetic moments inside the nanoparticle [26] where 0 is typically between 10 −12 and 10 −8 s. [27] The second mechanism is the Brownian rotation of the nanoparticle, de- scribed by the Brownian relaxation time [28] which depends strongly on the hydrodynamic radius R h of the nanoparticle. One can change this radius by modifying the nanoparticle surface or by clustering. In this paper, we focus on the latter. The hydrodynamic radius can be estimated using several different approaches. One of these is based on the Kirkwood formula, which relates the inverse of the hydrodynamic radius to the thermal averages of the inverse distances r ij between the beads of the particle where N is the number of beads of the conjugated nanoparticle, Figure 4. As expected, the hydrodynamic radius increases with the nanoparticle core radius r p and the polyethylene glycol (PEG) chain length n l , Figure 5. The influence of the latter is more visible at lower particle radii, due to the lower total number of nanoparticle beads for smaller nanoparticle sizes. An alternative way of calculating the hydrodynamic radius is through the Einstein-Stokes equation, which relates the particle self-diffusion constant D md and the solvent viscosity where r i (t) is the distance of the ith nanoparticle from its initial position r i (0) and N p the number of the nanoparticles.  Equation (6) holds for infinitely large systems. For a finite size cubic simulation box with periodic boundary conditions the correction reads [29,30] where D is the corrected self-diffusion constant, L is the box size and is the analogue of the Madelung constant and is equal to ≈ 2.837298. For L = 60r c this correction yields 0.0005 r 2 c ∕ . The dependence of MSD is shown in Figure 6. A system of N p = 30 nanoparticles in a box with dimensions 60r c × 60r c × 60r c was simulated for t = 1500 . Averaging and error estimation was performed on blocks of length 50 . Due to the small number of particles, a rather large error of the MSD curve at high lag times is observed, see Figure 6. Using Equations (6) and (7), we get for the hydrodynamic radii R h = (2.92 ± 0.12) r c , R h = (3.96 ± 0.16) r c , R h = (5.64 ± 0.24) r c , which are larger than what Equation (5) gives. This could be due to the hydrodynamic interactions between the nanoparticles. Assuming an average box size L∕N 1∕3 p for a single particle, the finite-size correction becomes larger and equals 0.0015 r 2 c ∕ . The corrected hydrodynamic radii are then a lot closer to those given by Equation (5):

Non-Equilibrium Conditions
To study the behavior of the nanoparticle under non-equilibrium conditions, we impose a constant flow along the x axis. We can achieve such a flow by applying, in addition to the ambient pressure force, a uniformly distributed force on the left buffer (Figure 1) pointing toward the region of interest. Throughout this numerical experiment we fix the position of the nanoparticle in the middle of the box to prevent it from leaving the simulation box. Experimentally one could achieve this by applying a sufficiently strong inhomogeneous magnetic field.
The flow is approximately constant, except in the center near the nanoparticle, Figure 7. The surfactant chains on the nanoparticle surface are bent toward the flow direction, Figure 8.
From the induced flow around the nanoparticle, one can in principle determine the hydrodynamic radius. [31] More precisely, the tangential flow velocity as a function of distance perpendicular to the flow direction is given by [32] v t (r) where v s is the slip velocity, v ∞ is the velocity far away from the nanoparticle and Re is the Reynolds number. Equation (8) includes an inertial correction given by Oseen in 1910. Another hallmark of the finite Re is a slightly asymmetric velocity profile as can be seen in Figure 7. All these parameters can be seen as fit parameters for the calculated tangential velocity, Figure 8. Shape of nanoparticle coating when exposed to flow along the x axis. The PEG chains (blue beads) are visibly bent in the flow direction, in contrast to dextran (yellow beads).  . We used the previously determined hydrodynamic radius Equation (5) to determine the slip velocity at the surface of the nanoparticle. As can be seen from this figure, there is a rather thick layer around the nanoparticle with zero tangential velocity. This corresponds to a depletion layer with a thickness of about r c with almost zero density of water beads (not shown) around the nanoparticle. This is caused by the relative high density of the nanoparticle beads as well as the dextran layer. One can also see that for our choice of model parameters the no-slip condition is satisfied.
In passing, we comment on the question of slip or no-slip condition at the boundary of the nanoparticle and its physical relevance. It has been shown that the interaction between the solvent and the nanostructures and the corresponding boundary conditions play an important role on the nanoscale, such as an enhancement of flow in carbon nanotubes [33,34] as well as tunability of the dielectric response of water. [35] The boundary condition at the nanoparticle surface can strongly influence the behavior of magnetic nanoparticle suspensions under certain conditions (e.g., frequency and amplitude of the external magnetic field). For a slip-corrected solution, [36] the Brownian relaxation time Equation (4) should read where is the slip length, which can vary from 0 for no-slip flow and up to infinity for perfect slip. For a large slip length, this expression therefore yields a very small relaxation time. Given that the overall relaxation dynamics of the magnetization is determined by the faster process (either Brownian or Néel), [37] this fact could be experimentally realized by, for example, modifying the nanoparticle surface.

Nanoparticle Cluster Formation
We now consider a suspension of multiple nanoparticles with complete magnetic interaction and magnetization dynamics treated by the Landau-Lifshitz-Gilbert equation (details in Methods). This empowers us to study the process of nanoparticle clustering, which is central to the detection of analyte amount through magnetic measurements.
The key feature of the model that enables clustering is the dynamic binding of biomarkers to the recognition ligands. Clustering, in turn, alters the rheological properties and the response to the external magnetic field. In this work, we set the valency of the biomarkers to 4. This means that a single biomarker can bind to four different recognition ligands, not necessarily from different nanoparticles. On the other hand, a single recognition ligand can bind to only one biomarker.
In the present study, we simulated N p = 20 coated nanoparticles with nanoparticle radius r p = 1r c in a simulation box of dimensions 20r c × 20r c × 20r c and periodic boundary conditions. This yields a volume concentration of about 1%.
The time dependences of the number of unbound biomarkers ( Figure 10) and nanoparticles not bound to other nanoparticles ( Figure 11) show a quick decrease at small times and a  slow decrease at large times. The binding dynamics and the resulting cluster configurations depend on the ratio c r = c bio ∕c lig of the concentrations of biomarkers c bio and recognition ligands c lig . For low c r , one expects mostly singlets (unbound nanoparticles) and doublets (pairs of bound nanoparticles). For larger c r , triplets, quadruplets, … start to form and the number of doublets decreases, Figure 12.
Moreover, we find that the concentration of singlets decreases with increasing c r until c r ≈ 1, where there are no more unbound nanoparticles, Figure 13. If c r is further increased, the number of singlets increases again. Namely, if there are too many biomarkers, by the time the nanoparticles have diffused to each other, the biomarkers will have occupied all available sites. Nanoparticles with all recognition ligands occupied cannot bind.
We studied the clustering of nanoparticles in absence of external magnetic fields. The simulations at two different values of the nanoparticle magnetic moment, corresponding to dipolar interaction strengths = 0.006 and = 60, Equation (34), do not show any major differences. The characteristic time of the clustering process can be roughly estimated by considering the times required for two nanoparticles to approach each other (t tr ) and the time for the nanoparticles to rotate sufficiently (t rot ) for binding of a biomarker, bound to one of the nanoparticles, to a free recognition ligand on another nanoparticle.
The characteristic time t tr is given by the translational diffusion and the concentration of the nanoparticles [38] 1∕3 . As can be seen from Equation (10), the clustering is faster (i.e., t r is smaller) when using solvent with lower viscosity . In addition to a different magnetic response, another physical consequence of lower viscosity is a larger translational self-diffusion coefficient of the nanoparticles as well as that of the biomarkers. We note that the viscosity of our DPD solvent can be matched to the experimental viscosity of water by using an extended DPD thermostat, see ref. [39], which will be considered in future work. The characteristic time t rot is much more difficult to estimate, since it depends on the number of ligands per nanoparticle N l , the distribution of the number of bound biomarkers per nanoparticle as well as the radius of the nanoparticle. In the limits c r N l ≲ 1 (multiple biomarker bonds per nanoparticle improbable) and c r ≳ 1 (more than one free ligand per nanoparticle improbable) it can be estimated from the rotational diffusion constant alone: For R h = 2 r c we get t tr ≈ 10 4 = 68 ns and t rot ≈ 1050 = 7 ns. In this regime of small nanoparticles the clustering is therefore determined by translational diffusion. As expected, the hydrodynamic radius of clusters also depends strongly on the concentration ratio c r , Figure 14. It increases with increasing c r until about c r ≈ 1, where it begins to decrease due to the over-saturation of biomarkers. This is a well-known  phenomenon in nanoparticle-based immunoassays and is commonly referred to as the hook effect [40] or the postzone (in the case of antigen excess) and prozone effect (in the case of recognition ligand excess). Such an effect may lead to false-negative results.
Increasing the magnetic moment does not seem to change the simulation results significantly, which could be due to the small number and size of the nanoparticles. The hydrodynamic radius at c r ≈ 1 might indicate that dipolar interactions do not favor a large cluster. More specifically, for = 0.006 we observe that all 20 nanoparticles in the system cluster together to form a complete graph of bonds, Figure 15. This is known as gelation, a phenomenon often observed in coagulation dynamics and described by the Smoluchowski equations. [41,42] A more thorough study of the gelation and the effects of the nanoparticle size and magnetic fields on the hydrodynamic radius will be discussed in our future work. We mention that the nanoparticle radius r = 1.0 r c used for simulating the clustering dynamics of the nanoparticles is rather small for considering Brownian dominating nanoparticles, which are typically used in immunomagnetic assay. We used a small radius simply to increase the total number of particles involved in the clustering dynamics. The fact that there exists an optimal value of the concentration ratio of the biomarkers and ligands for a large hydrodynamic radius, should not be affected significantly by the nanoparticle size. On the other hand, the nanoparticle size affects the timescales of the clustering dynamics and the behavior of the cluster in oscillatory magnetic fields. This will be covered in our subsequent study.

Conclusions
In this work we have combined the dissipative particle dynamics with the Landau-Lifshitz-Gilbert equation to simulate a dispersion of conjugated magnetic nanoparticles. To the best of our knowledge such a hybrid DPD-SD solver that takes into account also the orientation (easy axis) and the rigid body dynamics of nanoparticles as well as their binding in the presence of analytes has not been considered before.
With such a numerical set-up we were able to simulate the nanoparticle dispersion in presence and absence of biomarkers. We have demonstrated that the presence of biomarkers promotes the clustering of magnetic particles. For low biomarker concentrations, increasing the concentration increases the hydrodynamic radius. Above a certain critical biomarker concentration the hydrodynamic size starts to decrease again, which is also known as the hook effect in immunoassay studies.
As a future prospect we have developed a coupling between the lattice-Boltzmann method [43] (LBM) to describe the solvent, the molecular dynamics of the nanoparticle cores and conjugated surface and magnetic dynamics of nanoparticle cores. The nanoparticles are coupled to the LBM fluid through the so called immersed boundary method (IBM). [44,45] This allows one to simulate larger spatio-temporal scales, which is essential to simulate dispersions of realistic 50 nm sized conjugated nanoparticles. The simulation results will be compared to our experimental measurements of the magnetic response at different analyte concentrations.

Model of Conjugated Magnetic Nanoparticles
The simulated system consisted of coated spherical magnetic nanoparticles, suspended in water. The particles were surrounded by biomarkers, which may attach to the nanoparticle coating.
The construction of magnetic nanoparticles is first discussed. To properly capture the rotational behavior, more than one coupling point per nanoparticle was needed. [46] The nanoparticle cores were constructed by placing a large number of DPD beads at a fixed radius r p from the center. This is the so called raspberry model, [47,48] which was convenient for the attachment of Figure 16. The nanoparticle core is constructed by placing many beads at a fixed radius from the center. The core is magnetic and is within the macrospin approximation represented by a magnetic moment m (blue arrow). The preffered orientation of the magnetic moment of an isolated magnetic particle in absence of external fields points along the easy axis n (orange arrow).
the nanoparticle coating and ligands. These beads were approximately homogeneously distributed on the surface of the nanoparticle, see Figure 16. In addition, one bead was put at the center of the particle, which carried the magnetic moment m of the entire nanoparticle core. This is the so called macrospin approximation, which was relevant for single-domain nanoparticles smaller than 40 nm. [49,50] The nanoparticles were decorated with a surfactant, in the case a shell of dextran, Figure 4. The physical advantages of such a shell were biocompatibility with the aqueous environment, a decreased toxicity and increased stability of the nanoparticle dispersion. [51,52] The nanoparticles were in addition conjugated with multiple recognition ligands. Each of these ligands consisted of a longer chain of a certain polymer, in the case polyethylene glycol (PEG), and a recognition part at the end.
The surfactants and the ligands were modeled by constructing chains of beads, connected by harmonic springs where K b is the elastic constant, r 0 is the equilibrium bond distance, and r is the distance between neighboring beads within a chain. A harmonic potential was also added between a chosen nanoparticle bead and the chains. The nanoparticle core beads to which the dextran chains are attached were chosen randomly at the start of the simulation and were covering 75% of the nanoparticle core surface. The same bond constants were used for all types of beads. Note that the usual factor 1 2 is already included in the bond constant K b . The number of beads used for each chain depended on their length. In the case the dextran shell consisted of chains with two beads, while the ligands consisted of seven beads of PEG and one bead corresponding to the recognition part. In all of the simulations six such recognition ligands were used per nanoparticle, which were distributed on the sphere according to the Thomson problem of point charges [53] to prevent any effects of inhomogeneous distribution of ligands. Furthermore, a bending potential V is added between neighboring bonds within the chain V ( ) = K ( − 0 ) 2 (13) where is the bond angle and 0 is the equilibrium bond angle. The equilibrium shape of the chain was determined by the complicated interplay of the interactions between the rest of the particles in the systems, especially the interaction between the ligands and the solvent. For PEG, the bond parameters were K b = 8500 kJ mol −1 nm −2 , K = 42.5 kJ mol −1 , r 0 = 3.30 Å, and 0 = 130 • , which are taken from the MARTINI coarse-grained model. [54] In the DPD model, three monomer units of PEG are represented as one bead (HO-[CH 2 CH 2 O] 3 -H). This was to ensure that different types of beads (water, dextran, and PEG) have a similar bead volume. A similar model for PEG as a surfactant has also been used in DPD simulations of endocytosis of PE-Gylated nanoparticles. [55] For dextran, the so called M3B coarse grained model was often used, where each monomer unit of dextran was represented by three beads. [56] Since the dextran chains were rather small compared to the ligands and were mainly there to prevent overlap of nanoparticles and to some extent to prevent water leaking inside nanoparticles, the same bond parameters K b , K , r 0 and 0 were taken as for PEG.

Dissipative Particle Dynamics
The simulations were performed using DPD method, [57,58] which was a mesoscopic method often used for simulating polymers, [25] biological systems such as cell membranes [59,60] and colloids. [61] It was recently also used in modeling of ultrasound propagation in water. [62] The system was represented by DPD beads, which interacted with each other via a soft repulsive potential, also called a conservative force F C ij , where i and j are the indices of the ith and jth bead, respectively. The force, which acted on a given bead consisted also of a dissipative force F D ij and a random force F R ij , which acted together effectively as a thermostat for the system. DPD thermostat conserved linear and angular momentum, which was suitable for modeling hydrodynamic systems. The total non-bonded force on the ith particle, F i , is thus a sum of these three forces where where a ij is the interaction parameter between beads i and j, is the friction parameter, is the noise strength, R and D are weight functions, ‚ r ij denotes the normalized vector of the inter-particle axis r ij = r i − r j , v ij = v i − v j is the relative velocity of www.advancedsciencenews.com www.advtheorysimul.com the beads and Θ ij = Θ ji is a Gaussian random variable with zero mean, ⟨Θ ij (t)⟩ = 0, and the following second moment From the fluctuation-dissipation theorem one obtained the relation D (r) = 2 R (r) and 2 = 2k B T . Typically, R (r ij ) is chosen to be equal to (r ij ) DPD is a mesoscopic method, which involved representing multiple atoms or molecules as a single bead. To find the length scale of the system one matched the mass densities of a DPD coarse grained description and the real density of water w [59] where v m =  (20) can also be derived by matching the volume occupied by one DPD bead The DPD number density of beads DPD as well as the mapping N m were free parameters that determine the length scale r c . It should be mentioned that one can start the other way around by first setting the scale r c and then deriving the mapping number N m . The density was set to DPD = 3.0 r −3 c throughout the work. For densities above DPD = 2 r −3 c , the following equation of state is valid [25] p = DPD k B T + a WW 2 DPD (22) where = (0.101 ± 0.001) r 4 c is a numerical factor and a WW is the interaction parameter for the beads of the solvent (e.g., water).
To determine the interaction parameter a WW one typically matched the dimensionless isothermal compressibility to the experimental one for water ( −1 ≈ 16). Using Equation (22) and ) sim , with n the number density of water molecules, one can solve for a ii The time scale is then  For the dissipative parameter , The choice of the interaction parameter for water a WW was therefore determined by the experimental value of the isothermal compressibility and the level of coarse-graining, given by the mapping number N m . It should be mentioned that this mapping was limited to about N m = 10, above which one typically observed solidification effects. [63] A different approach was needed to estimate other interaction parameters, see Table 1. Following the model of Groot and Warren, interaction parameters between the beads of the same type were set to the same value, in this case a ii = a WW ≈ 210 k B T r c was taken. An alternative method estimated the interaction parameters for the like-beads directly from solubility parameters. [64] For the nanoparticle surface one must ensure that it was impermeable to the solvent and the surfactants. Since the surface number density s was rather large s = 4 r −2 c it was found that the value 0.5 a WW was sufficient. Using an interaction parameter that was too large leads to large density oscillations near the rigid walls. [65,66] The interaction between the ligand beads and the water beads was estimated using the Flory-Huggins parameter ij a ij = a ii + 3.497 ij k B T r c (25) which was valid for the value of the bead number density used in this work, DPD = 3. [25,59] The Flory-Huggins parameter can be estimated from Hildebrand solubility parameters i , [67,68] see Table 2 Adv. Theory Simul. 2023, 6, 2200796 www.advancedsciencenews.com www.advtheorysimul.com where V i is the volume of bead i.

Open Boundary Molecular Dynamics
To study the behavior of the conjugated nanoparticle under nonequilibrium conditions the OBMD method was employed. [71][72][73] This was an algorithm that allowed one to simulate grand-canonical ensembles, where the system exchanges energy and mass with the environment. The number of particles in the simulation box was therefore not constant.
In this paper open boundaries were employed along the x axis, while the simulation box had periodic boundary conditions in the other two directions. The simulation box was split into three regions: a region of interest in the middle surrounded by two buffers (left and right), where the insertion and the deletion of particles was performed, see Figure 1. At each time step the chemical potential was maintained by adding the particles to the buffers according to the following dynamic equation where B is the buffer relaxation time, which is typically B ≈ 100 t, [71] ΔN is the number of particles to be inserted or deleted, N 0 is the equilibrium number of particles in either of the buffers and is a user-defined parameter, which was set for stability reasons to 0.8 in this work. In short, particles were inserted if ΔN > 0. If ΔN < 0 the algorithm was not performed, but particles were allowed to diffuse over the open boundaries and were deleted as soon as they leave the simulation box. The insertion was performed using USHER algorithm. [74,75] In USHER, the particle was in the first step randomly inserted in the chosen buffer region. Then an energy minimization was performed until the total energy of the inserted particle was below a certain threshold. Since DPD potential was rather soft an implementation of adaptive resolution scheme (AdResS) [76,77] was not necessary.
One imposed boundary conditions by adding external forces f i to the particles in the buffer regions and were determined by the momentum balance where J is the momentum flux tensor, A is the surface area of the interface between a buffer and the region of interest, n is normal to this interface and points toward the region of interest, while Δ(m i ′ v i ′ ) represents the momentum change when the particle was inserted or deleted from the buffer. The momentum flux tensor was composed of the pressure contribution and the tangential shear force JA ⋅ n = p xx An + xy At (29) where p xx is the pressure, t points along the desired shear force (ê y ) and xy is the imposed shear stress.

Hybrid Solver
To model the dynamics of the magnetic moments of the nanoparticles the Landau-Lifshitz-Gilbert equation was employed. In short, this equation described the relaxation of the magnetization toward the local magnetic field and the precession around it. Since the magnetic interactions were of long range, the calculation of this local magnetic fields involved summation of the fields of all the surrounding magnetic particles and their periodic images. The simulations were carried out using LAMMPS simulation package. [78] A solver was implemented, which incorporated the translation of the particles, their rigid body dynamics as well as the dynamics of the magnetization. The equations for the positions and the velocities of all the particles in the system were updated using the well known velocity-Verlet algorithm [79] v (k) where m (k) , r (k) , v (k) , f (k) are the mass, the position, the velocity, and the total force acting on the kth bead in the system. In the last step, Equation (32), the force is calculated using positions at t + t and velocity at t . It should be emphasized that the magnetic force acting between the nanoparticles was already included in this step in the force f (k) where f (l) stands for the magnetic force acting on the center bead of the lth nanoparticle with m (l) the corresponding magnetic moment, 0 is the magnetic constant 0 = 4 × 10 −7 N A −2 and r = |r| is the distance between the centers of the ith and lth nanoparticle, with r = r (l) − r (i) . The magnetic moments were measured in The sum in the expression for the dipolar force, Equation (33), extended to periodic images as well. For this the particle-particle particle-mesh method was used, [80] which split the dipolar force contributions into short range, and long-range contributions, which were efficiently calculated in the Fourier space using fast Fourier transform.
It was noted that the strength of the dipolar interactions compared to the thermal energy is given by where 2r is the inter-particle distance of two touching spheres. can be estimated from the typical saturation magnetizations M s of the nanoparticle core materials, for example, for magnetite nanoparticles M s ≈ 400 kA m −1 . [81] For a nanoparticle with radius r in DPD units, this yields |m| = 0.22 r 3 in units of m * and = 0.006 r 3 . Dipolar effects become important when ≫ 1, which is in the case of monodomain magnetite particles equivalent to r ≫ 5.5 r c ≈ 5 nm. In the analysis of nanoparticle clustering rather small nanoparticles with radius r = 1 r c was used, which means = 0.006. To study also the effects of larger magnetic parameters the magnetization of the nanoparticle core was artificially increased by a factor of 100. This is equivalent to = 60. The study of more realistic larger particle sizes will be covered in a future work.
The nanoparticle cores were treated as rigid bodies with three translational and three rotational degrees of freedom. The dynamics of a given nanoparticle was then determined by the total force and the total torque on the center of mass, which yielded the updated coordinates, velocities, and angular velocities. The numerical algorithm for the angular velocity and the orientation is similar to the velocity-Verlet algorithm where L is the angular momentum, is the total torque on a nanoparticle and q is a quaternion which conveniently describes the orientation and evolves according to the angular velocity [82] q = 1 2 q (38) where on the right-hand side is represented by a quaternion ≡ (0, Ω x , Ω y , Ω z ). The central particle carried in addition the magnetic moment and an orientation, defined by the easy axis n, which determined the direction of the magnetic moment of an isolated nanoparticle in equilibrium. In general, due to the dipolar forces, external magnetic fields and the clustering m ∦ n.
At each step the additional orientational degree of freedom n (l) was taken into account, by updating it according to the angular velocities (l) of the corresponding nanoparticle [20] The rotation of the magnetic moment was determined by the magnetic relaxation and precession processes, described by the Landau-Lifshitz-Gilbert equation [83][84][85] (l) where (l) = m (l) |m (l) | is the unit vector representing the direction of the magnetic moment of the lth nanoparticle, X is the magnetization current In Equation (41) R represents the Gilbert gyromagnetic ratio, R = | e | 0 ≈ 2.21 × 10 5 m A −1 s −1 with e the electron gyromagnetic ratio e ≈ −1.76 × 10 11 T −1 s −1 . Expressed in DPD units, the gyromagnetic factor equals R = 11.42. The parameter was the dimensionless Gilbert damping constant, which is set to = 0.01 in this work.
The anisotropic contribution H ani i can be derived from a potential H ani,(l) i where K a is the anisotropy constant and V c is the volume of the nanoparticle magnetic core. In this work K a = 223 kJ m −3 was used, which corresponded to cobalt ferrite nanoparticles. [86] The dipole contribution was a summation of the magnetic fields of all the nanoparticles and their periodic images (including the periodic images of a given nanoparticle for which the force was calculated) where r = r (l) − r (j) Related to the dipolar contribution is the demagnetization field of a given nanoparticle H dem i . In general it can be written as H dem i = D ij M j where D ij is the demagnetization tensor. For the case of a sphere this tensor was isotropic, D ij = − 1 3 ij . Since the dynamic equation conserved the magnitude of the magnetic moment, this contribution did not have an effect on the dynamics of spherical nanoparticles.
The field H th i was a stochastic contribution [85,87,88] with zero mean ⟨H th i ⟩ = 0 and the following second moment It has been shown in ref. [89] that both Itô and Stratonovich interpretation of the stochastic term lead to the same result when applied to Equation (40). Therefore, forth order Runge-Kutta method was used to solve this equation. The last term in Equation (42) was needed to correctly account for the conservation of the total angular momentum, [21,90] J = L + S, where L = I is the angular momentum of the nanoparticle with I the moment of inertia tensor and S = − 0 M s V c R is the spin momentum. This contribution is also known as the Barnett Adv. Theory Simul. 2023, 6, 2200796 field, [91,92] and is related to the Barnett effect [93] where magnetization is induced by rotation. The moment of inertia tensor was in the case isotropic and equal to I ij = 3 5 m V c r 2 ij for a hollow sphere with radius r and nanoparticle core mass density m .
Since the particles are small and single domain, the exchange interaction term E ex ≈ A|∇M| 2 with A the spin stiffness constant was discarded.
The conservative part of the torque c m acting on the nanoparticle can be calculated by inducing an infinitesimal rotation ( n = × n). A nonzero torque appeared due to the misalignment of the easy axis and the magnetic moment To guarantee the conservation of the total angular momentum J there is also a dissipative torque, which is oppositely equal to the dissipative torque, second term of Equation (41) The overall rigid body dynamics of the nanoparticle core is then given bẏ where x (l) is the position of the lth nanoparticle, v (l) is its velocity, F tot is the total force acting on the nanoparticle core, w is the total torque of the solvent and other nanoparticles, analytes, and the surfactants acting on the nanoparticle surface. Since the solvent was treated explicitly and was a source of random uncorrelated kicks to the nanoparticle surface beads, adding a stochastic torque to w in the spirit of rotational Brownian motion ref. [94,95] was not needed. As a side note, it was mentioned that the rotational dynamics of the nanoparticles was characterized by low Reynolds numbers, which meant that the angular velocity of the nanoparticle could in principle be estimated from the balance of hydrodynamic and magnetic anisotropy torques, yielding with n the rotational viscosity, which is proportional to the viscosity of the solvent. In Equation (52) it was assumed that the magnetization has already relaxed toward the direction of the magnetic field. According to the Einstein-Stokes equations, n = 8 R 3 h . Since n depended on the hydrodynamic radius, which was in principle not known beforehandthis approach will not be used. Therefore, Equation (47) was used to calculate the torque and corresponding forces on the nanoparticle, which did not require knowledge of any additional mesoscopic parameters, such as n .

Dynamic Binding
To model the attachment of the biomarkers to the recognition parts, which were located at the end of the recognition ligands, the so called dynamic binding was employed. This is a model of the binding process, with which it was already able to be described the cluster formation of nanoparticles. In short, if the biomarker came sufficiently close (e.g., less than a critical distance d c ) to the recognition part of the recognition ligand, a harmonic bond Equation (12) was formed between the biomarker and the recognition part. d c = r c was set in this work. For the newly formed bond, the same parameters was used as for PEG. Due to the known high affinity between the analyte and the ligand, this bond was not allowed to break once it was formed. A single biomarker (analyte) can bind to several recognition ligands (e.g., 1, 2, 4, …) thereby modeling the analytes as multivalent. This last aspect was crucial for the cluster formation. In ref. [96] multivalent binding of DNA oligonucleotide probes attached to a surface was shown to increase the sensitivity of bacterial genome detection.

Magnetic Nanoparticles and Colloidal Characterization
The employed magnetic nanoparticle were iron oxide nanoparticles (IONPs) coated with dextran and carboxylic PEG (product code is 104-56-701) synthesized by micromod Partikeltechnologie GmbH (Germany). The IONP crystal size was (30 ± 4) nm, the diffusion coefficient was (8.2 × 10 6 ± 5.2 × 10 5 ) nm 2 s −1 , the hydrodynamic diameter was 64 nm whose PDI = 0.08. IONP size and shape of the studied IONPs were evaluated by transmission electron microscopy (TEM) (see Figure 17) in a JEM1400 Flash (Tokio, Japan) TEM operating at a 100 kV located at Servicio de Microscopía Electrónica del Centro de Biología Molecular "Severo Ochoa" CSIC-UAM. TEM images were examined through manual analysis of more than 150 particles randomly selected in different areas of TEM micrographs using Image-J software to obtain the mean size and size distribution. The hydrodynamic size (D h ) of IONPs was determined by dynamic light scattering (DLS) with a Zetasizer Nano ZS equipment (Malvern Instruments, USA). IONP suspensions were diluted in double distilled water (DDW) to a final concentration of 0.05 g Fe L −1 in a commercial cuvette. The energy source was a laser emitting at 633 nm, and the angle between sample and detector was 173 • . For the determination of the diffusion coefficient of IONPs per unit volume in magnetic suspensions with iron contents 1 g Fe L −1 , Nanosight NS300 (Malvern Instruments, USA) was employed. 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 was collected for 60 s and repeated three times per sample.

Magnetic Characterization
SQUID measurements under quasi-static conditions were performed in a range of temperature between 4 and 300 K in IONP suspensions with a magnetic mass concentration of 1 g Fe L −1 , in a commercial SQUID model MPMS-XL from Quantum Design at Servicio de Tecnicas Físicas en la Facultad de Físicas de la Universidad Complutense de Madrid, see Figure 18.