key: cord-0934477-34938tb0
authors: Ye, Huilin; Shen, Zhiqiang; Li, Ying
title: Adhesive rolling of nanoparticles in a lateral flow inspired from diagnostics of COVID-19
date: 2021-02-22
journal: Extreme Mech Lett
DOI: 10.1016/j.eml.2021.101239
sha: 7ea978bb4e8ec5f6d35f08ee0ce954373252c59f
doc_id: 934477
cord_uid: 34938tb0

Due to the lack of therapeutics and vaccines, diagnostics of COVID-19 emerges as one of the primary tools for controlling the spread of SARS-COV-2. Here we aim to develop a theoretical model to study the detection process of SARS-COV-2 in lateral flow device (LFD), which can achieve rapid antigen diagnostic tests. The LFD is modeled as the adhesion of a spherical nanoparticle (NP) coated with ligands on the surface, mimicking the SARS-COV-2, on an infinite substrate distributed with receptors under a simple shear flow. The adhesive behaviors of NPs in the LFD are governed by the ligand receptor binding (LRB) and local hydrodynamics. Through energy balance analysis, three types of motion are predicted: (i) firm-adhesion (FA); (ii) adhesive-rolling (AR); and (iii) free-rolling (FR), which correspond to LRB-dominated, LRB-hydrodynamics-competed, and hydrodynamics-dominated regimes, respectively. The transitions of FA-to-AR and AR-to-FR are found to be triggered by overcoming LRB barrier and saturation of LRB torque, respectively. Most importantly, in the AR regime, the smaller NPs can move faster than their larger counterparts, induced by the LRB effect that depends on the radius [Formula: see text] of NPs. In addition, a scaling law is found in the AR regime that [Formula: see text] (rolling velocity [Formula: see text] and shear rate [Formula: see text]), with an approximate scaling factor [Formula: see text] identified through fitting both theoretical and numerical results. The scaling factor emerges from the energy-based stochastic LRB model, and is confirmed to be universal by examining selections of different LRB model parameters. This size-dependent rolling behavior under the control of flow strength may provide the theoretical guidance for designing efficient LFD in detecting infectious disease.

Diagnostics of COVID-19 is the primary tool to control the spread of SARS-COV-2 after the report of the first case in Wuhan City (China) [1] , due to the lack of effective therapeutics and vaccines. Current diagnostic methods have three main approaches: nucleic acid amplification tests (NAATs), antigen-directed virus detection, and serological assays [2, 3, 4] . 5 Among them, rapid antigen diagnostic tests, usually based on lateral flow device (LFD), have expanded widely due to its inexpensive, fast, qualitative and point-of-care testing [5, 6] . Although LFD demonstrates fast detection, the accuracy and sensitivity are limited compared to the NAATs. This low reliability may in part due to the inaccurate control of the flow and inadequate understanding of the interaction between antibody and antigen under the flow 10 environment.

