OBJECT-IN-FLUID FRAMEWORK IN MODELING OF BLOOD FLOW IN MICROFLUIDIC CHANNELS OBJECT-IN-FLUID FRAMEWORK IN MODELING OF BLOOD FLOW IN MICROFLUIDIC CHANNELS

fluid and immersed elastic objects. The fluid is modeled using lattice - Boltzmann method [4] and the objects using immersed boundary method and spring network model of the membrane. This boundary, i.e. object’s surface, is covered with triangular mesh. Each mesh node is connected with several neighboring nodes and there are elastic forces corresponding to five elastic moduli, which are responsible for elastic behavior of the objects. They preserve the surface area (both locally and globally), volume and shape of the object. This model is implemented as Object-in-fluid framework in open-source We present a fully three-dimensional computational model of red blood cells and their flow in a fluid. This model includes all components necessary to capture important physical and biological aspects of the cell flow. It comprises descriptions of elasticity of the cell membrane, cell-cell interactions, two-way cell-fluid interaction, and adhesion of cell to surfaces. Using this model, we analyze different processes involving flow of cells. We present the results of ongoing research concerning the development of model for cell adhesion, the analysis of microfluidic devices with periodic obstacle arrays, the optimization of microfluidic connectors and biological process of red blood cells formation from reticulocytes.


Introduction
Currently, blood flow modeling is very useful in various applications. It may be highly abstract (e.g. using electromechanical analogies to model arterial trees [1]) or very realistic (e.g. comparing simulations and experiments of passage of cells through a sensor that measures their deformability [2]). In this work, we describe one such model, which includes a homogeneous fluid -blood plasma -and moving objects immersed in it -red blood cells (RBCs). RBCs are very elastic and therefore, in modeling on single-cell level, it is very important to capture their elastic properties. The model further includes adhesion, cell-cell interactions and two-way cell-fluid interaction. We use this model to simulate processes inside microfluidic devices. These lab-on-a-chip devices include various types of microchannels that may also contain arrays of obstacles. The blood flows through the micro-channels and it is possible to sort the cells or capture specific rare cells from the total cell populations. This kind of separation or capture is useful in diagnostics and monitoring of various diseases (e.g. cancer, sicklecell anemia, malaria). More details on modeling blood flow in microfluidic devices can be found in [3].
Our goal is to optimize the design of such micro-channels. We investigate the connectors that are used for inflow into microfluidic devices. In general, the cells in these connectors undergo significant shear stress that may lead to their damage. We are working on deriving a modified blood damage index that would allow future optimization of such connectors. In case of channels that include periodic obstacle arrays (POA), we focus on the capture rate of rare cells, should the surface of obstacles be coated with rare-cell antibodies. We explain the adhesion model that can be used for modeling the capture. And finally, we look at a practical problem of modeling the formation of RBCs that are used for further simulations. The combined progress in these partial topics will bring us closer towards development of improved microfluidic devices.

Model
Our model consists of two main parts: fluid and immersed elastic objects. The fluid is modeled using lattice -Boltzmann method [4] and the objects using immersed boundary method and spring network model of the membrane. This boundary, i.e. object's surface, is covered with triangular mesh. Each mesh node is connected with several neighboring nodes and there are elastic forces corresponding to five elastic moduli, which are responsible for elastic behavior of the objects. They preserve the surface area (both locally and globally), volume and shape of the object. This model is implemented as Object-in-fluid framework in open-source in the same simulation, we consider average values from several simulations with different number of cells.
Inlets and outlets of microfluidic devices, where the flow is redirected, may be a major source of shear stress which acts on the cells. The relative deviation of global area of a cell is computed as We randomly seeded 50 cells inside such a typical inlet of a commercially available analysis chamber (Fig. 1). We run 10 simulations, each with a different seeding of cells. For each cell we computed relative deviation of global area at any given time. The black line in Fig. 2 denotes CDI. We found a certain amount of cells with low damage, plenty of cells with mid-range damage and a few cells with high damage. Such histograms can indicate the quality of a microfluidic channel. Less damage is indicated by shift of the black line to the left. Further work will include different inlet geometries and we will compare simulation results of the standard continuum model (Navier-Stokes equations with BDI) and our Object-in-fluid model with newly developed CDI. scientific package ESPResSo [5]. More details on the model can be found in [6].
The model of red blood cell was calibrated by comparing the simulation results with results from stretching experiments of biological red blood cells [7]. We have also investigated proper mass distribution in the immersed boundary points and seeding of dense suspensions of cells, which are necessary for simulations with high hematocrit [8].

