A unified analysis of nano-to-microscale particle dispersion in tubular blood flow

Blood is a complex fluid suspended with multiple species, primarily including red blood cells (RBCs), platelets, white blood cells, and various biomolecules (such as von Willebrand factors, albumin, and fibrinogen) that covers length scales ranging from nanometers to micrometers.^1,2 In microvessels under physiological flow conditions, RBCs migrate toward the axis of the tube and leave a cell-free layer (CFL) near the wall.^3,4 Such a phenomenon, well known as the Fahraeus-Lindquist effect,⁵ contributes to the hemorheological heterogeneity of the blood flow. Unraveling the dispersion properties of solutes and cells of various sizes ranging from nanometer to micrometer in such heterogeneous blood flows under vascular confinement can potentially lead to the optimal design of drug carriers and better understanding, intervention, and control of vascular diseases.

As a relevant example of microscale particle transport in blood, platelets margination has shown to play an important role in affecting the rate of clot formation in hemostasis and thrombosis.⁶ Motivated by that, a plethora of studies over the past decades have been dedicated to unraveling the mechanistic mechanisms of margination or segregation of microscale particles/cells in blood(-like) flows through perfusion experiements,^7–10 continuum-level modeling,^11,12 and direct numerical simulations.^10,13–21 The platelet margination is found to be primarily driven by the cross-stream hydrodynamic fluctuation^14,15,22 or equivalently the RBC-enhanced shear-induced diffusion^12,18 in the RBC-laden region synergistically accompanied by the sinklike effect of the CFL.¹⁸ 

Nanoscale particle (NP) dispersion in blood flow, on the other end of the spectrum, has recently received considerable attention due to the fast development of nanodrug delivery techniques that have the potential to revolutionize the traditional therapeutics.²³ Although the effective diffusivity of nanoscale solutes in blood flow were measured decades ago,²⁴ it is not until the past several years multiscale particle-level simulation techniques^25–28 become feasible. Tan et al.²⁵ apply a coupled Brownian dynamics and immersed finite-element (FE) method to study the influence of RBCs on the NP dispersion in blood flows, showing substantial margination behavior for 100 nm particles. Through both in vivo and in silico techniques, Lee et al.²⁶ show that submicron particles (>500 nm) can marginate while NPs (∼100 nm) are mostly trapped in the RBC-laden region. Muller et al.²⁷ performed two-dimensional (2D) simulations and suggest that microscale particles compared to submicroscale particles show better margination propensity. Liu et al.²⁸ develop a multiscale complex blood solver and evaluate the role of Brownian diffusivity (BD) vs RBC-enhanced shear-induced diffusivity (RESID) in affecting the biodistribution of NPs. Recently, Liu et al.²⁹ characterize the complete 3-D diffusivity tensor of NP in blood under various shear rates and hematocrits, which can be employed to modeling large-scale NP biotransport applications.

Although the transport of both nanoscale and microscale particles in blood has been understood to a large extent, there is still a lack of a systematic interrogation of the particle dispersion behavior across nano-to-microscale sizes using a unified computational approach. Consequently, questions such as whether nanoscale particles exhibit margination qualitatively the same as microscale particle does still remain controversial. A recent effort by Cooley et al.³⁰ using in vitro experiment and 2D in silico simulation to understand the cross-length-scale particle margination and adhesion propensity has set an example for a unified understanding of the nano-to-microscale particle dispersion in blood flows. However, the general physical mechanisms behind the multiscale particle dispersion/margination phenomenon in blood are still not presented; besides, the 2D simulation could still overlook the 3D nature of the tubular blood flow phenomena.

In this work, we employ a recently developed 3D multiscale and multicomponent blood flow solver^2,28,29,31 to tackle the dispersive characteristics of spherical, rigid particles with sizes spanning nano-to-microscale in a tubular blood flow. Particle suspension dynamics in the presence of thermal fluctuation, RBC-particle direct and hydrodynamic interactions and wall-bounded confinement effect are captured under a unified 3D computational framework. The strong correlation between the nonuniform distribution of particle radial diffusivity and the equilibrium distribution of particle radial concentration is highlighted to gain mechanistic understanding of the occurrence of particle-size-induced dispersion-to-margination transition.