To improve the reliability of LFD, attentions are paid to focus on the ligand-receptor binding (LRB) process [2, 5, 6, 7, 8] . In LFD analysis, the viruses are usually modeled as nanoparticles (NPs) coated with ligands on the surface, as presented in Fig. 1(a) . The detection is conducted through the capture of ligands by the biological receptors on the substrate. 15 Typically, there are three types of motion of NPs on the substrate: initial contact; adhesive rolling on substrate; and final firm-adhesion. These kinetic processes are strongly determined by the size-dependent physiological (Stokes' drag) and biophysical (LRB) conditions. Extensive efforts have been made to study the adhesion behaviors of NPs [9, 10, 11, 12] . Nevertheless, the complex interplay between hydrodynamics, Brownian motion, and stochastic This bistability is also confirmed by in vitro experiments [29] and theoretical analysis [30] .

In the corresponding theoretical work, the models that describe the kinetics of cell rolling on the vessel wall under laminar shear flow have been established by torque and energy balance analyses. Although the target is cell instead of NP and the analysis is limited to two-dimensional circumstance, the theoretical approach such as energy balance may shed and the receptors are distributed on the substrate. The NP is immersed in a semi-infinite simple shear flow bounded by the substrate. The viscosity and shear rate of the flow are η andγ, respectively. We fix the environmental temperature of this system at T = 310 K. It is assumed that the distributions of ligands and receptors are uniform, and their densities are n l and n r , respectively. When the NP moves close to the substrate, the ligands and receptors 70 can interact with each other and then they might bind together to form biological bonds as depicted in Fig. 1(b) .

Under the shear flow, besides the LRB, the NP also experiences hydrodynamics and thermal fluctuation. The thermal fluctuation is attributed to the ambient fluid, acting as heat reservoir with temperature T . The Péclet number is defined as P e = 2Ru max /D ∞ , where u max =γR corresponds to the maximum flow velocity around the NP and D ∞ = k B T 6πηR is the Stokes-Einstein diffusivity. k B is the Boltzmann constant. Since the minimum size of NP we study here is R = 40 nm and the shear rate of flow has the orderγ ∼ 10 −4 ns −1 , the Péclet number corresponds to the order P e = R is the radius of the NP, v is the moving velocity of the NP, δeq represents the equilibrium length of the ligand receptor bond, x is the distance between ligand and receptor, r is the corresponding radius of the crossed circle and β is the angle with shear direction. don is the cutoff of the ligand receptor interaction, and A is the corresponding critical position at time t, which has a angle β0 with the shear direction. And at time t + dt, A moves to A . The ring on the plane is the projection of the moving area (denoted by red line) in the 3D configuration of spherical NP, which describes the movement of NP at time dt on the substrate. θ in the moving part is defined as the angle between the r and vertical direction.

fluctuation can be negligible. To confirm this, we still include the thermal fluctuation effect in the numerical simulations, and we find it contributes to the motion of NPs significantly only when the shear rate is extremely small (low Péclet number, cf. Appendix B). However, for simplicity, the thermal fluctuation is not considered in the theoretical analysis. Through removing the thermal fluctuation, the NP's motion is governed by the competition between hydrodynamics and LRB. Through this competition, the NP can present either rolling or static on the substrate. For the theoretical analysis, we assume the NP keeps rolling with velocity v on the substrate with no slipping, therefore, v = Rω with the angular velocity ω. The kinetic energy of the NP is

where I cm is the inertial moment, and J = mR 2 + I cm . The hydrodynamic work exerted on the NP done by the shear flow in the time interval (0, t) is

where F h and M h are hydrodynamic force and torque, respectively. The work consumed by the LRB can be described as

where M LR is the torque exerted on the NP due to LRB. We assume the NP is initially static, then we have the energy balance form:

which can be expressed as

Taking the derivation with respect to t, we have

Assuming the NP reaches a steady state, saying dv(t) dt = 0, then we have the governing equation of the steady motion for the NP 5   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 J o u r n a l P r e -p r o o f

In the following, we restrict our discussion to the steady motion of NP unless otherwise stated.

To gain insights into the motion of NP, the concrete expression of the Eq. (7) needs to be further deduced. Firstly, based on Stokes' law [32] , [31] derived the exact solution of Stokes' equations for a sphere moving near a plane wall under simple shear flow. According to [31] , the force and torque that the NP experiences in the flow are:

where F * s and F * t are hindered force coefficients induced by flow and translation of NP.

Similarly, T * s and T * r are hindered torque coefficients led by flow and rotation of NP. Furthermore, the empirical expressions of the coefficients are provided:

where Λ = R/h, and h is the distance between the center of NP and the substrate. Although the confinement of the substrate can affect the motion of NP, it is obvious that, from the force and torque expressions, the nature 80 of Stokes' prediction that larger particles will move faster than smaller ones in simple shear flow is maintained.

Secondly, the LRB is characterized by the biological bonds formed between ligands and receptors. In our model, the biological bond association and dissociation are governed by the probabilistic adhesion model developed by [33] based on the Bell model [34] . This model has been widely used in describing cell and NP adhesive behaviors and has performed well in terms of the accuracy and simplicity [35, 36, 37, 20] . When a ligand on the NP surface is close to a receptor on the substrate within a characteristic length d on , there is a probability P on to create a new biological bond between them. Reversely, an existing bond suffers a breakup probability P off within a critical length d off . The probabilities are mathematically defined as

where k on and k off are the association and dissociation rates, respectively. They have the forms:

, where σ on and σ off are the effective on and off strengths, denoting a decrease and increase of the corresponding rates 85 within the interaction lengths d on and d off , respectively. k B T is the energy unit. k 0 on and k 0 off are association and dissociation rates at the equilibrium length l = l 0 between ligand and receptor, respectively. l 0 represents the equilibrium length of the biological bond. Thus, the force of a biological bond is: F b = k s (l − l 0 ), which will exert on both receptor and ligand, but with different directions along with the bond. Here, k s is the bond strength which is 90 attributed to the rupture energy E r = − 1 2 k s ∆x 2 m that is required to break up a single bond. ∆x m is the maximum extension of the bond. ∆t in the expressions is the time step in both theoretical analysis and computer simulation. Note that only a single bond is allowed to form at each ligand or receptor. The selection of parameters in the LRB process and the corresponding references are listed in Table 1 .

References At time t, the NP is in torque-free state and the biological bonds are called 'closed bonds' under this state. When the NP starts to move, the LRB torque is induced by the imbalance between the association torque of leading edge led by newly formed bonds and the dissociation torque of trailing edge induced by stretching of closed bonds. In the following, we strive to find the analytical expressions of association torque and dissociation torque by geometrical analysis. The smallest distance between NP and substrate is defined as δ eq (cf. Fig. 1(c) ).

x denotes the distance of ligand from the substrate and r is the corresponding radius of the crossed circle. β is the angle with shear direction as shown in Fig. 1(c) . Similarly, β 0 is the angle between the critical position where the ligand has a distance d on from the substrate and shear direction. We track the critical point A (distance d on from te substrate ) at the leading edge of the NP from t to t + dt. At time t + dt, it will move to the position A with distance x from the substrate. We project this part (denoted as red line in Fig. 1(c) ) of the NP on the substrate, it is nearly a ring region, namely influence region, with intra radius r = R 2 − (R + δ eq − x) 2 and width vdt is the displacement of point A within time dt. Inside this region, there are new interactions between ligands and receptors. x can be expressed as R(1 − cos(β 0 − vdt R )) + δ eq . In the influence region, for a parcel with angle dθ, the number of association bonds is

Because the ligand density is about two times higher than that of the receptor, receptors in the influence region are expected to be fully saturated [10] , and this estimation is confirmed in our numerical studies through counting the number of active bonds (cf. Appendix A). The biological bond, in general, assumes an arbitrary angle with the substrate. For example, the closed bonds should demonstrate stochastic directions. However, this cannot apply to the bonds in the influence region. The angle of a bond will be around π/2 in the influence region at the leading edge. Because, for a specific ligand on the NP, the newly created bond forms in the vertical direction with high probability as the distance is minimum under this 7 J o u r n a l P r e -p r o o f circumstance, leading to a lower energy state of the LRB. Similarly, in the limited influence region at the trailing edge, the bonds will be stretched compared to their last state. And those bonds with angles larger or smaller than π/2 for specific ligands suffer a high probability to break up before stretching. Hence, the bonds in the influence regions are assumed to be in vertical direction (shear direction). To further validate this scenario, the resultant bond force in the simulation is calculated. It reveals that the bond force in the flow direction is nearly zero, though there is a small increase for the cases with a large moving velocity of NP (cf. Appendix A). Therefore, we only consider the bond force in the shear direction and it can be expressed as F b = k s (x − δ eq ). The arm of the bond force is L = r sin θ, then we have the torque on NP exerted by the association bonds dM on = π 0 dN b LF b . Ignoring the second order (dt) 2 , we have association torque

Similarly, we can obtain the dissociation torque

but here x = R(1 − cos(β 0 − β * + vdt R )) + δ eq , where β * = vts R corresponds to the angle of NP's rotation within time t s . The integration time t s should reflect the kinetics of interactions between ligands and receptors. Here, k 0 on and k 0 off characterize the LRB process. According to the experiment [41] , the association rate k 0 on is typically few orders of magnitude larger than the dissociation rate k 0 off . The fast association process is difficult to capture, and is not of significant when considering flow with small shear rate. Therefore, k 0 off is adopted to quantify the LRB process, and the integration time t s is chosen as 1/k 0 off . Eventually, we obtain the torque due to LRB as

In the numerical simulation, the Lattice Boltzmann method (LBM) is adopted to solve the fluid flow, which is an efficient and accurate method for fluid dynamics [42] . The linearized Boltzmann equation has the form of:

where f i (x, t) is the distribution function for fluid particles with velocity e i at position x and time t, f eq i (x, t) is the equilibrium distribution function and τ is the non-dimensional relaxation time, F i is an external forcing term. In this simulation scheme, D3Q19 model is 8   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 used [43] . The equilibrium distribution function f eq i (x, t) can be calculated as:

. The term c s represents the sound speed which equals ∆x/( √ 3∆t). The relaxation time is related to the kinematic viscosity in Navier-Stokes equation in the form of ν = (τ − 1 2 )c 2 s ∆t. The external forcing term can be discretized as [44] : (14) is solved by the integration algorithm proposed by [45] . Once the particle density distributions are known, the fluid density and momentum are calculated as

The NP is represented by a rigid spherical shell with surface uniformly discretized into 9   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 ξ is the damping coefficient. v is the velocity of the NP node, and u is the velocity of the fluid at the same location as the NP node. Because the fluid flow is described in Cartesian mesh and the fluid nodes are fixed, u should be interpolated by the surrounding fluid nodes. Based on the immersed boundary method [46], a smoothing kernel is adopted to interpolate the surrounding fluid nodes to the location of NP nodes.

where x and X denote fluid node and NP node, respectively. We use a 4 point approximation to the Dirac delta function, providing a support of 64 grid nodes

The dissipative force F d is directly applied to the NP nodes, and the motion of NP is driven by the resultant force and torque. In addition, to reflect the existence of a NP in the fluid flow, F d is also spread into the surrounding fluid nodes. And the same smoothing kernel as the interpolation of velocity should be used. Although this numerical scheme about 100 the particle-fluid coupling has been validated in [43] and our previous works [47, 48, 20] , we further validate it using the case of translation of a single NP in quiescent fluid. The NP is placed in the center of a cubic domain of side length 24R with initial velocity u 0 . As shown in Fig. 2 (a), at time t = 5 ns, a starting recirculation flow structure is observed, which is also reported in [49] . It demonstrates the strong momentum transfer between the flow and NP. At t = 20 ns, the momentum transfer becomes weak, and the radius of curvature of streamlines becomes larger. Eventually, at time t = 100 ns, the streamlines become parallel to each other, and there is no momentum transfer between the NP and surrounding fluid flow. We plot the evolution of NP's velocity in Fig. 2(b) . It indicates that, for a short time t < 40, the velocity of NP presents analytical exponential decay exp( ξt m ), where m is the mass of NP. In 110 the long time regime, we also find the power decay t −3/2 , which was pointed out in [49] . This long-time tail decay behavior is due to the hydrodynamic interaction between the NP and flow field. These good agreements demonstrate the accuracy of the hydrodynamic coupling between NP and fluid field. More details about our numerical scheme and implementation can be found in [48] .

The thermal fluctuation of the system is accounted for by introducing Gaussian white noise [50] . The fluctuation stress tensor s αβ satisfies a fluctuation dissipation relation with the form

where r is the position of the fluid node and δ ij is the Kroneckers symbol. With the same setup of the movement of NP in quiescent fluid, we apply thermal perturbation with temperature T = 310 K instead of initial velocity u 0 . With the thermal effect, the NP will make stochastic motion, and we calculate the time-averaged mean squared displacement (M SD) with the definition M SD = (r i − r 0 ) 2 . · represents the ensemble average of experiments for a single particle in the fluid flow. For a self-diffused Brownian particle in quiescent fluid, the diffusion coefficient D theo is predicted by the theoretical analysis in [51] that D theo = k B T /ξ. We also calculate the diffusion coefficient D sim in the simulation with D sim = M SD/6t, and the evolution of D sim is presented together with D theo in the Fig. 2(c) . It is found that, after a short transition, D sim converges to the theoretical prediction D theo . This confirms that the 125 thermal perturbation added into the fluid model is accurate to capture the Brownian motion of the NP. Other validations for the thermal fluctuation are referred to [43] .

For the LRB process, each node distributed on the NP surface denotes a ligand. And the receptors are illustrated by the nodes distributed with a square pattern on the substrate. The biological bonds formed between ligands and receptors are controlled by the stochastic 130 model (Eq. (9)). The specific validations of this model have detailed in the work of [47] .

According to Eq. (7), the motion of NP is governed by the balance between LRB torque and hydrodynamic torque. To quantify the torques with relation to NP's velocity v, we first 135 integrate the expression of M on and M off from analytical equations Eq. (11) and Eq. (12), respectively. Fig. 3(a) shows the results of a typical case of R = 80 nm, E r = −45 k B T and other parameters given in Table 1 . It is found that the association torque M on demonstrates three stages with the increment of moving velocity v: i) extremely small when v < 0.01 nm/ns; ii) continuous increase in the regime 0.01 nm/ns < v < 0.035 nm/ns; and iii) approximately 140 constant when v > 0.035 nm/ns. In the stage of extremely small, because v is very small, the influence region is narrow and the newly formed bonds are limited. Furthermore, part of the formed bonds will quickly become closed bonds due to the very large association rate k 0 on . Hence, the association torque M on exerted on the NP is extremely small. With the increment of v and enlargement of influence region in the stage of continuous increase, more 145 and more active bonds form in the leading edge, which makes the association torque M on become significant and continuously increase. However, a further increase of v will lead to the saturation of newly formed bonds, resulting in the plateau of M on (stage of approximately constant).

By contrast, only two stages are found for the dissociation torque M off . The first stage is 3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 In order to compare the LRB torque M LR and hydrodynamic torque M Ho = F h R + M h , the M LR is first obtained by Eq. (13) and presented in Fig. 3(a) . M LR is found to have three stages: i) dramatic increase; ii) small decrease; and iii) approximately constant with the 160 increase of v. Among them, the stage of dramatic increase is dominated by the dissociation torque M off and the stage of small decrease is attributed to the increase of association torque M on .

