Inflow/Outflow Boundary Conditions for Particle-Based Blood Flow Simulations: Application to Arterial Bifurcations and Trees

When blood flows through a bifurcation, red blood cells (RBCs) travel into side branches at different hematocrit levels, and it is even possible that all RBCs enter into one branch only, leading to a complete separation of plasma and RBCs. To quantify this phenomenon via particle-based mesoscopic simulations, we developed a general framework for open boundary conditions in multiphase flows that is effective even for high hematocrit levels. The inflow at the inlet is duplicated from a fully developed flow generated in a pilot simulation with periodic boundary conditions. The outflow is controlled by adaptive forces to maintain the flow rate and velocity gradient at fixed values, while the particles leaving the arteriole at the outlet are removed from the system. Upon validation of this approach, we performed systematic 3D simulations to study plasma skimming in arterioles of diameters 20 to 32 microns. For a flow rate ratio 6:1 at the branches, we observed the “all-or-nothing” phenomenon with plasma only entering the low flow rate branch. We then simulated blood-plasma separation in arteriolar bifurcations with different bifurcation angles and same diameter of the daughter branches. Our simulations predict a significant increase in RBC flux through the main daughter branch as the bifurcation angle is increased. Finally, we demonstrated the effectiveness of the new methodology in simulations of blood flow in vessels with multiple inlets and outlets, constructed using an angiogenesis model.


Introduction
Blood is a biological fluid that delivers nutrients and oxygen to living cells and removes their waste products. The two major components of whole blood are red blood cells (RBCs) and plasma, i.e., RBCs constitute approximately 40% of the total blood volume, plasma around 55%, while the rest is taken up by white blood cells (WBCs) and platelets. Lateral migration of RBCs takes place in blood flow, leading to the formation of two phases, i.e., a core consisting mainly of RBCs and a cell-free plasma layer adjacent to the vessel wall where the platelets tend to concentrate [1,2]. WBCs, which are larger and more rigid than RBCs, migrate toward the vessel wall through a process called margination [3]. The tendency of RBCs to concentrate at the vessel center also leads to plasma skimming-an asymmetric distribution of RBCs and plasma between two daughter branches when blood flows through a microvascular bifurcation. The RBCs prefer a daughter branch with higher flow rate leaving very few RBCs (even reaching zero) flowing into lower flow rate daughter branch [4].
Blood-plasma separation in bifurcations has been extensively investigated in the past few decades, and it is generally believed that three factors, feed hematocrit [5][6][7], size of parent channel [8,9] and flow rate ratio of daughter branches [1], mainly determine the degree of plasma skimming that will occur [8]. Studies of blood flow through bifurcations have revealed significant variability for a complete RBC separation from the whole blood (all-or-nothing phenomenon). The theoretical critical flow rate ratio between the daughter branches for predicting such phenomenon is approximately 2.5:1 [5]. However, more recent experimental measurements showed that for this flow rate ratio only 88.7% of RBCs enter into the higher flow rate daughter branches [10]. This raises the question as to what ratio value is more meaningful in determining total blood-plasma separation.
Computational modeling and simulations can help us to investigate this issue. In past decades, numerical modeling of blood flow in capillaries has attracted increasing attention [11,12]. For example, dynamic simulations can model how blood flow behaves in microfluidic channels [13][14][15][16][17][18] and predict human blood viscosity in silico [19]. Different cell models have also been employed for various qualitative and quantitative interpretations as well as predictions of biomechanical properties of RBCs with hematological diseases [20][21][22][23]. Examples include dynamic cell deformability for various stages of malaria-infected RBCs [24][25][26][27][28] and vaso-occlusion phenomena in sickle cell anemica [23]. However, most of these blood flow simulations were performed in systems with periodic boundary conditions (PBCs) along the flow direction, whereas very few studies so far have reported simulations of non-periodic flow [29][30][31]. In a previous study, we simulated the blood-plasma separation for healthy and diseased blood in microfluidic channels with geometrically symmetric bifurcation and confluence to satisfy the periodic flow assumption along the flow direction [32]. However, for a simulation study of plasma skimming in capillary bifurcations, the blood flow properties such as velocity and pressure fields differ drastically at the inlet and outlet regions. Therefore, the choice of PBCs is inappropriate for general cases, especially in arterial trees, and hence a new open (nonperiodic) boundary is required; this is a non-trivial issue, especially for particle-based Lagrangian methods.
For an open boundary system, the velocity profile at the inlet is generally specified, whereas the outflow profiles are rarely known. For a single-phase system, the inflow condition could be simply obtained by extending the inflow length so that the flow becomes fully developed at the inlet. However, for multiphase systems the inflow conditions even for a fully developed flow are unknown-a situation similar to turbulent inflow in single phase. For example, for a blood flow, the flow and viscous properties as well as the cell-free layer (CFL) distribution in arteries differ greatly with change in hematocrit level and shear rate. Thus, the inflow length should be long enough to generate inflow condition for blood flow. As a consequence, it is totally inefficient and perhaps impossible to perform blood flow simulations using this "brute-force" approach because of prohibitively expensive computations. Recent works have focused on the development of new methods to solve these problems. For example, an attempt to develop new boundary conditions has been presented by Flekkoy et al. [33], in which the simulation domain includes an auxiliary buffer domain for particle generation. However, the complexity of the flux control makes it difficult to perform flow simulations. Recently, a new method for such open systems has been developed by Lei et al. [34], where they generated particles at the inlet according to the local flux and introduced adaptive forces to control the flow rate at the outlet. This method has been successfully applied to single phase flow in straight channels and in bifurcations [34]. However, in multiphase systems, e.g., flows with colloids, polymer chains or RBCs, it is difficult to insert them at the inlet and remove them at the outlet. Thus, existing methods cannot be readily extended to the cases of complex flows such as blood flow.
In this paper we present a general framework for open boundary systems including the inflow and outflow boundaries for particle-based approaches targeting simulations of multiphase flows. We implemented this framework in the context of parallel computations. We show that the particles flowing in a complex computational domain can be treated as a system in contact with a simpler subsystem with a fully developed flow for the inflow combined with an osmotic membrane to control the outflow.