The remainder of the paper is organized as follows. In Sec. II, we describe the unified multiscale complex blood solver and layout the techniques for evaluation of particle radial concentration and diffusivity. In Sec. III, we present the simulation results, where the particle radial distribution at equilibrium state is discussed under various confinement ratio, hematocrit, and particle sizes. The mechanisms that drive the particle size-dependent dispersion-to-margination transition will be discussed. In Sec. IV, we conclude this systematic study.

The numerical method used to simulate the bidisperse particle-RBC suspensions confined in a tubular flow is through a multiscale and multicomponent complex blood flow method^28,29 that couples the lattice-Boltzmann/Spectrin-link (LB-SL) method³¹ with the lattice-Boltzmann/Langevin-dynamics (LB-LD) method.² This method leverages the off-lattice nature of the LB-LD approach and the efficiency of the course-grained SL RBC membrane method to concurrently simulate the dynamics of across nano-to-microscale particles and microscale deformable capsules with a fixed LB lattice resolution.^2,28 The hybrid LB-LD-SL method has previously been verified with theory^2,28 and validated against experiments.^29,31,32 Figure 1 demonstrates a nanoscale particle-RBC bidisperse suspension flow through a 40 μm vessel, where the computational methods for each module are denoted accordingly and presented in detail as follows.

Simulation of the suspending fluid is based on the Aidun-Lu-Ding (ALD) LB method.^33,36,37 The LB method solves the discretized Boltzmann transport equation in velocity space through the streaming-collision process. In streaming, the fictitious fluid particles propagate along discrete velocity vectors forming a lattice space. In collision, the fluid particles at each lattice site collide with each other, causing the relaxation of the particle distribution function (PDF) toward a local “Maxwellian” equilibrium PDF. The collision term is linearized based on the single-relaxation-time Bhatnagar, Gross, and Krook (BGK) operator.³⁸ The temporal evolution of the particle distribution function is given as

where f_i is the fluid PDF, $fi(0)$ is the equilibrium PDF, r is the lattice site, e_i is the discrete lattice velocity, t is time, τ is the single relaxation time, and $fiS$ is a forcing source term introduced to account for the discrete external force effect. The method has a pseudo speed of sound, $cs=Δr/(3Δt)$⁠, and a fluid kinematic viscosity, ν = $(τ−Δt/2)cs2$⁠, where Δt is the time step and Δr is the unit lattice distance. The positivity of ν requires τ > Δt/2. In the LB method, time and space are typically normalized by Δt and Δr, respectively, such that Δt_LB = Δr_LB = 1 are employed to advance Eq. (1). In the near incompressible limit (i.e., the Mach number, Ma = u/c_s ≪ 1), the LB equation recovers the Navier-Stokes equation³⁹ with the equilibrium PDF given in terms of local macroscopic variables as

where ω_i denotes the set of lattice weights defined by the LB stencil in use. The macroscopic properties such as the fluid density, ρ, velocity, u, and pressure, p, are obtained via moments of the equilibrium distribution functions as $ρ=∑i=1Qfi(0)(r,t)$⁠, $u=1ρ∑i=1Qfi(0)(r,t)ei$⁠, and $pI=∑i=1Qfi(0)(r,t)eiei−ρuu$⁠, respectively. Here, I is the identity tensor and pressure can be related to density and the speed of sound through p = $ρcs2$⁠. For the D3Q19 stencil adopted in the current study, Q is equal to 19. Along the rest, nondiagonal, and diagonal lattice directions, ω_i is equal to 1/3, 1/18, and 1/36, and |e_i| is equal to 0, Δr/Δt, and $2(Δr/Δt)$⁠, correspondingly.

The nano-to-microscale particle suspensions are resolved through a two-way coupled LB-LD method which has been verified^2,28 and validated against experiments.²⁹ This approach treats suspended particles in Stokesian regimes as point particles, while the volume exclusion effect of the particles are resolved through potential equations. The dynamics of LD particles is governed by the Langevin equation (LE)

where m_p is the mass of a single particle. The conservative force, C_p, specifying the interparticle and particle-surface interaction forces, is determined by calculating the directional derivatives of the total potential energy U_total as

where in this study U_total accounts for the particle-cell/wall short-distance interactions, as discussed in Sec. II F. The frictional force, F_p, is assumed to be proportional to the relative velocity of the particle with respect to the local viscous fluid velocity,^40,41