Optimization of device inlets
Cell sensitive applications like dynamic cell culturing rely on a predefined shear stress. Deviations from this optimal state lead to unwanted cell activation and even cell destruction. Microfluidic channels usually lack the information about cell activation or destruction. Complex geometries inside microfluidic devices increase the deviation from a blood cells' native condition. Recent efforts to quantify the blood cell damage concerned conventional blood pumping devices and a blood damage index (BDI). The BDI sums up the local shear stress along streamlines inside a microfluidic channel to estimate possible mechanical stress, leading to activation and/or hemolysis of RBCs. A conventional BDI calculation uses finite element solvers for Navier-Stokes equations to estimate the shear stress distribution of the whole blood volume. The blood damage index calculation typically shear stress x and empirical constants K, a and b (e.g. [9]). Publications [10] and [11] summarize several attempts to relate the blood damage index to hemolysis of cells (destruction of cell membrane). However BDI often fails to reproduce experimental hemolysis [11].
A certain flow geometry can be optimized in such a way that BDI is minimized. Using the Navier-Stokes equations, the blood is considered as a continuum, there are no mutual interactions of cells or cells with fluid. Since in our model we have implemented these types of interactions, we introduce a cell damage index (CDI), in which the single cell behavior is included. There are several possibilities, how to specify the CDI. The RBC membrane can withstand a finite strain, beyond which it ruptures. The strain is related to both shape and area change of the given object. In our model, we know the current deviation of global area of each cell at any given time. We use this for an approximation of total strain. Moreover, the structural damage of the membrane is influenced not only by the strain magnitude, but also by the duration of strain. In [12], the authors claim that the exposure time of areal strain plays a significant role in the cell rupture, therefore we consider cumulative value of global area deviation during the time, when the cell passes a fluidic channel. Since the values for any individual cell can vary significantly compared to other cells