Multiscale RBC model
We simulated the blood flow in arterioles with the help of a multiscale RBC (MS-RBC) model [35] based on the dissipative particle dynamics (DPD) approach [36][37][38]. For completeness, the MS-RBC model is briefly reviewed below, whereas details on the RBC model are available elsewhere [35,39].
In the MS-RBC model, the cell membrane is modeled by a 2D triangulated network with N v vertices connected by springs, where each vertex is represented by a DPD particle. The RBC membrane model takes into account the elastic energy, bending energy, and constraints of fixed surface area and enclosed volume, which is defined as where V s is the elastic energy that mimics the elastic spectrin network, given by where k B T is the energy unit, A α is the area of triangle α formed by three vertices. Also, x i = l i / l m , x 0 = l 0 /l m , where l i is the length of spring i, l 0 and l m are the equilibrium spring length and maximum spring extension, and p is the persistence length. The cell membrane viscoelasticity is imposed by introducing a viscous force on each spring, which has the form, where γ T and γ C are dissipative parameters; v ij is the relative velocity of spring ends, and is the traceless symmetric part of a random matrix representing the Wiener increment.
The bending resistance of the RBC membrane is modeled by where k b is the bending modulus constant, θ αβ is the instantaneous angle between two adjacent triangles having common edge, and θ 0 is the spontaneous angle. In addition, the RBC model includes the area and volume conservation constraints, which mimic the area-incompressibility of the lipid bilayer and the incompressibility of the interior fluid, respectively. The corresponding energy terms are given by where k a and k v are the area and volume constraint coefficients. Here A 0 and V 0 are the equilibrium area and volume of a cell, respectively.
The MS-RBC model is multiscale, as the RBC can be represented on the spectrin level, where each spring in the network corresponds to a single spectrin tetramer with the equilibrium distance between two neighboring actin connections of * 75 nm. On the other hand, for more efficient computations, the RBC network can also be highly coarse-grained with the equilibrium spring lengths of up to 500 * 600 nm. In most simulations, we use N v = 500, a highly coarse-grained RBC model which has been employed to conduct efficient simulations of RBCs in microcirculation [16,35,40,41]. For comparison, we also consider two finer DPD cases with N v = 2560 and N v = 5000. The internal and external fluids are modeled by free DPD particles.