where u_p denotes the particle velocity and u(r_p, t) is the interpolated LB fluid velocity at the center of the particle. The friction coefficient, ζ, is determined by Stokes’ drag law, ζ = 3πμd_p, where μ is the dynamic viscosity of the suspending fluid. The stochastic force, S_p, explicitly gives rise to the Brownian motion of the particle and satisfies the fluctuation-dissipation theorem (FDT)⁴² by

where i, j ∈ {x, y, z}, α and β run through all the particle indices, δ_ij and δ_αβ are Kronecker deltas, δ(t − t′) is the Dirac-delta function, k_B is the Boltzmann constant, and T is the absolute temperature of the suspending fluid. The angle brackets denote the ensemble average over all the realizations of the random variables. Since we are concerned with long-time scale phenomenon, the over-damped LE is adopted in the current study as suggested in the work of Liu et al.^2,28

The modeling of RBC dynamics and deformation is through the coarse-grained spectrin-link (SL) membrane method^43,44 coupled to the LB method.³¹ The hybrid LB-SL method has been extensively validated against experimental measurements and is capable of capturing both the deformation and dynamics of single RBC³¹ and the rheology of RBC suspensions at physiological hematocrit³² with good accuracy and efficiency.

In the LB-SL model, the RBC membrane is modeled as a triangulated network with a collection of vertices mimicking actin vertex coordinates. The Helmholtz free energy of the network system, E, including in-plane, bending, volume, and surface area energy components,⁴⁵ is given by

where the in-plane energy, E_IP, characterizes the membrane shear modulus through a wormlike chain (WLC) potential⁴⁶ coupled with a hydrostatic component;⁴⁴ the bending energy, E_B, specifies the membrane bending stiffness, which is essential in characterizing the equilibrium RBC biconcave morphology;^44,45 the volumetric constraint energy, E_Ω, and the area constraint energy, E_A, preserve the RBC volume and area conservation, respectively, when subject to external forces.

The dynamics of each vertice are updated according to Newton’s equations of motion,

where $v$_n is the velocity of the vertice at the position x_n and M is taken as the fictitious mass of the RBC that is evaluated as the total mass of the cell divided by the number of vertices, N_v.^17,29 The number of vertices used to discretize the RBC membrane is N_v = 613, which has shown to yield adequate resolution to resolve the hydrodynamic forces³⁴ and capture single RBC dynamics³¹ and concentrated RBC suspension rheology³² when coupled with the LB method. $fnLB$ specifies the forces on the vertex due to the fluid-solid coupling. $fnCC$ are the forces due to cell-cell interactions. The forces due to the Helmholtz free energy based on the SL model is determined by

The SL method is solved by integrating equations (8) at each LB time step using a first-order-accurate forward Euler scheme in consistency with the LB evolution equation to avoid excessive computational expense.^28,31

The coupling between fluid and RBC is accomplished through the ALD fluid-solid interaction scheme.³⁷ In this method, the momentum transfer at the fluid-solid interface is accounted for by applying the bounce-back operation along lattice links that cross solid surfaces. As a result, the no-slip condition is enforced by adjusting the PDFs of the fluid nodes at the end point of a link along the i direction through

where i′ is the direction opposite to i, f_i(r, t⁺) is the postcollision distribution, and u_b is the solid velocity at the intersection point with the link. The fluid force exerted on the vertex on the RBC membrane mesh can be determined by

which is applied to the advancement of the RBC dynamic equation through Eq. (8).

The LD particles with Brownian effect are coupled to the nonfluctuating LB fluid in a two-way fashion using spatial extra/interpolation schemes,^41,47 through which treatment the long-distance many-body hydrodynamic interactions and the correct temperature scale can be captured simultaneously without empirical renormalization.^2,28,48 Specifically, the hydrodynamic force exerted on the particle, $FpH$⁠, is systematically decomposed into frictional and stochastic components as

where the fluid velocity at the particle site, u(r_p, t), is interpolated based on surrounding LB velocities and applied to update the LD particle dynamics through Eq. (3). The weighting functions, w(r, r_p), for interpolation is constructed using a trilinear stencil.^28,40 Since F_p and S_p are both originated from the “collision” between NP and liquid molecules, $FpH$ (instead of F_p) is assigned back to the fluid phase to satisfy momentum conservation. The same weighting function is then applied to constructing the local forcing source term as