Estimating capture rates using simulations of periodic obstacle arrays
There are several applications in biomedicine that rely on isolation of rare cells (e.g. circulating tumor cells -CTCs) from large populations of other cells. This isolation can be achieved using periodic microfluidic obstacle arrays, which use various physical or biological properties of cells: some of them work on the principle of transverse displacement, others use antibody coated surfaces that capture the rare cells if these come into contact with the coated surface for sufficiently long time.
One of the assumptions sometimes used in modeling of periodic obstacle arrays (POAs) is that it is sufficient to know the rare cell's radius and its offset from obstacles that identifies the future trajectory, in order to determine the cell's collision mode and collision rate, e.g. using Lagrangian tracers in [13].
We studied the interplay of column radius and hematocrit and their combined influence on the number and duration of the rare cell -column contacts. We found that the interactions of rare cell with other cells may have impact on the contacts with obstacles, especially in dense suspensions, where the interactions may result in displacement of rare cell from the average expected trajectories [14].
We simulated a chamber periodic in x-and y-directions with cylindrical obstacles, Fig. 3. Due to the large computational complexity of the problem, we simulated one passage of the rare cell through the chamber over the range of parameters (varied column radius and number of RBCs, different initial CTC positions at the left side of the chamber, different random seedings of RBCs). We traced the distance of one CTC from the nearest obstacle and the duration of contact(s), if any. We then analyzed the CTC trajectories an determined that regardless their initial position, in RBC-sparse simulations they tend to fall into a small number of typical trajectories, Fig. 5. While we observe the same general pattern in simulations with large number of RBCs, Fig. 6, the trajectories are perturbed by the CTC-RBC interactions, which may influence the capture rate, as evident in Fig. 4. The capture rate was calculated using where is the cell-obstacle contact area, is shear stress with respect to the closest obstacle. Constants a and b are lumped parameters determined using [15]. They include several biological quantities such as ligand and receptor densities, association constants, characteristic length of the ligand-receptor bond and others.
Dissociation and association dynamics are represented by probability of bond formation and bond breaking. At the moments when the receptor and ligand are close enough, the bond can be created with a specific rate kon and the probability that this happens during time interval t d is given by Pon . On the other hand, any bond can rupture with the rate koff with dissociation probability Poff over time.
, exp e xp P k t P k t 1 1 In [18], we have introduced the simplification of previously known adhesion model from [17]. The model from [17] includes the dependence of koff on the magnitude of the force exerted on the bond. It is a natural assumption, however, the dependence of koff on the bond prolongation is quite complicated, includes parameters such as bond detachment force and reactive compliance that are quite difficult to measure, and are not available for many receptor-ligand couples. Our simplification is based on the fact that koff is a fixed constant.   [19]. This leads us to conclude that it is important to include actual cells and their cell-cell interactions when computing the capture rates in microfluidic channels with periodic obstacle arrays. In the future, this approach will be verified by comparisons with model that includes adhesion as described in the following section.

Simplified adhesion model
Modeling of the adhesion and rolling of biological cells on functionalized surfaces gives valuable insights into rare cell isolation. We study cells moving in shear flow above a wall to which they can adhere via specific receptor-ligand bonds based on receptors from selectin as well as integrin family. The adhesion mechanism is modeled by adhesive bonds including stochastic rules for their formation and rupture. Once the cell comes close enough to such a wall, receptors and ligands start forming bonds, slowing the motion of the cell down [16]. We explore a simplified model with dissociation rate independent of the length of the bonds.
Before a receptor and a ligand can form a bond, their physical locations need to be brought in close vicinity and a bond can be created with certain probability. Reversely, any bond can rupture, again with certain probability. Finally, whenever receptor and ligand positions move apart more than a certain threshold, the bond is broken with probability one, see [17].
The computational model of a cell described in previous sections can be readily extended with any other phenomenon based on forces acting on the boundary mesh points of the elastic object. The adhesion mechanism is based on bonds between the receptor site on a cell and the ligand site on a wall. This bond can be modeled as a harmonic spring. Once the stiffness of such spring is known, the bond exerts repulsive or attractive force F l aff $ l = on the corresponding mesh points.
The model we use, allows us to simulate and explore the properties of the process, using both the inflation and deflation scenarios. We based the first one on an elastic model of a sphere with the volume of a real RBC and with a geometrically corresponding, smaller surface. Using the global area conservation forces, we relaxed this triangulation into a shape with the surface of RBC. With the second scenario, we deflated a sphere with the RBC surface and larger volume to smaller RBC volume, using the global volume elastic modulus. In both cases, we achieved shapes that visually correspond to real RBCs (Figs. 9 and 10).
The simulation accuracy can be firstly verified by looking at basic characteristics such as surface and volume. Shapes obtained by inflation and deflation scenarios intrinsically have the surface and volume consistent with biological RBC. Another accuracy check of the obtained cell is verification of its curvature. Since the cell surface is discretized by a triangular mesh, there is a natural opportunity to characterize its curvature using the distribution of inner angles formed by the incidental triangles of the triangulation. We looked more closely at statistics of this distributions.