Inflow and outflow boundary conditions
For simple computational domains where PBCs can be applied, our framework is shown to exhibit exactly the same flow characteristics as those obtained by imposing PBCs. Of course, what is more important is that it can be applied to domains where PBCs can not be employed. Unlike previous works, the proposed approach has a great advantage of being applicable to both simple fluid flows and suspensions. In the proposed scheme, the computational domain is divided into three regions as illustrated in Fig 1. Here, the simulation is performed in region B, while regions A and C are auxiliary. In order to generate an inflow in the main computational domain, we use the generating region as a source of new particles. At the same time, when a particle leaves the main simulation domain, it enters into the region C and is removed from the system.
To illustrate the approach in detail, let us consider a microtube flow of RBCs and plasma suspension as shown in Fig 2. In order to have a fully developed inflow in the main computational domain, the generating region is expected to mimic the flow in an infinite tube and be independent of the simulation in the main simulation domain. To reach these requirements, the PBCs in the generating region were implemented by ghost particles along with the particles shifting from one side of the domain to another side when they leave the periodic box. For all particles from zone A2 we create ghost particles and place them in zone A4, and a similar procedure is used for zones A3 and A1. The width of zones A1-A4 is the maximum of the cutoff radii of the force interaction used in simulations. Ghost particles are created after integration in the velocity-Verlet algorithm, but before computer processors exchange forces. The independence of the generating region from the main simulation domain is achieved by turning off forces acting from particles in the main simulation domain on the particles in the generating region. At the same time, the interactions in the opposite direction are preserved because it is desirable to prevent penetration of created particles and their topological structures into the generating region. To connect the aforementioned regions, we design a procedure to duplicate particles from the generating region to the main simulation domain. That is, when a particle in the generating region crosses the copy border, a duplication of the particle is created in the main simulation domain. This duplicated particle is created in the main simulation domain, and hence the inflow is fully developed.
The extension of the proposed inflow boundary conditions for blood flow requires extra care for the cell topology. Specifically, for the ghost interaction implementation, a ghost particle corresponding to a cell vertex needs to keep bonded potential terms such as bonds, angles, and dihedral angles. In addition, the size effect of the periodic domain for computing the corresponding interactions needs to be considered. The procedure of RBC generation at the inlet is also different from the one for a single particle. First, the RBC cannot be duplicated particle by particle, instead, the whole RBC is expected to be duplicated once its center-of-mass (COM) crosses the copy border between the generating domain and main simulation domain. Second, when a new RBC has been generated in the main simulation domain, some vertices of the new duplicated RBC may still be in the generating domain. Thus, we have to pay more attention to the RBC particles and their duplications in the generating domain because they may be in close contact and cause artificial strong repulsive interactions. To avoid this, we enforce the duplicated RBC to move exactly like the original RBC in the generating domain until the entire RBC lies fully inside the main simulation domain. Without extra effort, the proposed procedure allows us to add RBCs and fluid naturally to the main simulation domain. In particular, there are no artificial interactions because of the one way interaction between particles in the generating region and particles in the main simulation domain. It is worth mentioning that the proposed method can achieve a "seamless" connection between the generating region and main simulation domain, so the RBCs flowing in a complex computational domain can be treated as a system in contact with a simpler subsystem with a fully developed flow for the inflow. The periodic pattern observed in regions A2 and A3 is not unexpected since we run the simulation in this region periodically to generate the full developed inflow. An increase of the length of the generating region leads to a better time-averaged flow properties at the inlet and may eliminate or reduce such periodic pattern, but its effect to the simulation results is insignificant since the flow has already been fully developed at the inlet. Also, it requires extra computation time.
In order to impose the outflow boundary conditions for the simple fluid flow, we employ a method similar to Lei et al. [34]. Specifically, those particles leaving the simulation domain and entering into particle deletion region are reflected back to the main simulation domain with a probability (1 − P) depending on the particle number density ρ we want to preserve. They are computed at each iteration using the following algorithm: The generating region is divided into zones: zones A2 and A3 are sources of ghost particles while zones A1 and A4 are used for placing ghost particles. As soon as a particle crosses the copy border, its copy is created. There is one way interaction between particles in the generating region and particles in the main simulation domain. , where h is a weighting factor and it is set at 0.05 in this study.
3. If ρ ρ target , P is updated to (P + dP); otherwise it is P = (P − dP). 4. For a particle crossing the outflow plane, reflect the particle back with probability (1 − P).
We note that identical outflow boundary conditions might be implemented by applying adaptive forces in the vicinity of the outflow plane [34]; however, the formulation proposed here is much simpler to implement. A similar reflecting membrane has been used to generate a fluid flow in previous molecular dynamics simulations [42]. The outflow boundary conditions for RBCs are implemented in a different way. When the whole RBC is inside the region for cell deletion, we destroy the cell topology but leave particles in place and change their properties from the cell-like particles to fluid-like particles, which are removed downstream as described above. We found that this method outperforms an alternative implementation, where the whole cell is deleted because the removal may create density artifacts in the region of cell deletion.