which is adopted by Eq. (1) to update the local hydrodynamics. The coupled LB-LD method, similar to the external boundary force (EBF) method,⁴⁹ modifies the conventional LB evolution equation into Eq. (1) by adding the forcing distribution function $fiS(r,t)$⁠, which is shown to approximate the Navier-Stokes equation in the macroscopic scale.⁵⁰ 

The short-distance interactions between particle and RBC or between particle and the vessel wall is through Morse potential that forbids particles from penetrating the RBC membrane or the vascular wall. This contact model has previously been used in the characterization of the NP long-time diffusion tensor in an unbounded sheared blood, where the calculated NP diffusivity compares favorably with experimental measurements.²⁹ The Morse potential function is given as

where r is the normal distance between the particle center to the RBC surface, r₀ is a cut-off distance in which no interaction forces are present, D_e is the potential well depth, and β is a scaling factor. The Morse potential is imposed when r ≤ r₀ to preserve the repulsive effect. Model parameters are adjusted to match the measured intercell potential energy, as discussed in the work of Liu et al.^28,29 Specifically, the scaling factor is set to β = 2 μm⁻¹, the surface energy has a value of D_e = 10⁷k_BT, and the equilibrium distance is set to r₀ = d_p/2 + 10 nm. This simple contact model, bridging the LB-LD approach^2,28 and the LB-SL method,^31,32 can capture the margination phenomenon of microscale particles comparably well as the DNS approach does.^18,31Figure 2 presents the temporal evolution of the ensemble average of the radial displacement of microscale particles, 2⟨r_p⟩/d_p, where the particle margination process through the LB-LD-SL approach and that via DNS compares favorably well especially when approaching the equilibrium stage (⁠$tγ̇w≥$ 2000).

The particle number concentration at specific radial location, C_n(r, t), can be evaluated as

where N denotes all LD particles in the simulation and L_v is the length of the tube. The radial bin width, Δr, is set to one tenth of the tube radius to accurately resolve the radial profiles of the particle concentration distribution.¹⁷ The bulk ensemble-averaged particle number concentration can be calculated as ⟨C_n⟩ = $4N/πdv2Lv$⁠, which is later used to normalize the particle local concentration.

The particle radial diffusivity is evaluated through a moving time-origin measurement⁵¹ of the particle mean squared displacement (MSD) based on a fixed sampling time interval (STI). The STI is properly chosen to exclude the short-time ballistic regime.^2,28 By measuring the radial MSD of particles at a radial location r, the local instantaneous particle radial diffusivity can be evaluated according to

where N denotes all LD particles in the simulation and Δt is chosen to be 1000 in lattice units.²⁹ Same technique can be applied to measure the radial distribution of RESID, $DrrRBC(r,t)$⁠, where the BD is excluded by setting S_p = 0. The bulk ensemble-averaged particle radial diffusivity is calculated as ⟨D_rr(t)⟩ = $∑α∈N{[rpα(t+Δt)−rpα(t)]2}2Δt$⁠; similarly, the bulk ensemble-averaged RESID can be obtained as $⟨DrrRBC(t)⟩$ = $∑α∈N{[rpα(t+Δt)−rpα(t)]2}2Δt|Sp=0$⁠. The equilibrium counterparts of the particle radial diffusivity are denoted as ⟨D_rr⟩ and $⟨DrrRBC⟩$ without time dependence.

The physical problem of particle-RBC suspension flow through a straight tube can be defined by the vessel diameter, d_v, the systemic hematocrit, ϕ, the particle diameter, d_p, the wall shear rate, $γ̇w$⁠, and temperature, T, given fixed RBC properties (hydrodynamic radius, a_RBC, and membrane shear modulus, G). Apart from the hematocrit, the corresponding nondimensional parameters are the confinement ratio, $dv*=aRBCdv$⁠, which determines the severity of the RBC finite size effect; the particle-cell size ratio, $dp*=dpaRBC$⁠, that quantifies the length-scale discrepancy between the two species suspended; the Peclet number, $Pe=3μπγ̇wdpaRBC2kBT$⁠, which describes the competition between the shear-induced diffusion and the Brownian diffusion; and the capillary number, $Ca=μγ̇waRBCG$⁠, which defines the deformability of the RBC capsule.