Furthermore, the hydrodynamic torque M Ho for cases with the same sized NP, but different flow strengths (shear rateγ) are demonstrated in the same figure of M LR (cf. Fig. 3(a) ). 165 Here, the shear rateγ is normalized by k 0 off that characterizes the time scale of LRB process. According to Eq. (8), M Ho monotonically decreases with the increase of v as shown in Fig. 3(a) , and this makes the M LR and M Ho curves intersect at a specific point. The intersection means the balance of torques exerted on the NP and represents the steady state. Also, to investigate how the moving velocity evolves, we recall the time-dependent governing equation 170 Eq. (6) and numerically integrate it. The evolution of moving velocity v is shown in Fig. 3(b) for corresponding cases. From Fig. 3(a) , we find that when the flow is weak (γ = 0.4 k 0 off ), the moving velocity v of NP is extremely low, and this steady state is called firm-adhesion (FA), in which the NP can hardly roll on the substrate. Under this circumstance, an initially high velocity of NP will continuously decrease and finally, it reaches the static state (FA, see 175 Fig. 3(b) ). When the shear rate increases to moderate (γ = 2.2 k 0 off ), the flow strength is comparable with the strength of LRB, and the movement of NP becomes significant. The NP can always reach a steady rolling state, which is irrelevant to the initial velocity (high (P − ) or low (P + ), cf. Fig. 3(b) ). We name this type of motion as adhesive-rolling (AR). With the further increase of flow strength, the intersection enters the regime where the M LR is 180 approximately constant and becomes irrelevant to the moving velocity v. Therefore, stronger flow (γ = 3.2 k 0 off ) will lead to the hydrodynamics-dominated regime where the NP motion is determined by M Ho . In this regime, even though the initial velocity is small, the strong flow will drive the NP to move and accelerate until the NP reaches a steady state ( Fig. 3(b) ). This freely moving motion is named free-rolling (FR). 

