Heterogeneous partition of cellular blood-borne nanoparticles through microvascular bifurcations
Zixiang L. Liu, Jonathan R. Clausen, Justin L. Wagner, Kimberly S. Butler, Dan S. Bolintineanu, Jeremy B. Lechman, Rekha R. Rao, Cyrus K. Aidun
HHeterogeneous partition of cellular blood-borne nanoparticlesthrough microvascular bifurcations
Zixiang L. Liu,
1, 2, ∗ Jonathan R. Clausen, Justin L. Wagner, Kimberly S. Butler, DanS. Bolintineanu, Jeremy B. Lechman, Rekha R. Rao, † and Cyrus K. Aidun
1, 2, ‡ The George W. Woodruff School of Mechanical Engineering,Georgia Institute of Technology, Atlanta, GA, 30332 USA The Parker H. Petit Institute for Bioengineering and Bioscience,Georgia Institute of Technology, Atlanta, GA, 30332 USA Thermal and Fluid Processes, Sandia National Laboratories, Albuquerque, NM, 87185, USA Aerosciences Department, Sandia National Laboratories, Albuquerque, NM, 87185, USA Molecular and Microbiology, Sandia National Laboratories, Albuquerque, NM, 87185, USA Fluid and Reactive Processes, Sandia National Laboratories, Albuquerque, NM, 87185, USA (Dated: June 19, 2020) a r X i v : . [ phy s i c s . c o m p - ph ] J un bstract Blood flowing through microvascular bifurcations has been an active research topic for manydecades, while the partitioning pattern of nanoscale solutes in the blood remains relatively unex-plored. Here, we demonstrate a multiscale computational framework for direct numerical simulationof the nanoparticle (NP) partitioning through physiologically-relevant vascular bifurcations in thepresence of red blood cells (RBCs). The computational framework is established by embedding anewly-developed particulate suspension inflow/outflow boundary condition into a multiscale bloodflow solver. The computational framework is verified by recovering a tubular blood flow without abifurcation and validated against the experimental measurement of an intravital bifurcation flow.The classic Zweifach-Fung (ZF) effect is shown to be well captured by the method. Moreover,we observe that NPs exhibit a ZF-like heterogeneous partition in response to the heterogeneouspartition of the RBC phase. The NP partitioning prioritizes the high-flow-rate daughter branchexcept for extreme (large or small) suspension flow partition ratios under which the complete phaseseparation tends to occur. By analyzing the flow field and the particle trajectories, we show thatthe ZF-like heterogeneity in NP partition can be explained by the RBC-entrainment effect causedby the deviation of the flow separatrix preceded by the tank-treading of RBCs near the bifurcationjunction. The recovery of homogeneity in the NP partition under extreme flow partition ratios isdue to the plasma skimming of NPs in the cell free layer. These findings, based on the multiscalecomputational framework, provide biophysical insights to the heterogeneous distribution of NPs inmicrovascular beds that are observed pathophysiologically. ∗ [email protected] † [email protected] ‡ [email protected] . INTRODUCTION Blood is a biological fluid that conveys nutrients and wastes throughout the body. Theprimary constituents of blood are red blood cells (RBCs) that comprise ∼
45% of the systemicblood volume and plasma that contains a variety of biomolecules including clotting factorsand nutrients, etc. The exchange of materials between the blood and tissues primarily occursin the microvasculatures where blood vessels feature complex tree-like bifurcating structures.Such microvascular bifurcations allow the blood flow to be partitioned and disseminated tothe bulk of tissues for efficient mass exchanges [1].Given its significance in physiology and fluid mechanics, blood flowing through microvas-cular bifurcations has been a topic of research for many decades [2–4]. Among all the studies,substantial efforts have been dedicated to the understanding of the pattern and biomecha-nistic mechanisms of RBC partitioning through single or multi-generation bifurcations usingin vivo [5, 6], in vitro [7, 8] and theoretical [2, 9, 10] approaches. It has been well acceptedthat the RBC volume concentration (i.e., hematocrit) in the two daughter branches are of-tentimes not equally partitioned due to the particulate nature of the blood cells. Rather,the hematocrit in the high-flow-rate daughter branch often gets elevated while the otherdaughter branch receives a diluted or even depleted hematocrit level. This heterogeneouspartition of RBCs, often referred to as the ZweifachFung (ZF) effect [5, 8], is physiologicallyessential since it biomechanistically explains the reduction and heterogeneous distributionof the hematocrit in microvasculatures that has been observed physiologically [11, 12]. Uti-lizing the ZF effect, novel microfluidic devices have also been developed for cell or bioproteinseperation and other relevant biomedical applications [13, 14].Besides the large body of research focusing on the RBC partition through microvas-cular bifurcations, a few studies investigate the distribution and adhesion of other bloodspecies through bifurcating structures. The effect of size and shape of microscale particles,mimicking white blood cells, on their adhesion propensity at bifurcations has been studiedin vitro using synthetic microvascular networks [15, 16]. The model platelet distributionin blood flow through idealized, vascular bifurcations or confluences has been simulateddirectly, where the platelets shows anti-margination behaviors at the downstream of a con-fluence [17]. The distribution and adhesion of nanoparticles (NPs) in vascular bifurcationshas been interrogated in silico without explicit consideration of RBCs [18, 19].3otwithstanding current understanding of the RBC partitioning characteristics throughmicrovascular bifurcations, there is clearly a lack of knowledge regarding the partitioning be-haviors of nanoscale solutes through capillary beds subject to the influence of RBCs, partlybecause of the difficulty imposed by the multiscale nature of the problem [20]. Nonethe-less, studying the partition of nanoscale solutes through microvascular bifurcations holdsconsiderable significance biophysically and clinically. For example, the undesired heteroge-neous distribution of nanomedicine in tumor microvasculature remains to be a major hurdlethat limits the overall therapeutic efficiency [21]. Understanding the partitioning of solutesthrough microvascular bifurcations may provide possible mechanisms that could be appliedto regulate the heterogeneity of solute distribution in microcirculation. Accordingly, noveldesign of NP-based drugs with improved drug delivery efficacy could be inspired [22].The goal of this paper is to develop a multiscale computational framework for direct nu-merical simulation of NP partition through microvascular bifurcations, while quantitativelyelucidating the NP partitioning in the presence of RBCs. Towards that end, an generalindex-based manipulation strategy for the particle inflow/outflow boundary conditions em-bedded in a verified and validated multiscale blood flow solver [20, 23, 24] has been developedand validated. Using the multiscale computational framework, the partition of the blood-borne NPs through microvascular bifurcations in response to the ZF effect is found to beheterogeneous as well. The heterogeneity in the NP partition is conditional depending on theratio of the flow rates between the daughter branches. The physical mechanisms responsiblefor the heterogeneity in the NP partition is found to be associated with the flow separatrixdeviation at the bifurcation junction caused by the RBC motion.The remainder of the paper is organized as follows. In § II, the detailed computationalmethods are presented. In § III, the in vivo experimental setup is described. The validityof the computational methods is demonstrated in § IV A and § IV B. § IV C discusses thepattern and the mechanisms for NP partitioning through microvascular bifurcations. In § V,we summarize the study.
II. COMPUTATIONAL METHODS
The computational framework proposed for this study is a 3D lattice-Boltzmann (LB)based multiscale complex blood flow solver [20, 23, 25–27], augmented with a particulate4nflow and outflow boundary condition. The LB method is a well-established numericalmodel for hydrodynamics and proves to be a highly scalable method for direct numericalsimulation (DNS) of dense particle suspensions [26, 28]. Modeling of the RBC dynamicsand deformation is via a coupled spectrin-link (SL)/LB method [27, 29]. The NP suspensiondynamics are resolved via a two-way coupled lattice-Boltzmann Langevin-dynamics (LB-LD)method with both particle Brownian motion and long-range hydrodynamic interactions (HI)directly resolved and validated [20, 23]. The solver has been applied to a variety of problemspertaining to particle and biomolecule transport in blood flow[20, 24, 29–33]. The newlydeveloped open boundary conditions (BCs) are especially suitable for complex geometriessuch as microvasculatures. All modules are two-way coupled and embedded in a massivelyparallelized framework, as shown in Fig 1.
FIG. 1. The multiscale computational framework for simulating blood-borne nanoparticle parti-tioning through microvasculatures. The method directly and concurrently resolves the dynamics ofmicroscale cells and nanoscale particles and molecules suspended in blood plasma through complexmicrovascular networks, where open boundary conditions (BCs) are prescribed.
A. Lattice-Boltzmann method for fluid motion
The fluid phase, including the blood plasma and the hemoglobin within RBCs, is simu-lated using the LB method particularly designed for suspension hydrodynamics [25, 26, 34].5he method solves the discretized Boltzmann transport equation in velocity space through astreaming-and-collision procedure. In the streaming step, the fictitious fluid particles propa-gate along discrete velocity vectors in a lattice space. In the collision step, the local particledistribution function (PDF) relaxes into a local ‘Maxwellian’ equilibrium distribution. Thecollision term is linearized based on the Bhatnagar, Gross, and Krook (BGK) operator [35]for simplicity. The governing equation of the PDF reads f i ( r + ∆ t e i , t + ∆ t ) = f i ( r , t ) − ∆ tτ [ f i ( r , t ) − f ( eq ) i ( r , t )] + f Si ( r , t ) , (1)where f i is the fluid PDF, f ( eq ) i is the equilibrium PDF, r is the lattice site, e i is the discretelattice velocity, t is time, τ is the single relaxation time and f Si is a forcing source termintroduced to account for the discrete external force effect that leads to corrections to themacroscopic variables. This method has a pseudo speed of sound, c s = ∆ r/ ( √ t ), anda fluid kinematic viscosity, ν =( τ − ∆ t/ c s , where ∆ t is the time step and ∆ r is the unitlattice distance. The positivity of ν requires τ > ∆ t/
2. In the LB method, time and space aretypically normalized by ∆ t and ∆ r , respectively, such that ∆ t LB =∆ r LB =1 are employed toadvance equation 1. In the near incompressible limit (i.e., the Mach number, M a = u/c s (cid:28) f eq ) i ( r , t ) = ω i ρ [1 + 1 c s ( e i · u ) + 12 c s ( e i · u ) − c s ( u · u )] , (2)where ω i denotes the set of lattice weights defined by the LB stencil in use. The macro-scopic properties such as the fluid density, ρ , velocity, u , and pressure, p , can be recov-ered via moments of the PDFs as ρ = (cid:80) Qi =1 f i ( r , t ), ρ u = (cid:80) Qi =1 f i ( r , t ) e i and ρ uu = (cid:80) Qi =1 f i ( r , t ) e i e i − p I , respectively. Here, I is the identity tensor and pressure can be relatedto density and the speed of sound through p = ρc s . For the D3Q19 stencil adopted in thecurrent study, Q is equal to 19. Along the rest, non-diagonal, and diagonal lattice direc-tions, ω i is equal to 1/3, 1/18, and 1/36, and | e i | is equal to 0, ∆ r/ ∆ t , and √ r/ ∆ t ),correspondingly. B. Spectrin-link method for RBC membrane dynamics
RBC deformation and dynamics are through a coarse-grained spectrin-link (SL) method[37, 38] coupled with the LB method. This method has been previously validated against6xperimental results and is proved to be a superior method compared to finite-elementmethod [27]. In the SL model, the RBC membrane is modeled as a triangulated network witha collection of vertices mimicking actin vertex coordinates, denoted by { x n , n ∈ , ..., N } . TheHelmholtz free energy of the network system, E ( x n ), including in-plane, bending, volumeand surface area energy components [39], is given by E ( x n ) = E IP + E B + E Ω + E A , (3)where the in-plane energy, E IP , prescribes the membrane shear modulus through a worm-like chain (WLC) potential [40] coupled with a hydrostatic component [37]. The bendingenergy, E B , specifies the membrane bending stiffness, which is essential in characterizingthe equilibrium RBC biconcave morphology [37, 39]. The constraint energies in volume, E Ω , and area, E A , ensure the conservation of the RBC volume and area, respectively, in thepresence of external forces.The dynamics of each vertex advance according to the Newton’s equations of motion, d x n dt = v n ; M d v n dt = f SLn + f LBn + f CCn (4)where v n is the velocity of the vertex, n , and M is the fictitious mass equal to the total massof the cell divided by the number of vertices, N . The number of vertices used to discretizethe RBC membrane is N =613. This corresponds to a RBC link length of 1.5 lattice units( ∼ nm ), which has shown to be a cost effective resolution to capture both the singleRBC dynamics [27] and concentrated RBC suspension rheology [29] when coupled with theLB method. f LBn specifies the forces on the vertex due to the fluid-solid coupling. f CCn arethe forces due to cell-cell interactions. The forces due to the Helmholtz free energy basedon the SL model is determined by f SLn = − ∂E ( x n ) ∂ x n . (5)The SL method is solved by integrating equations 4 at each LB time step using a first-order-accurate forward Euler scheme in consistency with the LB evolution equation to avoidexcessive computational expense [20, 27]. C. Langevin dynamics for nanoparticle suspension dynamics
The nanoscale particle suspensions are resolved through a two-way coupled LB-LDmethod [20, 23]. This method introduces the thermal fluctuation from the Brownian par-7icles and performs momentum exchange between the particle phase and non-fluctuatingLB fluid. Different from approaches coupled with fluctuating LB method [41, 42], the cur-rent method captures the theoretical Brownian diffusivity and the long-range many-bodyHI simultaneously without introducing empirical rescaling [23]. The dynamics of the LDparticles is described through Langevin equation (LE), m d u p dt = C p + F p + S p , (6)where m is the mass of a single particle. The conservative force, C p , is determined bycalculating the directional derivatives of the total potential energy U total as C p = − dU total d r p , (7)where the U total specifies the interparticle and particle-surface interaction forces. The fric-tional force, F p , is assumed to be proportional to the relative velocity of the particle withrespect to the local viscous fluid velocity [41, 42], F p = − ζ [ u p ( t ) − u ( r p , t )] , (8)where u p denotes the particle velocity, and u ( r p , t ) is the interpolated LB fluid velocity atthe center of the particle. The friction coefficient, ζ , is determined by the Stokes drag law, ζ = 3 πµd p , where µ is the dynamic viscosity of the suspending fluid. The stochastic force, S p , satisfies the fluctuation-dissipation theorem (FDT) [43] by (cid:104) S αp,i ( t ) (cid:105) = 0; (cid:104) S αp,i ( t ) S βp,j ( t ) (cid:105) = 2 k B T ζδ ij δ αβ δ ( t − t (cid:48) ) , (9)where i, j ∈ { x, y, z } , α and β run through all the particle indices, δ ij and δ αβ are Kroneckerdeltas, δ ( t − t (cid:48) ) is the Dirac-delta function, k B is the Boltzmann constant and T is the absolutetemperature of the suspending fluid. The angle brackets denote the ensemble average overall the realizations of the random variables. Since only the time scales equal to and greaterthan the Brownian diffusion time scale is of interest, this study solves the over-dampeddiscretized LE as suggested in [20, 23]. D. Open boundary strategies for particulate suspension inflow and outflow
Periodic boundary condition (P-BC), due to its computational efficiency, has been widelyused for simulations of blood flows through straight tubes/channels [20, 31, 33, 44] or artifi-cial bifurcation structures [10, 45]. Physiologically, microvascular bifurcations often involve8ne inlet but multiple outlets where flow conditions are non-periodic. To overcome themethodological limitation, several open boundary strategies have been developed for inflowand outflow of particle or cell suspensions in the context of dissipative particle dynamics(DPD) [46, 47] or immersed-boundary (IB)/LB method [48]. Here we propose a general par-ticulate suspension inflow and outflow boundary condition (PSIO-BC), which in principle isapplicable to any index-based particle suspension methods. Applying the PSIO-BC to boththe LB-SL method [27] and the LB-LD method [23], we simulate the RBC-NP suspensionthrough a model microvascular bifurcation (Fig 2a) where the number of particles initiallyincreases and eventually plateaus when the system reaches equilibrium (Fig 2b).
FIG. 2. Demonstration of the particulate suspension inflow and outflow boundary condition (PSIO-BC). The schematic illustration of the PSIO-BC applied to a NP-RBC bidisperse suspension flowingthrough a two-generation microvascular bifurcation (a). The time course of the particle numbersimulated in the bifurcation (b), where at equilibrium the particle number becomes quasi-conserved. . Fluid phase The Dirichlet BCs are applied to prescribe the velocity or pressure of the flow fieldat the inlet and outlet. In the context of LB method, we employ the regularization BCprocedure [49, 50] to construct the incoming unknown PDFs based on prescribed macroscopicproperties at the inlet and outlet as f Ri = f ( eq ) i + f ( neq ) i = f ( eq ) i ( ρ, u ) + ω i c s Q i : Π (1) , (10)where the tensor Q i is defined as Q i = e i e i − c s I and the stress tensor Π (1) is evaluatedbased on the non-equilibrium components of the outward-pointing PDFs (in the oppositedirections of the unknown PDFs) [50]. Single phase regions are arranged to the upstream ofthe particle introduction region as well as to the downstream of the particle removal regions,as shown in Fig 2a. This way, the velocity or pressure can be prescribed within the singlephase region such that simple analytical solutions can be used. The aspect ratio of the singlephase regions is chosen to be roughly one to allow smooth transitions of the properties fromthe single-phase to multi-phase regions or vice versa.
2. Particle phase
The introduction of particles at inlet is realized by inserting a particle introduction sec-tion downstream proximal to the single phase region, as depicted in Fig 2a. Inside theseeding section, the number of particles is conserved to control the particle concentrationand number flow rate introduced to the downstream. Similar to other approaches [46–48],particles in this section periodically leave and re-enter the domain, while at the exit newparticles are introduced to the downstream by duplicating the particle (and its propertiessuch as mesh coordinates, velocity and mechanical properties, etc) that are wrapped to theentrance (upstream end) of the particle introduction region. Besides, random motions canbe introduced to the particles within the introduction region to reduce the periodic pat-terning of the particle distribution. Through these operations, new particles at controlledconcentration and flow rate are introduced to the downstream with the surrounding flowfield being well-perturbed.Similarly, the particle removals are realized by embedding a particle removal region up-stream proximal to a single phase region where the outlet fluid BCs are imposed. At the10
IG. 3. Schematics for the general particulate inflow and outflow boundary condition (PSIO-BC).The particle array represents the array storing the particle information according to the particleindex, which increases from 1 to N max from left to right. Here for example N max = 5 at time t particle removal region, the particle information will be deleted at the local processor (as-suming multi-processor parallelization is employed). Meanwhile, the particle index and thetotal number of removed particles at each time step will be recorded and communicated tothe processors in charge of the particle introduction. To minimize the memory allocation for11he particle array and perform minimal particle index manipulation, the indices of newly in-troduced particles are assigned conditionally according to the number of inflow and outflowparticles ( N pin and N pout ). The schematics for the PSIO-BC operations are shown in Fig 3. Algorithm 1
The Particulate Suspension Inflow and Outflow (PSIO) Algorithm
Require: i, j, i p , j p , k p ∈ Z + ∪ { } ,1: Count the number of inflow particles N pin
2: Store properties (location, velocity, etc) of inflow particles P pin ( i ) | i ∈{ ,...,N pin − }
3: Count the number of outflow particles, N pout
4: Store indices of outflow particles, I pout ( i ) | i ∈{ ,...,N pout − }
5: Update total particle number: N pt ← N pt − + N pin − N pout if N pin = N pout for i ← N pout −
18: Assign P pin ( i ) to particle I pout ( i )9: end for else if N pin > N pout for i ← N pout − P pin ( i ) to particle I pout ( i )13: end for
14: Allocate memory for particles N pt − , ..., N pt − for i ← N pt − to N pt − i p ← i − N pt − + N pout
17: Assign P pin ( i p ) to particle i end for else if N pin < N pout j p ← for i ← N pin − i p ← I pout ( i )23: while i p ≥ N pt j p ← j p + 125: i p ← I pout ( j p )26: end while j p ← j p + 128: Assign P pin ( i ) to particle i p end for k p ← for i ← N pt to N pt − − for j ← N pout − if i (cid:54) = I pout ( j )34: I ps ( k p ) ← i k p ← k p + 136: end if end for end for for i ← j p to N pout − k p ← i − j p
41: Assign P pin ( I ps ( k p )) to particle I pout ( i )42: Free the memory of particle I ps ( k p )43: end for end if Specifically, if N pin = N pout , the inflow particle properties are assigned to the indices12orresponding to the outflow particles, as indicated in Fig 3a. If N pin > N pout , in addition tothe same procedure as above, the data structure for storing all the particle information in thecomputational domain (referred to as particle array) will be extended to accomodate for theincremental particles, shown in Fig 3b. If N pin < N pout , the particle array will be truncated toreflect the reduction of the total particle number. When the truncated particle indices areall associated with outflow particles (denoted as Scenario I in Fig 3c), the outflow particlesand their memory allocation will be simply eliminated, and the inflow particles will take theindices of the outflow particles starting from the lowest index. When the truncated particleindices contain particles remaining in the computational domain (denoted as Scenario IIin Fig 3c), in addition to the operations in Scenario I, the remained particle informationwill be assigned to the unoccupied outflow particle indices that are within the range of theupdated particle array. The detailed procedure for the PSIO-BC algorithm is listed in thepseudo-code shown in Algorithm 1. E. Reconstruction of physiologic microvascular bifurcations
Microvascular networks are comprised of many generations of single bifurcations connect-ing all the way from arterioles or venules to capillaries [51]. Each single bifurcation followscertain physiologic laws that constrain the geometric relation between parent and daughterbranches. The radii of parent and daughter branches have been found to obey the Murray’slaw [52] stated as, R = R + R , (11)where R is the radius of the parent branch; R and R are the radii of the daughter branches.Complying to this simple geometric relationship leads to a more uniform distribution of vesselwall shear stress and consequently minimizes the power consumption [53]. The bifurcatingangles of two daughter branches also satisfy certain physiologic optimality relating to thevessel radius based on Zamir’s law [54],cos α = R + R − R R R ; cos α = R + R − R R R , (12)where α and α are the bifurcating angles of the two daughter branches with respect to theaxial direction of the parent branch. Obeying Zamir’s law minimizes the lumen volume ofa single bifurcation [54]. In addition, we also assume the axes of three branches of a single13ifurcation are in-plane as suggested previously [55]. The vessel length is assumed to beproportional to the vessel diameter with certain randomness [56]. The cross-section shapeof the branches are assumed to be circular.Based on the above physiologic constraints, microvascular bifurcations of multiple gen-erations of branches can be reconstructed as shown in Fig 4, provided the radius of theparental branch, r , and the size ratio of two daughter branches, η D = R R , being prescribedrandomly within the physiological range [55]. FIG. 4. Demonstration of three reconstructed microvascular bifurcations of physiologic relevance.
F. Coupling schemes
The coupling between fluid and RBC is accomplished through the Aidun-Lu-Ding (ALD)fluid-solid interaction scheme [25–27]. For the suspension dynamics of NP subjected tohydrodynamics and long-ranged many-body HI, the LD particle and the non-fluctuating LBfluid phase are two-way coupled, as discussed and verified in previous studies [20, 23].The RBC-RBC and RBC-vessel wall interactions that specify f CCn are based on the sub-grid contact functions originally formulated in Ding and Aidun [57] and later improved byMacMeccan et al. [58] and Clausen et al. [59]. It prevents the RBC from overlapping whencell-cell/wall separation is below one LB lattice spacing. In this model, the lubrication termis replaced with an exponential contact function to avoid numerical instability driven by thesingular nature of the lubrication hydrodynamics and the discrete nature of the interparticleseparation calculation, as explained in detail in [58, 60].The NP-RBC contact model that provides U total is based on the Morse potential that canbe calibrated to match experimentally measured surface interaction forces [61, 62]. Due tothe lack of statistics for actual NP-RBC short-distance interaction, this study employs the14easured cell-cell interaction potential [61] for the NP-RBC interactions. The model pa-rameters are specified according to [62] but with a cut-off distance to exclude the attraction.The short-distance NP-NP interaction is neglected due to the extreme dilution of the NPconcentration ( (cid:28) III. IN VIVO EXPERIMENT
To validate the numerical methods, the blood flow through a microvascular bifurcationwithin the chicken chorioallantoic membrane (CAM) in vivo model has been characterizedthrough microscale particle image velocimetry ( µ PIV). The CAM model is a highly vascu-larized extraembryonic membrane that is readily accessible for intravital imaging. In-housesynthesized mesoporous silica nanoparticles (MS-NP) were used as tracking particles.The µ PIV system is shown in Fig 5 (left), which mainly consists an upright microscope(Zeiss Axio Examiner Z.1) modified with a heated stage operated at 63 magnification toimage the CAM model (10 days post fertilization). The microscope is fitted with a waterimmersion lens (LDC-APO Chromat 63 1.15W) to facilitate imaging of deep vessels. Imagingwas performed using a high-speed camera (Hamamatsu Orca Flash 4.0) operated at 142frames/sec, enabling a resolution of 0.1 microns/pixel. A broadband light source was usedalong with careful image and PIV processing methods through the LaVision PIV softwarepackage DaVis v8.4 to mitigate image blur effects especially at peak velocities.The CAM model is depicted in Fig 5 (middle). The ex-ovo avian embryos were handledaccording to published methods [63] with all experiments conducted following an institu-tional approval protocol (11-100652-T-HSC). Embryos were utilized for experimentation 7days post removal from shells (day 10 post fertilization). A dosage of 50 µ g MS-NPs withsize ∼
130 nm were injected into tertiary veins of the CAM. After injection, embryos wereplaced in a customized avian embryo chamber for imaging.Applying the preceding methodology, a greyscale image of a vascular bifurcation is shownin Fig 5 (right). The image location was chosen such that the vessel branches remainrelatively flat over the entire field-of-view. To obtain the dynamic velocity field, 5000 imageswere acquired in total. A sliding sum-of-correlation method was used to calculate vectorfields [64]. In § IV B, the mean velocity field obtained based on the NP flow is comparedwith the computed counterparts through the current multiscale computational framework.15
IG. 5. The benchtop microscale particle image velocimetry ( µ PIV) system (left), the chickchorioallantoic membrane microvasculature model (middle) and the grey-scale visualization of theflow field of a selected region of interest (right).
IV. RESULTS AND DISCUSSION
In this section, we first verify the PSIO-BC by recovering a tubular blood flow andcomparing the results with that obtained from the typical periodic boundary condition.Furthermore, the multiscale computational framework is further validated by reconstructingthe CAM bifurcation flow in silico and comparing simulation against the experimental mea-surement. As follows, the NP partitioning pattern in the presence of RBCs is interrogatedquantitatively by simulating a bidisperse RBC-NP suspension through single physiologically-relevant microvascular bifurcation under different flow partition ratios with or without geo-metric asymmetry. Physical mechanisms leading to specific NP partitioning characteristicsare discussed in detail.
A. Recovering tubular flows
Flow through a straight tube with or without RBCs is simulated with both PSIO-BC andP-BC. The tube has a diameter of 2 R =20 µm and a length of L =80 µm . The fluid viscosityis selected as µ =1.2 cP. A hematocrit of φ =30% is selected for the multiphase simulations.The RBC has a shear modulus of G =0.0063 dynes/cm . For the P-BC simulations, a constantpressure gradient, − dpdz =240 P a/mm , is applied to drive the flow (from left to right according16o Fig 6a) such that the wall shear rate is ˙ γ w =1000 s − . For the single phase simulationusing the PSIO-BC, parabolic velocity profiles, u ( r )= | dpdz | ( R − r )4 µ , are prescribed at both theinlet and outlet. For RBC case through PSIO-BC, inlet and outlet pressures are prescribedat the outmost ends of the single phase regions. The inlet pressure is p i = p + | dpdz | ( L + L BC )and the outlet pressure is p o = p , where p is chosen as the atmospheric pressure and L BC =60 µm that accounts for the extra tube sections for PSIO-BC. FIG. 6. Normalized axial flow velocity plotted against vessel radial locations for a tubular flowwith or without red blood cells, computed with the P-BC or the PSIO-BC.
The velocity profiles for the tubular flow simulations using both PSIO-BC and P-BCare plotted in Fig 6b. The velocity profiles for the cases with φ =30% are averaged over 1017ifferent axial locations at one time point during equilibrium. It is shown that the results viaPSIO-BC agrees well with that via P-BC, showing the validity of the PSIO-BC algorithm. B. Comparing with in vivo measurements
To further validate the current multiscale computational framework, we simulate a NP-RBC suspension flow through a CAM microvascular bifurcation (see supplemental video 1)and compared the results against experimental measurements through µ PIV. An arterialbranch in the CAM model is selected for the study due to its relatively big size and flatorientation, which are suitable for measurement. The parent branch has a radius of R =20 µm and the two daughter branches have radii of R =7.5 µm and R =17.5 µm respectively,shown in Fig 7a. Slight pulsatility of the flow is observed with a frequency of f = 3 Hz .The Womersley number for this problem can be estimated as α = (cid:112) ρf R /µ ≈ . (cid:28) µ PIV measurement are demonstratedin Fig 7a on the left. The mean streamlines are overlaid, indicating the flow direction.We simulate the NP-RBC suspension flow through the same reconstructed geometry usingthe multiscale computational method, where Dirichlet velocity BCs based on the in vivomeasurements are applied at the inlet and outlet using the PSIO-BC. The hematocrit isassumed to be 20% which is an average physiological value for CAM blood flows [65]. TheNP number concentration is set to 10 /ml for statistic purposes. The RBC size is selected tobe the same as human RBC size, while the rigidity (shear modulus) of the RBC membrane iselevated by three folds to be consistence with low deformability observed in avian RBCs [66].As shown in Fig 7, the fluid velocity distribution reproduced by simulation qualitativelymatches well with the PIV results. 18 IG. 7. (a) Comparisons of the mean velocity distribution between the experiment and simulationthrough a CAM vascular bifurcation. The velocity contour for the simulation is an snapshotafter the system reaches equilibrium (i.e., the number of particles reach plateau). (b) The radialdistribution of the mean axial velocity at the A-A cross-section annotated in (a). The asymmetryof the experimental curve may be related to the non-circular shape of the actual cross-section area.
The velocity profiles across the diameter of the larger daughter branch at the A-A location(denoted) are further plotted in Fig 7b. The velocity profile in the simulation is obtainedby tracking the streamwise velocity of the NP flowing through the A-A cross-section at19quilibrium to mimic the µ PIV process. The mean velocity near the wall has a finite numberfor both the simulation and experiment as a result of the sliding nature of the NP adjacentto the wall. The velocity profile in the RBC-laden region shows overall good agreementsbetween experiment and simulation.
C. NP partition through microvascular bifurcations
To understand how NPs partition at the junction of a microvascular bifurcation, a tradestudy of NP-RBC suspension flow through single microvascular bifurcations is performed,taking into account the asymmetry in both the flow partitioning and vascular geometry.To control the partitioning of the bulk fluid, a flow partition ratio (FPR) is defined as η FQ,i = Q i /Q where Q i ( i =1 or 2) denotes the suspension flow rate at one daughter branchand Q is the suspension flow rate at the parental branch. Similarly, a particle partition ratio(PPR) is defined as η pQ,i = N i /N , where N denotes the number of particles entering thecorresponding branch. The PPR is used to quantitatively describe the partition pattern ofthe particle phase (RBC and NP) under different FPR. For the following studies, the parentalbranch has a size of 2 R =20 µm , a wall shear rate of ˙ γ =2000 s − and a hematocrit of φ =20%,reflecting physiologic conditions of arterioles [3, 51]. The Dirichlet velocity boundaries areapplied to both the inlet and outlet, where the Poiseuille flow velocity profiles are imposed atthe single phase region pertaining to specific η FQ,i . The vascular geometric effect is controlledby the daughter branch size ratio, η D . Provided R and η D , the exact size and bifurcatingangle of the two daugher branches are uniquely constrained by the physiologic laws discussedin II E.
1. Symmetric bifurcations
We first study the NP-RBC suspension flow partitioning through a geometrically sym-metric bifurcation by setting η D = 1, as shown in Fig 8a. As a control, when there is noRBC in the bifurcation, the η pQ of NP changes linearly as the η FQ changing from 0.1 to 0.9,as shown in the inset of Fig 8b. This partitioning pattern is termed homogeneous parti-tion, as the flow of the particle phase follows the fluid phase while branching towards thedownstream. 20 IG. 8. (a) Snapshots of RBC-NP suspension flow through symmetric bifurcations under varioussuspension flow flux partition ratio, η FQ . (b) Particle flow partition ratios, η pQ , plotted againstsuspension flow partition ratio, η FQ , for symmetric microvascular bifurcations. The inset shows theresults without the presence of RBCs. η pQ of the RBCphase exhibits a sigmodal trend with respect to η FQ , as shown in Fig 8b. This trend agreeswell with both the experimental measurement and theoretical correlation by Pries et al. [6],indicating the classic ZF effect is well captured by the current computational methods. Suchheterogeneous partition of RBC is most prominent at η FQ = 0 . η pQ for NP also deviates from the linear line, suggestingthe NP partitioning also becomes heterogeneous. However, such heterogeneity in the NPpartitioning seems to be conditional as the η pQ falls back on the linear line at about η FQ =0.1or 0.9 while complete phase separation of RBC occurs.
2. Asymmetric bifurcations
To introduce the effect of geometric asymmetry, we set η D = 0 . et al. [55]. The resulting geometry of the asymmetricbifurcation is shown in Fig 9a. Similar to the symmetric bifurcation case, NP shows mainlyhomogeneous partition when RBCs are not present, as shown in the inset of Fig 9b. However,due to geometric asymmetry, the η pQ of NP is deviated from the linear curve at η FQ = 0 . et al. [6]. Different from the symmetric bifurcation, the heterogeneity of RBC partitioneven occurs during equal flow partitioning ( η FQ = 0 .
5) as another manifestation of the effectof the geometric asymmetry. The enhanced heterogeneity of the RBC partitioning furtherelevates the heterogeneity in the NP partitioning compared to the symmetric bifurcationcase. Again, the partitioning of NP recovers homogeneity at extreme flow partition ratioswhen severe phase separate occurs. 22
IG. 9. (a) Snapshots of RBC-NP suspension flow through asymmetric bifurcations under varioussuspension flow flux partition ratio, η FQ . (b) Particle partition ratio, η pQ , plotted against suspensionflow partition ratio, η FQ , for asymmetric microvascular bifurcations. The inset shows the resultswithout the presence of RBCs. . Mechanisms for conditional heterogeneity in NP partition Next, we take a detailed look at the particle trajectory and flow structure of the asym-metric bifurcation case to understand the underlying mechanisms governing the conditionalheterogeneity in the NP partition through microvascular bifurcations. Specifically, we lookat two flow partition ratios, η FQ = 0 . η FQ = 0 . η pQ = 0 .
83 being higher than η FQ . The RBC tank-treading motionat the bifurcation junction will further alter the local flow field. As shown in Fig 10c,at the condition of η FQ = 0 .
7, without RBCs the flow separatrix is well aligned with thegeometric dividing line at the junction; with RBCs lingering and tank-treading towards theupper branch, the flow separatrix tends to locally shift towards the lower branch, whichlocally leads to the recruitment of more NPs towards the upper branch. This RBC-inducedentrainment effect seems to be the mechanism that causes the heterogeneity in the NPpartition when RBCs are present.When the recovery of homogeneous partitioning occurs at η FQ = 0 .
9, a certain amountof NPs can still flow into the RBC-free branch, as shown in Fig 10b. As indicated by thetrajectories of NP throughout the bifurcation, the majority of the NPs that enter the lowerbranch come from the cell-free layer (CFL) of the parental branch. Unlike the deformableRBCs that tend to concentrate to the center of the tube as a result of the lift force exertedby the vessel wall [67], NPs tend to disperse throughout all the radial locations due to thesevere Brownian effect [20, 33]. In other words, the NP phase can always enter the cell-free branch so long as the plasma skimming occurs, which will reasonably contribute to thehomogeneity in the NP partitioning. As further shown in Fig 10b, the flow separatrix at η FQ = 0 . IG. 10. (a-b) The cell distribution, NP trajectories and a cross-section of interest within anasymmetric microvascular bifurcation. (c) The contours of the Y velocity component at the cross-section adjacent to the bifurcation junction. The left column shows the case without the presenceof RBCs. The other six columns corresponds to six time points when a RBC sliding towards theupper branch. The black dashed line denotes the zero velocity line, which shows the flow separatrixthat divides the flow entering two daughter branches.
V. SUMMARY
In this work, we have developed a multiscale computational framework to directly simu-late the NP distribution and partition through the microvascular bifurcations in the presenceof RBCs. A general particulate suspension inflow and outflow boundary condition has beendeveloped and embedded in the computational framework, and later applied to both theRBC and NP suspension dynamics. To validate the model, in vivo measurements of the NPflow through a CAM microvascular bifurcation have been performed. The current computa-tional framework is able to quantitatively reconstruct the detailed physiological bifurcationflow by directly simulating both the RBC and NP suspensions. The velocity distribution forthe NP phase through simulation matches well with the experimental measurement, showingthe accuracy and validity of the present multiscale computational framework.Applying the multiscale computational framework, we further interrogate the RBC-NPsuspension flow through idealized single bifurcations under physiologically relevant geom-etry and flow conditions. The classic Zweifach-Fung effect is found to be well capturedby the method. Moreover, we find the partition of the blood-borne NP through microvas-cular bifurcations in response to the ZF effect also shows a substantial heterogeneity. Theheterogeneity in the NP partition is however conditional depending on the level of dispropor-tionality in the flow partitioning between the daughter branches. The physical mechanismsresponsible for the heterogeneity in the NP partition has been analyzed in terms of the flowstructure and the particle trajectories. It is found that the heterogeneity in the NP partitionin the presence of RBCs seems to be associated with the flow separatrix deviation at thebifurcation junction caused by the RBC tank-treading motion. The recovery of homogeneity26n the NP partitioning however has to do with the plasma skimming of the NPs in the CFLaccompanied with the normalization of flow separatrix at the extreme flow partition.This study presents a numerical tool that can be used to directly analyze the blood-borneNP distribution within microcirculation that often comprises complex branching structures(see the supplemental video 2). The results imply the importance of the presence of theRBC phase in directing the blood solute distributions and provides possible physical mecha-nisms that contribute to the heterogeneous distribution of blood solute in microvasculaturesobserved clinically [21]. Such insights can provide valuable information for a upscaling toa network model of the full vasculature [56] and may help the design of nanocarriers withimproved drug delivery efficacy [22]. It is also noteworthy that the proposed particulatesuspension inflow and outflow boundary condition is general and can be applied to otherapplications such as modeling arterial thrombosis [68] or other reacting flows that ofteninvolve consumption of particles or species.
ACKNOWLEDGEMENTS
The authors acknowledge the partial funding from Sandia National Laboratories madeavailable through contract number 2506X36. The computational resource granted byXSEDE under contract number TG-CT100012 is acknowledged. This work is supported bythe Laboratory Directed Research Development Program at Sandia National Laboratories.Sandia National Laboratories is a multimission laboratory managed and operated by Na-tional Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary ofHoneywell International Inc. for the U.S. Department of Energys National Nuclear SecurityAdministration under contract de-na0003525. This paper describes objective technical re-sults and analysis. Any subjective views or opinions that might be expressed in the paper donot necessarily represent the views of the U.S. Department of Energy or the United StatesGovernment. [1] L. Sherwood,
Human physiology: from cells to systems (Cengage learning, 2015).[2] A. R. Pries, T. W. Secomb, P. Gaehtgens, and J. Gross, Circulation research , 826 (1990).[3] A. R. Pries and T. W. Secomb, in Microcirculation (Elsevier, 2008) pp. 3–36.
4] T. W. Secomb, Annual Review of Fluid Mechanics , 443 (2017).[5] K. Svanes and B. Zweifach, Microvascular Research , 210 (1968).[6] A. Pries, K. Ley, M. Claassen, and P. Gaehtgens, Microvascular research , 81 (1989).[7] G. Bugliarello and G. C. Hsiao, Science , 469 (1964).[8] Y.-C. Fung, Microvascular research , 34 (1973).[9] Z.-Y. Yan, A. Acrivos, and S. Weinbaum, Microvascular research , 17 (1991).[10] P. Balogh and P. Bagchi, Physics of Fluids , 051902 (2018).[11] P. C. Johnson, American Journal of Physiology-Legacy Content , 99 (1971).[12] G. Kanzow, A. Pries, and P. Gaehtgens, International journal of microcirculation, clinicaland experimental , 67 (1982).[13] S. Yang, A. ¨Undar, and J. D. Zahn, Lab on a Chip , 871 (2006).[14] R. Fan, O. Vermesh, A. Srivastava, B. K. Yen, L. Qin, H. Ahmad, G. A. Kwong, C.-C. Liu,J. Gould, L. Hood, et al. , Nature biotechnology , 1373 (2008).[15] N. Tousi, B. Wang, K. Pant, M. F. Kiani, and B. Prabhakarpandian, Microvascular research , 384 (2010).[16] N. Doshi, B. Prabhakarpandian, A. Rea-Ramsey, K. Pant, S. Sundaram, and S. Mitragotri,Journal of Controlled Release , 196 (2010).[17] C. B¨acher, A. Kihm, L. Schrack, L. Kaestner, M. W. Laschke, C. Wagner, and S. Gekle,Biophysical journal , 411 (2018).[18] J. Tan, S. Shah, A. Thomas, H. D. Ou-Yang, and Y. Liu, Microfluidics and nanofluidics ,77 (2013).[19] S. Sohrabi, S. Wang, J. Tan, J. Xu, J. Yang, and Y. Liu, Journal of biomechanics , 240(2017).[20] Z. Liu, Y. Zhu, R. R. Rao, J. R. Clausen, and C. K. Aidun, Comput. Fluids , 609 (2018).[21] R. K. Jain and T. Stylianopoulos, Nature reviews Clinical oncology , 653 (2010).[22] S. Wilhelm, A. J. Tavares, Q. Dai, S. Ohta, J. Audet, H. F. Dvorak, and W. C. Chan, Naturereviews materials , 1 (2016).[23] Z. Liu, Y. Zhu, J. R. Clausen, J. B. Lechman, R. R. Rao, and C. K. Aidun, InternationalJournal for Numerical Methods in Fluids , 228 (2019).[24] Z. Liu, J. R. Clausen, R. R. Rao, and C. K. Aidun, J. Fluid Mech. , 636 (2019).[25] C. K. Aidun, Y. N. Lu, and E. J. Ding, J. Fluid Mech. , 287 (1998).
26] C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. , 439 (2010).[27] D. A. Reasor, J. R. Clausen, and C. K. Aidun, Int. J. Numer. Methods Fluids , 767 (2012).[28] J. R. Clausen, D. A. Reasor, and C. K. Aidun, Comput. Phys. Commun. , 1013 (2010).[29] D. A. Reasor, J. R. Clausen, and C. K. Aidun, J. Fluid Mech. , 497 (2013).[30] J. Reasor, D. A., M. Mehrabadi, D. N. Ku, and C. K. Aidun, Ann. Biomed. Eng. , 238(2013).[31] M. Mehrabadi, D. N. Ku, and C. K. Aidun, Phys. Rev. E , 023109 (2016).[32] F. Ahmed, M. Mehrabadi, Z. Liu, G. A. Barabino, and C. K. Aidun, J. Biomech. Eng. ,061013 (2018).[33] Z. Liu, J. R. Clausen, R. R. Rao, and C. K. Aidun, Physics of Fluids , 081903 (2019).[34] C. K. Aidun and Y. Lu, J. Stat. Phys. , 49 (1995).[35] P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. , 511 (1954).[36] M. Junk and W.-a. Yong, Asymp. Anal. , 165 (2003).[37] D. A. Fedosov, B. Caswell, and G. E. Karniadakis, Biophys. J. , 2215 (2010).[38] I. V. Pivkin and G. E. Karniadakis, Phys. Rev. Lett. , 118105 (2008).[39] M. Dao, J. Li, and S. Suresh, Mater. Sci. Eng. C , 1232 (2006).[40] C. Bustamante, Z. Bryant, and S. B. Smith, Nature , 423 (2003).[41] P. Ahlrichs and B. D¨unweg, Int. J. Mod. Phys. C , 1429 (1998).[42] P. Ahlrichs and B. D¨unweg, J. Chem. Phys , 8225 (1999).[43] R. Kubo, Rep. Prog. Phys. , 306 (1966).[44] H. Zhao and E. S. G. Shaqfeh, Phys. Rev. E , 061924 (2011).[45] P. Balogh and P. Bagchi, Journal of Computational Physics , 280 (2017).[46] V. W. A. Tarksalooyeh, G. Z´avodszky, B. J. van Rooij, and A. G. Hoekstra, Computers &fluids , 312 (2018).[47] G. Li, T. Ye, and X. Li, Journal of Computational Physics , 109031 (2020).[48] K. Lykov, X. Li, H. Lei, I. V. Pivkin, and G. E. Karniadakis, PLoS computational biology (2015).[49] J. Latt and B. Chopard, Mathematics and Computers in Simulation , 165 (2006).[50] J. Latt, B. Chopard, O. Malaspinas, M. Deville, and A. Michler, Physical Review E ,056703 (2008).[51] Y.-c. Fung, Biomechanics: circulation (Springer Science & Business Media, 2013).
52] C. D. Murray, Proceedings of the National Academy of Sciences of the United States ofAmerica , 207 (1926).[53] T. F. Sherman, The Journal of general physiology , 431 (1981).[54] M. Zamir, The Journal of general physiology , 837 (1978).[55] M. Zamir, S. Wrigley, and B. Langille, The Journal of general physiology , 325 (1983).[56] T.-R. Lee, J.-A. Hong, S. S. Yoo, and D. W. Kim, International journal for numerical methodsin biomedical engineering , e2981 (2018).[57] E. J. Ding and C. K. Aidun, J. Stat. Phys. , 685 (2003).[58] R. M. MacMeccan, J. R. Clausen, G. P. Neitzel, and C. K. Aidun, J. Fluid Mech. , 13(2009).[59] J. R. Clausen, D. A. Reasor, and C. K. Aidun, J. Fluid Mech. , 202 (2011).[60] J. R. Clausen and C. K. Aidun, Phys. Fluids , 123302 (2010).[61] B. Neu and H. J. Meiselman, Biophys. J. , 2482 (2002).[62] Y. Liu, L. Zhang, X. Wang, and W. K. Liu, Int. J. Numer. Methods Fluids , 1237 (2004).[63] H. S. Leong, N. F. Steinmetz, A. Ablack, G. Destito, A. Zijlstra, H. Stuhlmann, M. Manch-ester, and J. D. Lewis, Nature protocols , 1406 (2010).[64] J. Y. Lee, H. S. Ji, and S. J. Lee, Physiological measurement , 1149 (2007).[65] M. Maibier, B. Reglin, B. Nitzsche, W. Xiang, W. W. Rong, B. Hoffmann, V. Djonov, T. W.Secomb, and A. R. Pries, American journal of physiology-heart and circulatory physiology , H913 (2016).[66] P. Gaehtgens, F. Schmidt, and G. Will, Pfl¨ugers Archiv , 278 (1981).[67] M. Abkarian, C. Lartigue, and A. Viallat, Physical review letters , 068103 (2002).[68] D. Kim, C. Bresette, Z. Liu, and D. N. Ku, APL bioengineering , 041502 (2019)., 041502 (2019).