In this work, we consider $dv*$ ranging from 0.07 to 0.29, corresponding to typical diameters of arterioles.⁵² The particle-cell size ratio considered ranges from $dp*$ = 0.003–0.86, covering typical size of biomolecules and cells (such as von Willebrand factor, vWF, and platelet) in blood flows. Given the low sensitivity of platelet margination to shear rate,¹⁸ a physiologically relevant wall shear rate, $γ̇w$ = 1000 s⁻¹, typical in arterioles or capillaries is considered for all cases. The fluid viscosity is set to the same as blood plasma, μ = 1.2 cp. The temperature is set to the body temperature, T = 310 K. The RBC membrane has a shear modulus of G = 0.0063 dynes/cm. The effective hydrodynamic radius of RBC is a_RBC = 2.9 μm. As a result, the dependence on Pe is determined by $dp*$⁠. The deformability of RBC is fixed with Ca_G = 0.55.

All simulations are initialized with the particles and RBCs uniformly and randomly mixed in the tube, except the particles are only seeded at 2r/d_v ≤ 0.6. Periodic boundary conditions are imposed on the two ends of the tube. The tube has a length of L_v/a_RBC ≥ 10 to ensure the periodic boundary treatment exerting negligible effect on the particle/cell transport. This is paper focuses on the dispersive characteristics at equilibrium albeit the transient effects may play a significant role in the particle distribution in microvascular bifurcating structures.^53–57 The equilibrium conditions are determined by tracking the particle accumulation at each radial location until it plateaus. As an example in Fig. 3, we present the temporal change of particle number percentage and radial diffusivity at different radial locations, where the particle number percentage, n(r)/N, is defined as the number of particles within certain radial layer, $n(r)=∑α∈N{δ[rpα(t)−r]}$⁠, normalized by the total particle number, N, within the simulation domain. The simulation is performed with d_p = 100 nm, d_v = 20 μm, and ϕ = 0.2. The equilibrium state is arrived at $tγ̇w∼2000$⁠, when the mean values of both n(r)/N and D_rr(r)/D_B remain unchanged with respect to time.

We first interrogate the dispersion characteristics of NPs under different confinement ratios controlled by adjusting vessel diameters in the range of d_v = 10 ∼40 μm (corresponding to $dv*$ = 0.29–0.073), which corresponds to typical size of arterioles or capillaries in human.⁵² The particle size is fixed to d_p = 100 nm. The wall shear rate is set to $γ̇w$ = 1000 s⁻¹, and the systemic hematocrit is set to ϕ = 0.2, which are within the range of physiological hemorheological ranges in human arterioles or capillaries.⁵² The number of particles simulated in the microvessels are N = 4000, 1000, and 250 from large to small vessels, respectively, to conserve the particle volume concentration.

Figure 4 presents the simulation snapshots of NP and RBC equilibrium distribution in microvessels under various confinement conditions. Qualitatively, the RBC dynamic mode changes from tank-treading/tumbling dominant to parachuting dominant^31,58,59 as the vessel diameter decreases from 40 to 10 μm. Such an increase in confinement does not alter the radial distribution of shear rate significantly but does change the radial distribution of local hematocrit to a large extent, as shown in Figs. 5(a) and 5(b). Specifically, for the case with d_v = 40 μm, the RBC-laden region shows a relatively uniform distribution except near the axis of the tube where the shear rate is close to zero. Consequently, the NP concentration, C_n(r), at 0 < 2r/d_v < 0.2 appears twice the bulk average NP concentration, ⟨C_n⟩, while C_n(r) near the wall exhibits slightly lower values than ⟨C_n⟩. As the vessel diameter decreases to 20 μm, the dimensionless CFL thickness δ_CFL/d_v increases from ∼0.2 to ∼0.4, i.e., the RBC-laden region becomes relatively more focused. Moreover, the local hematocrits get intensified especially at the inner boundary of the CFL and at the tube axis. These hemorheological changes substantially affect the equilibrium radial distribution of C_n(r)/⟨C_n⟩. As a result, the location of peak NP concentration shifts toward the CFL region, as shown in Fig. 5(c). Further confining the system to d_v = 10 μm, allowing only one train of RBCs parachuting through the vessel, appears to slightly enhance the peak concentration of NP at the CFL region while decreasing the NP concentration at the RBC-laden region. Previous microfluidic experiments by Nott, Guazzelli, and Pouliquen⁶⁰ also show the enhancement of the number percentage of particles adhered to the wall when the width of a ∼40 μm channel is reduced by half.