where a 1 = 6F * s + 8T * s and a 2 = 6F * t + 8T * r . In Eq. (19), we defineγ =γ l +γ h , wherė 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 We integrate the Eq. (19) , and the analytical v −γ curve is presented in Fig. 4 . To 190 demonstrateγ h andγ l , we also plot the v −γ curve for free-rolling of NP (black dashed line) under the circumstance that the NP is driven to freely move on the substrate only by shear flow. Then we can identifyγ h andγ l as shown in Fig. 4 . Furthermore, the simulation results (symbols) are presented in the same figure, which are in good agreement with theoretical predictions. The v −γ curve can be separated into three parts: i) v is very small; ii) v 195 continuously increases; and iii) v dramatically increases with the increment of flow strengtḣ γ. These three parts exactly correspond to the three types of motion: FA, AR, and FR respectively, which will be discussed below.

It should be noted that, in the FA regime, the simulation result demonstrates that the NP maintains its averaged static state under a low shear rate, which is slightly inconsistent with 200 the theoretical analysis. This discrepancy is caused by the assumption that, in the theory, the NP is always expected to move on the substrate. However, the NP can start to move only when the flow strength exceeds a critical value (γ >γ on ). It can be explained as follows. Without shear flow, the torques on the trailing and leading edges of NP are balanced in the static state. However, when the NP tends to move with the increment of the flow strength, 205 the biological bonds on the trailing edge are strongly stretched while those on the leading edge are compressed, in comparison with their equilibrium state. The torque exerted on the NP induced by the biological bonds on the trailing edge is larger than that on the leading edge. This imbalance is named LRB barrier that hinders the start of the NP. Therefore, before NP can move, the LRB barrier should be overcome by the compensation from hydrodynamics 210 through continuously increasing the flow strength, which is like mechanical static friction. Furthermore, we find that the critical shear rateγ on will decrease if the thermal fluctuation is considered in the simulation, which indicates the thermal fluctuation can lower the LRB barrier through introducing thermal torque (cf. Appendix B). 14 1 2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 In the second part of the v −γ curve (AR regime), the moving velocity v continuously 215 increases with the increment ofγ. Comparingγ h withγ l , we can seeγ l dramatically increases, whileγ h increases linearly, which subjects to free rolling motion. This means that the hydrodynamic torque increases to be comparable with the LRB torque. At this stage, only a small portion of the hydrodynamic torque is consumed to maintain the free-rolling-like motion, while the majority is used to balance the LRB torque. 220 However, the curve in the third part (FR regime) is found to be parallel to the freerolling curve (black dashed line). In this part, theγ l will be constant, which corresponds to the plateau stage of the LRB torque. On the contrary, theγ h will continuously, though linearly, increase and exceed theγ l . Therefore, hydrodynamics will dominate the motion of free rolling.