Results/Discussion
To validate the proposed open boundary conditions (OBCs), first a single phase flow (without RBCs) in straight microtubes is simulated and compared with an analytical solution. Numerical simulations are carried out in a 3D geometry representing the microtube used in DPD simulations. In all simulations, the solid walls are modeled by freezing layers of particles with bounce-back reflection to satisfy the no-slip boundary condition [43,44]. Here, for simple fluid and, later, for plasma in blood suspension, the following DPD parameters are employed [34]: a = 4.0, γ = 30.0, r c = 1.5, k B T = 0.0945, n = 2.96. A generalized weight function, w(r) = (1 − r/r c ) s , for dissipative force with s = 0.5 is also used in order to increase the viscosity of the DPD fluid [45,46]. An external body force with magnitude of g = 0.1 is exerted on each fluid particle to generate a Poiseuille flow in the microtube. The microtube diameter is d = 10.0 μm. Fig 3a shows the average velocity profiles obtained from the simulations with the OBCs. The velocity profile is parabolic, which agrees well with the analytical prediction and proves the correctness of the scheme. For a more quantitative analysis, we compute the pressure profile along the flow (z-axis) direction (see Fig 3b). We find that the DPD simulation results are in good agreement with the analytical prediction given by dP/dz = 16v max η/d 2 = ng, where η is the viscosity of the DPD fluid.
Having verified the single phase flow, we simulate the motion of RBC suspension through a straight tube. Fig 4 shows the average velocity profiles for blood flow at two different hematocrit levels (H t ), i.e., H t = 15.0% and 30.0%. In this plot, quasi-parabolic (or flat plug-like) shapes of the typical velocity profiles of blood flow are shown and compared against results obtained with PBCs. These results indicate that the blood flow in microtubes can be simply and accurately implemented by the described scheme.
Next, we apply the proposed OBCs to model the blood flow in microvascular bifurcations. We used flow rate ratios, ϕ d , between branches 1.0, 1.2, 2.5, 4.0 and 6.0 and studied the particle recovery efficiency (the proportion of parent RBC flux entering each daughter branch) with respect to the flow rate ratios between two daughter branches. We perform simulations with three hematocrit levels (H t = 15.0%, 30.0% and 45.0%) and find that the higher flow rate ratio yields the higher particle recovery efficiency. A 100% recovery efficiency is achieved for the flow rate ratios starting at 6:1. We also demonstrate that the higher the hematocrit is the higher the probability for a RBC is to follow a lower flow rate branch (see Fig 5a). This trend is similar to our previous observation from simulations of blood flow with the PBCs [32]. We also simulate the blood flow and study the particle recovery efficiency at different levels of coarse graining with N v = 500, 2560 and 5000, using the MS-RBC model, see Fig 5b. We find that the particle recovery efficiency value increases somewhat with finer DPD resolution. For example, for ϕ d = 1.5, the value of particle recovery efficiency shifts from 66.7% with N v = 500 to 69.9%  with N v = 2560 and further to 70.9% with N v = 5000; for ϕ d = 2.0, it rises from 79.1% with N v = 500 up to 81.0% with N v = 2560 and then to 83.3% with N v = 5000. At ϕ d = 6:1, we find that 100% recovery efficiency can be achieved for all three of these cases.
The blood flow and stress characteristics in human arteriolar bifurcations are affected by the branch location and bifurcation angle variation [47,48]. The proposed approach offers an effective method in investigating these effects on the behavior of RBCs flowing through a microvascular bifurcation. We then control the flow rate ratio by changing the bifurcation angle, θ, between the two daughter branches (see Fig 6a) at a fixed flow rate in the parent branch. We find that the particle recovery efficiency is clearly different in the small angle and in the large angle bifurcation, see Fig 6b. More RBCs move into the side branch at a smaller angle bifurcation, while more RBCs move into the main branch at a larger θ. A critical angle value can be estimated at θ % 78 o when nearly half of the RBCs move into each branch. Therefore, our 3D DPD simulations demonstrated that the bifurcation angle influences the RBC flux  to the daughter branches so that its effect on the RBC flux distributions in microvascular bifurcation cannot be neglected.
Finally, we demonstrate that the proposed approach can be employed to target blood flow simulations for multiple inlets and outlets. Here, for illustration purposes, we present a simulation of blood flow in a complex arterial network, which was constructed using the angiogenesis model [49]. As shown in Fig 7 and the S1 Video, the network has three inlets and multiple outlets. An interesting observation is that the number of fluid particles and the number of RBCs are almost constant during the simulations. Thus, we can simulate arterial blood flow and study the effect of combined different outflow and Dirichlet boundaries on the flow pattern.
In summary, in this paper, we have developed a general parallel framework for open boundary conditions including the inflow and outflow boundaries for particle-based methods. The approach presented above offers a straightforward way for an open system simulations such as blood flow in arteriolar bifurcations, which provides the possibility to simulate particulate flows for various systems with open boundaries. It was implemented as an extension to the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [50] and extensively tested in simulations for different domains in the High Performance Computing environment. To the best of our knowledge, this is the first simulation of blood flow in an arterial network with the MS-RBC model. Due to the nature of the proposed technique, there are some disadvantages to consider. First, the inflow at the inlet is duplicated from a fully developed flow generated in a pilot simulation with PBCs, thus, the proposed technique can not be used if the recirculation region is present at the inflow. Second, as the implementation of the proposed technique requires significant increase in communication among computational nodes in comparison to that of normal PBC systems. It is efficient for flow of bodies whose size is much smaller than that of the computational subdomains assigned to a process; however, it may be inappropriate when modeling an immersed body with a long chain structure (such as long Fig 7. Snaphot for simulation of the blood flow in part of arterial network with three inlets and multiple outlets. The complex microvascular network was constructed using an angiogenesis model [49]. The walls of the microvascular network were formed by stationary DPD particles. An extra bounce-back reflection was applied to fluid particles to prevent them from entering the solid wall domain. Supporting Information S1 Video. Simulation of blood flow at H t = 15.0% in a 3D arterial network. The arterial network has three inlets and multiple outlets, and each of them has an internal diameter ranged from 28.0 μm to 40.0 μm. The network is modeled in a simulation box of size 192.0 μm × 250.0 μm × 311.0 μm with a total of 3,254,000 plasma particles in the system. (GIF)