The distributions of the NP radial concentration can be better understood by evaluating the NP radial diffusivity, as depicted in Fig. 5(d). The NP diffusivity in the velocity-gradient direction based on the unconfined linearly sheared blood flow²⁹ is also plotted for comparison using shear rates and hematocrits of the d_v = 40 μm case. In general, increasing the confinement reduces the magnitude of NP radial diffusivity, where the unconfined case shows up to two folds the radial diffusivity, D_rr(r)/D_B, of the d_v = 40 μm case. Besides changing the magnitude of D_rr(r)/D_B, adjusting confinement ratio also alters the radial distribution of D_rr(r)/D_B. For the d_v = 40 μm case, the NP radial diffusivity shows high values near the CFL inner boundary and low value in the RBC-laden region, which is similar to the D_rr(r)/D_B distribution in the unconfined case. This distribution of D_rr(r)/D_B seems to be the cause of the low concentration of NPs near the CFL region and the high concentration at the RBC-core region. The increase in confinement (decrease in vessel diameter to 20 or 10 μm) renders the radial location of high D_rr(r)/D_B to move toward the RBC-laden region, which appears to be responsible for the shift of the high NP concentration region toward the CFL as observed in the high confinement cases (d_v = 10 and 20 μm).

Overall, the increase in the confinement ratio enhances the NP near-wall concentration by inhibiting the NP diffusion near the wall. This however does not warrant the margination of NPs, given no excessive NP concentration [C_n(r) < 1.5⟨C_n⟩] is observed near the wall as the vessel confinement increases to capillary scale. It is noted that when vessel size decreases to capillary scale, retention of microscale particles in the RBC-induced recirculation is reported,⁶¹ which however is not observed in the current study with NPs under the studied hemorheological conditions.

Changing hematocrit significantly modifies the apparent viscosity of blood^32,62 and could drastically influence the RBC-enhanced shear-induced diffusivity of NPs in sheared blood flow.²⁹ To understand how the variation of systemic hematocrit in microvessels changes the local hemorheology and hence the NP radial distribution, we investigate the hematocrit dependence of the NP radial dispersion behavior under various systemic hematocrits in the range of ϕ = 0–0.3. For the cases considered here, we select a fixed vessel diameter of d_v = 20 μm and a NP size of d_p = 100 nm. The wall shear rate is set to $γ̇w$ = 1000 s⁻¹. The number of NPs are set to N = 1000.

Figure 6 plots two snapshots of NP-RBC distribution in a 20 μm vessel, where the high hematocrit case (ϕ = 0.3) exhibits a thinner CFL compared to the low hematocrit case (ϕ = 0.15) as expected. Figures 7(a) and 7(b) present quantitative analysis of the hemorheological response to the change of systemic hematocrit, where the increase in systemic hematocrit affects the radial distribution of both the local shear rate, $γ̇(r)$⁠, and the local hematocrit, ϕ(r), which are two competing drivers for the particle cross-stream migration. On the one hand, it alters the flow structure from a Poissuelle-type flow toward a plug-type flow; as a result, the local shear rate in the RBC-laden region decreases, which drives the Brownian particles toward the tube axis.⁶³ On the other hand, it increases the local hematocrit in the RBC-laden region, which drives the particles to migrate to the wall.