In the transition from AR to FR, there should exist a critical velocity point (γ c , v c ). After this point, the moving of NP becomes irrelevant to LRB, which is due to the saturation of LRB torque and continuous increase of hydrodynamic torque with the increment of moving velocity shown in Fig. 3(a) . Therefore, the critical velocity v c corresponds to the maximum LRB torque. We have known that the dissociation torque M off increases with the increment of v, and eventually reaches plateau regime. While association torque M on slightly increases with the increment of v only after v reaches the relatively large value, and then also becomes approximately constant. Hence, the maximum of LRB torque occurs when M off reaches maximum and M on slightly increases. To determine this critical point, we consider the essential parts of association and dissociation probabilities P on and P off . Fig. 5 presents the probabilities that are dependent on the length of biological bond. We define x o as the length of bond where P on = P off and x min that corresponds to β * which points to the minimum length of biological bond in the influence region on trailing edge. At velocity v = v c , we let x min = x o . When v increases, x min will be smaller than x o . Although P off is 15   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 lower, the influence region will enlarge a little and dissociation torque M off will have a small increase. However, the association probability P on has a substantial increase, which leads to the significance of association torque M on . Therefore, LRB torque will decrease under this circumstance. When v decreases, x min will be larger than x o . Although P off increases, the influence region will shrink and less closed bonds will be stretched, which results in the decrease of dissociation torque M off . And the decrease of P on makes the association torque M on negligible. Hence, LRB torque will also decrease in this case. Combining the above two aspects, we can determine that the bond length at the intersection point between P on and P off should be the minimum length of biological bond (x min ) in the influence region, which corresponds to the maximum LRB torque. Furthermore, according to the relation between β * and x min , we can obtain the critical velocity

It is found that the critical velocity v c is strongly dependent on the radius R of NP, which means different sized NPs have different critical velocity points for the transition from AR to FR. In the next section, we will further study the effect of NP size in both AR and FR regimes.

To quantify the motion of different sized NPs, we conduct three more cases with the radii R = 40, 60, and 100 nm. The corresponding v −γ curves are presented in Fig. 6(a) .

Except for the FA regime, they have different performances in both AR and FR regimes. In the regime of AR, we find that the smaller NP moves faster than the larger ones under the same flow strength. This is counter-intuitive as Stokes' prediction says the larger particle will move faster than the smaller one in the same flow. However, the inclusion of LRB process breaks this prediction and predicts the opposite result in the AR regime. Furthermore, the critical points of the transition from AR to FR vary with the radius of NP, which has been confirmed by Eq. (20) . The corresponding critical shear ratesγ c are also calculated based on the expression of v c (Eq. (20)). It is plotted against the radius of NP in Fig. 6(b) . We 240 find that the larger the NP is, the higher the critical shear rate will be, which indicates the stronger flow is needed to drive the larger NP to move with the same velocity as the smaller one.