Verification of the simplified model
To show that this model is capable of capturing the mesoscopic properties, such as velocity of rolling cells, we perform two computational experiments. In the first experiment, we test whether our model can capture statistical properties of velocity profile. We track the velocity of individual cells during the adhesion movement. The movement of adherent cell is much slower than that of non-adherent cell due to bonds between ligands and receptors. Since the bonds on a cell in the downstream brake and bonds in the upstream are continuously created, the cell rolls on the surface. The velocity of the cell itself randomly oscillates according to how often bonds brake and how often they are created again. Because of the stochastic nature of bond formation, every simulation run gives different velocity profile. However, one can look at the mean value and the standard deviation and compare the simulated movement with the measured data for real cells.
The second experiment reveals that our model can capture the dependence of cell velocity on the different fluid shear stresses. Namely, the question is whether one set of parameters gives the same values of averaged velocity of the cell for different fluid shear stresses. Detailed description of the computational experiments is given in [18]. We present the result in Figs. 7 and 8. In [20], the authors performed laboratory experiments and obtained the rolling velocity of the cells. Our model is capable of reconstructing these velocities. In Fig. 7, we can see the measured velocities depicted by green symbol cross and the velocities obtained from our model depicted with red symbols plus. Note that we were able to set the model parameters so that both velocity profiles have the same mean value.
In [19], the authors investigated the mechanics of leukocyte adhesion to endothelial cells using in vitro side-view flow assay. The authors showed how the averaged velocity of the cell depends on shear stress. We demonstrate that we are able to capture this dependence. The experiment was performed for shear stresses ranging from 0 to 2 Pa. For each value of stress, we tracked the movement of the cell and we computed the averaged velocity over the time 1ms. The results are presented in Fig. 8.

Modelling blood cell formation
Overall, RBCs represent up to 95% of the solid blood component volume. Therefore, correct modeling of RBC behavior is the key to the correct blood flow simulations. A natural validation of the RBC model is based on its specific properties associated with RBC formation process. It is known that the shape of RBC with a given volume maximizes its surface. We also know that RBC is formed from the spherical shape of the proerytroblast or reticulocyte, with the aim to maximize the surface-to-volume ratio. The second was a comparison of the curvature histogram of the RBC created from sphere to its known biological shape. However, here we are confronted with the problem of the existing triangulations of this surface. They exhibit considerable regularity, which is reflected in the histograms and they are difficult to compare to our significantly stochastic triangulations (Fig. 12,   Figures 9 and 10 show the histograms of these angles for the inflation and deflation scenarios, and Table 1 shows the maximum, minimum, average and standard deviation of those angles. Note that this concerns the interior angles: the values above 180° correspond to the concave parts of the cell surface and the values below 180° correspond to the convex parts of the cell surface. The obtained data indicate that the two processes, where different elastic moduli were dominant, result in remarkably similar distributions. This result strengthens the confidence in consistency of the model. Table 1 With regard to the stochastic nature of the simulations, we also confirm this result by the statistical tests of the selection from identical empirical distribution. The necessary precondition of the results comparison is using comparable regular meshes with almost the same number of nodes and adjacent triangles. Since in our model we can easily use various regular meshes, we do not have to deal with the complicated problem of normalizing the obtained curvature spectrum for significantly different triangulations.

Statistical values of angle distribution
We tested the consistency of the RBC curvature characterization by looking at the angles distribution statistics in two more ways. The first one involved a simulation experiment and examined the RBC deformation during its flow through a narrow opening (Fig. 11, Table 1). As expected, the values of the convex and concave angles in the central deformed part of the RBC have larger range. The histogram shape and the increase of the standard deviation correspond to that. Further computational experiments from section 3.4 confirmed the consistency of the model. Using two different scenarios for formation of a red blood cell, we have obtain objects with the same statistical properties. Thus inflation and deflation scenario give comparable results. This justifies the use of inflation scenario to seed cells into dense suspensions: we have successfully created dense suspension by starting with smaller objects and gradually increasing their surface to real values.