The adjustment of the two competing effects leads to certain variation of the NP radial distribution, as presented in Fig. 7(c). At ϕ = 0, the NP dispersion is purely driven by the Brownian diffusivity and the shear-gradient driven dispersion. The former is isotropic, while the latter tends to drive the particle toward the low shear region.⁶⁰ As a result, the NP distribution shows a high NP concentration in the core and a low concentration near the wall. Increasing the systemic hematocrit generally alters the NP distribution such that the high NP concentration region shifts to the CFL. Interestingly, a slight increase in ϕ from 0 to 0.05 appears to be enough to shift this paradigm of NP distribution, leading to about 3-fold decrease in NP concentration at the core and ∼1.5 folds increase in NP near wall concentration. Further increasing ϕ slightly increases the near wall NP concentration but also increases the NP concentration at the RBC-laden region. Correspondingly, in Fig. 7(d), the D_rr(r) value (especially in the RBC-laden region) shows a nonmonotonic change with respect to ϕ, which first increases by up to 3 folds as ϕ rises to 0.05 and gradually gets inhibited to be close to the theoretical Brownian diffusivity as ϕ further increases to 0.3. The inhibition of NP radial diffusivity at high hematocrit can be explained by the excessive local ϕ(r) and low $γ̇(r)$⁠, as shown in Figs. 7(a) and 7(b).

Therefore, low systemic hematocrits appear to be optimal to enhance the NP near-wall concentration in microvessels, owing to the relatively high local shear rates and moderate local hematocrits that does not inhibits the NP dispersion in the tubular core. Nevertheless, changing hematocrit does not lead to the margination of NP.

So far we have focused on the long-time dispersion behavior of NPs in microvessels under various confinement ratios and hematocrit conditions. In these cases, particles do not show margination behavior. Instead, a nonuniform radial distribution of particles is observed with the particle concentration near the wall being less than 1.5 times its bulk average concentration. Besides, the near-wall concentration of NP is dynamically conserved at equilibrium, accompanied by the cross-migration of NPs between the CFL and the RBC-laden region due to severe Brownian effect.²⁸ 

In this section, we consider the size-dependent dispersion behavior of nano-to-microscale particles in microvessels. Particles with sizes ranging from d_p = 10–2500 nm are considered, covering particles ranging from nanoscale biomolecules such as vWFs in globular conformation to microscale cells such as platelets. The vessel diameter is fixed to d_v = 20 μm. The wall shear rate is set to 1000 s⁻¹. The systemic hematocrit is kept at ϕ = 0.2. To maintain the volumetric concentration of the particle phase in the dilution limit (i.e., ≪1%), the number of large particles considered in the system is reduced accordingly but kept above 50 to ensure statistical significance as consistent with our previous margination study.¹⁷ 

Figure 8 presents the equilibrium distributions of RBCs and particles. Qualitatively, NPs show nonuniformly dispersed distribution across the vessel, where NPs at any radial position can disperse to a random radial location given enough time owing to severe Brownian effect. These features are qualitatively different from the margination behavior of microscale particles, where excessive concentration and retention of particles in the CFL can be observed.^17,18 Figure 9(a) further depicts the radial distribution of C_n(r)/⟨C_n⟩ for various particle sizes. As the particle size increases above 1 μm, the CFL region exhibits a prominently high particle concentration. Specifically, for particles with a diameter d_v = 2.5 μm, a five times bulk average particle concentration, C_n(r) ≈ 5⟨C_n⟩, can be observed at the CFL region. These observations are consistent with particle margination study using a microfluidic perfusion system by Namdee et al.,⁶⁴ where the number percentage of particles adhered to the wall gets increased by 5–7 times as the particle size changes from nanoscale to microscale. The change of the particle number concentration at the CFL, $CnCFL$⁠, vs that at the RBC-laden region, $⟨Cn⟩non−CFL$⁠, as a function of particle size is further plotted in Fig. 9(b). When the particle size is below 1 μm, the $CnCFL/⟨Cn⟩non−CFL$ value shows weak dependence on the particle size with only a slight increase from 1.0 to about 1.5 as d_p changing from 10 nm to 1000 nm. As particle size exceeds 1 μm, the $CnCFL/⟨Cn⟩non−CFL$ value increases abruptly (up to ∼17) and margination occurs. Particle size d_p = 1 μm seems to be a critical watershed that divides the dispersion state and the margination state, as denoted in Fig. 9(b).