The critical points for different sized NPs together consist of the critical v −γ curve. We put this curve in Fig. 6 (a) and it illustrates the boundary between the AR and FR regimes.

The v −γ plane is further split into two main regions: AR (green) and FR (grey) plus a narrow FA region (red).

In the FR region, the motion of NP is dominated by the hydrodynamics, which recalls the Stokes' prediction that larger NP moves faster than the smaller one. However, it is not true in the whole FR region. Since the critical point depends on the size of NP, larger NP 250 has higher moving velocity when entering FR region compared to the smaller one. Therefore, at the initial stage of the FR region where different sized NPs enter FR regime, the smaller one can still move faster than the larger one. For example, we zoom into the initial stage of FR region for the cases R = 40 nm and R = 100 nm. It is clear that, for the same flow strength, the smaller NP (R = 40 nm) moves faster than the larger one (R = 100 nm).

However, further increment of the flow strength will make the larger NP exceed the smaller one. Therefore, there is a closed region between the critical v −γ curve, the v −γ curve for R = 40 nm and the v −γ curve for R = 100 nm, as denoted in the inserted zoom-in figure.

We name this region as memory zone where LRB torque still exists and contributes to the motion of NP, though it keeps constant. It is found that the larger size difference between 260 NPs is, the larger the memory zone will be. It means that, compared to a specific small NP (R = 40 nm), if the size of the other NP is larger, higher hydrodynamic torque is required to compensate for the LRB torque and eliminate the memory.

However, the memory of LRB will become negligible when the hydrodynamic torque further increases and Stokes' predication is fully recovered. It is characterized by the intersection of two v −γ curves of different sized NPs. Assuming two NPs with radii R 1 and R 2 , their shear rates can be expressed aṡ

where A = a 2 /a 1 . At the intersection point,γ R 1 =γ R 2 , then we have the velocity Finally, we can have theγ in based on the v −γ relation. Furthermore, we adopt the case of R 2 = 500 nm as a reference, and plot theγ in against the radius of NP in Fig. 6(b) . It is 265 found thatγ in increases with the increment of NP radius, which illustrates that the smaller NP will be easier to exceed the referenced one (R 2 = 500 nm) in the FR regime.

The NP size has demonstrated different influences in the AR and FR regions. Compared to the relation v ∝γR in the FR region, the v −γ relation in the AR region remains to 270 be further explored due to the complex expression of LRB torque with respect to v. In the v −γ plane, the v −γ curves of different sized NPs are scaled with respect to the reference (R = 100 nm). By fitting both theoretical and numerical results, we find that there exists a scaling factor α ≈ −0.2 ± 0.04 that can make all v −γ curves approximately collapse into that of the referenced NP. Furthermore, to have a consistent view, the v −γ curves are also scaled 275 in the FR region according to scaling ∼ R. And we find all the v −γ curves will be parallel to each other. The vertical shift between the scaled curves is caused by the different critical points (v c ,γ c ) that depend on the radii of NPs. Nevertheless, they can collapse together through vertical shifting with the distances that are determined by the critical points. And the final results are presented in Fig. 7 with simulations results represented by symbols. 280 We have known that R scaling in the FR region stems from the hydrodynamics based on Stokes' prediction. To gain insights into R α scaling in the AR region, we recall the expression of M LR . It mainly consists of three parts that are relevant to the radius of NP: (i) area of influence region Σ; (ii) length of bond x; and (iii) association and dissociation probabilities P on and P off . Firstly, the area of influence region Σ scales to r 2 , hence, we have Σ ∝ R 2 .

Secondly, the length of bond x ∝ R. According to Eq. (19) , these two parts cancel R 3 and 18   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 the scaling R α should be mainly dependent on the last part -contributions of association and dissociation probabilities P on and P off . In this work, the probabilities are established based on the energy of the bond E b , which is related to (l − l 0 ) 2 . The dependence of R on the length of bond l − l 0 in these probabilities leads to the scaling factor α. Based on this 290 analysis, the scaling can be also extended to the force-based adhesion model, in which the force is defined as F b = k s (l − l 0 ). Although the scaling factor α is an approximation, it can still predict that smaller NP should move faster than the larger counterpart. This universal scaling law should direct all the stochastic model-based adhesive motion of NPs.

According to above discussions, the scaling law is determined to be relevant to the stochas-295 tic model of LRB. We further investigate the influence of important parameters in the probability model of LRB. Because the dissociation torque on the trailing edge dominates the LRB process, the dissociation rate k 0 off at the equilibrium length and effective off strength σ off are systemically varied. The corresponding v −γ curves are presented in Fig. 8(a) and (b), respectively. In Fig. 8(a) , it is found that k 0 off has negligible effect in the AR regime, 300 while it affects the critical points of AR-to-FR transition. In Eq. (12) , with the increase of dissociation rate k 0 off , the dissociation probability P off will increase. However, the integration time t s = 1/k off that characterizes the LRB timescale will decrease. Therefore, the increase of k 0 off cannot significantly affect dissociation torque. For the association torque on the leading edge, when NP's moving velocity is not large, the association probability P on is small and 305 the association torque can be negligible (cf. Fig. 3 ). Hence, when v is not large enough, k 0 off has negligible effect on the v −γ curve in the AR regime. However, when the moving velocity v increases to be close to the critical point, the association probability becomes significant. And the association torque decreases due to the decrease of the integration time t s . Hence, the maximum LRB torque enhances with the increase of k 0 off , which leads to the raise of the 310 critical point. In Fig. 8(b) , the effective off strength only affects the dissociation probability. When the moving velocity v is small, the term exp( σ off (l−l 0 ) 2 2k B T ) is small. The change of σ off has negligible effect on P off and thus has minimal effect on v −γ curve in the AR regime. However, when the moving velocity v is close to that at the critical point, the increase of σ off will result in the significant increase of dissociation probability and the maximum LRB 315 torque, which makes the right shift of the critical point as shown in the figure.