To shed light on the size-dependent dispersion behavior of particles in tubular blood flows, the distribution of particle radial diffusivity, D_rr(r)/D_B, is plotted in Fig. 9(c). For NPs, the D_rr(r)/D_B distribution tends to be uniform due to the dominance of isotropic Brownian diffusivity. As the particle size increases above one micronmetre, the RBC-laden region shows prominent enhancement of D_rr(r)/D_B compared to the CFL. Moreover, both at the CFL edge (2r/d_v ∼ 0.8) and close to the tube axis (2r/d_v ∼ 0.2), the magnitude of D_rr(r)/D_B peaks and the inner peak is more pronounced than the peak close to the CFL. For the 2500 nm particles, the inner peak shows more than ten times D_rr(r)/D_B values of that at the CFL. The radial distribution of D_rr(r)/D_B seems to be inversely correlated with the radial distribution of C_n(r)/⟨C_n⟩ in terms of the radial location, suggesting that the margination of microscale particles is probably due to the large magnitude difference in D_rr(r)/D_B between the RBC-laden region and the CFL region.

In Fig. 9(d), we plot the ensemble-averaged radial diffusivity, ⟨D_rr⟩, as a function of the particle size. The ensemble average is performed among all particles located at various radial locations at the equilibrium state. For small d_p, the ⟨D_rr⟩ value asymptotically matches the Stokes-Einstein relation due to the dominance of Brownian diffusion. Increasing the particle size decreases the effect of thermal fluctuation and leads to the deviation of ⟨D_rr⟩ from D_B. The value of ⟨D_rr⟩ eventually plateaus at the microscale size regime, where the RESID is dominant over BD. The bulk ensemble-averaged RESID, $⟨DrrRBC⟩$⁠, seems to be weakly dependent on the particle size, similar to the particle diffusivity observed in a unbounded sheared blood flow.²⁹ Subtracting ⟨D_rr⟩ with $⟨DrrRBC⟩$ shows an overlap of the dataset $⟨Drr⟩−⟨DrrRBC⟩$ with the theoretical Brownian diffusivity, confirming the RESID is linearly superimposed with the Browanian diffusivity.²⁸ More interestingly, the increase in $⟨DrrRBC⟩/DB$ is strongly correlated with the increase in $CnCFL/⟨Cn⟩non−CFL$⁠, as shown in Figs. 9(b) and 9(d). This strong linkage between the dispersion-to-margination transition and the particle-size relevant change of the RESID-to-Brownian diffusivity ratio is consistent with the margination characterization based on a kinetic theory-based analysis,⁶⁵ where the margination of microscale particles is shown to be weakened as strong Brownian effect comes into play.

Using a three-dimensional multiscale complex blood flow solver,^2,28,31 we have interrogated the long-time dispersive characteristics of rigid spherical particles with sizes across nano-to-micrometers in blood flow through microvessels. The role of the confinement ratio and the systemic hematocrit in altering the nanoparticle radial dispersion is quantitatively analyzed in terms of the radial distribution of particle concentration and particle radial diffusion rate. The effect of changing particle size on the alteration of the particle dispersive characteristics is highlighted.

In the range of parameters considered here, it is found that nanoscale particles do not marginate under various confinement effects or hematocrit levels in the same way as microscale particles do, but rather show a nonuniform radial distribution across the vessel. Increasing the confinement effect by decreasing the vessel diameter hinders the particle radial diffusivity but also enhances the equilibrium concentration of nanoscale particles in the cell free layer. Low hematocrit level (ϕ ∼ 5%) in the microvessel appears to be optimal to the radial dispersion of nanoscale particles, leading to high radial diffusion rate and near-wall concentrations being higher than the average concentration. High hematocrits (ϕ = 30%) slightly increases the near-wall concentration but meantime inhibits the dispersion of nanoscale particles in the RBC-laden region.

Microscale particles exhibit pronounced margination behavior, where at equilibrium, the microscale particles get concentrated in the cell-free layer at up to 5 times the particle average concentration in the bulk (or more than 10 times the particle concentration in the RBC-laden region). The margination propensity seems to be enhanced with the particle size. For microscale particles, the RBC-enhanced shear-induced diffusivity is dominant over the Brownian diffusivity, where the RBC-laden area shows more than 10 times higher diffusivity compared to that in the RBC-free layer. The particle-size induced alteration of particle radial diffusivity in both distribution and magnitude gives rise to margination of microscale particles in confined tubular blood flows.

The authors acknowledge the partial support from Sandia National Laboratories under Grant No. 2506X36 and the computational resource provided by the National Science Foundation under Grant No. TG-CT100012. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under Contract No. DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.