To further exclude the dependence of parameter selections, we study effects of receptor density n r and rupture energy E r as shown in Fig. 8(c) and (d), respectively. It is found that the receptor density n r plays an important role in the AR region. For the same flow strength, the moving velocity v decreases with the increase of the receptor density. The 320 increase of n r leads to more bonds formed on the leading edge and more bonds are stretched on the trailing edge, which enlarges the LRB torque that inhibits the moving of NP on the substrate. Also, in Fig. 8(d) , the rupture energy E r denotes the strength of a single bond and higher rupture energy means stronger biological bond. The stronger bond will result in the larger LRB torque with the same number of bonds. Therefore, the moving velocity v of 325 the NP will decrease with the increase of the rupture energy E r for the same flow strength.

As σ off and k 0 off have the same effects on the motion of NP in the AR region, we adopt σ off to examine its dependence on the size scaling. Fig. 9 (a) presents the size scaling result for σ off = 0.15 pN/nm. It is found that the scaling factor α also works well, which illustrates that the scaling law is irrelevant to the parameter selection of the probability model. Furthermore, we investigate the dependence of rupture energy E r as a typical example of parameter selection in the model of torque M LR . The scaling law is shown in Fig. 9 (b) for E r = −30 k B T . We find that the v −γ curves also collapse together with the factor α in the AR region. The examinations of parameter selections further confirm the universal size scaling law of NPs in the AR region as v ∝γR α and α ∼ −0.2 ± 0.05. 335 

In the present work, both theoretical and numerical analyses are utilized to study the adhesive motion of a spherical NP on the substrate under LRB and shear flow. In the theoretical formulation, the motion of NP is governed by the interplay of hydrodynamics and LRB. The governing equation is subsequently derived through geometrical analysis and 340 energy balance. It predicts three types of motion: FA, AR, and FR with the increase of flow strength. When the flow is weak (lowγ), the hydrodynamics is small compared to the LRB and the flow cannot drive the moving of NP. Moderate flow strength will make the hydrodynamics comparable with the LRB, which results in adhesive-rolling of NP. The FR regime is characterized by the strong flow strength that dominates the motion of NP. In this 345 regime, the NP can freely move, which is like the motion of NP only under simple shear flow.

Between these three types of motion, two transitions are identified. The first transition is FA-to-AR. The reason why weak flow cannot drive the NP to move is the existence of LRB barrier. Only when the hydrodynamics overcomes the LRB barrier, the NP can start to move. The second transition is AR-to-FR. Due to the saturation of LRB torque, a contin-350 uous increase of flow strength will make hydrodynamics overwhelm the LRB, triggering the transition from AR to FR.

The motion of NP is also found to be strongly dependent on its radius. In particular, in the AR regime, the smaller NPs move faster than their larger counterparts. This conflicts with the Stokes' prediction that larger NP will move faster than the smaller one in the same 355 shear flow. Furthermore, this contradiction is found to extend to the initial stage of the FR regime. Since the critical point that characterizes the transition of AR-to-FR depends on the size of NP, different sized NPs have different critical velocities. Therefore, there exists LRB memory zone between two different sized NPs. Stronger flow strength is required to eliminate the LRB memory for larger NP to pursue and exceed the smaller one. 

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

In the setup of parameters related to the LRB process, the density of receptors on the substrate is two times lower than that of ligands on the NP surface. Therefore, we assume receptors on the substrate are saturated, which means that, inside the region on the substrate where receptors can interact with ligands, all the receptors will form bonds with corresponding ligands. Through projection of NP on the substrate, we can calculate the projection area 395 A 0 , hence, the number of active biological bonds should be N active = n r A 0 in the steady adhesive rolling state according to the theoretical assumption. To confirm this aspect, in the numerical simulations, we count the number of active biological bonds when the NP reaches steady state, and collect the results for different sized NPs, which is shown in Fig. A.10(a) . From the comparison, the numerical results are consistent with the theoretical predictions, 400 which confirms our theoretical assumption on saturation of receptors in the steady state regime. 22   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 The formation of biological bond between ligand and receptor is a stochastic process, therefore, the direction of the bond should be, in principle, assumed arbitrary. However, this has limited application to the bonds in influence region. The angles of these bonds should be 405 around π/2 in the influence region on the leading edge, because those with other angles will break up due to larger length than the critical length d on . This is also true for the bonds in the influence region on the trailing edge. To validate this aspect, in the simulation, we calculate the resultant bond force and show the force in the flow direction f x in the Fig. A.10(b) . It is found that, for different sized NPs, the force f x is nearly zero, despite of the limited force 410 for NPs moving with large velocity. For comparison, we also plot the Stokes' drag F x of the NP with size R = 40 nm moving in the flow. It is obvious that the bond force that the NP experiences in the flow direction is negligible compared to the Stokes' drag. Therefore, it validates the assumption that the ligand-receptor bonds in the influence region is nearly vertical.

Appendix B.

In the simulation, we find that when the shear rate is lower than a critical valueγ * on , the NP demonstrates averaged static state as shown in Fig. B.11 . This is not predicted in the theoretical analysis due to the assumption that NP is moving with velocity v. Wheṅ γ >γ * on , the numerical results are consistent with the theoretical predictions. Furthermore, 420 we consider the thermal fluctuation of the fluid. It is found that only when the moving velocity of NP is small, the thermal fluctuation has a significant effect on the motion of NP. This is because the Péclet number is small for low shear rate (slow moving particle). In addition, we find the critical shear rateγ on is smaller than that without thermal fluctuation (γ on <γ * on ). It 23   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 means the NP is easier to move under thermal fluctuation, which presents that the thermal 425 fluctuation can help NP overcome the LRB barrier. With the increase of shear rate, the hydrodynamics is much stronger than thermal fluctuation. The thermal fluctuation effect cannot contribute to the regular motion of NP. Therefore, there is no significant difference between those with and without consideration of thermal fluctuation after NP starts to move. 3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 J o u r n a l P r e -p r o o f  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 

A pneumonia outbreak associated with a new coronavirus of probable bat origin

The sars-cov-2 spike 435 protein binds sialic acids and enables rapid detection in a lateral flow point of care diagnostic device

Current best practices for respiratory virus testing

How nanophotonic label-free biosensors can contribute to rapid and massive diagnostics of respiratory virus infections: Covid-19 case

Signal amplification and quantification on lateral 445 flow assays by laser excitation of plasmonic nanomaterials

Assay techniques and test development for covid-19 diagnosis

Point-of-care diagnostics for respiratory viral infections

Quantifying nanoparticle adhesion mediated by specific molecular interactions

The adhesive strength of non-spherical particles mediated by 455 specific interactions

The role of specific and non-specific interactions in receptormediated endocytosis of nanoparticles

Modeling particle shape-dependent dynamics in nanomedicine

Influence of red blood cells on nanoparticle targeted delivery in microcirculation

Multiscale modeling and uncertainty quantification in nanoparticle-mediated drug/gene delivery

Microparticle shape effects on margination, near-wall dynamics and adhesion in a three-dimensional simulation of red blood cell suspension

Nanoparticle diffusion in sheared cellular blood flow

Motion of a nano-spheroid in a cylindrical vessel flow: Brownian and hydrodynamic 475 interactions

Self-propelled colloidal particle near a planar wall: A brownian dynamics study

Interplay of deformability and adhesion on localization of elastic micro-particles in blood flow

Shape-dependent transport of microparticles in blood flow: From margination to adhesion

Traction dynamics of filopodia on compliant substrates

Cell movement is guided by the rigidity 485 of the substrate

Matrix elasticity directs stem cell lineage specification

Neutrophil recruitment and function in health and inflammation

The state diagram for cell adhesion under flow: leukocyte rolling and firm adhesion

The kinetics of l-selectin tethers and the mechanics of selectin-mediated rolling

Bistability of cell adhesion in shear flow

Probing the cytoadherence of malaria infected red blood cells under flow

Rolling adhesion of cell in shear flow: a theoretical model

Slow viscous motion of a sphere parallel to a plane wall-ii couette flow

Course of theoretical physics

Simulation of cell rolling and adhesion on surfaces in shear flow: general results and analysis of selectin-mediated neutrophil adhesion

Models for the specific adhesion of cells to cells

Lifetime and strength of adhesive molecular bond clusters between elastic media

Wall shear stress-based model for adhesive dynamics of red blood cells in malaria

Quantifying the biophysical characteristics of plasmodium-falciparum-parasitized red blood cells in microcirculation

Ligand-receptor interactions between surfaces: the role of binary polymer spacers

Multiscale modeling of blood flow and soft matter

The reaction-limited kinetics of membrane-to-surface adhesion and detachment

Binding constants of membrane-anchored receptors and 525 ligands depend strongly on the nanoscale roughness of membranes

Lattice boltzmann method for fluid flows

Hydrodynamic forces implemented into lammps 530 through a lattice-boltzmann fluid

Discrete lattice effects on the forcing term in the lattice boltzmann method

Fluctuating lattice-boltzmann model for complex fluids

The immersed boundary method

Cell stiffness governs its adhesion dynamics on substrate under shear flow

Computational modeling of magnetic particle margination within blood flow through lammps

Decay of the velocity autocorrelation function

The immersed molecular finite element method

On the movement of small particles suspended in stationary liquids required

This work was mainly supported by the National Science Foundation (OAC-1755779). The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources (Frontera project and the National Science Foundation award 1818253) that have contributed to the research results reported within this paper. This work was partially supported by a fellowship grant (to H.Y. and Z.S. ) from