Chapter 1

Introduction

This thesis concerns blood flow in the microvasculature. Section 1.1 outlines the basic functioning of the circulatory system. Section 1.2 reviews the main components of whole blood, including potential pathological conditions of red blood cells. Section 1.3 is dedicated to the endothelium, whose varied interactions with blood flow serve as the subject of this thesis. Section 1.4 describes blood flow conditions in the microvasculature, including a brief overview of two important phenomena in the study of blood flow: the Fahraeus and the Fahraeus-Lindqvist effects. Section 1.5 briefly introduces the main contributions of this thesis.

Circulatory System

The main function of the cardiovascular system is to deliver necessary substances such as nutrients and oxygen to cells and to transport waste products such as carbon dioxide away. The human cardiovascular system is a closed system whereby the blood is pumped from the heart through a ramified network of blood vessels ranging in size from aorta (2 -3cm in diameter) to capillary (5 -8µm in diameter) which traverse the body and bring the blood back to the heart again. The cycle begins as oxygenated blood is pumped by the left ventricle of the heart into the aorta, the largest vessel in the arterial system. From the aorta, the blood branches off into arteries and then flows into the microcirculation through the arterioles and capillaries. Microcirculation has three major components: pre-capillary, capillary and post-capillary. In the pre-capillary sector, arterioles are surrounded by smooth muscle cells. Their function is to regulate blood flow before it enters the capillaries and venules by the contraction and relaxation of the smooth muscle found on their walls. Then, the blood flows into the capillaries, the smallest vessels of the blood vessel network. The nutrient and oxygen transfer takes place in the fine network of capillaries, referred to as a capillary bed, where the blood flow rate is reduced, facilitating transfer. Blood flows out of the capillaries into the venules then into the veins which bring the deoxygenated blood back to the heart. The blood is brought back to the heart through the right atrium and right ventricle whereby it is pumped into pulmonary circulation, releasing carbon dioxide and waste and reoxygenating the blood. This newly oxygenated blood is brought back to the heart through the left atrium and is pumped through the left ventricle to the aorta and so the cycle continues. In its journey throughout the circulatory system, blood experiences a large range of spatial scales. Between the aorta and capillaries, there is about a 3,000-fold difference in diameter. The size of the blood vessel has important ramifications for the blood flow properties: we will see later that blood flow properties depend on the size of the conduit through which it flows.

Blood Components

Blood is a suspension composed of plasma (∼ 55% by volume) and suspended blood elements: red blood cells (RBCs) (∼ 44% ), white blood cells (WBCs), and platelets (together ∼ 1%). Fig 1 .2 shows a scanning electron microscope image of the larger blood components. Blood accounts for 7% of human body weight with an average density slightly higher than that of pure water, 1060kg/m 3 . The average adult has a blood volume of roughly 5 litres. Plasma is essentially an aqueous solution containing ∼ 92% water, ∼ 8% blood plasma proteins, and trace amounts of other materials. Plasma circulates dissolved nutrients, such as glucose, amino acids, and fatty acids (dissolved in the blood or bound to plasma proteins), and removes waste products, such as carbon dioxide, urea, and lactic acid.

WBCs (also called leukocytes) are part of the immune system, with one microliter of human blood containing ∼ 4, 000 -11, 000. In general, WBCs have a spherical shape with a diameter of ∼ 12µm. The role of WBCs as part of the immune system is to destroy and remove old or aberrant cells and cellular debris, as well as attack infectious agents (pathogens) and foreign substances.

Platelets (also called thrombocytes) number between 200, 000 -500, 000 in one microliter of blood and take part in blood clotting to stop bleeding at the site of wounded endothelium. Circulating unactivated, they take the form of a bioconvex discoid 2 -3µm in diameter. At the site of injury to the blood vessel wall, they adhere, activate whereby they change shape, turn on receptors, and secrete chemical messengers, and aggregate, connecting to each other through receptor bridges, forming a "platelet plug" preventing blood loss from the site of injury and providing a barrier to contaminants from outside.

RBCs (also called erythrocytes) are the most abundant cells in blood and represent the largest volume fraction among blood constituents. One microliter of blood contains about 5 million RBCs. RBCs are partly made up of an internal solution containing hemoglobin (abbreviated Hb or Hgb), an iron-containing protein which facilitates oxygen transport by binding to dissolved O 2 in the blood to form oxyhemoglobin, HbO 2 . The oxygen from the lungs diffuses into the blood and subsequently into the RBCs following a down gradient (the alveoli of the lungs are richer in oxygen than the blood entering into the capillaries of the lungs; therefore, a diffusion of O 2 toward the poorer region will occur). Hemoglobin has an oxygen binding capacity between 1.37 -1.39 ml O 2 per gram hemoglobin, (Dominguez de [START_REF] Dominguez De Villota | Equality of the in vivo and in vitro oxygen-binding capacity of hemoglobin in patients with severe respiratory disease[END_REF]) which increases the total blood oxygen capacity sixty-fold, compared to if oxygen were solely carried by its solubility of 0.03 ml O 2 per liter blood per mm Hg partial pressure of oxygen (assuming 1.3 g/L Hg in blood and about 100 mm Hg in arteries) [START_REF] Costanzo | Physiology[END_REF]. Still, if all human hemoglobin were free in the plasma rather than being contained in RBCs, the plasma fluid would become extremely viscous and thus the cardiovascular system less efficient. As it is, oxygen is carried by the RBCs in the blood until released in the capillaries, making mass transfer involving RBCs crucial to the proper functioning of the cardiovascular system. RBCs also carry the majority of carbon dioxide in circulation, either as the bicarbonate ion (HCO -

3 ) as part of the ph-regulating bicarbonate buffer system or reversibly bound to hemoglobin forming carbaminohemoglobin, while the remaining is dissolved in the plasma.

RBC Pathological Conditions

The typical shape of a healthy RBC is a biconcave disk with diameter of ∼ 8µm and thickness of ∼ 2.5µm. The biconcave shape is characterized by a small volume to surface area ratio (reduced volume of a normal, healthy RBC is said to be 0.65), allowing a large deformation while maintaining constant membrane surface area. An essential feature of RBCs is their ability to deform. As capillaries have a comparable size to an individual RBC, RBCs must deform in order to pass through them. The high deformability can be attributed to the biconcave shape exhibited by healthy RBCs and to the flexible membrane. Their major features include a permeable membrane made up of phosphobilipid bilayer and transmembrane proteins connected to intracellular network of spectrin proteins forming a 2-D skeletal shell, known as the cytoskeleton, and an internal hemoglobin solution. Fig 1 .3 shows a schematic of the major features of a healthy RBC. Defects in proteins linking the RBC membrane to the cytoskeleton can result in changes to RBC morphology. In the case of hereditary spherocytosis, genetic defects in the proteins involved in vertical interactions, those connecting the membrane to the cytoskeleton (cf Fig 1 .4), lead to an uncoupling of the membrane from the cytoskeleton. As the RBC travels through the blood system and through the spleen, where it must deform significantly in order to pass through narrow passages, tiny parts of the membrane are broken off in a process known as vesiculation. Vesiculation results in a decreased surface area to volume ratio, driving the morphological change from a biconcave disc into a sphere. Spherocytic red blood cells show decreased deformability compared to their healthy biconcave disc counterparts and they are easily destroyed, lasting perhaps only 30 days in the body, compared to healthy RBCs which survive on average 120 days. Often the patient's bone marrow cannot produce enough RBCs to compensate for the high destruction rate and the person may become anemic. Figure 1.4 shows a schematic of the differences at the scale of the RBC structural proteins in healthy and pathological RBCs which leads to the pathological condition of spherocytotic RBCs. Other pathologies affecting RBC morphology and mechanics exist. Similar to hereditary spherocytosis, hereditary elliptocytosis involves disruptions of horizontal cytoskeleton interactions (i.e. binding within cytoskeleton) which result in decreased membrane mechanical stability leading to RBC vesiculation resulting in reduced surface area to volume ratio and reduced deformability leading to greater risk of anemia in affected patients. Hereditary elliptocytosis is characterized by RBC morphologies with an elliptical or oval morphology. Sickle cell disease arises from a mutation resulting in an abnormal form of normal hemoglobin (HbA) called sickle cell hemoglobin (HbS). RBCs cycle between sickled and unsickled state as reoxygenation causes HbS to depolymerize and the cell to thus unsickle. After many cycles, RBCs lose ability to recover the discocyte shape upon reoxygenation and become permanently sickled. Sickled cells exhibit decreased deformability and increased intracellular viscosity, affecting passage through small blood vessels and leading to decreased blood flow rates [START_REF] Diez-Silva | Shape and Biomechanical Characteristics of Human Red Blood Cells in Health and Disease[END_REF]). As we will see later, the morphology of individual RBCs has profound effects on blood flow properties.

In our investigation of blood flow, we assume that blood is a suspension consisting of deformable particles. The plasma is treated as a simple Newtonian fluid and the RBCs are modeled as a flexible membrane enclosing an internal solution. The membrane is responsible for conferring bending resistance and maintaining cell surface area, while the spectrin cytoskeleton forms the inner shell of the RBC and confers in-plane shear elasticity, which for our 2D simulations is disregarded [START_REF] Tahiri | On the problem of slipper shapes of red blood cells in the microvasculature[END_REF]). As we have seen previously, defects in the linkage between the membrane and cytoskeleton network result in abnormal RBC morphologies with increased volume to surface area ratios. The lipid bilayer being ∼ 4nm in thickness, compared to the size of a RBC, the thickness of the membrane can be ignored. The membrane resists bending, with a bending modulus of 2 × 10 -9 Nm [START_REF] Scheffer | Atomic Force Pulling : Probing the Local Elasticity of the Cell Membrane Atomic force pulling : probing the local elasticity of the cell membrane[END_REF]. The dynamics of a single RBC, modeled as an inextensible, flexible membrane enclosing a fluid in Poiseuille flow have been previously studied [START_REF] Tahiri | On the problem of slipper shapes of red blood cells in the microvasculature[END_REF][START_REF] Kaoui | Two-dimensional vesicle dynamics under shear flow: Effect of confinement[END_REF][START_REF] Aouane | Modeling and simulation of the motion of deformable interfaces in a confined geometry: application to the study of the flow of red blood cells in microcirculation[END_REF]. In studying the dynamics, phase diagrams are often employed in order to delineate the different morphologies and dynamics observed under different conditions. A phase diagram is presented in Fig 1 .5 taken from [START_REF] Kaoui | Two-dimensional vesicle dynamics under shear flow: Effect of confinement[END_REF], showing different morphologies and dynamics as a function of the strength of confinement (x-axis) and strength of flow (y-axis).

Endothelium

The monolayer of cells lining blood vessels known as the endothelium also plays an important role in the functioning of the cardiovascular system. The endothelium, which is composed of a single layer of endothelial cells, lines the entire circulatory system, from the heart (where it is referred to as the endocardium) to the smallest capillaries. The endothelium serves many functions in the maintenance of a healthy cardiovascular system. It serves as a semi-selective permeable barrier between the vessel lumen and surrounding tissues, permitting flux into and out of the bloodstream of oxygen, nutrients, and glucose to nourish tissue and WBCs to fight infection in underlying tissues. The endothelium also plays an important role in vascular flow regulation such as vasoconstriction and vasodilation events. Vasodilation refers to the widening of blood vessels resulting from relaxation of smooth muscle cells which results in increased blood flow and decreased blood pressure. This process is observed during periods of heavy cardiac activity such as exercise whereby the elevated blood flow increases the rate of transport of O 2 , nutrients, and glucose to cells as well as increases the rate of removal of cellular waste in order to accommodate the heightened cellular demands during strenuous activity. The regulatory mechanisms Hemodynamic shear stress has also been linked to endothelial dysfunction whereby regions of low and oscillatory shear stress tend to be the sites of atherosclerotic plaques [START_REF] Frangos | Localization of atherosclerosis: role of hemodynamics[END_REF]. Hemodynamic forces are also said to play an important role in the formation of new blood vessels, a process known as angiogenesis [START_REF] Garcia | Vascular development and hemodynamic force in the mouse yolk sac[END_REF]. The endothelium is covered by a dense, fuzzy layer of glycoproteins and glycolipids known as the glycocalyx, or endothelial surface layer, which is bound to the endothelium and extends into the lumen. It is said to increase in thickness with increasing diameter of the blood vessel, with figures ranging from tens of nanometers to over a micron thick. It is said to be important in the increased hydrodynamic resistance observed in vivo compared to experiments in glass tubes of the same diameter [START_REF] Pries | The endothelial surface layer[END_REF]. It is also thought to play an important role in mechanotransduction events involving Chapter 1. Introduction the endothelium as it is directly exposed to the hydrodynamic forces generated by the flowing blood.

Blood Flow Characteristics

When we consider blood flow in the microcirculation, certain dimensionless numbers can inform us of the dominant dynamics at play. In arterioles, the mean blood flow velocities range between 7 -12mm/s with diameter sizes 15 -60µm, leading to channel Reynolds numbers (ratio of inertial to viscous effects) between 0.03-0.2 and Womersley numbers (measure of importance of pulsatile effects) between 0.2-0.8 [START_REF] Popel | Microcirculation and Hemorheology[END_REF]. This shows that inertial and pulsatile effects cannot be neglected under all circumstances when studying arteriole blood flow.

The mean blood velocity in capillaries is quite small, about 0.2mm/s. The Reynolds number and Womersley number are about 1.4e -3 and 0.07, respectively [START_REF] Popel | Microcirculation and Hemorheology[END_REF]. Thus, one characteristic of the blood flow in capillaries is that both of the Reynolds number and Womersley number are small. Throughout this work both inertial and the pulsatile effects will be neglected.

In this section, we also introduce two important phenomena regarding blood flow in vessels smaller than 300µm.

The Fahraeus effect

In Poiseuille flow, RBCs and blood do not travel at the same average velocity; RBCs exhibit migration towards the center of the channel, leaving behind a cell-free layer (CFL). Therefore, the average RBC velocity is greater than the average velocity of the suspension. We first define tube hematocrit as the volume fraction of RBCs to the total volume of the vessel measured before entrance into the channel and discharge hematocrit as the ratio of the RBC flux to the wholeblood flux exiting the channel. The Fahraeus effect is the observation that the tube hematocrit is consistently smaller than the discharge hematocrit. This effect is directly correlated with the cell migration to the centerline, forming the CFL and RBC-dense core resulting in the mean RBC velocity being higher than the mean blood velocity. A mass balance analysis reveals the Fahraeus Effect is a consequence of these features (RBC migration to center of channel and development of CFL next to channel walls). The inverse effect is observed for white cells and platelets: the concentration increases with decreasing diameter of the tube and is known as the inverse Fahraeus effect [START_REF] Goldsmith | Effect of red blood cells and their aggregates on platelets and white cells in flowing blood[END_REF].

Fahraeus-Lindqvist effect

The second effect is also a consequence of RBC migration to the center of the channel and the development of the CFL. The Fahraeus-Lindqvist effect refers to the decrease of the apparent 1.4. Blood Flow Characteristics viscosity with decreasing vessel diameter until the vessel diameter is of comparable size to the RBC. As illustrated in Fig 1 .7, the apparent viscosity falls until reaching a minimum at diameters 5 -7µm, then as the tube diameter decreases further, the apparent viscosity spikes as the size of the tube becomes considerably smaller than the size of the RBC. The development of the CFL facilitates the passage of RBCs in smaller tubes by reducing flow resistance, with the net effect that the apparent viscosity is less than that for whole blood. It has been shown recently that the organization of the cells affects the relative apparent viscosity; a relation between the confinement, the hematocrit and the viscosity in linear shear flow has been proposed in [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]). We will explore the effect of these parameters on the apparent viscosity in pressure-driven flow in Chapter 4.

Figure 1.7: Consolidated experimental data of flowing RBC suspensions (hematocrit of 0.45) in glass tubes illustrating the Fahraeus-Lindqvist effect whereby the relative apparent viscosity of a suspension of RBCs decreases with decreasing tube diameter until the diameter of the tube is comparable to the size of an individual RBC [START_REF] Pries | Blood viscosity in tube flow: dependence on diameter and hematocrit[END_REF]).

These effects show that taking into consideration the discrete nature of blood flow in the microcirculation is necessary in order to correctly and completely describe its dynamics.

Original Contributions of this Thesis

This thesis is dedicated to the investigation of the interactions between blood flow and the endothelium. These interactions occur through mechanical as well as chemical signaling. Chapter 3 is dedicated to a single RBC under normal or potential pathological conditions and its effect on the wall shear stress. Also in Chapter 3, we take into account the waviness of the endothelium. Chapter 4 deals with suspensions of RBCs, ranging from dilute to dense, in straight and wavy channels. Finally, Chapter 5 addresses the interaction between solute agents released by RBCs and their reaction at the endothelium, for which we develop a numerical method capable of solving coupled membrane-fluid-solute problems. Below we summarize the original contributions of each chapter.

RBC morphology and endothelium waviness

In Chapter 3, we investigate how changes in RBC morphologies affect shear stress at the endothelium. Pathological conditions of RBCs are known to present challenges in the proper functioning of the circulatory system. Our results reveal morphological alterations greatly affect the values of wall shear stress, the RBC creating a local "footprint" disturbance of the shear stress along the wall. Within the microvasculature, we find the inverse minimum distance of the RBC-wall gap normalized to the channel width is predictive of the maximum wall shear stress. When investigating the extremal values of the wall shear stress, we find that the maximum wall shear stress shows more variation over the parameter space than the minimum wall shear stress. In addition, we consider the waviness of the endothelium. We find that the time-averaged behavior of extremal wall shear stress is not predicted by the contiuum description of blood and that subcellular gradients in wall shear stress are significant. This work has given rise to a publication in Biomechanics and Modeling in Mechanobiology (BMMB) [START_REF] Hogan | Shear Stress in the Microvasculature: Influence of Red Blood Cell Morphology and Endothelial Wall Undulation[END_REF]. Also in this chapter, we investigate how deformability mediates an RBC's trajectory in an undulated channel, which could be potentially useful in the design of microfluidic devices. We also find that with large enough amplitude, recirculation zones can appear between the valleys of the undulated walls and these can be used to trap cells.

Suspensions

In Chapter 4, we study suspensions of RBCs in straight and undulating channels. We make an unexpected finding when investigating the apparent viscosity as a function of hematocrit. We observe a decrease in apparent viscosity with hematocrit in small channels, an observation previously made in linear shear flow [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]. Upon investigation, we are able to link the rheological behavoir to the spatial organization of the RBCs in the channel. The RBCs organize themselves in such a way which reduces the increase of viscosity in the face of increasing 1.5. Original Contributions of this Thesis 13 hematocrit. These results have important consequences for the efficiency of blood flow in the capillaries. Our results indicate a decrease in suspension flux in wavy channels. We observe an increase in discharge hematocrit in suspensions flowing in wavy channels versus straight. We developed a 2-phase 2D analytical model to explore the effect of the RBC flux as a function of hematocrit. Our numerical results indicate that waviness decreases the RBC flux. We see large variation in the cell-free layer for the smallest channel. For small hematocrit, we find that the cell-free layer is larger in the wavy case compared to the straight, while for large hematocrit, the cell-free layer is smaller in the wavy case than the straight due to spacial exclusion effects. In investigating the individual trajectories of the center of masses of the RBCs in the suspensions, we observe a robust two-file organization in wavy channels, while in straight channels, we observe "cell flipping" between rows. We also observe significant variation and periodicity in the elongation index of the cells while they pass through the wavy channel undulations. Lastly, we find significant WSS perturbations at high hematocrits even for large channels. The effect of waviness on perturbation WSS is greatest in small channels with low to moderate hematocrits. In larger channels, our results suggest waviness could even work to decrease the wall shear stress compared to the straight channel.

Solute Release

In Chapter 5, we developed a 2D advection-diffusion solver for moving, deformable particles within the LBM framework. We used our model to explore two modes of release from the particles. In the first mode of release, the concentration is initialized inside the particle and diffuses into the domain through the concentration gradient, and in the second, the concentration is located on the membrane and is released to the domain as a source term. We use our model to perform simulations in the microcirculation, investigating the effect of RBC morphology on the solute dynamics. Related to the first mode of release, we find several interesting conclusions. In strongly convection-dominated flows, RBCs act to accelerate clearance of the solute species from the channel, resulting in smaller residence times. A marked reduction in the mean residence time is observed for strongly convection-dominated flows. We also generally find stronger "hot spots" of the concentration along the wall when the RBC is present. Morphological differences in the RBCs can lead to important differences in residence time and dilution upon contact with the wall. In the second mode of release, we find that the residence time and dilution upon contact with the wall are strongly moderated by the exact morphology of the RBC.

Chapter 2

Methods

In this chapter, we present an introduction to the basic numerical methods which are employed and which will serve as a basis for this thesis work. Sections 2.1 and 2.2 introduce the numerical methods used in this thesis to simulate the coupled dynamics of red blood cells (RBCs) in fluid flow: the lattice Boltzmann and immersed boundary methods. The model used to represent RBCs is presented in Section 2.3. A summary of the basic numerical procedure, coupling the methods and models described in Sections 2.1, 2.2, and 2.3, is briefly outlined in Section 2.4. We present the relevant dimensionless numbers in Section 2.5, which represent the link between our simulations and the physical problem, and in Section 2.6, we describe the procedure for relating simulation and physical units.

Lattice Boltzmann Method for Fluid Dynamics

The lattice Boltzmann method (LBM) [START_REF] Guo | Lattice BGK Model for Incompressible Navier-Stokes Equation[END_REF][START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF][START_REF] Krüger | The Lattice Boltzmann Method -Principles and Practice[END_REF]) has been employed to solve various equations; however, in this introduction, we will focus on the use of LBM for solving hydrodynamic equations. In LBM, we consider the evolution of populations of particles moving on a lattice. LBM is thus a mesoscopic approach whereby we are able to recover macroscopic dynamics. The link between the scales is accomplished through a Chapman-Enskog expansion [START_REF] Guo | Lattice BGK Model for Incompressible Navier-Stokes Equation[END_REF][START_REF] Krüger | Shear stress in lattice Boltzmann simulations[END_REF]). The Lattice Boltzmann equation is given by:

∂ f k ∂t + c k ∇ f k = Ω k + F k (2.1)
where c k is the mesoscopic speed in the lattice such that c k = ∆x/∆t, f k is the distribution function in the direction k, Ω k is the collision operator, and F k is the bulk force term. Bhatnagar, Gross, and Krook noticed that the main effect of the collision term is to bring the distribution functions closer to the equilibrium distributions. In the Bhantnagar-Gross and Krook (BGK) approach, we approximate the collision operator as a relaxation of the population to the equilibrium distribution such that:

Ω k = - 1 τ ( f k -f eq k ) (2.2)
Chapter 2. Methods with the equilibrium distribution being given by a 4th order expansion of the Maxwell distribution such that:

f eq k = w k ρ 1 + 1 c 2 s ( c k • u) + 1 2c 4 s ( c k • u) 2 - 1 c 2 s ( u • u) (2.3)
where w k represents the weight factors in each of the k lattice directions, c s represents the speed of sound in the lattice, ρ is the density of the fluid, and u is the macroscopic fluid velocity.

Completing the necessary terms to solve the equation, [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF] give the form of the bulk force term, F k , as:

F k = 1 - 1 2τ w k c k -u c 2 s + c k • u c 4 s c k • F (2.4)
The discretized form of Eq. 2.1 is written as:

f k ( x, t + ∆t) -f k ( x, t) ∆t + c k f k ( x + ∆x, t + ∆t) -f k ( x, t + ∆t) ∆x = - 1 τ ( f k ( x, t) -f eq k ( x, t)) + F k (2.5)
Using c k = ∆x/∆t, we can rewrite the equation as:

f k ( x, t + ∆t) -f k ( x, t) + f k ( x + c k ∆t, t + ∆t) -f k ( x, t + ∆t) = - ∆t τ ( f k ( x, t) -f eq k ( x, t)) + ∆tF k (2.6)
We define ω = ∆t/τ so that we can write the final discretized form of the LB equation (LBE):

f k ( x + c k ∆t, t + ∆t) = f k ( x, t)(1 -ω) + ω f eq k ( x, t) + ∆tF k (2.7)
The routine for solving Eq. 2.7 at each time step consists of the following two steps: local collision and streaming to the new positions on the lattice such that:

Collision step: fk ( x, t) = f k ( x, t) -ω f k ( x, t) -f eq k ( x, t) -∆tF k Streaming step: f k ( x + c k ∆t, t + ∆t) = fk ( x, t)
where f k and fk denote the pre-and post-collision states of the distribution function, respectively. It is the moments of the distribution function that provide the macroscopic fluid properties, such that the fluid density is given by:

ρ = ∑ k f k (2.8)
and the fluid velocity is given by [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF]:

u = 1 ρ ∑ k c k f k + 1 2 ∆t F (2.9)
and the fluid pressure is obtained as an equation of state such that P = c 2 s ρ. The LBE does not exactly describe an incompressible flow. However, in the limit of small Mach (Ma = u 0 /c s where u 0 is the macroscopic velocity of the flow) and Knudsen number (Kn = l/L, where l is the molecular mean-free path and L is the characteristic length scale of the system), through the Chapman-Enskog expansion, the LBE recovers Navier Stokes equations (NSEs) for incompressible flows, written as:

∇ • u( x) = 0 ρ( x) ∂ u( x) ∂t + u( x) • ∇ u( x) = -∇P( x) + µ( x)∇ 2 u( x) + F( x) (2.10)
where µ is the dynamic viscosity and F represents the membrane force applied on the fluid. We employ the D2Q9 for two-dimensional flows with nine discrete velocities c k (k = 0, 1, ..., 8). In a D2Q9 square lattice, each node has eight neighbors as shown in Populations residing on a node move along their respective discretized velocities when streaming occurs. Other lattices are possible where DmQk refers to a lattice in m dimensions with k discretized velocity directions. The speed of sound in the lattice is set for a chosen lattice, such that for the D2Q9 model, c s = (1/ √ 3)∆x/∆t. Also for the D2Q9 model, the relaxation parameter, τ, is related to the kinematic viscosity, ν by ν = (τ -1/2)c 2 s ∆t. When the LBE is used to solve the NSE, τ = 1 is chosen for the best stability in results. In order to have a displacement of unity of the populations at each time step, we set ∆x = 1 and ∆t = 1. For the D2Q9 model, the discrete velocity set c k is written as: 6, 7, 8 (2.11) and the lattice weights are 6, 7, 8 (2.12)

c k =          (0, 0), k = 0 (cos((k -1) π 2 ), sin((k -1) π 2 )), k = 1, 2, 3, 4 (cos(2k -9) π 2 , sin(2k -9) π 2 ), k = 5,
w k =          4/9, k = 0 1/9, k = 1, 2, 3, 4 1/36, k = 5,
The LBM used as a NSE solver has several advantages over competing solvers, including a relatively straight-forward implementation, relatively easy handling of complex geometries, and its amenability to parallelization which means it can be applied to large simulations.

Boundary Conditions

There have been many papers addressing the question of implementing various boundary conditions in the Lattice Boltzmann framework. The trouble arises since, unlike the NSE where we solve for macroscopic fluid quantities such as the velocity and pressure and give boundary conditions in terms of the macroscopic quantities, the LBE solves the evolution of distribution functions. Thus, the various boundary conditions we know at a macroscopic scale must be translated in terms of the distribution functions. Periodic boundary conditions The periodic boundary condition states that when the populations arrive at the outlet of the computational domain, they reenter the domain at the inlet. Given the length of the domain, L, the periodic boundary condition can be written as

Inlet f 1 (1, y, t + 1) = f 1 (L, y, t) f 5 (1, y, t + 1) = f 5 (L, y, t) f 8 (1, y, t + 1) = f 8 (L, y, t) Outlet f 3 (L, y, t + 1) = f 3 (1, y, t) f 7 (L, y, t + 1) = f 7 (1, y, t) f 6 (L, y, t + 1) = f 6 (1, y, t) (2.

13)

No-slip boundary condition The Dirichlet, or no-slip, boundary condition imposes that the fluid velocity at the boundary matches that of the boundary, so that in the case of a non-moving boundary, the fluid velocity is zero. The implementation in LBM is referred to as the "bounce-back" scheme. In the "bounce-back" scheme, when a population reaches a wall node, the population will be reflected back to the nearest fluid node. Therefore, where the populations cancel each other out lies the effective no-slip wall, which is located half-way between the wall node and the nearest fluid node. The half-way bounce-back can thus be written as,

f k * ( x, t + 1) = fk ( x, t) (2.14)
for post-collisional populations fk streaming into the wall where x is the position of the nearest fluid nodes to the wall and k * represents the direction opposite k. For example, as illustrated in Fig 2 .3, in the case of a flat wall, the incoming post-collisional populations to the solid boundary f 3 , f 6 , f 7 are reflected back into the domain. The Dirichlet boundary condition is applied as

f 1 ( x, t + 1) = f3 ( x, t) f 8 ( x, t + 1) = f6 ( x, t) f 5 ( x, t + 1) = f7 ( x, t) (2.15)
While the half-way bounce-back is second-order accurate when the boundary is aligned parallel to the flow, the case of curved boundaries requires special treatment. For curved boundaries, a refined bounce-back scheme combined with interpolation has been proposed by Bouzidi et al. in order to more accurately handle curved boundaries [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF]).

Immersed Boundary Method

The immersed boundary method (IBM) is used to couple the dynamics of the RBCs and the fluid flow. The IBM was developed by Peskin where he proposed the method as a way to couple the dynamics of flexible membranes in fluid flow [START_REF] Peskin | The immersed boundary method[END_REF]. While the LBE is solved on an Eulerian scheme, the RBC membrane is described by a Lagrangian scheme. The idea behind IBM is to spread the forces of membrane nodes to nearby fluid nodes and to update the velocity of membrane nodes through interpolation from nearby fluid nodes. The forces from the membrane become bulk forces in the LBE for the nearby fluid nodes. The translation between schemes is accomplished by use of a smeared discretized delta function ∆( r) where r = X n (t)x such that X n (t) are the Lagrangian coordinates of the n membrane nodes and x are the Eulerian coordinates of the fluid. We write the smeared delta function in 2D as:

∆( r) = φ(x)φ(y) (2.16)
Although it is possible to use various delta functions, Yang et al. developed several smoothed discrete delta functions to suppress non-physical oscillations of the body forces [START_REF] Yang | A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations[END_REF]).

We adopt their smoothed 3-point delta function given as:

φ(x) =                      17 48 + √ 3π 108 + |x| 4 -x 2 4 + 1-2|x| 16 -12x 2 + 12|x| + 1 - √ 3 12 arcsin( √ 3|x| - √ 3 2 |x| ≤ 1 55 48 - √ 3π 108 -13|x| 12 + x 2 4 + 2|x|-3 48 -12x 2 + 36|x| -23 + √ 3 36 arcsin( √ 3|x| -3 √ 3 2 ) 1 ≤ |x| ≤ 2 0 |x| ≥ 2 (2.17)
They also provide a smoothed 4-point delta function given by:

φ(x) =                      3 8 + π 32 -x 2 4 |x| ≤ 0.5 1 4 + 1-|x| 8 -4x 2 + 8|x| -2 -1 8 arcsin( √ 2|x| - √ 2) 0.5 ≤ |x| ≤ 1.5 17 16 -π 64 -3|x| 4 + x 2 8 + |x|-2 16 -4x 2 + 16|x| -14 + 1 16 arcsin( √ 2|x| -2 √ 2) 1.5 ≤ |x| ≤ 2 0 |x| ≥ 2 (2.18)
The discrete 3-point and 4-point delta functions are plotted in Fig. 2.4. The velocity interpolation and force spreading are then accomplished using one of the discrete delta function schemes. A schematic of the velocity interpolations and force spreading Chapter 2. Methods operations is shown in Fig. 2.5. The velocity at membrane node n can be interpolated from the lattice nodes falling in the shaded region as:

˙ X n (t + 1) = ∑ x u( x, t + 1)∆( X n (t) -x) (2.19)
where ˙ X n is the velocity of the membrane node n and u( x, t) is the fluid velocity. The force spreading at point x is calculated as contributions from the membrane nodes (red) falling in the shaded region by:

F( x, t) = ∑ n F n (t)∆( X n (t) -x) (2.20)
where F( x, t) is the fluid force density acting at the fluid node x and F n (t) is the membrane force at node n. The IBM is a way of achieving the coupling of the fluid and membrane dynamics. As we have shown, the principle involves a translation between Lagrangian and Eulerian frames, which represent the RBC and fluid respectively. The fluid-structure interaction is accomplished in two steps: membrane force spreading and velocity interpolation. The force on each membrane node is spread to nearby fluid nodes through a smeared delta function. This force enters the NSE as a body force term. The resulting fluid velocities of the fluid nodes surrounding the membrane are then interpolated to the membrane nodes and the membrane positions are updated accordingly. We note that the IBM can be used with any NSE solver with the presence of a body force term (ex, finite element, finite difference, etc). Zhang et al. combined LBM and IBM for the simulation of RBCs [START_REF] Zhang | An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows[END_REF]).

RBC Model

A simple description of an RBC depicts an elastic membrane enclosing a viscous solution. RBCs also contain a cytoskeleton which confers upon it a shear elasticity; however, in our 2D simulations, this aspect has no meaning. The RBCs are modelled as a collection of springs, as shown in Fig. 2.6 [START_REF] Tsubota | Particle method for computer simulation of red blood cell motion in blood flow[END_REF][START_REF] Tsubota | Effect of the natural state of an elastic cellular membrane on tank-treading and tumbling motions of a single red blood cell[END_REF]. The RBC membrane is modelled by N membrane points. We define three sources of energy potentials associated to the RBC: stretching/compression, area incompressiblity, and bending. The stretching/compression potential exists between neighboring nodes such that:

E l = 1 2 k l N ∑ n=1 l n -l 0 l 0 2 (2.21)
where k l is the spring constant and l 0 is the equilibrium spring length. Lipid bilayers are characterized by a high surface incompressibility due to the presence of strong bonds between lipid molecules. Additionally, the membrane is impermeable to water exchange and the solution inside the RBC, hemoglobin, is incompressible. These two features imply that both the cell surface and volume are conserved in normal conditions. Therefore, we set an inextensibility condition in our model which is achieved by setting k l large, thereby imposing a constant perimeter. In addition, in order to maintain the area of the RBC constant (no water exchange and incompressibility of internal fluid), we use the following potential:

E s = 1 2 k s s -s 0 s 0 2 (2.22)
where k s is the area constant set to be large ensuring the internal fluid incompressibility. Finally, the bending potential is given by

E b = 1 2 k b N ∑ n=1 tan 2 θ n -θ 0n 2 (2.23)
where k b is the bending constant and θ n is the angle between two adjacent springs with θ 0n set to 0. With these potentials, we can calculate the force acting on membrane points by:

F n = - ∂(E l + E s + E b ) ∂ X n (2.24)
The force arising from each of the three potentials can be written in terms of the membrane coordinates directly. The derivation of the forces used in the spring model (applying Eq. 2.24 to Eqs. 2.21-2.23) can be found in Shen's thesis [START_REF] Shen | Blood flow in microfluidic architectures: collective behavoirs of deformable particles in confined flow[END_REF]). E l , Es, and E b represent elastic energies [N • m] per unit thickness h = 1m, and therefore the units of

E l , E s , E b , k L , k b are Newton [N].
In order to avoid numerical complications, a repulsive force is added to keep the distance between two RBC membranes greater than the interpolation distance. The force is given by the Morse potential given by [START_REF] Liu | Coupling of Navier-Stokes equations with protein molecular dynamics and its application to hemodynamics[END_REF])

Φ(x) = D e e 2β(x 0 -x) -2e β(x 0 -x) (2.25)
where x 0 denotes the repulsion range which is set to 2∆x in our simulations (when using the 3-pt discrete delta function). Beyond the repulsion range, the effect is cut off in order to avoid imposing an attractive force. In the interest of investigating RBC attraction, the range of interaction can be extended beyond the range of repulsion where the interaction changes sign.

Numerical Procedure

We briefly outline the numerical procedure for the coupled LB-IBM scheme. Each step can be thought of as a subroutine. 1. Initialization. Define the domain (fluid/solid nodes), initialize the distribution functions, f k ( x, t = 0), to their respective equilibrium distributions at each fluid node x in the domain, and calculate the initial positions of the membrane nodes X n (t = 0). 2. Calculate membrane forces. As the positions of the membrane nodes X n (t) are known at each time step t, one can calculate the forces on the membrane nodes F n (t).

3. Force spreading. Spread the membrane forces F n (t) to nearby fluid nodes. 4. LB collision. Execute local collisions of the distribution function such that f k ( x, t) → fk ( x, t). 5. LB streaming. Stream the distributions along their respective discretized velocities such that fk ( x, t) → f k ( x + c k , t + 1).

6. LB boundary conditions. Implement boundary conditions. 7. Macro values. Calculate the macroscopic values, including the fluid velocity u( x, t + 1). 8. Velocity interpolation. Interpolate the velocity on the RBC membrane nodes, ˙ X n (t + 1), from nearby fluid nodes and update the membrane positions X n (t + 1). 9. Repeat starting at step 2. Now all the variables in t + 1 become variables in t.

The original code for the coupled LB-IBM simulating RBC dynamics in fluid flow was developed by Zaiyi Shen [START_REF] Shen | Interaction and rheology of vesicle suspensions in confined shear flow[END_REF]).

Dimensionless Numbers

The NSEs given in Eq 2.10 can be written in dimensionless form by introducing dimensionless variables. Dropping for the moment (we will come back to this shortly), the furthest term on the RHS in Eq 2.10, F n , which represents an external force applied on the fluid at point x generated by the force on the membrane at point X n , we can then nondimensionalize all of the variables appearing in the equation by: ˆ

x = x/R 0 ˆ u = u/U t = tT P = PR 0 /µ ext U
where we introduce µ ext as the dynamic viscosity of the plasma and we choose the characteristic length scale as the effective radius of the RBC such that R 0 ≡ (S/π) where S is the area of the RBC, the characteristic time scale as T, and the characteristic velocity as U. We can substitute the new dimensionless variables back into Eq. 2.10 and multiply it by the quantity (R 2 0 /µU) so that the equation becomes

ρR 2 0 µT ∂ ˆ u ∂ t + ρUR 0 µ ˆ u • ∇ ˆ u = -∇P + µ( ˆ x) µ ext ∇2 ˆ u (2.26)
which we rewrite in the form

β ∂ ˆ u ∂ t + Re( ˆ u • ∇ ˆ u) = -∇P + λ ∇2 ˆ u (2.27)
where β is the frequency number, Re is the Reynolds number, and λ is the viscosity ratio. We can define the characteristic time scale of the flow as T = R 0 /U and thus we see that β = Re. We can now go back to the force term, F. Retaking Eq. 2.20 and 2.24, we can rewrite the force term as

F = - ∂(E l + E s + E b ) ∂ X n ∆( X n -x) (2.28) we introduce a new dimensionless variable ∆( X n -x) = ∆(R 0 ( ˆ X n -ˆ x)) = ∆( ˆ X n -ˆ x)R 0
using the following property of the delta function: ∆(ax) = 1/|a| ∆(x).

We will briefly explain the significance of these dimensionless numbers we have found as they are used throughout the work.

• Reynolds number represents the relative importance of inertial versus viscous forces.

Re = ρUR 0 µ ext (2.29)
In the case that the Re << 1, we can neglect inertial terms and the NSEs simplify to the Stokes equations which we will briefly discuss shortly.

• Viscosity ratio is defined as the ratio between the viscosity of the fluid at position x and the external (i.e., plasma) viscosity, µ ext . In the case that x ∈ Ω + where we define Ω + and Ω - as the domains inside and outside the RBC, respectively, the viscosity ratio takes on the value

λ = µ( x) µ ext (2.30)
while if the node is located in the outer fluid, such that x ∈ Ω -, λ = 1.

• Capillary number is defined as the ratio between the bending resistance and the strength of the flow.

Ca ≡ τ c /τ f = µ ext R 3 0 κ γ (2.31)
Therefore, for a RBC where Ca < 1, the bending forces dominate, and the cell can be considered as effectively "stiff," while for the case where Ca > 1, the flow forces dominate, and the cell can be considered effectively "soft." The Ca number can also be presented as the ratio between the time relaxation of curvature of the cell τ c = µ ext R 3 0 /κ to the time scale of the flow τ f = γ-1 . We define the characteristic shear rate γ as the imposed velocity gradient evaluated at y = R 0 /2, which we calculate as γ = R 0 U/(W/2) 2 .

• Confinement relates the size of the RBC to the width of the channel, W in which it is passing

Cn = 2R 0 W (2.32)
• Reduced area is another dimensionless parameter that denotes the extent of roundness of the RBC and is defined as the area of the RBC divided by the area of a circle having the same perimeter, p, written below as:

α = S π(p/2π) 2 (2.33)
The equilibrium shapes obtained for various reduced areas, 0.5 A healthy RBC is typically considered to have α = 0.6 -0.7. When we study the effect of the reduced area, we vary the reduced area by varying the perimeter.

≤ α ≤ 1,
In the microcirculation, typical flow velocity values vary from 10 -4 -10 -3 m/s, the size of the RBC is ∼ 10 -5 m and we take the properties of water for the density of plasma and viscosity of the plasma, 10 3 kg/m 3 and 10 -3 Pa • s, respectively. Therefore, Re ∼ 10 -3 -10 -2 and can be disregarded in our simulations, such that the minimum necessary dimensionless parameters to describe our system are (Ca, Cn, λ, α). The parameter space can be systematically studied for these dimensionless parameters and the resulting rich dynamics explored.

Unit Conversion

Conversions between physical units and lattice units is necessary when carrying out LBM simulations. The speed of sound is fixed for a given lattice; in the case of the D2Q9 lattice, it is fixed to c s,LBM = 1/ √ 3∆x/∆t. As stated previously, in order to have displacements of unity of the populations at each time step, we will take ∆x = 1 and ∆t = 1. The fluid is said to have the same properties as water such that the kinematic viscosity and speed of sound are ν = 10 -6 m 2 /s and c s = 1500m/s, respectively. In order to calibrate the units in LBM, we use the law of similarity to say that the Reynolds numbers associated to the physical problem, Re, and the LBM simulation, Re LBM , are equal. We define the Reynolds numbers in both cases such that

Re = Re LBM = ul/ν = u LBM l LBM /ν LBM (2.34)
where l LBM is the characteristic lengthscale for LBM simulations and ν LBM is the kinematic viscosity in LBM simulations, which can be found as ν LBM = 1 2 c 2 s,LBM (2τ -1). Similarly, we can Chapter 2. Methods write for the Mach number, for the physical and simulated cases that 2.35) where u and u LBM are, respectively, the characteristic velocities associated to the physical and simulated systems. Using these two relations, we find that l = (νc s,LBM /ν LBM c s )l LBM . Plugging in the values, we find that l = 2.31 × 10 -9 l LBM , meaning that one spatial unit in LBM is equal to 2.31nm. At this scale, the characteristic length of a single RBC would need to be represented by 10 4 nodes in LBM, such that a computational domain of around 10 4 × 10 4 lattice nodes would need to be resolved, which exceeds computer memory capability. However, we can find a way to circumvent this problem. We do not impose that the simulated Mach number correspond to the physical Mach number. The physical Mach number corresponding to the flow in the microcirculation is very small (on the order of 10 -7 ); however, just a small simulated Mach number (Ma LBM < 0.3) will satisfy the incompressibility condition in LBM. Our simulated Mach numbers will therefore be much larger than the physical Mach number. As we already discussed in the previous section, the physical Reynolds number in the microcirculation is small (Re ∼ 10 -3 -10 -2 ). We will not calibrate the simulated Reynolds number to reality either, but we will set it Re LBM < 1 in our simulations. In fact, the most important dimensionless number concerning the dynamics of RBCs in the microcirculation is the capillary number. As explained in the section above, the capillary number is the ratio between the resistance of the membrane, or the membrane relaxation time, to the strength of the flow, or the inverse shear rate. We will calibrate the capillary number to reality by selecting a reasonable resolution of the mesh for RBC flow as allowed by accuracy and efficiency, such that l = 2 × 10 -7 l LBM . Accordingly, we find a value of the radius of an RBC as R 0 = 12. We then equate Eq. 2.29, which represents the simulated capillary number, with Eq. 2.32. Once the capillary number of the problem is fixed, it can be shown from Eq. 2.29 that the value of k b in LBM should be satisfied with Re = Cak b 0.45 < 1 in accordance with the small Re number assumption.

Ma = Ma LBM = u/c s = u LBM /c s,LBM (

Benchmarking

Planar Poiseuille Flow

We apply our LBM-BGK model previously outlined to solve steady planar Poiseuille flow. The flow is driven by a pressure gradient, periodic boundary conditions are implemented at the inlet and outlet, and no-slip on the walls, as shown in the schematic in The steady Poiseuille flow has the familiar exact solution for velocity in the horizontal direction:

u(y) = ∆P 2µL y(y -H) (2.36)
where u is the velocity in the horizontal direction, ∆P is the pressure gradient applied across the inlet and outlet, L is the length of the channel, and H is the height of the channel. The criterion of steady state is given as

∑ ij u n+1 ij -u n ij ∑ ij u n ij ≤ 1 × 10 -9 (2.37)
where u n ij is the horizontal velocity component at (x i , y j ) at the nth time step.

Planar Flow with RBC

We can compare the deformation of our RBC with those previously shown in literature [START_REF] Kaoui | Two-dimensional vesicle dynamics under shear flow: Effect of confinement[END_REF][START_REF] Tahiri | On the problem of slipper shapes of red blood cells in the microvasculature[END_REF]). Previous works have described the complex dynamics of RBCs in microvessels under pressure-driven flows. As it would be time consuming, we do not benchmark our code against all dynamics previously shown in literature, but rather, we compare against three steady-state solutions: the unconfined slipper, the parachute, and the confined slipper. Since we are expected to have stationary solutions for these three cases, we define the steady state criterion by the y-position of the center of mass of the RBC, such that:

y cm = ∑ i y RBCi N i N total (2.38)
To look for the steady state solution, we observe the center of mass in the y-direction. Fig. 2.10 shows the stationary morphologies of the three cases. In this chapter, the effect of red blood cells and the undulation of the endothelium on the shear stress in the microvasculature is studied numerically using the LB-IBM framework. We introduce the importance of the waviness of the endothelium and the attention it has received in literature in Section 3.1. In Section 3.2, we describe the necessary Bouzidi boundary condition which will provide us with second-order accurate results when using LBM with wavy walls. Section 3.3 is dedicated to the derivation of the wall shear stress on a curved boundary, and Section 3.4 is devoted to the implementation in LBM. Finally, in Section 3.5, we present and discuss the results, including single cells flowing in flat and wavy channels. The results demonstrate a significant effect of both the undulation of the endothelium and red blood cells on wall shear stress. Our results also reveal that morphological alterations of red blood cells, as occur in certain pathologies, can significantly affect the values of wall shear stress. The resulting fluctuations in wall shear stress greatly exceed the nominal values, emphasizing the importance of the particulate nature of blood as well as a more realistic description of vessel wall geometry in the study of hemodynamic forces. We find that within the channel widths investigated, which correspond to those found in the microvasculature, the inverse minimum distance normalized to the channel width between the red blood cell and the wall is predictive of the maximum wall shear stress observed in straight channels with a flowing red blood cell. We find that the maximum wall shear stress varies several factors more over a range of capillary numbers (dimensionless number relating strength of flow to membrane elasticity) and reduced areas (measure of deflation of the red blood cell) than the minimum wall shear stress. We see that waviness reduces variation in minimum and maximum shear stresses among different capillary and reduced areas. Finally, in Section 3.6, we investigate undulated channels for which there exist solutions in which a cell can be "trapped." Inspired by the idea behind deterministic lateral displacement devices (DLD), we also explore how the trajectory of an RBC is mediated by its deformability.

Introduction

The fluid shear stress due to viscous blood flow is an important quantity in the study of the circulatory system. In the endothelial cells, shear stress regulates important mechanotransduction events that control vessel development during embryogenesis [START_REF] Culver | The effects of hemodynamic force on embryonic development[END_REF][START_REF] García-Cardeña | Hemodynamic Control of Endothelial Cell Fates in Development[END_REF][START_REF] Roman | Mechanotransduction in Embryonic Vascular Development[END_REF], regulates vascular tone [START_REF] Davies | Flow-Mediated Endothelial Mechanotransduction[END_REF][START_REF] Hahn | Mechanotransduction in vascular physiology and atherogenesis[END_REF], and plays a role in the localization of atherosclerotic lesions [START_REF] Barakat | Responsiveness of vascular endothelium to shear stress: Potential role of ion channels and cellular cytoskeleton (review)[END_REF]. Endothelial cells show sensitivity to the details of the shear stress applied upon them. In vitro experiments have shown that temporal changes in applied shear stress patterns have an important influence on endothelial cell phenotypic expression [START_REF] Uzarski | Adaptation of Endothelial Cells to Physiologically-Modeled , Variable Shear Stress[END_REF]. While the transient and space-varying nature of shear stress is important to endothelial cell function, a study of blood from a continuum point of view will not suffice to bring out these important time and spatial fluctuations due to RBCs. Thus, studies are necessary in order to determine the shear stresses present in the vessel taking into account the particulate nature of the flowing blood.

Due to the fact that RBCs constitute 40-45% of the total volume of blood under normal conditions, they play an important role in hemodynamics. Consequently the presence of RBCs would be expected to have a significant impact on the shear stress at the endothelial cell surface. Furthermore, pathologies of RBCs can affect their deformability and therefore their morphologies and dynamics in blood flow. To date, few studies have addressed the interplay among RBC morphological changes due to pathology, RBC dynamics, and the stresses experienced at the vascular wall. Two particular cases of pathological RBC morphology have received some attention, including the "schizont stage" of malaria, known to lead to reduced RBC rigidity and morphological changes whereby the RBC assumes a spherical shape. In this stage, RBCs infected by the Plasmodium falciparum (Pf-RBCs) parasite exhibit enhanced cytoadherence to vascular endothelium and other healthy and infected RBCs, which may lead to capillary occlusions [START_REF] Dondorp | Reduced microcirculatory flow in severe falciparum malaria: Pathophysiology and electron-microscopic pathology[END_REF]). Fedosov et al. (2011) developed an adhesion model of Pf-RBCs as a function of wall shear stress (WSS) and other parameters, studying how WSS affects the propensity of Pf-RBC to adhere to the vessel wall. Similarly, [START_REF] Barabino | Endothelial Cell Interactions with Sickle Cell, Sickle Trait, Mechanically Injured, and Normal Erythrocytes Under Controlled Flow[END_REF] studied in vitro how wall shear rate, calculated from the WSS, mediated adherence to the vessel wall in sickle cells and normal RBCs and found that adherence of sickle cells is enhanced at low shear rates but is comparable to that of normal RBCs at high shear rates. The wall shear rate calculated in the study is calculated considering a homogeneous, Newtonian fluid, and what our study suggests is that the RBCs themselves and the details of their shape and rigidity will induce significant alterations in the WSS, especially in the microvasculature.

Other RBC pathologies include spherocytosis and elliptocytosis which are caused by genetic defects in intracellular proteins responsible for the structural connections either between the lipid bilayer membrane and the cytoskeleton or within the cytoskeleton, respectively, and lead to altered cell shape. These genetic defects result in morphologies with reduced surface-to-volume ratios, or higher reduced volumes, through a shear stress-driven mechanism known as vesiculation [START_REF] Diez-Silva | Shape and Biomechanical Characteristics of Human Red Blood Cells in Health and Disease[END_REF]. Even below the hemolytic threshold and independent of genetic abnormalities, shear stress has been shown to induce alterations in RBC properties, most notably a deterioration in cellular deformability [START_REF] Baskurt | Red Blood Cell Mechanical Stability[END_REF]. Reduced RBC deformability has also been observed in RBCs stored for blood transfusions [START_REF] Xu | Stiffness increase of red blood cells during storage[END_REF]. RBC deformability, in addition to hematocrit and RBC size, has been shown to influence platelet adherence to vessel wall, an early step in thrombus formation. This increase is said to be due to increased platelet transport towards the vessel wall induced by the red blood cells [START_REF] Aarts | Red Blood Cell Deformability Influences Platelets-Vessel Wall Interaction in Flowing Blood[END_REF]. Platelets thus being marginalized to the layer of fluid closest to the wall will be exposed to shear stress gradients generated at the vessel wall [START_REF] Kroll | Platelets and Shear Stress[END_REF]. Our results show that variation among extremal wall shear stresses experienced among different reduced areas and RBC rigidities lead to important fluctuations in wall shear stress values.

In wishing to understand how different RBC morphologies may affect the wall shear stress at the endothelial surface, we acknowledge an additional complexity in that the endothelium is not a flat surface. Atomic force microscopy (AFM) has been used to image the endothelial surface topography in vitro [START_REF] Barbee | Shear Stress-Induced Reorganization of the Surface Topography of Living Endothelial Cells Imaged by Atomic Force Microscopy[END_REF][START_REF] Oberleithner | Human Endothelium: Target for Aldosterone[END_REF]). [START_REF] Barbee | Shear Stress-Induced Reorganization of the Surface Topography of Living Endothelial Cells Imaged by Atomic Force Microscopy[END_REF] reported that unsheared cells had an aspect ratio (length/width) of 1.12 ± 0.31 and a height of 3.39 ± 0.70 µm. After exposure to flow, the cells elongated and shrank in height so that the aspect ratio increased to 2.16 ± 0.53 and the height decreased to 1.77 ± 0.52 µm. [START_REF] Oberleithner | Human Endothelium: Target for Aldosterone[END_REF] Chapter 3. Results: Red Blood Cells under Channel Flow: Effect of Endothelium and Red Blood Cell Morphology Although these experimental results clearly show the non-flat nature of the endothelium, there have been a limited number of theoretical or numerical studies taking into account the waviness of the endothelial cell geometry. [START_REF] Satcher | The Distribution of Fluid Forces on Model Arterial Endothelium Using Computational Fluid Dynamics[END_REF] modeled the wavy surface as a two-dimensional sinusoidal surface and showed through a linear perturbation analysis that the uneven endothelial surface leads to a nonuniform shear stress distribution at the cellular level. Their results showed that the perturbation shear stress due to the wavy surface could reach values as much as a third greater than the stresses experienced in the flat case, with the peak perturbations occurring at the crests of the wavy surface, as would be expected. [START_REF] Barbee | Subcellular distribution of shear stress at the surface of flow-aligned and nonaligned endothelial monolayers[END_REF] used their AFM measurements to generate three-dimensional geometries of the wavy endothelial surface that were subsequently used in numerical simulations of the wall shear stress. Similar to the findings of [START_REF] Satcher | The Distribution of Fluid Forces on Model Arterial Endothelium Using Computational Fluid Dynamics[END_REF], their results showed large subcellular gradients of wall shear stress. Furthermore, they demonstrated that flow-induced elongation of the endothelial cells had the effect of significantly reducing these gradients. [START_REF] Pontrelli | Modelling wall shear stress in small arteries using the Lattice Boltzmann method: Influence of the endothelial wall profile[END_REF] investigated the wall shear stress in wavy-walled channels using 2D simulations. Although the effect of the particulate nature of blood on wall shear stress was only slightly touched upon, the authors reported that the mean effect of erythrocytes matches well with continuum hydrodynamics, which logically follows since the simulations involved a pressure-driven flow, thereby imposing the average shear stress along the walls. Interestingly, they indicated that they did find that the instantaneous wall shear stress distribution in the presence of RBCs in the flow could not be recovered by a consideration of blood as a continuous medium alone.

In the present work, we use a lattice Boltzmann-immersed boundary method (LB-IBM) code where RBCs are simulated as discrete, deformable structures interacting with the flow to investigate the interplay between flow in microvessels, RBC dynamics, and the wall shear stress on both a flat and wavy endothelial surface. There have been numerous previous computational studies which focused on different aspects of normal and pathological RBCs in flow, in addition to those mentioned above, [START_REF] Ye | Particle-based simulations of red blood cells-A review[END_REF][START_REF] Ye | Numerical modelling of a healthy/malariainfected erythrocyte in shear flow using dissipative particle dynamics method[END_REF][START_REF] Wu | Simulation of malaria-infected red blood cells in microfluidic channels: Passage and blockage[END_REF]Fedosov et al. 2011a;[START_REF] Bagchi | Mesoscale simulation of blood flow in small vessels[END_REF]; however, the interaction of pathological RBCs and the role of RBC rigidity on vessel shear stress has yet to be widely investigated. We show that RBC morphology has a profound effect on the shear stress at the endothelial surface and that endothelial waviness further amplifies these effects.

Bouzidi boundary condition

In the Methods Chapter, we discussed the implementation of the "half-way bounceback scheme" in LBM in order to give the Dirichlet boundary condition on flat walls. While this works well when the wall is aligned with the Cartesian mesh and located halfway between the first solid node and its neighboring fluid node, it is insufficient on a curved boundary. Since the fluid is represented on a regular Cartesian grid, curved boundaries must be approximated by staircase profiles which are aligned with the coordinate grid. The resulting inaccuracies are particularly worrisome when high local resolution is desired, which is the case when searching for the local shear stress at the wall. In order to increase the accuracy on a curved boundary, [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] proposed a treatment of walls which involves the bounce-back scheme and interpolations, depending on the distance between the wall and the fluid node adjacent to the wall. A schematic of the Bouzidi boundary condition is shown in Fig 3 .2. Fluid nodes are shown on the left: nodes F, E, A, such that x A = x, x E = xc i , and x F = x -2 c i ; the wall is shown as the vertical dotted line, and the wall node is shown as B. The location of the wall is given by q = |AC| /|AB|. The Bouzidi boundary condition consists of two conditional cases, depending on the value of q. Below we present the quadratic interpolation formulation: The procedure involves identifying all of the links of the lattice that cross the solid boundary so that x is a fluid node and x + c i is a wall node (i.e., direction i is streaming into the wall). Then the distance is calculated and the condition is satisfied by the appropriate equation. As indicated, the RHS of Eq 3.1 uses post-collisional distributions, and the distribution on the LHS is the postcollision, post-streaming value (i.e., the value after one complete time step). The reader may note for completeness that the interpolation formula reduces to the half-way bounceback when q = 1/2. We implemented this treatment for the cases involving wavy walls.

f i * ( x A , t + 1) = q(2q + 1) fi ( x A , t) + (1 + 2q)(1 -2q) fi ( x E , t) -q(1 -2q) fi ( x F , t), q < 1 2 f i * ( x A , t + 1) = 1 q(2q + 1) fi ( x A , t) + (2q -1) q f i * ( x A , t) + (1 -2q) (1 + 2q) f i * ( x E , t), q ≥ 1 2 (3.1)
Chapter 3. Results: Red Blood Cells under Channel Flow: Effect of Endothelium and Red Blood Cell Morphology

Wall shear stress on a curved boundary

When fluid flows over a rigid surface, the velocity vanishes at the wall due to the no-slip condition. In the following, we derive the wall shear stress on a curved boundary. The stress tensor in fluids represents the sum of pressure and viscous terms, given by

σ ij = -pδ ij + σ ij (3.2)
where δ ij is the Kronecker delta and

σ ij = µ(∂ i u j + ∂ j u i )
is the deviatoric stress tensor such that i = j and where Einstein summation convention is applied and i, j = (x, y). The Cauchy formula gives the full stress vector, T

i , on the surface with a normal vector, n, such that

T (n) i = σ ij n j (3.3)
The tangential shear stress is given by

f t = t i σ ij n j (3.4)
Applying the Einstein summation convention for the indices in two dimensions leads us to

f t = t x σ xx n x + t x σ xy n y + t y σ yx n x + t y σ yy n y (3.5)
Substituting Eq. 3.2 into Eq. 3.5, we find

f t = t x (-p)n x + t x σ xy n y + t y σ yx n x + t y (-p)n y (3.6)
In terms of velocities, Eq. 3.6 becomes

f t = t x (-p)n x + t x [µ(∂ x u y + ∂ y u x )]n y + t y [µ(∂ y u x + ∂ x u y )]n x + t y (-p)n y (3.7)
The pressure term disappears and using the commutative property, Eq. 3.7 can be rewritten as

f t = µ[n y (∂ x u y + ∂ y u x )t x + n x (∂ y u x + ∂ x u y )t y ] (3.8)
In the numerical method, we use the approximation that the tangent vector to the surface does not vary much to the point where we compute the interpolated velocity from the coordinate nodes, such that (∂ i u j )t j = ∂ i (u j t j )u j (∂ i t j )where u j (∂ i t j ) ∼ 0 so that (∂ i u j )t j ∼ = ∂ i (u j t j ). We can therefore rewrite Eq. 3.8 as

f t = µ[n y (∂ y u x t x ) + n x (∂ x u y t y )] (3.9)
which can be rewritten in vector form as

f t = µ ∂( u • t) ∂n (3.10)
where u is the velocity vector at an off-grid position close to the wall and t is the tangent to the curved surface. This interpretation involves only the interpolation of the velocity vector at the off-grid position. Since we have modeled the wavy surface as a sine function, the tangent can be found analytically. It has been shown that the stress tensor can be calculated at each lattice node [START_REF] Krüger | Shear stress in lattice Boltzmann simulations[END_REF])Matyka2013, such that

σ ij = 1 - 1 2τ ∑ k c ki c kj - δ ij 2 c k • c k f neq k (3.11)
where f neq k

represents the nonequilibrium populations given by f neq k

= f k -f eq k .
Applying this equation properly would have led to the interpolation of the 9 distribution functions. I instead chose to use the macroscopic equation given by Eq. 3.10 implemented with a finite-difference scheme.

Implementation in LBM

As we have discussed, the presence of endothelial cells gives rise to a wavy or corrugated wall surface. For our model, we use a sine function to give a continuous, undulated wall layer, as shown in Fig. 3.3, where we can define the height of the endothelial cell, or the amplitude, A, of the sine function, the length of a single endothelial cell, or the wavelength, λ endo , and we can define the alignment of the the top and bottom curves by the phase shift, n, where 0 ≤ n ≤ 1 so that we can define the wavy boundary as Therefore, position x (represented by the blue star) can be identified as a fluid node adjacent to the wall having a lattice link intersecting the boundary, x F . If we were to draw a line along link 7 (red line), we would find the following equation of the line:

y boundary (x) = A 2 sin 2π x λ endo + n + A 2 (3.12)
y l7 (x) = x + (y F -x F ) (3.13)
We search for the point of intersection such that

y boundary (x) = y l7 (x) (3.14)
so that we need to solve

A 2 sin 2π x λ endo + n + A 2 = x + (y F -x F ) (3.15)
This yields a transcendental equation with no closed-form solution and thus numerical methods are required to solve for x. We use the Newton-Raphson method, defining the function, f (x), for which we would like to find the roots. In this case,

f (x) = A 2 sin 2π x λ endo + n + A 2 -x -(y F -x F ) (3.16)
To apply the Newton-Raphson method, we begin with an initial guess, x 1 = (x F + x S )/2. The second approximation involves the derivative of f (x) such that x 2 = x 1f (x 1 )/ f (x 1 ). Successive approximations are made until f (x) reaches a value below the tolerance, tol; that is, f (x) < tol.

For the calculation of the wall shear stress, we apply Eq. 3.10. We will briefly discuss the details of implementation. For this equation, we discretize the "real-valued" sine function at points (x s , y s ), and we define a distance d n from the real-valued sine function such that it lies on the normal to the curve at (x s , y s ) and then we search for the coordinates of this point, (x d , y d ), such that d n = (x sx d ) 2 + (y sy d ) 2 where y d = y boundary (x s ) -1 x d + (y sy boundary (x s ) -1 x s ). We can use this equation to solve for the coordinates of this point using the quadratic equation. We then complete a bilinear interpolation to find the value of the velocity at this coordinate position u(x d , y d ). We can then apply Eq. 3.10, using a finite difference scheme such that

f t (x s , y s ) = µ (n y (x s , y s )u x (x d n , y d n )t x (x s , y s ) + n x (x s , y s )u y (x d n , y d n )t y (x s , y s ) d n y boundary (x s ) 2 + 1 2
(3.17)

Note that we want to choose the value of d n not too far away from the curve so that the approximation to the tangent vector is still valid, but we need to choose d n such that it falls between 4 points on the grid in order to perform the bilinear interpolation, the principle of which is illustrated in Fig. 3.5.

Results

Simulations were carried out with a D2Q9 lattice grid with periodic boundary conditions in the direction of flow induced by imposing a pressure gradient across the inlet and outlet of the channel and with the Dirichlet no-slip boundary condition at the walls. The microvessel is first represented by a rigid, flat-walled channel and later on by wavy walls given by sinusoidal functions. The RBC is modeled in 2D as an elastic membrane filled with a viscous fluid as described previously. In our simulations, the fluids separated by the membrane were taken to have the same density and viscosity. In LBM units, R 0 = 3µm is represented by 15 discrete meshes in consideration of both efficiency and accuracy. During all simulations in this study, the Reynolds number was held constant such that Re = UR 0 /ν = 0.1 to approximate the Stokes flow in microcirculation. The channel length is chosen as L = 50R 0 in our simulations. We have validated that this length is long enough to eliminate artifacts due to periodic boundary conditions.

Single RBC in Flat-Walled Channel

We begin by considering an individual RBC flowing through a 12 µm-wide flat-walled channel and explore the impact of the RBC on the wall shear stress. It has been shown that in terms of the morphology and dynamics of a single RBC in Poiseuille flow, rich dynamics can be seen [START_REF] Tahiri | On the problem of slipper shapes of red blood cells in the microvasculature[END_REF]. As has been previously reported, the morphology of RBCs in confined Poiseuille flow is a function of the three dimensionless parameters Ca, C n , and α defined above. Among the stable RBC morphologies, several are permanent including the parachute. We begin by analyzing the effect on the flow of a single RBC with this morphology. Figure 3.6a demonstrates the perturbation or induced flow due to the presence of the parachuteshaped RBC. In our simulations, the shapes of the RBCs are the result of an interplay between the bending and tension forces of the RBC and hydrodynamic forces. For a symmetric parachute shape, there is no tank-treading or normal displacement. Thus, when we plot the velocity in the co-moving frame with the RBC in Fig. 3.6b, we see zero flow field inside, like a rigid particle having the same shape. We investigate the sign of the induced flow for information about how the flow is perturbed due to the presence of the RBC. The total flow (TF) is the sum of the 3.5. Results

45

imposed flow (ImF) and the induced flow, or the perturbation flow (PF), such that

TF = ImF + PF (3.18)
The sign of PF (positive in the direction of the ImF which is left to right; negative if against the ImF) indicates how the flow is being perturbed by the presence of the RBC. We see in Fig. 1a that the PF generally acts to accelerate the flow near the wall, except at two recirculation zones where the perturbation has resulted in local reversal of the induced flow at the wall. These two zones correspond to the point where the cell is closest to the wall, an observation which was also previously reported in [START_REF] Freund | The wall-stress footprint of blood cells flowing in microvessels[END_REF][START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF][START_REF] Xiong | Shear stress variation induced by red blood cell motion in microvessel[END_REF]. These zones of flow reversal are not observed in the unbounded geometry [START_REF] Ghigliotti | Why and how does collective red blood cells motion occur in the blood microcirculation?[END_REF].

Previous work has reported how changes in the degree of confinement (value of C n ) in Poiseuille flow alter RBC morphology [START_REF] Tahiri | On the problem of slipper shapes of red blood cells in the microvasculature[END_REF] .

We varied C n in the range of 0.2 -0.7 and examined the effect of three stable RBC morphologies, namely the unconfined slipper, the parachute, and the confined slipper, on the wall shear stress. Figure 3.8 shows how the RBC morphology affects the wall shear stress footprint normalized to the Poiseuille solution of flow through parallel plates, such that:

τ τ Pois = τ (dP/dx) im W/2 (3.19)
where τ is the wall shear stress measured along the wall obtained from simulations, τ Pois is the wall shear stress obtained from the Poiseuille solution of flow through parallel plates, and (dP/dx) im is the imposed pressure drop along the channel axis per unit length. Figure 3.8 demonstrates that a "peak-valley-peak" form of the wall shear stress generally emerges. This form of the wall shear stress was previously observed in [START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF][START_REF] Xiong | Shear stress variation induced by red blood cell motion in microvessel[END_REF]. [START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF] perform a lubrication analysis for a cylinder moving close to a wall in order to investigate the wall shear stress. In obtaining an expression for the pressure along the wall, they find that it is the reverse pressure gradient generated between the moving cell and wall which results in the curious waveform observed in the wall shear stress. In light of their analysis, we extracted the pressure gradient along the wall for each of the three morphologies investigated. In Fig. 3.9, we see the features of the WSS are reflected in the plots of the pressure gradient along the wall. Similarly to [START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF], we observe the largest pressure gradient opposite to the RBC motion at the narrowest part between the RBC and the wall. This reverse pressure gradient generated between the RBC and the wall impairs the flow velocity and the WSS, consistent with the flow reversal at the narrowest gap observed in Fig. 1a of the induced flow. We find the behavior of the WSS captured in the profile of the pressure gradient at the wall for the different morphologies. We see that the pressure gradient and WSS variations are more profound as the size of the channel is decreased, also in accord with observations of [START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF].

A summary of the characteristics of the wall shear stress variation is given in Table 1 for the different results presented in Fig. 3 wall shear stress observed along either channel wall, while max(τ/τ Pois ) corresponds to the maximum wall shear stress observed in the channel, IR represents the influencing range, or the width over which the WSS perturbation to the Poiseuille value persists, such that |(τ(x)τ Pois )/τ Pois | ≥ 0.02 from the trailing and leading edges of the cell, corresponding to a minimum of 2% perturbation from the nominal shear stress, and T is the time duration of the perturbation for a fixed point on the endothelium calculated by T = IR/V RBC where V RBC is the velocity of the RBC. We see that channel size has an inverse effect on the amplitude of WSS perturbation.

For the largest channel (case a), we see the smallest variation 0.98 -1.01, while for the smallest channel (case c), we see the largest variation 0.83 -1.26. Considering that the Poiseuille solution for shear stress in the small channel is higher than that in a wider channel with the same mean flow velocity, the absolute variation of the shear stress in the smaller channel is even more significant. The influencing range IR cannot be defined in the case of the unconfined slipper, case a, because the perturbation is always below the threshold. It is largest for the parachute morphology, since we see that its range extends beyond the length of the cell (the second peak is found well after the trailing edge of the cell in Fig. 3.8b)). The influencing range of the confined slipper is different for the top and bottom walls, since the morphology is not axisymmetric. Lastly, the RBCs travel at roughly the same velocity so that the time duration results largely reflect the differences in influencing range.

In order to understand the effect of the morphology of RBCs on the shear stress at the wall, we varied the reduced area, which represents the degree of "swollenness" of the cell. By varying the reduced area, we can simulate the morphologies of different potential pathological conditions of RBCs. A relatively low confinement vessel, C n = 0.42, representing a vessel with a width of 12 µm is investigated first. As illustrated in Fig. 3.10a, the RBC shapes obtained are generally asymmetric. Blood Cell Morphology the centerline determines whether the values of the shear stress on the top and bottom walls are equal or not. While the maximum shear stress values can range from several to ∼ 30% greater than the Poiseuille solution in the absence of the RBC; the minimum wall shear stress values vary less than ∼ 11% over the same range of reduced areas. We hypothesized that the dependence of the maximum wall shear stress on reduced area is determined by the extent of the cell in the channel in the perpendicular streamwise direction. To test this hypothesis, we determined the distance between the closest point of the cell and the channel wall. In Fig. 3.10c, we plot a quantity related to the instantaneous cell-free layer, or the RBC-free layer of plasma next to vessel walls, such as that previously reported in [START_REF] Oulaid | Temporal and Spatial Variations of Wall Shear Stress in the Entrance Region of Microvessels[END_REF]. We calculate from the instantaneous cell-free layer, the minimum distance between the RBC membrane and the wall, and plot the inverse of this normalized quantity in Fig. 3.10c. We observe that the inverse normalized minimum distance exhibits exactly the same dependence on reduced area as the maximum wall shear stress (cf: Fig. 3.10b).

Results

51

To investigate how robust the results shown in Fig. 3.10 are, we increased the value of Ca by a factor of 100 (from 0.1 to 10), effectively creating a "softer cell" and repeated the simulations for the same confinement. of RBC rigidity. We see that the point closest to the wall is further from the wall for larger reduced areas, shown in Fig. 3.11c. We again observe a larger variation in maximum wall shear stress, ∼ 17%, while the minimum wall shear stress varies only at most ∼ 10% from the no-cell Poiseuille solution.

Comparing Figs. 3.10b and 3.11b, we observe that for the same confinement, the wall shear stress is lower for a softer RBC (except for α ≈ 1 in which case the morphologies are the same and therefore the distance to the wall equal). [START_REF] Freund | The wall-stress footprint of blood cells flowing in microvessels[END_REF] also observe that rigid cells induce greater perturbations to the wall shear stress. In our analysis, we show that the morphologies adopted for a larger capillary number have minimum distances that are further from the wall than those observed for the morphologies adopted in the case of a smaller capillary number. The cells become more elongated in the perpendicular flow-wise direction in the case of the smaller capillary number and therefore the distances between the cell and the wall are smaller.

Finally, we examine the effect of confinement by considering a smaller vessel of 7 µm in width.

55

Figure 3.12a demonstrates that an asymmetric slipper shape is observed for α < 0.80. At α = 0.80, the slipper shape becomes axisymmetric, and for all reduced areas beyond this value, symmetric morphologies are obtained. The difference in RBC morphology behavior above and below α = 0.80 leads to the branches found in the plots of the maximum and minimum wall shear stresses observed in Fig. 3.12b, where the maximum shear stress along the top wall is identical to that along the bottom wall for α > 0.80 but not so for reduced areas below this cutoff. Again, the maximum shear stress shows more variation, greater than 50%, while the minimum wall shear stress varies at most ∼ 21% for the same range of reduced areas. Finally, Fig. 3.12c demonstrates that here again, (d min /W) -1 is the physical determinant of the wall shear stress behavior observed in Fig. 3.12b. Comparing Fig. 3.12b to Figs. 3.10b and 3.11b, we see that the maximum shear stresses are higher in a smaller vessel than in a larger vessel. We also see again that RBCs are closer to the wall in the case of the smaller channel. Interestingly, we see that among the reduced areas, that corresponds the most to a circle exhibits the smallest maximal wall shear stress fluctuation in the case of the 12µm channel; however, in the case of the smaller 7µm channel, it exhibits a much larger wall shear stress than all the other reduced areas due to its greater extension into the channel than the other RBCs which instead elongate in the streamwise direction. This behavior could be of interest in the study of the schizont stage in malaria.

We recast the results in Figs. 3.10, 3.11, and 3.12 onto a single graph plotting τ max /τ Pois vs (d min /W) -1 . The inverse minimum distance normalized to the width of channel appears to be an accurate predictor of the normalized maximal wall shear stress in the channel, irregardless of cell rigidity and channel size. We note one possible cause of variation from a strictly linear relationship could be from the different resistance produced by different morphologies of RBCs which would result in different nominal shear stresses in the channel for a given imposed pressure gradient for which our model does not take into account.

Effect of a Wavy Wall on Wall Shear Stress

At the scale of the microcirculation, the endothelium's wavy topography becomes important, and a pertinent question then becomes what effect the vessel waviness has on the flow field. Sinusoidal curves are implemented as boundaries. The magnitude of shear stress can be extracted from simulations using bilinear interpolation of velocities located on grid points close to the wall, as explained in the Modeling Methodology section. For the time being, we will keep the sine curves on the top and bottom boundaries in phase. The effect of dephasing will be studied later. We simulate a channel with wavy walls, with effective width of 12 µm and amplitude and wavelength of undulation A = 2 µm and λ = 30 µm, respectively. Figure 3.14 depicts the magnitude of wall shear stress at the wavy wall normalized to the Poiseuille solution for a straight channel with the same width as the effective width of the wavy channel, where effective width is the mean width over one period. In agreement with previous studies [START_REF] Satcher | The Distribution of Fluid Forces on Model Arterial Endothelium Using Computational Fluid Dynamics[END_REF][START_REF] Barbee | Subcellular distribution of shear stress at the surface of flow-aligned and nonaligned endothelial monolayers[END_REF], the simulations demonstrate that the maximal shear stress can be > ±30% that of the Poiseuille value for a straight channel with the same effective width, the extremes corresponding to the crests and valleys of the wall.

Results

57

RBC in a Wavy-Walled Channel

We examined the effect of confinement for two Ca numbers in wavy-walled channels where the waviness of the top wall was in-phase with that of the bottom wall. As the RBC travels across one wavy period, changes in the local flow field result in changes to its morphology. We examined the morphologies obtained and the extremal wall shear stresses. Fig. 3.15 presents for the lowest and highest reduced areas investigated, the wall shear stress and its extremal values along one period of the wall undulation such that τ max = max(τ(x)/τ w,nocell (x)) and τ min = min(τ(x)/τ w,nocell (x)) where τ w,nocell is the wall shear stress obtained in the wavy-walled channel in the absence of RBC. We see that the average WSS is roughly equal to the steady-state, continuum solution in the absence of the RBC. However, the time-averaged minimum and maximum shear stresses reveal important differences in subcellular extremal WSS values. In the case of α ≈ 1, the variation over one period is greatest, the maximum varying greater than 20% over one period, while in the case of α = 0.6, there is less than a 10% difference in extremal values over one period. It

Results

59

has been speculated that subcellular gradients are potentially important in mechanotransduction processes, and our results show that different reduced areas of cells lead to different timeaveraged subcellular gradients in extremal WSS values [START_REF] Satcher | The Distribution of Fluid Forces on Model Arterial Endothelium Using Computational Fluid Dynamics[END_REF][START_REF] Barbee | Shear Stress-Induced Reorganization of the Surface Topography of Living Endothelial Cells Imaged by Atomic Force Microscopy[END_REF]).

As we did in the case of the flat-walled channel (cf: Figs. 3.10-3.12), we examine the behavior of the extremal shear stress as the reduced area is varied. This is performed for the case of a relatively rigid RBC (Ca = 0.1) and a relatively soft RBC (Ca = 10). The extremal wall shear stresses are plotted as a function of reduced area in Fig. 3.16a. .16a shows that the maximum wall shear stress in the wavy wall case is less sensitive to changes in reduced area than in the flat wall case. In extreme, the maximum wall shear stress only varies by ∼ 10%, whereas in the flat wall case it can vary by ∼ 30% among reduced areas. Waviness reduces the variation in maximum shear stresses observed for different reduced areas; however, the values are in general higher than the equivalent in the straight case. The inverse of the minimum distance normalized to the effective width of the channel is not predictive of the maximum wall shear stress as is suggested by Fig. 3.5. Results

61

Figure 3.17 shows snapshots of one RBC as it flows over one period of wall undulation. We see that what would be a stationary morphology in the case of a flat-walled channel is no longer so due to the variation of the flow field along the undulation. Because the RBC is effectively "rigid," it does not undergo a noticeable deformation in either the streamwise or perpendicular streamwise directions as it passes through the constriction and expansion of one wall undulation. Instead, we see that the tilt angle of the RBC changes over one period. As the cell approaches the maximal constriction of the two walls, the RBC tilts towards the horizontal axis, and as the cell exits the constriction, the tilt moves toward the vertical axis. The results above were for a rigid RBC (Ca=0.1). For Ca=10, representing a softer cell, all morphologies obtained are axisymmetric, so the maximum shear stress along the top and bottom walls are identical, as shown in Fig. 3

63

The morphologies of the RBC as it passes over one wall undulation are always axisymmetric as shown in Fig. 3.18. We can analyze the passage of the cell over one wall undulation. Before entering the constriction, the cell is stretched. The cell is "squashed" in the middle of the constriction. Upon exiting, the cell experiences a rebound, whereby the rear-end of the RBC moves faster than the front end so that we see stretching in the perpendicular streamwise direction. The RBC elongates along the streamwise direction in the constriction and expands in the perpendicular streamwise direction in the valley. These same dynamics were also observed for capsules passing through constrictions [START_REF] Rorai | Motion of an elastic capsule in a constricted microchannel[END_REF].

All the wavy wall results presented thus far were for the case where the undulations in the top and bottom walls are in phase. In this case, RBCs pass periodically through areas of strong confinement (at peaks) and weak confinement (at valleys). This scenario is unlikely to be the most representative configuration in real situations. This has naturally led us to study the effect of dephasing between the two walls. We dephased the sine curves at the top and the bottom walls of the channel and examined the maximum shear stresses experienced as a function of the dephasing angle, n, in Eq 3.12. The morphologies of one cell as it passes over one wall undulation are shown in Fig. We see the effect of dephasing on the morphology of the RBCs in this case where Ca = 10. For this case, when the walls were in-phase, the RBC maintained axisymmetric morphology at all times while passing over the length of one period; however, now the RBC shows asymmetric morphologies when the top and bottom walls are out-of-phase. The trailing edge of the parachute elongates in the direction of the further wall. In order to see what effect the dephasing has on the wall shear stress, we plot the maximum wall shear stress as a function of the dephasing angle in Fig. 3.20. The maximum wall shear stress takes its largest value when the top and bottom walls are inphase and its smallest value when the walls are half a wavelength out-of-phase. We observe that the effect of dephasing can decrease the maximum wall shear stress in the channel by nearly 30%. Although the properties of the RBCs are held constant as is the applied pressure gradient and the effective width of the channel, the local changes to the flow field resulting from undulation dephasing have a significant effect on the observed wall shear stress. Thus, even small changes to the geometry of the endothelial surface have important ramifications on RBC morphology and the resulting wall shear stress.

Cell Trapping

In this section, we investigate undulated channels for which there exist solutions in which a cell is "trapped." Cell trajectories in general are the defining element behind deterministic lateral flow displacement devices. The idea behind the deterministic lateral flow displacement device is to Chapter 3. Results: Red Blood Cells under Channel Flow: Effect of Endothelium and Red Blood Cell Morphology passively sort cells based upon their characteristics, such as size [START_REF] Davis | Deterministic hydrodynamics : Taking blood apart[END_REF]. Our results further show that differences in particles' deformability will lead to different trajectories. Krüger et al. 2014 previously explored this finding with simulations of healthy and malaria-infected RBCs which tend to be stiffer than their healthy counterparts. Microfluidic devices which can detect diseased particles based on the their trajectories has the potential to save expensive and time-consuming biochemical preparation steps.

The channel we investigate here has A/W = 0.4 and λ/A = 2. Here we define the Reynolds number as Re = U A/ν. We find that for large enough amplitude, recirculation zones appear, even at low Reynolds number. Varying, the Re number in the channel, we observe for a fixed amplitude, the recirculation zones increase with increasing Re number in a linear fashion, as illustrated in Fig. 3.21f. We fix Re = 1.5 and now consider varying initial positions of a rigid RBC in the channel (Ca = 0.25) and the trajectories of their center of mass in the y-direction. We see the RBC trapped in the valley at y0 = 25 in Fig. 3.22d. In Fig. 3.22a, we see its center of mass remains stationary in time. When the RBC is placed just slightly higher in the valley at y0 = 30, it eventually escapes, and when it is placed at mid-amplitude height in the valley, y0 = 50, escapes more quickly. In Fig. 3.22b, we see that the center of mass of the RBC is affected by the undulations of the wall. When the RBC is placed just level with the peak of the undulation at y0 = 100, its oscillations are large in the beginning and then dampen. For the cell initially placed just off-center of the middle of the channel at y0 = 145, we see that its oscillations are stable in time as it transits past the peaks the valleys of the channel walls. For an RBC placed with its major axis along the centerline (y0 = 150, 0 • ), the center of mass is stationary. When the RBC is placed at the center of the channel with a 45 • angle to the horizontal, it oscillates but remains in the vicinity of the center; however, a cell placed with a 90 deg appears to move off-center.

We then increase the Ca number by twenty-fold, creating an effectively "softer" cell (Ca = 5). , we see striking differences in the trajectories of the "rigid" and "soft" cells. When a "soft" RBC is placed at the centerline with no angle to the horiztonal, its center of mass will remain stationary in the beginning; however, eventually, it will move off-center. When the RBC is placed with a 90 • angle, it will exhibit steady oscillations of its center of mass, this is due to the periodic morphological changes it goes through as it traverses one period of the wall undulation. These morphological changes are shown in Fig. 3.24. We see that the trajectory of the center of mass of an RBC will be mediated by its initial position and upon its deformability.

Conclusion

We have studied the effect of morphological changes in RBCs, confinement, capillary number, reduced area, and vessel wall geometry on the wall shear stresses in microvessels. We have computed the possible shear stress fluctuations that each of these physical parameters might induce using a simulation model for flow in the microvasculature with the idea that these fluctuations could be important in mechanotransduction events. Our results have demonstrated that the presence of RBCs has a significant effect on the wall shear stress. Not surprisingly, the effect of RBCs on wall shear stress becomes more pronounced with increased confinement (i.e. smaller vessels), and the effect becomes more dynamic under conditions that induce fluctuations in RBC morphology. Interestingly, we have determined that the maximum wall shear stress is correlated to the minimum distance between the cell and the adjacent wall in small confinements. This observation runs counter to previous notions that the wall shear stress is directly proportional to the local distance between the RBC and the wall (i.e. τ = µV RBC /δ where δ is the local distance between the RBC and the wall) [START_REF] Namgung | Effect of cell-free layer variation on arteriolar wall shear stress[END_REF].

Another important result of the present work stems from accounting for the waviness of the endothelial surface. We determined that in the case of wavy walls, the exact geometry of the endothelium changes RBC dynamics and hence has a significant impact on the shear stresses experienced along the wall, highlighting the importance of experimental studies focusing on the topography of confluent endothelial monolayers [START_REF] Chtcheglova | AFM Functional Imaging on Vascular Endothelial Cells[END_REF]) so that a more realistic geometry could be incorporated into future modeling. We have seen that in the cases studied, the effect of waviness is to reduce the variation in maximal wall shear stress among different reduced areas when compared to the flat wall counterpart. In the final section of the chapter, we briefly investigate how deformability mediates an RBC's trajectory in an undulated channel, which could be potentially useful in the design of microfluidic devices. We find that with large enough amplitude, recirculation zones appear between the valleys of the undulated wall and these can be used to trap cells. work of Fahreus and Lindqvist in which a striking decrease in apparent blood viscosity in glass tubes with diameters ranging between ∼ 500 -50 µm was observed [START_REF] Fahraeus | The viscosity of the blood in narrow capillary tubes[END_REF]. The phenomenon was confirmed in vitro by many investigators who demonstrated that the decrease continued down to diameters of ∼ 10 µm [START_REF] Pries | Blood viscosity in tube flow: dependence on diameter and hematocrit[END_REF][START_REF] Goldsmith | Robin Fåhraeus: evolution of his concepts in cardiovascular physiology[END_REF][START_REF] Reinke | Blood viscosity in small tubes: effect of shear rate, aggregation, and sedimentation[END_REF]). The CFL is also related to another important characteristic of blood flow, known as the Fahraeus Effect [START_REF] Fahraeus | The suspension stability of the blood[END_REF]. In Poiseuille flow, RBCs and plasma do not travel at the same average velocity, which results in important differences in tube and discharge hematocrits, whereby tube hematocrit is defined as the volume fraction of RBCs to the total volume of the vessel measured before entrance into the channel and discharge hematocrit is the ratio of the RBC flux to the whole-blood flux (or suspension flux) exiting the channel. The Fahraeus Effect refers to the observation that in blood flow, the tube hematocrit is consistently smaller than the discharge hematocrit. This effect is directly correlated with cell migration to the centerline, forming the CFL and RBC-dense core which results in the mean RBC velocity being higher than the mean blood velocity. The CFL, Fahraeus-Lindqvist, and Fahraeus effects have been welldemonstrated for straight channels since the seminal work of Fahraeus and Linqvist [START_REF] Popel | Microcirculation and Hemorheology[END_REF].

Because bulk blood rheology is dependent on the interactions at the level of a single RBC, we were interested in the effect of the local geometry of the endothelium on hemodynamic and hemorheological properties, which have been demonstrated to depend heavily on the geometry of the conduit through which blood flows [START_REF] Popel | Microcirculation and Hemorheology[END_REF]. As discussed in Chapter 3, the endothelium has been shown to have a height profile varying up to several microns along the length of an individual cell, a feature which could become especially significant in the microvasculature. The influence of the non-flat geometry of the vessel walls on rheological properties of bloof has not currently received any attention in the literature. In Chapter 3, we focused on the effect of single cells specifically on the shear stress at the endothelial surface. In this chapter, we will study suspensions of RBCs in the microvasculature and the effect of wavy-walled channels on rheological and dynamical properties. In particular, we investigate their dependence on channel width and hematocrit, for which the straight case has received much attention. We assess the effects of RBC concentration and channel size on blood flow in wavy-walled channels. Our results show the formation of a CFL next to undulating channel walls and a decrease in viscosity as a function of channel diameter, corresponding to the well-known Fahraeus-Lindqvist effect, as well as an increase in discharge hematocrit corresponding to the Fahraeus Effect. This represents a first step toward incorporating key elements, including RBCs and a more realistic endothelial geometry, into a comprehensive model of blood flow.

Results

We present the results of simulations of blood flow in microchannels with non-flat and flat walls for comparison. In the case of non-flat walls, the top and bottom walls were modelled as sine functions in phase with an amplitude of 3 µm and a wavelength of 30 µm. The channels varied in width ranging from 12 µm to 57 µm. The properties of the cell suspensions were evaluated for hematocrits varying from 5 to 40%. The channel hematocrit, or in our 2D simulations, the cell area fraction, φ, in the channel can be calculated by:

φ = N c A rbc A ch (4.1)
where N c is the number of RBCs each with an area A rbc in a channel whose area is A ch = A domain -A walls where A domain is the area of the entire computational domain and A walls is the area occupied by the wavy walls if present.

Fahraeus-Lindqvist Effect

The non-Newtonian nature of blood flow makes the investigation of its rheological properties a subject of intense study. The components of blood, including RBCs, white blood cells, and platelets are suspended in plasma, an aqueous solution that is generally considered to follow Newtonian dynamics. It has been shown that the consideration of RBCs in blood flow leads to important consequences for hemodynamic and hemorheological properties. One of the most important properties of blood is the ease with which it flows. We quantify flow resistance in terms of apparent viscosity and relative viscosity, which relate blood flow to Newtonian flow. To find the apparent viscosity of blood containing RBCs, we assume a Poiseuille relation between the imposed pressure gradient ∆P/L and the average volume flux per unit time over the channel length (i.e. volume flow rate per unit width in the z-direction), Q, such that :

µ app ∼ ∆PW 3 12QL (4.2)
where ∆P is the pressure drop in length L and W is the width of the channel. The dimensionless relative apparent viscosity is found as:

µ rel = µ app µ ext (4.3)
where µ ext is the viscosity of the plasma (which is taken to be the kinematic viscosity of water at room temperature, or 10 -6 m 2 /s). All else being equal, we can compute the relative viscosity as:

µ rel = Q p Q s (4.4)
where Q p represents the average volume flux in the channel in the absence of RBCs and Q s is the average volume flux of the suspension in time found by:

Q s = 1 T T 1 L T T L W u(x, y)dydxdt (4.5)
where u(x, y) is the velocity at point (x, y), T T is the total time over which the volume flux is averaged, L is the length of the channel, and W is the width of the channel. For straightwalled channels, Q p could be found analytically from the Poiseuille solution given by Eq. 4.2.

For wavy-walled channels, we calculate it from simulations. Apparent and relative viscosity are not intrinsic properties of blood but rather vary with RBC properties, hematocrit, and vessel geometry. Using Eq. 4.2, apparent viscosity can be measured in any flow in which the Newtonian solution is known; we need only measure flow rate and pressure with and without RBCs to calculate the relative viscosity.

The seminal work of Fahraeus and Lindqvist involved flow of blood in in vitro glass tubes. It was observed that in tubes <∼ 200 µm in diameter, the apparent viscosity strikingly falls with decreasing diameter until reaching a minimum at diameters ∼ 5 -7 µm [START_REF] Fahraeus | The viscosity of the blood in narrow capillary tubes[END_REF]. As tube diameter increases further, the apparent viscosity spikes as the size of the tube becomes considerably smaller than the size of the RBC. The drop in apparent viscosity is related to the migration of RBCs toward the tube centerline, creating a layer depleted of RBCs next to the channel walls, known as the CFL. In vivo experiments of blood flow in microvessels have shown the existence of the Fahraeus-Lindqvist effect arising from the CFL formed next to the vessel walls [START_REF] Pries | Resistance to blood flow in microvessels in vivo[END_REF][START_REF] Kim | Temporal and spatial variations of cell-free layer width in arterioles[END_REF][START_REF] Maeda | Erythrocyte flow and elasticity of microvessels evaluated by marginal cell-free layer and flow resistance[END_REF]). However, it was observed that in vivo, blood flow resistance was strikingly higher than that in microtubes [START_REF] Pries | Resistance to blood flow in microvessels in vivo[END_REF]. In contrast to glass microtubes, microvessels are elastic, relatively short, may be irregular in shape, and are lined by endothelial cells and the glycocalyx, the thin layer of proteins 60 -90 nm thick coating the luminal surface of the endothelium. These differences may account at least in part for the observed increased effective viscosity of blood in vivo measurements by attenuating RBC migration to the center of the vessel. How each of these factors plays into the dynamics of the RBCs is important for a fundamental understanding of in vivo blood flow. [START_REF] Pries | Blood viscosity in tube flow: dependence on diameter and hematocrit[END_REF] have proposed the following empirical formula for the variation of relative apparent blood viscosity in tubes with hematocrit φ d and diameter D:

µ rel = 1 + (µ rel,0.45 -1) (1 -φ d ) C -1 (1 -0.45) C -1 (4.6)
where µ rel,0.45 = 220e -1.3D + 3.2 -2.44e -0.06D (4.8) Equation 4.6 was obtained from an extensive analysis of the available data in literature as well as their own experimental data for blood flow through tubes ranging in diameter from 3.3 µm to 1978 µm at high shear rates and at hematocrits up to 90%. Figure 4.1 illustrates the dependence of the relative apparent viscosity on hematocrit and channel width for flat-walled channels as calculated from our simulations and from the empirical formula provided by Pries et al. The results demonstrate a generally monotonic increase in relative apparent viscosity with hematocrit as does the empirical relationship of Pries et al. (dashed lines), except for W = 12 µm which we will come back to later. At equivalent hematocrits, the relative apparent viscosity is Chapter 4. Suspensions in wavy-walled channels higher in wider channels than in smaller channels, in accordance with the Fahraeus-Lindqvist Effect. Our two-dimensional simulations actually simulate flow between two parallel plates and not that in a cylindrical tube, and in order to apply the Pries et al. relation given by Eq. 4.6, we have interchanged D for W. For this reason, we do not expect an exact match between the simulations and empirical relations; however, the non-linear effect of hematocrit on apparent viscosity is apparent from Fig. 4.1.

One may remark a striking feature of the 12 µm channel: the apparent viscosity increases from φ = 0.05 to φ = 0.15 and then decreases from φ = 0.15 to φ = 0.20 before recovering the monotonically increasing trend exhibited by the Pries et al. relation. In fact, when we look at the placement of the cells in the channel, we see that the shapes adopted by the RBCs are not identical for different hematocrits. The effect of increasing the hematocrit and therefore adding new cells changes the local flow field around the cells, and thus the local Ca number, so that different morphologies are attained. In Fig. 4.2, we show the difference in morphologies observed for φ = 0.15 and φ = 0.20. For φ = 0.15, we see that the cells are more extended in the perpendicular streamwise direction, while for φ = 0.20, the cells form a more compact core, remaining closer to the channel centerline. This fact is also reflected in the plot of the CFL as we will show later (Fig. 4.9), where we observe an increase in CFL thickness between φ = 0.15 and φ = 0.20. We see the expected blunted velocity profiles due to the RBC-rich core in the central region of the channel for the suspension cases with φ = 0.15 and φ = 0.20. The Poiseuille solution (shown in red) represents the pure plasma flow under the same pressure gradient. Q s is smaller in the case of φ = 0.15 than φ = 0.20, and the relative apparent viscosities found are µ rel = 1.074 and µ rel = 1.044, respectively, indicating a decrease in flow resistance between φ = 0.15 to φ = 0.20 which is due to the differences in the morphologies of the individual RBCs between the two cases (cf: Fig. 4.2). Previously, a plateau and then decrease in apparent viscosity with increasing hematocrit was reported under linear shear for low hematocrit suspensions. [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF].In linear shear flow, this is due to the RBC organization in parallel files. However, as we have just seen in channel flow, it results from the cells being less extended in the perpendicular streamwise direction at higher hematocrit. This feature is a consequence of the confinement and the arrangement of cells.

We then investigated to see over which widths the decrease was observed. We find that it persists for widths below 12µm, but not above. For larger widths, we see the intuitive behavior; however, when we enter into the capillaries, we see this striking behavior. From Fig. 4.4, we see that this phenomenon persists within the range of channels with size comparable to capillaries. These results have important implications for transport in the capillaries. RBCs appear to organize themselves in a way to reduce as much as possible the increase of viscosity due to an increase of hematocrit within certain hematocrit ranges. The result of this is enhanced efficiency of oxygen transport. It would be interesting to see how in-vivo hematocrits at each of these levels compare to those where decreased apparent viscosity is observed in our simulations. We next move to study wavy channels. Figure 4.5a depicts the dependence of the relative apparent viscosity on hematocrit and channel size in channels with undulating walls, as predicted by our simulations.

Chapter 4. Suspensions in wavy-walled channels

The relative apparent viscosity again shows the non-linear, monotonic increase with hematocrit, reminiscent of the trend observed in the case of straight-walled channels. In order to better understand the effect of waviness on the apparent viscosity, we look at the ratio between the two. Figure 4.5b shows the ratio of the suspension flow rate between the straight and wavy cases. The results demonstrate that the presence of waviness always decreases the flow rate of the suspension, as the ratio is always greater than unity, with the effect being the greatest for channels of smaller width. In the case of the smallest channel, W = 12 µm, the effect of waviness on the flow rate of the suspension is dramatic with the flow rate in the wavy channel decreasing significantly as the hematocrit increases. For the largest hematocrit, the flow rate of the suspension in the wavy channel is nearly 50% of that in the straight channel. For W = 22 µm, the effect of waviness is always greater than 10%; however, the effect of hematocrit is less striking (only a few percent difference between the highest and lowest hematocrits). For channels of larger width, W = 32 -57µm, the effect becomes largely independent of channel width, as can perhaps be expected as the perturbation of the endothelial cell height becomes less important with respect to the width of the channel. The effect of hematocrit is present but slight.

Similar 2D systems have previously been investigated for the straight-walled case [START_REF] Fedosov | Blood Flow and Cell-Free Layer in Microvessels[END_REF][START_REF] Zhang | Effects of Erythrocyte Deformability and Aggregation on the Cell Free Layer and Apparent Viscosity of Microscopic Blood Flows[END_REF][START_REF] Bagchi | Mesoscale simulation of blood flow in small vessels[END_REF][START_REF] Sun | Particulate nature of blood determines macroscopic rheology: A 2-D lattice Boltzmann analysis[END_REF]. For a suspension of rigid RBCs in a 20 µm channel with φ = 0.305, [START_REF] Sun | Particulate nature of blood determines macroscopic rheology: A 2-D lattice Boltzmann analysis[END_REF] found a relative apparent viscosity of µ rel = 1.50. [START_REF] Bagchi | Mesoscale simulation of blood flow in small vessels[END_REF] also studied a 20µm channel with a channel hematocrit of φ = 0.33 and reported a relative apparent viscosity of µ rel = 1.60, which is much larger than our result of µ rel = 1.25. Bagchi's result is even higher than that obtained by Sun and Munn in the rigid cell limit; one possibility for the discrepancy is the different mean flow viscosities as discussed by Bagchi. Our result agrees well with that of [START_REF] Zhang | Effects of Erythrocyte Deformability and Aggregation on the Cell Free Layer and Apparent Viscosity of Microscopic Blood Flows[END_REF], who found for a W = 20 µm channel with φ = 0.322 a relative apparent viscosity of µ rel = 1.27.

Fahraeus Effect

The Fahraeus effect is the observation that in blood flow, the channel hematocrit φ (given by Eq. 4.1) is smaller than the discharge hematocrit φ d , defined as the ratio of volume flux of RBCs to the volume flux of whole blood. This effect is directly correlated with RBC migration to the centerline, thus resulting in the mean RBC velocity being higher than the mean plasma velocity. We find the ratio of the channel to discharge hematocrit by:

φ/φ d = U/U RBC (4.9)
where U is the mean velocity of the suspension (plamsa+RBCs) calculated by U = Q/A and U RBC is the mean RBC velocity calculated by averaging over all cells present in the channel and over time. [START_REF] Pries | Blood viscosity in tube flow: dependence on diameter and hematocrit[END_REF] compiled experiments on blood flow through glass tubes of different

Results

83

diameters at various φ d values to obtain an empirical relationship between φ and φ d given by:

φ φ d = φ d + (1 -φ d )(1 + 1.7e -0.35D -0.6e -0.01D ) (4.10)
where the tube diameter D, which we replace with the width W of the channels in our 2D simulations. Figure 4.6a shows the hematocrit ratio obtained in the straight-wall case calculated from our simulations and from the empirical relation given by Pries et al. (dashed lines). We see that our results follow the same general trend as the empirical relation given by Pries et al. For smaller hematocrits, RBCs are concentrated near the centerline and thus their mean speed is higher. As hematocrit increases, the cells are pushed further toward the walls, and thus the mean speed decreases. Also, for a given hematocrit, the ratio increases with the width of the channel. As the channel width for a given hematocrit increases, the number of cells increases and therefore so does the size of the RBC core. One will remark that our simulations predict higher hematocrits, and thus lower discharge hematocrits, than is predicted by the Pries empirical formula. This discrepancy is attributable to the fact that the Pries relation was derived for flow in cylindrical tubes, while our simulations depict two-dimensional flow between parallel plates. In fact, 2D simulations are expected to underpredict the CFL. The CFL cross-sectional area ratios for channel and tubular flow are 2δ/W and 4δ(Dδ)/D 2 , respectively, implying a thinner CFL in 2D simulations compared to the corresponding 3D systems (for example, δ = 2 µm and W = D = 20 µm, so that 2δ/W = 0.2, whereas 4δ(Dδ)/D 2 = 0.36.

We also performed numerical simulations in channels with wavy walls. We are again able to observe the Fahraeus effect, whereby all the ratios of the tube to discharge hematocrits fall below unity, as shown in Fig. 4.6b. We recover a similar trend for the hematocrit ratio of wavy-walled channels as a function of hematocrit and width as seen in the straight case. [START_REF] Zhang | Effects of Erythrocyte Deformability and Aggregation on the Cell Free Layer and Apparent Viscosity of Microscopic Blood Flows[END_REF] report for a W = 20 µm channel with φ = 0.322 (for whom we mentioned earlier the good agreement for µ rel ), a discharge hematocrit of φ d = 0.388, which is in good agreement with our value of φ d = 0.371 for the straight case. For this case, in the presence of wavy walls, we calculate a discharge hematocrit of φ d = 0.380, representing an increase in discarge hematocrit in the presence of undulating walls. As seen in Fig. 4.6c, with only a few exceptions, it is the case that the discharge hematocrit in the wavy case is larger than in that observed in the straight case. Furthermore, the results of the hematocrit ratio found in [START_REF] Sun | Particulate nature of blood determines macroscopic rheology: A 2-D lattice Boltzmann analysis[END_REF] for rigid particles are φ d ∼ 0.355 for the W = 20 µm channel and φ d ∼ 0.335 for the W = 40 µm channel at channel hematocrit φ = 0.305, while our results are φ d = 0.371 and φ d = 0.359, respectively. These results are coherent with the trend we would expect: because RBC migration towards the centerline should be attenuated for more rigid particles and thus would have lower discharge hematocrits. For the wavy case where W = 40 µm, we find φ d = 0.362, which is slightly larger than the discharge hematocrit found in the straight case. Since RBC deformability has been shown to increase CFL thickness, generating a denser RBC core region, we would expect that in our cases where we model deformable RBCs, that the CFL would be larger and therefore the RBC volume flux to be greater and therefore the discharge hematocrit to be than for corresponding suspensions of rigid particles.

Analytical Model

We develop an analytical model for the dependence of the RBC flux Q RBC on the channel hematocrit φ. We adopt a continuous two-fluid model for suspension flow. This approach has previously been used to explore other problems related to blood flows [START_REF] Fung | Biomechanics: Mechanical Properties of Living Tissues[END_REF], and a 3D equivalent of this analytical model was developed in [START_REF] Farutin | Optimal cell transport in straight channels and networks[END_REF]). We consider an inner core containing a homogeneous suspension of area fraction φ c surrounded by a CFL of thickness δ, flowing in a channel of width W. For a given channel hematocrit φ, we find the hematocrit of the inner core by

φ c = φ( A ch A core ) = φ W W-2δ
where A ch is the total area of the channel and A core is the area of the inner core. The RBC flux is found by

Q RBC (φ) = φ W/2-δ -W/2+δ u(φ c , y)dy (4.11)
For a given imposed pressure gradient, we can calculate the Stokes velocity profiles, u(φ c , y), in the inner core and outer cell-free layers. We assume the inner core is a homogeneous, dense suspension of concentration φ c , and its relative viscosity can be estimated from the Krieger-Dougherty law:

µ rel = (1 -φ c /φ m ) -[µ]φ m
, where φ m is the maximum fraction, [µ] is the intrinsic viscosity, and µ ext is the viscosity of the cell-free fluid defined previously above [START_REF] Krieger | A Mechanism for NonNewtonian Flow in Suspensions of Rigid Spheres[END_REF]. We assume that the thickness of the outer cell-free layer decreases linearly as hematocrit increases, setting δ = δ 0 (1φ/φ m ) where δ 0 is a constant characterizing the cellfree layer at low hematocrit. We set the intrinsic viscosity [µ] = 2, the well-known value for suspensions of discs [START_REF] Brady | The Einstein Viscosity Correction in n Dimensions[END_REF], φ m = 1, and δ 0 = 2µm. Solving, we find

Q RBC = 1 2 dP dx 2 3µ rel W 2 -δ 3 -- W 2 + δ 3 -Wδ -δ 2 W -2δ (4.12)
Figure 4.7a illustrates the prediction of this reduced model. We present our simulation results for straight channels in Fig. 4.9b. We would expect in the case of deformable particles that the The effect of waviness, shown in Fig. 4.7d, is to decrease the RBC flux, drastically for the smallest channel, W e f f = 12µm.

Cell-Free Layer

Through our simulations of suspensions in wavy-walled channels, we have observed the Fahraeus-Lindqvist effect whereby the apparent viscosity of the blood flow decreases with decreasing diameter of the tube as well as the Fahraeus effect where the discharge hematocrit exceeds the channel hematocrit. The primary explanation to describe these phenomena is the formation of a layer of reduced RBC concentration near the walls, known as the CFL. The CFL refers to the layer of blood adjacent to the confining walls which is absent of RBCs, since they are subject to migration towards the channel center. The fluid viscosity of the CFL region is much smaller than that of the core region populated with RBCs, providing an effective lubrication for the core to flow. An example depicting the edges of the CFL in the case where W e f f = 57 µm and φ = 0.40 is shown in The instantaneous edges of the CFL are shown in black. An average CFL thickness, δ, for the channel, defined as the closest cell-to-wall distance for both the top and bottom walls, can be calculated as:

δ = δtop + δbottom 2 (4.13)
where δtop and δbottom are respectively the average closest cell-to-wall distances in time and space between the RBC membranes and the top and bottom walls. The average CFL for suspensions Chapter 4. Suspensions in wavy-walled channels in straight channels as a function of hematocrit for different channel widths is reported in Fig. 4.9. For channel widths W = 22 -57 µm, we see that the CFL ratio is greatest for smaller hematocrits and decreases as the hematocrit increases. In fact, for smaller hematocrits, it appears that the thickness of the CFL as a proportion to the total width of the channel is fairly independent of channel width. However, as we increase the hematocrit, we observe greater dispersion of the CFL ratio among the different sizes of channels. The clustering at low hematocrits can be explained as follows. For low hematocrits, the cells are aligned single-file in the channel so that in a channel of width W, 2δ + 2R 0 = W so that δ W = 1 2 -R 0 W . Thus, for channels of small width, the ratio R 0 /W is important; however, as we increase the width of the channel and W >> R 0 , δ/W remains more or less constant. This is seen in the φ = 0.05 case. Since for smaller hematocrits, the cells are single-file in the center of the channel, the thickness is governed by the size 4.2. Results 89 of a single RBC; however, as we increase hematocrit, we spread the cells to the periphery of the channel and this relation is no longer valid. For the smallest channel, W = 12 µm, we see that from φ = 0.05 -0.20, the ratio increases but then decreases with hematocrit after φ = 0.20. In fact, we touched upon this point earlier in Fig. 4.2 where we showed that the morphology of the cells in this single-file configuration changes thereby affecting the elongation of the cells along the perpendicular streamwise direction and thus increasing the CFL from φ = 0.15 to φ = 0.20. In Fig. 4.10, we display the morphologies of the RBCs for W = 12 µm for hematocrits in the range φ = 0.05 -0.25.

When comparing our CFLs to previous values found by numerical studies in the literature, we find some discrepancies, which we think are mostly due to the difference in channel length used and a difference in initial RBC distribution in the channel. For example, [START_REF] Zhang | Effects of Erythrocyte Deformability and Aggregation on the Cell Free Layer and Apparent Viscosity of Microscopic Blood Flows[END_REF][START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF] simulated a suspension of RBCs in a 19.5 µm × 58.5 µm channel with 27 biconcave RBCs with normal membrane elasticity. Our simulations are done using a longer channel length so that the initial RBC distributions greatly differ. In Zhang et al.'s (2009) Therefore, while the average hematocrit in the channel is similar to that in Zhang et al., the differing RBC distributions can lead to large differences in the CFL thickness. However, our results agree well with those of [START_REF] Sharan | A two-phase model for flow of blood in narrow tubes with increased effective viscosity near the wall[END_REF]. We recast the data in Fig. 4.9 in terms of dimensional CFL thickness all widths at three different hematocrits.

Figure 4.12: Dependence of the CFL width on the channel width and hematocrit for a straight channel [START_REF] Sharan | A two-phase model for flow of blood in narrow tubes with increased effective viscosity near the wall[END_REF] similarly found an increase in dimensional CFL thickness with channel width for this range of channel widths (∼ 10 -60 µm). They also found that the thickness of the CFL decreases with increasing hematocrit.

The results above were obtained for flat-walled channels. In vivo, vascular walls constitute a wavy surface due to the undulating endothelial cells that line the inner vascular walls. Figure 4.13a shows the ratio of the width of the CFL to the effective width of each of the channels studied with undulating walls having wall waviness similar to that in vivo. The effective width W e f f is the average width over one period of wall undulation and is set such that W e f f = W. We see largely the same general trend for the wavy case as for the straight-walled case, i.e. the ratio decreases with increasing hematocrit with the exception of the 12 µm channel. To compare both cases, we plot the ratio of the wavy to straight cases in Fig. 4.13b. For the smallest hematocrit examined, φ = 0.05, the effect of wall waviness is negligible for all channel widths except the smallest (W = 12 µm). Even at larger hematocrits for all but W = 12 µm, a trend appears indecipherable, with the ratio fluctuating within ∼ 5% for most channel widths among the different hematocrits. For the W = 12 µm case, the ratio remains near unity except at the two extremes of hematocrit values. At the smallest hematocrit φ = 0.05, the CFL obtained in the wavy-walled case is larger than that for the straight-walled case, whereas the opposite occurs at the high hematocrit limit φ = 0.40. If we think of the transit of a single RBC (representative of Chapter 4. Suspensions in wavy-walled channels the low hematocrit limit) over one period in the wavy-walled channel with in-phase waviness between the top and bottom wals, the RBC, which flows largely along the channel centerline, is closer to the wall than in the flat-wall case in the peak region but is further from the wall in the valley region. However, for higher hematocrits, because of space constraints, this no longer holds true, and the cells are pushed outward upon exiting the peak region into the valley region, thus acting to decrease the CFL layer in the wavy-walled case. This exclusion effect is not acting in the straight case, and it is especially prominent for the smallest diameter at higher hematocrits because the space is particularly limited.

Results

93

Trajectories of Cells

We look at the effect of waviness on the evolution of the center of mass of the RBCs for the smallest and largest channel sizes studied and for the lowest and highest hematocrits. In the case of the small-width channel where W = 12µm and with a low hematocrit, corresponding to 4 RBCs in the flow, the cells align themselves in a single-file at the channel centerline in both the straight and wavy cases, as shown in Fig. 4.14a and b, respectively. The RBCs form two rows in the channel. In the straight case, we see that some RBCs switch rows, creating a "zipper effect" evidenced in Fig. 4.15a. In the wavy case, the cells also form two rows. However, we see no instances of row flipping. The evolution of the center of mass follows the periodicity of the boundary, changing ∼ 10% the width of the channel while traversing one period.

Increasing the width of the channel, and looking to the low hematocrit case, the effect of waviness on the trajecories of the center of mass of the cells is negligible. When we increase the hematocrit eight-fold, and look at the effect on the center of mass, we find that the effect is still negligible. We then look at the evolution of the elongation index, EI, of the RBCs, calculated by:

EI = L RBC W RBC = max(x N ) -min(x N ) max(y N ) -min(y N ) (4.14)
where L RBC is the length of the RBC in the streamwise direction and W RBC is the width of the RBC in the perpendicular streamwise direction. We find no change in the EI in the straight case; however, we see that the EI for the smallest channel studied follows the periodicity of the channel boundary such that it is the most extended in the streamwise direction (i.e., EI is the largest) at the narrowest part of the channel and is the most extended in the perpendicular streamwise direction at the widest part of the channel (i.e., EI is the smallest), as shown in Fig. 4.18b. Now increasing the hematocrit, in Fig. 4.19, we see that the in the widest part of the channel, the EI reaches a value slightly lower than that observed in the straight case, while the value observed in the narrowest part of the channel is greater than twice the value in the straight case. When we increase hematocrit, the EI does not show a discernible regular effect in the straight case; however, in the wavy case, the EI shows a behavior consistent with the periodicity of the channel walls. The maximal value of the EI is obtained in the narrowest part of the channel (i.e., greatest extension in the streamwise direction) and the minimal value in the widest part

Chapter 4. Suspensions in wavy-walled channels of the channel (i.e., greatest extension in the perpendicular streamwise direction). The values of the EI are higher than what we observe in the lower hematocrit case. As we saw in Fig. 4.15b, the RBCs form two rows in the channel such that the cells are more limited the available space in the perpendicular streamwise direction. Since the cells are more limited to extend in the perpendicular streamwise direction, the values of the elongation index are higher than what we observe in the low-hematocrit case.

In the wider channel with low hematocrit, there is a range of EI which the majority of cells fall between which appears to stay constant between both cases. However, there are more cells which exhibit EI outside of this range in the straight case than the wavy case, as shown in Fig. 4.20. There is no discernible effect for the wider channel with a high hematocrit between the straight and wavy cases, as shown in Fig. 4.21.

103

The linear shear model proposed in Eq. 4.15 thus is not accurate in describing the flow in the gap or the resulting WSS. When moving from x = 29 µm to x = 31 µm, we see that the cell-wall gap expands (i.e., CFL increases) from 2.8 µm to 4.8 µm. From the linear velocity profile proposed by Eq. 4.15 and shown in Fig. 4.24b (dashed lines), we would predict the WSS to thus decrease (as evidenced by the increasing slope of the dashed lines). However, from the local velocity gradient, (∂u/∂y) y=0 as well as the WSS shown in Fig. 4.23, we see that the WSS increases. From x = 31 µm to x = 32 µm, the CFL thickness decreases from 4.8 µm to 4.0 µm, and from the slope of the dashed lines in Fig. 4.24b, we would predict the WSS to increase. However, our simulations show that the WSS drops.

In order to examine the relationship between perturbation WSS (a measure of the fluctuations in WSS) and hematocrit and channel width, we compute the average of the WSS fluctuations from the nominal wall shear stress in the channel. The mean WSS values are not reported since they are always close to the nominal WSS, τ 0 , found by averaging the wall shear stress in the channel in the absence of cells, with relative errors < 2%. The average perturbation is then found by applying the following equation:

τ * = 1 T f 1 2L T f ∑ n=1 L ∑ x=1 |τ bottom (n, x) -τ bottom,empty (x)| + L ∑ x=1 |τ top (n, x) -τ top,empty (x)| (4.16)
where T f is the simulation time, L is the length of the channel, τ bottom (n, x) and τ top (n, x) are the WSS along the bottom and top walls, respectively, at time n and position x, τ bottom,empty (x) and τ top,empty (x) are the WSS along the bottom and top walls, respectively, for an empty channel (i.e., in the absence of RBCs) at position x (a constant given by the Poiseuille solution in the case of the straight channel, otherwise a function of streamwise channel coordinate in the wavywalled case). The average perturbation WSS represents the average change in WSS that a given point on the endothelium might experience over time as a suspension of RBCs flows over it. For straight-walled channels, the results are shown in Fig. 4.25. With the exception of the smallest channel, W = 12 µm, we see that for a given channel width, the average perturbation shear stress increases with increasing hematocrit. It has been shown in previous works as well as in Chapter 3 of this work that the local WSS variation is inversely related to the local thickness of the CFL. If we look at the trend of the smallest channel, W = 12 µm, we see that the nonmonotonicity we see between φ = 0.15 and φ = 0.20 whereby the perturbation WSS decreases, follows the increase in average CFL we observed earlier in Fig. 4.9. However, it is the larger channels with the larger CFLs which experience the most change with hematocrit. At lower hematocrits, the average perturbation WSS is only a few percent, while at the largest hematocrit, the effect can be as large as 15 -20%. This is perhaps a bit surprising, since the absolute CFL, δ, is greater for larger channel widths than for smaller channel widths; however, the nominal WSS, τ 0 , changes among channel widths, since the Reynolds number in the channel was held constant. Therefore, τ 0 becomes smaller as we increase the channel width. The results show, however, that the hematocrit can have a large effect on the nominal WSS, even in larger channels. We depict the flow in a large channel at a high hematocrit in Fig. 4.26a.

We see from Fig. 4.26b again the same features between the CFL and WSS plots we identified in Chapter 3 as well as in Fig. 4.23 (peaks corresponding to peaks, valleys to valleys). We again look at a subset of the simulation domain in order to examine the velocity within the gap in Fig. 4.27. Chapter 4. Suspensions in wavy-walled channels These results underscore the fact that the WSS variation can still be important even when the overall CFL is large and that it is inadequate to estimate WSS merely from the local CFL thickness. From x = 72 µm to x = 75 µm, we see that the CFL thickness increases. Equation 4.15 would predict a decreasing WSS; however, our simulations show the WSS increasing. In this case, the CFL varies from 4.6 µm to 8.0 µm, nearly double that of the magnitude of the CFL examined in the W = 12 µm channel. Previous studies concluded that shear stress variation is less pronounced in situations with wider CFLs; however, these studies were conducted for a single width of channel. Keeping the channel width constant and varying flow and RBC parameters, [START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF] find that a larger CFL corresponds to a smaller WSS variation from the nominal WSS. They explain this by saying that for a wider CFL, the reverse pressure difference generated over the cell length is less significant. However, in our example, where we examine two different sizes of channel, we achieve a constant flow rate between the two cases by imposing a smaller pressure drop for the larger-width channel. The dimensional CFL in the larger channel is greater than the dimensional CFL in the smaller channel; however, the perturbation WSS varies greatly, more than 50% of the nominal value. This example shows that perturbations to nominal WSS can still be important for larger channels with larger dimensional CFLs. These findings underscore the notion that flow conditions, including the imposed pressure gradient, play an important role and that the conclusion previously reached that the WSS variation becomes less significant with wider CFLs cannot be applied if the width of the channel is varied.

Results

Figure 4.28a depicts the perturbation WSS obtained for suspensions in wavy-walled channels. We see that for a given channel width, except for the smallest channel, W = 12 µm, increasing the hematocrit results in a monotonically increasing average perturbation WSS. We can see that Generally, as the average CFL thickness decreases for a given channel width with increasing hematocrit, the average perturbation WSS increases; however, we see again, from W = 12 µm φ = 0.30, the decrease in the average WSS perturbation is not predicted by the average CFL. The trend of lower WSS values generally corresponding to locations with narrower gaps between cell and wall is generally observed; however, we cannot predict the magnitude or entire dynamic of instantaneous WSS from the instantaneous CFL. We see from the plot of the ratios in the wavy vs straight cases that waviness has the greatest effect in smaller channels and at lower hematocrits. Surprisingly, for the larger channels, W = 47µm and W = 57µm, we see that waviness may even reduce the perturbation WSS observed in the straight channel.

Conclusion

Firstly, we observe a striking decrease in apparent viscosity with hematocrit in small channels, an observation previously made in linear shear flow [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]. Upon investigation, we are able to link the rheological behavoir and the spatial organization of the RBCs in the channel. Our results indicate a decrease in suspension flux in wavy channels. We observe an increase in discharge hematocrit in suspensions flowing in wavy channels. We developed a 2phase 2D analytical model to explore the effect of the RBC flux as a function of hematocrit. Our numerical results indicate that waviness decreases the RBC flux. We see large variation in the CFL for the smallest channel. For small hematocrit, we find that the CFL is larger in the wavy case compared to the straight, while for large hematocrit, the CFL is smaller in the wavy case than the straight due to spacial exclusion effects. In investigating the individual trajectories of the center of masses of the RBCs in the suspensions, we observe a robust two-file organization in wavy channels, while in straight channels, we observe cell flipping between the two rows. We also observe significant variation and periodicity in the elongation index of the cells while they pass through the channel undulations. Lastly, we find significant WSS perturbations at high hematocrits even for large channels. The effect of waviness on perturbation WSS is greatest in small channels with low to moderate hematocrits, and in larger channels, waviness may even work to reduce perturbation shear stress.

Introduction

Coupled mass transfer-fluid flow problems generate interest both from a fundamental point of view as well as for a host of practical applications. For the biomedical community, natural phenomena occurring in the bloodstream including oxygen and ATP release from RBCs generate interest into processes involving coupled fluid-interface-solute dynamics [START_REF] Gong | An immersed boundary method for mass transfer across permeable moving interfaces[END_REF][START_REF] Zhang | ATP Release by Red Blood Cells under Flow : Model and Simulations[END_REF]). In addition, the engineering of controlled release systems, such as targeted drug delivery systems whereby capsules or particles release drugs into the bloodstream, involve studies of coupled mass transfer-fluid flow systems [START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF][START_REF] Gekle | Dispersion of solute released from a sphere flowing in a microchannel[END_REF].

Chapter 5. Advection-Diffusion of Passive Scalar Released by Red Blood Cells in Flow

The underlying ingredients to conduct a numerical study of these diverse problems are similar. One would need a two-way coupling between the dynamics of a moving, deformable particle and the flow, a one-way coupling between the solute dynamics and the flow (assuming the passive scalar quantity does not affect the flow) and perhaps a coupling between the solid mechanics of the particle and the solute dynamics. The model presented in this Chapter represents a general numerical model for solving coupled membrane-fluid-solute problems and can be adapted to suit the individual needs of a variety of problems involving mass transfer from particles in a flow. Two modes of solute release are considered: one in which the solute concentration is initialized inside of the particle with a perfectly permeable membrane and the other where the solute is localized at the membrane and released as a source term.

Lattice Boltzmann for Advection-Diffusion

In order to study the release of solute species from RBCs, we will solve a coupled fluid-membranesolute problem using LBM. We will show how the advection-diffusion equation (ADE) can be solved by the LB equation in very much the same way that it can be used to solve the Navier-Stokes equations. This allows us to have a framework where we can solve coupled fluid-solute problems.

The macroscopic advection-diffusion equation (ADE) for a scalar field C with an isotropic and homogenous diffusion coefficient D and source term q reads:

∂C ∂t + ∇ • (C u) = D∇ 2 C + q (5.1)
In fact, the Navier-Stokes equations can be thought of as an ADE for the fluid momentum density vector, ρ u. Therefore, we can easily adapt LBM for advection-diffusion problems. The LB equation (LBE) for solving the advection-diffusion equation is given by:

g k ( x + c k ∆t, t + ∆t) -g k ( x, t) = Ω k ( x, t) + Q k ( x, t) (5.2)
where Q k represents source terms and Ω k the collision operator. In the familiar BGK model, or a single-time relaxation to equilibrium, we can write:

Ω k ( x, t) = - 1 τ g g k ( x, t) -g eq k ( x, t) (5.3)
where we now define τ g as the relaxation time associated with the distribution functions solving for the passive scalar dynamics. The main difference between the different formulations of LBM for solving the NSE and the ADE is that the momentum of the distribution functions in the latter case is not conserved. Instead, the speed is imposed externally (taken from the NS solver).

Lattice Boltzmann for Advection-Diffusion

113 Therefore, the quadriatic equilibrium distribution takes the form:

g eq k = w k C 1 + c k • u c 2 s + ( c k • u) 2 2c 4 s - u • u 2c 2 s (5.4)
It can be shown that the diffusion coefficient, D, can be found by:

D = c 2 s τ g - ∆t 2 (5.5)
which is analagous to the kinematic viscosity, ν, in the LBM formulation for the NSE. However, the BGK model is limited in accuracy for large τ g ; it has been shown that the BGK approximation is optimum for τ = [1, 0.5). Thus, in order to be able to simulate a larger range of diffusivities, we will employ another collision operator with two relaxation times, namely TRT. The TRT results in better accuracy for higher values of τ. As the name of the model implies, it involves two relaxation rates, ω + and ω -. The so-called magic parameter, Λ, relates the two and characterizes the truncation error and stability properties. The expression relating the two relaxation times for the TRT model is given by:

Λ = 1 ω + ∆t - 1 2 1 ω -∆t - 1 2 (5.6)
Specific values of Λ which ensure accuracy and stability under certain conditions have been found. [START_REF] Ginzburg | Optimal Stability of Advection-Diffusion Lattice Boltzmann Models with Two Relaxation Times for Positive / Negative Equilibrium[END_REF] proposed that Λ = 1/4 is optimal and works well for any Peclet number. One might simulate a large range of diffusion coefficients by changing ω -according to Eq. 5.4 and then adjust the free parameter ω + . The TRT model is based on the decomposition of populations into symmetric and anti-symmetric components, g k = g + k + g - k where g + k = (g k + g k)/2 and g - k = (g kg k)/2; k denotes the opposite direction of the kth discrete velocity c k and for the rest particle, g + 0 = g 0 and g - 0 = 0 and the same applies for the equilibrium distribution functions. The evolution of the particle distribution function is described as:

g k ( x + c k , t + 1) = g k ( x, t) -ω + (g + k -g eq+ k ) -ω -(g - k -g eq- k ) + Q k (5.7)
We have included the source term so that we can account for any creation or destruction of the solute species in the domain. The source term takes the form,

Q k = 1 - 1 2τ g w k q
(5.8)

Chapter 5. Advection-Diffusion of Passive Scalar Released by Red Blood Cells in Flow which can be found by a Champman-Enskog analysis for which there is an unwanted term which is a discrete lattice artefact. We then can find the concentration field using:

C = ∑ k g k + Q k ∆t 2 (5.9)
This is analagous to the procedure to remove the discrete lattice artefact when solving the LBE for the NSE with an external forcing term. Further details on the use of LB to simulate the ADE equation can be found in [START_REF] Krüger | The Lattice Boltzmann Method[END_REF].

Benchmark Testing of Code

In order to evaluate the performance of our ADE solver, we evalute our results against the analytical 2D solution for advection-diffusion of an initial Gaussian concentration distribution in the domain. Accordingly, the initial concentration profile is given by:

C( x, t 0 ) = C 0 exp - ( x -x 0 ) 2 2σ 2 0
(5.10) with initial position x 0 = (x 0 , y 0 ) T ∆x and width σ 0 . In the case of a uniform advection velocity, u, the evolution of the concentration field in time is given as:

C( x, t) = σ 2 0 σ 2 0 + σ 2 D C 0 exp - ( x -x 0 -ut) 2 2(σ 2 0 + σ 2 D ) , σ D = √ 2Dt
(5.11) Equation 5.11 represents the analytical solution to an initial Gaussian concentration profile in a uniform velocity field. We choose a large 2D domain so that we will avoid special boundary condition considerations with L x × L y = 500∆x × 500∆x on the familiar D2Q9 grid. We initialize the center of the Gaussian hill at x 0 = (200, 200) T ∆x with an initial width of σ 0 = 10∆x, and the initial peak concentration as C 0 = 1 (in lattice units). In the first case, we consider where the dynamics are diffusion-dominated, setting Pe = 0 by setting u = (0, 0) T ∆x/∆t and we give a non-zero diffusion coefficient such that D = 1.5∆x 2 /∆t (i.e., with the BGK collision model τ = 5∆t and with the TRT model, τ -= 5∆t). In Fig. 5.1, we show the concentration contour plot and the concentration profile at y = 200∆x at t = 200∆t. As expected, a large relaxation time results in significant errors for the BGK collision operator; the TRT outperforms in terms of accuracy whereby we choose Λ = 1/6 as the "magic parameter" that best minimizes the TRT error in the case of pure diffusion. We then test the case for Pe = 0 such that u = (0.1, 0.1) T ∆x/∆t for a small diffusion coefficient, D = 0.0043∆x 2 /∆t. In this advection-dominated case, we expect the Gaussian hill to maintain its shape while being advected. The "magic parameter" in TRT Λ = 1/12 for best results in advection-dominated cases [START_REF] Krüger | The Lattice Boltzmann Method[END_REF]. The results of this section validate the accuracy of our LBM solver for advection-diffusion computations.

A non-dimensionalization of Eq. 5.1 will yield the characteristic dimensionless number for advection-diffusion problems, the Peclet number. It is a measure of the strength of the flow to the diffusiveness of the scalar field. In our simulations, we define it as Pe = u max W 2D (5.17) Therefore, advection-dominated problems are defined by a large Peclet number, and diffusivedominated problems are defined by small Peclet numbers. In the microvasculature, with channel widths between 5-15 µm, flow speeds between 0.1-5 mm/s, and molecular diffusion coefficients between 10 -11 and 10 -10 (Popel and Johnson 2005), we would expect a range of Peclet numbers between 5-7500.

Diffusion of Internal Solute

We examine the first mode of release whereby an initial concentration is prescribed inside the RBC while everywhere else in the domain the concentration is zero. In order to accomplish this task, the internal nodes of the RBCs must be identified. We use a seed algorithm to locate the internal nodes of the RBC, the details of which can be found in [START_REF] Shen | Blood flow in microfluidic architectures: collective behavoirs of deformable particles in confined flow[END_REF]. In summary, the algorithm identifies the nodes adjacent to the coordinate points of the membrane which fall inside the membrane contour and then uses these nodes as seeds to search neighboring nodes along the x and y axes. Once an unidentified node is identified, it is chosen as the new seed and the process continues until all unidentified nodes within the membrane are found. By identifying the internal nodes of an RBC, we can study the effect of different shapes of RBCs on the diffusion of solute initialized inside the RBC. We accomplish the study of different shapes of RBCs by varying the reduced area whereby we fix the area of the RBC and vary the perimeter such that one has a reduced area of α = 0.7 and the other α ≈ 1. At the time of introduction of the solute in the domain, the solute is located within the RBC, such that

C( x, t = 0) =    C 0 , if x ∈ Ω in 0, if x ∈ Ω ex (5.18)
where Ω in and Ω ex denote the RBC internal domain and external plasma domain, respectively. The initial solute concentration for the three cases investigated, pure Poiseuille, RBC with α ≈ 1, and RBC with α = 0.7, are illustrated in Fig. 5.3. We first investigate a channel with C n = 0.51, where the effective radius of the RBC is defined as R = 2.6 µm and corresponding to a channel width of W = 10.2 µm. The residence time that the substance remains in the flow before it reaches the endothelium where it is instantaneously absorbed can be calculated from the simulations for the three different initial conditions presented in Fig. 5.3 and for different Pe numbers. The RBCs move in the channel, perturbing the flow, as illustrated in Fig. 5.4. These changes in the flowfield will mediate the differences we see in residence time and reaction position. In Fig 5 .5a, we see that increasing the Pe number to 100 leads to a large reduction in the residence time probabilities at all times and even more of a reduction is observed for Pe = 1000. Because the initial solute concentrations are initialized in the center of the channel, where the velocity is the greatest, flow has a significant effect on the solute residence time. Comparing the effect of the presence of an RBC at different Pe, we see from Fig. 5.5b that there is no significant effect of the RBC on the residence time distribution for Pe = 10; however, Figs. 5.5c and 5.5d show that for higher Pe, RBCs increase the peak of the distribution from the no RBC case. The effect of the presence of an RBC in strongly convection-dominated flows is to increase the distribution of smaller residence times. The effect on the mean residence time, t m (which is scaled by the diffusion time in the y-direction, (W/2) 2 /D), is shown in Fig 5 .6. We see that in the absence of an RBC (therefore in pure Poiseuille flow), t m is constant and independent of Pe. For smaller Pe, the effect of the presence of an RBC on the mean residence time is very small; however, for larger Pe, the mean residence time is significantly reduced, up to ∼ 13% for Pe = 1000. We then consider the total remaining concentration, C(t), in the channel after t = 0 in Fig. 5.7. The results show that stronger convection considerably slows down the clearance process; however, the presence of the RBC accelerates the clearance of the solute from the channel for larger Pe. We then turn our attention to the reaction position on the channel walls, shown in Fig. 5.8a. We see that lower Pe exhibit greater peak heights than higher Pe, suggesting that the concentration is more concentrated at the wall in certain places than in the case of weak convection. The effect of higher Pe is to decrease concentration "hot spots" on the wall and spread the solute in the direction of flow. Interestingly, we see for higher Pe, while the peak heights are significantly smaller than those in the low Pe cases, the tail ends of the distributions reach higher values than in the case of lower Pe. Intuitively, we understand this as such: in the case of lower Pe, the solute is absorbed closer to the point of release, so less of it is available to be absorbed downstream; whereas, for higher Pe, the solute is quickly swept downstream and it finally arrives at the walls further downstream. These results underscore the importance of downstream sites in the case of strong convection. The effect of RBC is negligible at lower Pe, as illustrated in Fig. 5.8b, but for higher Pe, the presence of the RBC increases the peak height for Pe = 100 and Pe = 1000 (cf. Fig. 5.8c and d), meaning that reactions with the wall are stronger closer to the point of release when the RBC is present. The mean reaction positions from the distributions shown in Fig. 5.8 are displayed in Fig. 5.9a. We see that the mean reaction position, x m , decreases for higher Pe from the no RBC case, meaning that upon reaching the wall, the solute concentration is larger than in the no RBC case closer to the point of release x c . However, the RBC does not have a huge effect on the standard deviation, displayed in Fig. 5.9b. Reflective of the greater tail distributions we saw in Figs. 8.8c and d, the dispersion is greatly increased after Pe = 10, with and without RBC. We then investigate the influence of the size of the channel on the results by decreasing the width of the channel. This smaller channel has a confinement ratio of C n = 0.70. We first begin by again examining the effect on the residence time probability distribution for a more confined channel. We see that for this larger confinement, the residence time distributions correspond to larger values; however, the qualitative behavior is much the same, as seen in Fig 5 .10a. However, we see a striking difference in the distributions obtained in the Pe = 1000 case for the different RBC morphologies . Both RBCs increase the peak height; however, α ≈ 1 shows a greater increase in peak height in the residence time distribution. The α ≈ 1 RBC clears solute from the channel much faster than the α = 0.7 RBC. The plot of the mean residence time as a function of Pe shows that the RBC corresponding to α ≈ 1 reduces mean residence time ∼ 8% from the α = 0.7 case and nearly 23% from the no RBC case. We then consider the total remaining concentration in the channel, as shown in Fig. 5.7. It largely follows the same behavior as the smaller confinement case, although the clearance process is slightly more enhanced for each case. However, we do observe a much greater effect of RBC on the total remaining concentration for larger Pe. The effect of the RBC is to accelerate clearance of the solute from the channel, with RBC α ≈ 1 having the greatest effect. Finally we look at the reaction position distribution functions. In the Pe = 100 case, we see that α ≈ 1 results in a slightly larger peak height than the pure Poiseuille case, while the RBC case α = 0.7 actually shows a small reduction in peak height from the Poiseuille case. These results suggest that α = 0.7 carries more of the solute downstream to react, while α ≈ 1 creates more of a "hot spot" on the wall. We turn to the mean reaction position and the degree of dispersion, shown in Fig 5 .15. We find that the mean reaction position for α = 0.7 occurs much earlier than in the pure Poiseuille and α ≈ 1 cases when Pe = 1000. The distribution for α = 0.7 appears "flatter" than the no RBC case and the α ≈ 1 cases. This could be due to the fact that the RBC is more extended in the streamwise direction, therefore spreading the solute over a larger area to the wall.

Diffusion of Internal Solute

Solute Release at Membrane

Here we consider that the solute is released for a given amount of time at a constant rate from the membrane. In this case, the release of the species from the membrane acts as a source term in Eq. 5.7 which we solve over the entire domain. We use the discrete delta function to interpolate the solute concentration from the membrane nodes to the Eulerian nodes. This mode of release can be used to probe the effect of release position on the residence time and reaction positions.

We investigate for angles varying between the front of the RBC to the tail, θ 0 = 0 and θ 0 = π, respectively. Figure 5.16 shows the position of the angles of release, θ r , for the two different RBC morphologies studied. We see for Pe = 0, α ≈ 1 is symmetric around θ r = π/2 with the smallest t m occurring at the north pole, θ r = π/2, corresponding to the closest release position to the wall. As Pe is increased, we see that the mean residence times become slightly longer, but are still symmetric with respect to the north pole. For α = 0.7, we see the smallest mean residence time occur at θ r = 2π/3, corresponding to the closest release position to the channel wall.

The influence of the release position on the reaction position on the wall show similar results. Interestingly, in the Pe = 1000 case, α = 0.7 shows increasing dispersion for small release angles towards the front of the RBC, suggesting that the solute is carried away downstream faster, making contact with the wall only later.

Both models we have described are extremely simple and there remains much work to be done going forward. In the case of solute release from RBCs, there should be a boundary condition on the membrane. We have considered only perfectly permeable membranes in these two models. This is especially important to caution in the interpretation of the results from the second mode of release. Because the RBC is permeable, we expect some of the solute to be released into the RBC itself. In modelling RBCs, we should instead seek a selectively permeable membrane or perfectly impermeable membrane (i.e., perfectly reflecting), and we expect very different results to those of a perfectly permeable membrane presented here.

Conclusion

In this chapter, we developed a 2D advection-diffusion solver for moving, deformable particles within the LBM framework. We explored two modes of solute release from the particle: one where a solute concentration is placed inside particles and the other where the solute is released at the membranes as a source term. Related to the first mode of release, we find several interesting conclusions. The presence of RBCs in strongly convection-dominated flows acts to accelerate the clearance of the solute species from the channel, resulting in smaller residence times. We see a marked reduction in the mean residence time for strongly convected flows. For a channel with strong confinement, Cn = 0.70, α ≈ 1 clears solute faster than α = 0.7 (i.e., smaller mean residence times observed for α ≈ 1). We also find stronger "hot spots" when the RBC is present, except in the case of Cn = 0.7, Pe = 100, where we find a decrease in peak activity; otherwise, results suggest that the presence of RBCs enhances stronger reactions in the channel. We also find that the residence time and dilution upon wall contact of the solute concentration released on the membrane are strongly moderated by the exact morphology of the RBC. We observe that the distance of the RBC-wall gap at the point of release corresponds to both the smallest mean residence times and the smallest dispersion values. In future works, we plan to implement more complex boundary conditions on the surface of the RBC to simulate specific biophysical mass transfer phenomenon such as oxygen transfer or ATP release, and we expect these boundary conditions to dramatically affect the dynamics of the system.

Chapter 6

Summary and Perspectives Chapter 3: RBC morphology and endothelium waviness

We see with our simulations of pathological RBCs the significant effect of rigidity and morphological alterations resulting from disease in exerting perturbation shear stress on the wall. When investigating the extremal values of the wall shear stress, we find that the maximum wall shear stress shows more variation over the parameter space than the minimum wall shear stress. Within the microvasculature, we find that the inverse minimum distance of the RBC-wall gap normalized to the channel width is predictive of the maximum wall shear stress. This relationship could be of value when considering the challenges of experimental procedures to measure the wall shear stress. Perhaps relationships could be found based upon distance between the RBC and wall. With our 2D model, we have confirmed the mechanism behind the variation seen in the wall shear stress originally proposed in the analytical solution by [START_REF] Yin | Cell-free layer and wall shear stress variation in microvessels[END_REF]. In addition, we consider the waviness of the endothelium. We find that the time-averaged behavior of extremal wall shear stress is not predicted by the contiuum description of blood and that subcellular gradients in wall shear stress are significant. Simulations in 3D would be an obvious next step for this work. With more realistic 3D simulations of blood flow in tubes, we would also be able to consider wall-normal and circumferential components of stress on the vessel wall in addition to the most often considered streamwise component. That being said, there is plenty more complexity to add within 3D models. The glycocalyx, the fuzzy layer of proteins covering the luminal size of endothelial surface, will itself affect the flowfield and the transmission of forces to the endothelium. Also in this chapter, we investigate how deformability mediates an RBC's trajectory in an undulated channel, which could be potentially useful in the design of microfluidic devices. We find that with large enough amplitude, recirculation zones can appear between the valleys of the undulated walls and these can be used to trap cells.

Chapter 4: Suspensions

In studying suspensions of RBCs, we find that for larger widths, we observe the monotonic increase of apparent viscosity with hematocrit as we expect; however, when we enter into the capillaries, we observe a different behavior. We observe a decrease in apparent viscosity with hematocrit in small channels, an observation previously made in linear shear flow [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]. Upon investigation, we are able to link the rheological behavoir to the spatial organization of the RBCs in the channel: the spatial organization of the cells affects the rheology of the suspension and leads to a decrease of apparent viscosity with increasing hematocrit. These results have important consequences for the efficiency of blood flow in the capillaries. As the rheological properties are dependent on the spatial organization of the cells, investing suspensions of pathological RBCs would be of great interest, as would spatial organisations of suspensions in 3D. We find wavy channels change the rheological and dynamic behavior of suspensions of RBCs. Our results indicate a decrease in suspension flux in wavy channels. We observe an increase in discharge hematocrit in suspensions flowing in wavy channels versus straight. We developed a 2-phase 2D analytical model to explore the effect of the RBC flux as a function of hematocrit. Our numerical results indicate that waviness decreases the RBC flux. We see large variation in the cell-free layer for the smallest channel. For small hematocrit, we find that the cell-free layer is larger in the wavy case compared to the straight, while for large hematocrit, the cell-free layer is smaller in the wavy case than the straight due to spatial exclusion effects. In investigating the individual trajectories of the center of masses of the RBCs in the suspensions, we observe a robust two-file organization in wavy channels, while in straight channels, we observe "cell flipping" between rows. We also observe significant variation and periodicity in the elongation index of the cells while they pass through the wavy channel undulations. Lastly, we find significant WSS perturbations at high hematocrits even for large channels. The effect of waviness on perturbation WSS is greatest in small channels with low to moderate hematocrits. In larger channels, our results suggest waviness could even work to decrease the wall shear stress compared to straight channels.

Chapter 5: Solute Release

We developed a 2D advection-diffusion solver for moving, deformable particles within the LBM framework. We used our model to explore two modes of release from the particles. In the first mode of release, the concentration is initialized inside the particle and diffuses into the domain through the concentration gradient, and in the second, the concentration is located on the membrane and is released to the domain as a source term. We use our model to perform simulations in the microcirculation, investigating the effect of RBC morphology on the solute dynamics. Related to the first mode of release, we find several interesting conclusions. In strongly convection-dominated flows, RBCs act to accelerate clearance of the solute species from the channel, resulting in smaller residence times. A marked reduction in the mean residence time is observed for strongly convection-dominated flows. We also generally find stronger "hot spots" of the concentration along the wall when the RBC is present. Morphological differences in the RBCs can lead to important differences in residence time and dilution upon contact with the wall. In the second mode of release, we find that the residence time and dilution upon contact with the wall are strongly moderated by the exact morphology of the RBC. Going forward, a study of solute release with suspensions of RBCs would provide more insight into the mechanisms at play in a more physiological condition. We recognize that the boundary conditions on the RBC will greatly affect the results (i.e., semi-permeable, impermeable, reacting). Also of interest to the field would be the implementation of boundary conditions at the endothelium to simulate the reaction of the bioactive substance at the endothelium. ATP for instance reacts at the surface of the endothelium following Michaelis-Menten kinetics while releasing ATP as a source term dependent on shear stress, both of which could be built into our model in the future.

Chapter 7

Résumé Chapitre 3: Morphologie des globules rouges et ondulation de l'endothélium

Nous voyons avec nos simulations de globules rouges pathologiques l'effet significatif de la rigidité et des altérations morphologiques résultant d'une maladie lorsqu'elles exercent une contrainte de cisaillement perturbée sur le mur. Lors de l'étude des valeurs extrêmes de la contrainte de cisaillement au mur, nous trouvons que la contrainte de cisaillement maximale sur la paroi présente plus de variation sur l'espace des paramètres que la contrainte de cisaillement minimale sur la paroi. Au sein de la microvascularisation, nous constatons que la distance minimum inverse de l'espace entre les parois de RBC et normalisé à la largeur du canal est prédictive de la contrainte de cisaillement maximale de la paroi. Cette relation pourrait être utile lorsque l'on examinera les difficultés que posent les procédures expérimentales pour mesurer la contrainte de cisaillement des murs. Des relations pourraient peut-être être établies en fonction de la distance entre le RBC et le mur. Avec notre modèle 2D, nous avons confirmé le mécanisme à l'origine de la variation observée dans la contrainte de cisaillement du mur proposée à l'origine dans la solution analytique de Yin et al. (2012). De plus, nous considérons l'ondulation de l'endothélium. Nous trouvons que la description continue du sang ne prédit pas le comportement moyenné dans le temps de la contrainte de cisaillement de la paroi extrême et que les gradients subcellulaires de la contrainte de cisaillement de la paroi sont significatifs. Les simulations en 3D seraient une étape évidente de ce travail. Avec des simulations 3D plus réalistes du débit sanguin dans les tubes, nous serions également en mesure de prendre en compte les composantes normales et circonférentielles de la paroi de la paroi vasculaire, en plus du composant le plus souvent considéré au niveau du flux. Cela étant dit, il y a beaucoup plus de complexité à ajouter dans les modèles 3D. Le glycocalyx, la couche floue de protéines recouvrant la surface endothéliale, affectera lui-même le champ d'écoulement et la transmission des forces à l'endothélium. Également dans ce chapitre, nous examinons comment la déformabilité influence la trajectoire d'une RBC dans un canal ondulé, ce qui pourrait être potentiellement utile pour la conception de dispositifs microfluidiques. Nous trouvons qu'avec une amplitude suffisante, des zones de recirculation peuvent apparaître entre les creux des parois ondulées et peuvent être utilisées pour piéger les cellules.

Chapitre 4: Suspensions

En étudiant les suspensions de globules rouges, nous trouvons que, pour les plus grandes largeurs, nous observons l'augmentation monotone de la viscosité apparente avec l'hématocrite, comme nous l'espérons; cependant, lorsque nous entrons dans les capillaires, nous observons un comportement différent. Nous observons une diminution de la viscosité apparente avec l'hématocrite dans les petits canaux, une observation faite précédemment dans le flux de cisaillement linéaire [START_REF] Thiébaud | Prediction of anomalous blood viscosity in confined shear flow[END_REF]. Après enquête, nous pouvons relier le comportement rhéologique à l'organisation spatiale des globules rouges dans le canal: l'organisation spatiale des cellules affecte la rhéologie de la suspension et conduit à une diminution de la viscosité apparente avec l'augmentation de l'hématocrite. Ces résultats ont des conséquences importantes sur l'efficacité du flux sanguin dans les capillaires. Les propriétés rhéologiques dépendant de l'organisation spatiale des cellules, il serait très intéressant d'investir dans des suspensions de globules rouges pathologiques, tout comme dans l'organisation en 3D de suspensions spatiales. Nous trouvons des canaux ondulés changent le comportement rhéologique et dynamique des suspensions de globules rouges. Nos résultats indiquent une diminution du flux de suspension dans les canaux ondulés. Nous observons une augmentation de l'hématocrite de décharge dans les suspensions s'écoulant dans les canaux ondulés par rapport aux droites. Nous avons développé un modèle analytique 2D à 2 phases pour explorer l'effet du flux de RBC en fonction de l'hématocrite. Nos résultats numériques indiquent que l'ondulation diminue le flux de RBC. Nous constatons de grandes variations dans la couche sans cellules pour le plus petit canal. Pour le petit hématocrite, nous constatons que la couche acellulaire est plus grande dans le cas ondulé que dans le segment droit, alors que pour un grand hématocrite, la couche acellulaire est plus petite dans le cas ondulé que dans le segment droit en raison des effets d'exclusion spatiale. En examinant les trajectoires individuelles du centre des masses des RBC dans les suspensions, nous observons une organisation robuste à deux fichiers dans les canaux ondulés, tandis que dans les canaux droits, nous observons un "basculement de cellule" entre les lignes. Nous observons également une variation et une périodicité significatives de l'indice d'élongation des cellules lorsqu'elles traversent les ondulations du canal ondulé. Enfin, nous trouvons des perturbations WSS significatives à des hématocrites élevés, même pour les grands canaux. L'effet de l'ondulation sur la perturbation WSS est plus important dans les petits canaux présentant des hématocrites faibles à modérés. Dans les grands canaux, nos résultats suggèrent que l'ondulation peut même contribuer à réduire la contrainte de cisaillement du mur par rapport aux canaux droits.

Chapitre 5: Dynamique des solutés

Nous avons développé un solveur 2D d'advection-diffusion pour les particules en mouvement et déformables dans le cadre du LBM. Nous avons utilisé notre modèle pour explorer deux modes de libération des particules. Dans le premier mode de libération, la concentration est initialisée à l'intérieur de la particule et diffuse dans le domaine à travers le gradient de concentration; dans le second, la concentration est localisée sur la membrane et est relâchée dans le domaine en tant que terme source. Nous utilisons notre modèle pour effectuer des simulations dans la microcirculation, en étudiant l'effet de la morphologie des GR sur la dynamique des solutés. En relation avec le premier mode de libération, nous trouvons plusieurs conclusions intéressantes. Dans les flux fortement dominés par convection, les globules rouges agissent pour accélérer la clairance des espèces de solutés du canal, ce qui réduit les temps de séjour. Une réduction marquée du temps de séjour moyen est observée pour les écoulements fortement dominés par la convection. Nous trouvons aussi généralement des "points chauds" plus forts de la concentration le long du mur lorsque le RBC est présent. Les différences morphologiques dans les globules rouges peuvent entraîner des différences importantes dans le temps de résidence et la dilution au contact du mur. Dans le deuxième mode de libération, nous trouvons que le temps de résidence et la dilution au contact du mur sont fortement modérés par la morphologie exacte du RBC. À l'avenir, une étude sur la libération de soluté avec des suspensions de globules rouges permettrait de mieux comprendre les mécanismes en jeu dans un état plus physiologique. Nous reconnaissons que les conditions aux limites sur le RBC auront une grande incidence sur les résultats (c'est-à-dire, semi-perméables, imperméables, réagissant). La mise en oeuvre de conditions aux limites au niveau de l'endothélium afin de simuler la réaction de la substance bioactive au niveau de l'endothélium serait également intéressante pour le domaine. L'ATP, par exemple, réagit à la surface de l'endothélium selon la cinétique de Michaelis-Menten en libérant l'ATP en tant que terme source dépendant de la contrainte de cisaillement, ces deux éléments pouvant être intégrés dans notre modèle à l'avenir.

Titre: Etude numérique de la microcirculation sanguine et de ses interactions avec l'endothélium Mots clés: écoulement sanguin, globules rouges, endothelium, méthode Boltzmann sur réseau, rhéologie, advection-diffusion Résumé: Cette thèse porte sur l'étude des interactions entre les globules rouges (GRs) et l'endothelium, la couche des cellules qui delimite les vaisseaux sanguins. Il a été demontré que l'endothelium et les GRs jouent des rôles actifs dans divers processus du système vasculaire, et leurs interactions produisent un signal bio chimique grace à des moyens à la fois chimiques (molécules de signalisation) et mécaniques (taux de cisaillement sur les parois). D'abord, nous étudions le rôle des GRs, y compris dans des conditions pathologiques, dans la creation de contraintes de cisaillement spatiallement et temporellement dynamiques sur l'endothelium. Il a ete montré que les contraintes de cisaillement constituaient un element critique dans le déclenchement d'un signal bio mécanique dépuis l'endothelium. Par ailleurs, étant donné qu'il a été montré que les parois des vaisseaux sanguins ondulent en raison des cellules endothéliales individuelles qui le composent, nous avons integré à notre modelisation cette géometrie. On trouve que cette undulation affecte la dynamique des GRs dans l'écoulement ainsi que le taux de cisaillement sur les parois. Nous étudions rapidement dans quelle mésure la déformabilité d'un GR affecte sa trajectoire dans un vaisseau ondulé.

Pour céla, nous nous inspirons du processus de fonctionnement un appareil de deplacement lateral deterministe (DLD) qui utilise les variations de trajectoires des particules en fonction de leur taille pour les séparer dans l'écoulement. Nous etudions par ailleurs l'éffet des suspensions de GRs sur les caracteristiques rhéologiques et les contraintes de cisaillement sur la paroi du vaisseau. Finalement, nous nous adressons à les interaction chimiques en développons un modèle numérique avec la méthode de Boltzmann sur réseaux-limite immergée (LB-IBM) pour résoudre la diffusion et l'advection d'un soluté libéré par un particule en mouvement et déformable. L'oxygène et l'adénosine triphosphate (ATP) sont toutes les deux libérées par les GRs, se diffusent dans lécoulement, et sont absorbées par l'endothelium. Ils représentent des facteurs de signalisation critiques pour les processus de l'inflammation et vasodilation. Nous montrons que la morphologie des GRs affectera le temps de résidence et la dilution des espèces chimiques lorsqu'elles rentreront en contact avec la paroi du vaisseau. Ensemble, ces éléments nous conduisent vers la développement d'un modèle capable de simuler des processus vitaux du système vasculaire qui résulent d'événements locaux de composants individuels.

Title: Numerical study of blood microcirculation and its interactions with the endothelium Keywords: blood flow, red blood cells, endothelium, Lattice Boltzmann method, rheology, advection-diffusion Abstract: This thesis is devoted to the study of the interactions between red blood cells (RBCs) and the endothelium, the monolayer of cells lining blood vessels. The endothelium and RBCs have been shown to be active participants in various processes in the vascular system, and their interactions trigger biochemical signalling by mechanical (wall shear stress) and chemical (signalling molecules) means. We first investigate the role of RBCs, including pathological conditions, in creating time-and space-varying shear stress on the endothelium. Shear stress has been shown to be a critical element in biochemical signalling from the endothelium. In addition, as it has been shown that the endothelium is undulating due to the individual endothelial cells comprising it, we take this into account in our model of the geometry of the vessel wall. We find that this undulation affects the dynamics of the RBCs in the flow and the wall shear stress. We briefly explore how the deformability of a single RBC affects its trajectory in undulating channels, inspired by the idea behind deterministic lateral displacement devices (DLDs) which exploit the differing trajectories of particles based upon their sizes to separate them in flow. We also investigate the effect of suspensions of RBCs in undulating channels on rheological properties and wall shear stress. Finally, we address the chemical interactions by building a numerical model with the lattice Boltzmann-immersed boundary method (LB-IBM) to solve advection-diffusion of solute released from moving, deformable particles. Oxygen and adenosine triphosphate (ATP) are both released by RBCs and are advected and diffused in the flow and uptaken by the endothlium and serve as critical signalling factors in inflammation and vasodilation. We find that the morphology of RBCs will affect the residence time and dilution of the chemical species upon contact with the wall. Together, these elements lead us towards the development of a model capable of simulating vital processes in the vascular system which result from local interactions of individual components.

Université Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

v

  List of Figures1.1 Schematic of the ramified cardiovascular system showing artery branching into arterioles and capillaries followed by venules and veins which bring the blood back to the heart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Scanning electron microscope showing the main components of blood: RBCs, WBCs, and platelets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 (a) Biconcave RBC morphology; (b) Staining electron microscopy image of spectrin network; (c) Schematic of RBC membrane and cytoskeleton network . . . . . . 1.4 Schematic showing the development of a spherocytotic RBC due to the deficiencies in the structural connections in the RBC. The deficiency in these structural connections results in a spheric morphology with reduced surface area to volume ratio. (Perrotta et al. 2008) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 A phase diagram taken from Kaoui et al. (2011) showing five distinct regions depending on flow (represented by the dimensionless number measuring the flow strength to the cell deformability) and the degree of confinement (represented by the vesicle diameter over the channel width) . . . . . . . . . . . . . . . . . . . . . . 1.6 ATP release pathways from a RBC. Pannexin 1 (Px 1) and Piezo1 can be activated by local shear stress; CFTR is activated by freed actin and upregulates Px1 (Zhang et al. 2018) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Consolidated experimental data of flowing RBC suspensions (hematocrit of 0.45)in glass tubes illustrating the Fahraeus-Lindqvist effect whereby the relative apparent viscosity of a suspension of RBCs decreases with decreasing tube diameter until the diameter of the tube is comparable to the size of an individual RBC[START_REF] Pries | Blood viscosity in tube flow: dependence on diameter and hematocrit[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Graphical illustration of streaming step in LBM (Bao and Meskas 2014). . . . . . . 2.2 The two-dimensional, nine-velocity (D2Q9) lattice model. There are 9 discrete velocities used to stream the distribution functions at each time step. c 0 is stationary; c 1-4 point in the east, north, west, and south directions, respectively; c 5-8 point in the north-east, north-west, south-west, and south-east directions, respectively. . . . vi 2.3 Schematic of the bounce-back boundary condition used in order to recover the no-slip boundary condition. The wall is located halfway between the solid nodes (depicted as open circles) and the adjacent fluid nodes (filled circles). The incoming post-collisional populations f 3 , f 6 , f 7 are reflected back equal and opposite into the fluid domain. (Bao and Meskas 2014). . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The discrete delta functions for the 3-point and the 4-point schemes (Shen 2016) . . 2.5 Depiction of velocity interpolation on membrane node X n (t) (left) and force spreading to fluid node x (right). The hollow points denote the Eulerian fluid nodes, while the red dots depict the membrane nodes. The velocity is interpolated on node n from the velocity of the fluid nodes falling in the shaded region, while the force at x is determined from the forces of the membrane nodes falling in the shaded region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 The membrane of the RBC is modelled as a collection of springs. (Tsubota et al. 2006) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Equilibrium shapes obtained for various reduced areas where ν represents the value of reduced area of the RBC (Shen 2016) . . . . . . . . . . . . . . . . . . . . . . 2.8 Schematic of conditions for planar Poiseuille flow . . . . . . . . . . . . . . . . . . .2.9 Parabolic velocity profile of planar Poiseuille flow showing close agreement of analytical and computed solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 (a) Unconfined slipper morphology (Ca = 10; Cn = 0.2); (b) Parachute morphology (Ca = 10; Cn = 0.45); (c) Confined slipper morphology (Ca = 10; Cn = 0.7) . . 3.1 AFM images and cell profiles of unsheared human endothelial cells (HUVECs) taken from Oberleithner, 2004. The image insets show lines corresponding to the height profile displayed at the bottom of each image. The study looked at the effect of the hormone aldosterone on HUVECs. . . . . . . . . . . . . . . . . . . . . . 3.2 1D projection of regular lattice and curved boundary illustrating the Bouzidi interpolation bounceback scheme for curved boundaries . . . . . . . . . . . . . . . . 3.3 An example of a domain with wavy boundaries . . . . . . . . . . . . . . . . . . . . 3.4 A schematic of the domain showing part of the curved boundary in green. Nodes falling underneath the boundary are defined as solid nodes and are shown in green, and fluid nodes are shown in blue. The fluid node of interest is shown with a blue star. The line associated to the link 7 from this fluid node is shown as a solid red line, which passes through the bordering solid node denoted by a green star. The coordinates of the point of intersection must be known in order to calculate the distance necessary for applying the Bouzidi boundary condition. . . . . . . . . 3.5 Bilinear interpolation. The value at the black spot is the sum of the value at each colored node multiplied by the area of the rectangle of the same color, divided by the total area of all four rectangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii 3.6 (a) Streamlines and velocity flow field corresponding to induced flow by a single RBC with a parachute morphology (C n = 0.45, Ca = 10, α = 0.7). Colorbar corresponds to relative magnitude of streamwise velocity to maximal velocity, U, in the channel; (b) Velocity vectors in the co-moving frame with the RBC showing zero velocity inside the RBC. Colorbar corresponds to relative magnitude of streamwise velocity in the co-moving frame to maximal velocity. . . . . . . . . . . 44 3.7 Phase diagram showing different morphologies adopted as a function of the confinement ratio, C n , and the flow strength, C k W/R 0 , where C k is the capillary number, Ca, defined previously where the shear rate is defined as γ = R 0 U/(W/2) 2 . Figure taken from Aouane thesis (2015). . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.8 (a) Large channel (small confinement) and cell exhibiting unconfined slipper morphology and resulting wall shear stress on the top and bottom channel walls (C n = 0.2, CaW/R 0 = 15); (b) Increased confinement and cell exhibiting parachute morphology. Since the parachute is an axisymmetric shape, the wall shear stress curves for the top and bottom walls lie on top of each other (C n = 0.45, CaW/R 0 = 15); (c) Cell very confined and exhibiting confined slipper morphology (C n = 0.7, CaW/R 0 = 15). Note the different scales for shear stress for optimal visibility of the effect. x cm denotes the x-position of the center of mass of the RBC. . . . 48 3.9 Pressure gradient along the top and bottom walls of the channel. (a) Large channel (small confinement) with RBC exhibiting the unconfined slipper morphology (C n = 0.2, CaW/R 0 = 15); (b) Intermediate confinement with RBC exhibiting the parachute morphology. Since the parachute is an axisymmetric shape, the curves for the top and bottom walls lie on top of each other (C n = 0.45, CaW/R 0 = 15); (c) Small channel (large confinement) with RBC exhibiting the confined slipper morphology (C n = 0.7, CaW/R 0 = 15) . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.10 Effect of reduced area on RBC morphology and wall shear stress for C n = 0.42 (W = 12 µm) and Ca = 0.1. (a) Stationary shapes obtained for different values of reduced area; (b) Maximum, τ max , and minimum, τ min , shear stress experienced in the channel with a single RBC normalized by the Poiseuille solution. The values for both the top and bottom of the channel are shown; (c) Plot of the normalized inverse minimum distance d -1 min W of the cell to the respective closest wall . . . . . 51 3.11 Effect of reduced area on RBC morphology and wall shear stress for C n = 0.42 (W = 12 µm) and Ca = 10. (a) Stationary shapes obtained for different values of reduced area; (b) Maximum and minimum shear stress experienced in the channel with a single RBC normalized by the Poiseuille solution. The values for both the top and bottom of the channel are shown; (c) Plot of the inverse minimum distance d -1 min W of the cell to the respective closest wall . . . . . . . . . . . . . . . . . . . . . . . . . . 53 viii 3.12 Effect of reduced area on RBC morphology and wall shear stress for C n = 0.72 (W = 7 µm) and Ca = 10. (a) Stationary shapes obtained for different values of reduced area; (b) Maximum and minimum shear stress experienced in channel with one RBC normalized by the Poiseuille solution. The values for both the top and bottom of the channel are shown; (c) Plot of the inverse of the normalized minimum distance, d -1 min W , of the cell to the respective closest wall . . . . . . . . . . . . . . . 3.13 Recasting results from Figs. 3.10c, 3.11c, and 3.12c such that τ max /τ Pois vs d -1 min W . 3.14 Normalized wall shear stress along a channel with a boundary given by a sinusoidal function where λ = 30 µm and A = 2 µm . . . . . . . . . . . . . . . . . . . . . 3.15 (a) One period of the wall geometry; (b) Average, maximum, and minimum WSS averaged in time a for channel containing a flowing RBC with α = 0.6; (c) RBC with α ≈ 1 and Ca = 0.1 and Ca = 10. The corresponding morphologies as the RBC transverses the undulating wall are shown in Fig. 10a and 11a and Fig. 10e. . 3.16 (a) Maximum and minimum wall shear stress along the top and bottom walls of a wavy-walled channel as a function of reduced area with C n = 0.42 (W e f f = 12 µm ) and Ca = 0.1 and Ca = 10; (b) The inverse of the normalized minimum distance as a function of reduced area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Six snapshots of the morphological changes of a single RBC as it passes over one wall undulation for different values of reduced area (a) α = 0.6; (b) α = 0.7; (c) α = 0.8; (d) α = 0.9; (e) α ≈ 1 with C n = 0.42 (W e f f = 12 µm) and Ca = 0.1 . . . . . 3.18 Snapshots at 6 instants in time showing the morphological changes of a single RBC as it passes over one wall undulation for different values of reduced area (a) α = 0.6; (b) α = 0.7; (c) α = 0.8; (d) α = 0.9; (e) α ≈ 1 with C n = 0.42 (W e f f = 12 µm) and Ca = 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Five snapshots of the morphology of a single RBC as it passes over one wall undulation for different dephasing angles between the top and bottom walls (a) n = 0; (b) n = 1/5; (c) n = 2/5; (d) n = 3/5 . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Maximum wall shear stress normalized to the Poiseuille solution (for a flat-walled channel having the same effective width) along the channel walls as a function of the degree of dephasing, n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 (a) Channel where amplitude is divided in half (i.e., A/W = 0.2). We no longer see the recirculation zones for Re = [1, 10]; (b) Re = 0.5; (c) Re = 5; (d) Re = 50; (e) Re = 125; (f) Height of the center of the recirculation zone, h, normalized to the channel amplitude, A, as a function of Re . . . . . . . . . . . . . . . . . . . . . . 3.22 Trajectories of the center of mass of a "rigid" RBC (Ca = 0.25) initially placed at (a) y0 = 25, y0 = 30, y0 = 50, y0 = 100 in the valley; (b) y0 = 100, y0 = 145; (c) y0 = 150 with the RBC positioned with a 0 • , 45 • , or 90 • angle to the horizontal; (d) Snapshot of trapped RBC y0 = 25 . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 3.23 Trajectories of the center of mass of a "soft" RBC (Ca = 5) initially placed at (a) y0 = 25, y0 = 30, y0 = 50, y0 = 100; (b) y0 = 100, y0 = 145; (c) y0 = 150 with the RBC positioned with a 0 • , 45 • , or 90 • angle to the horizontal . . . . . . . . . . . . . 3.24 (a) Morphology exhibited in valley region; (b) Morphology exhibited in peak region 4.1 Relative apparent viscosity as a function of channel width , W, and channel hematocrit, φ, for a channel with straight walls. Dashed lines represent the empirical relation of Pries et al. (1992) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Top: Placement of cells in the channel for a hematocrit of φ = 0.15; Bottom: Placement of cells in the channel for a larger hematocrit of φ = 0.20. The RBCs are less extended in the perpendicular streamwise direction for φ = 0.20. . . . . . . . . . . 4.3 Streamwise averaged velocity profiles . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Relative apparent viscosity as a function of hematocrit in capillaries. . . . . . . . . 4.5 (a) Relative apparent viscosity as a function of channel width and hematocrit for channels with wavy walls; (b) Ratio of volume flux for suspensions in straight vs. wavy channels as a function of hematocrit and channel width . . . . . . . . . . . . 4.6 (a) Ratio of channel hematocrit, φ, to discharge hematocrit, φ d , versus channel hematocrit for the straight case; (b) Ratio of tube hematocrit, φ, to discharge hematocrit, φ d , versus tube hematocrit for the wavy-walled case. Dashed lines in (a) and (b) represent empirical relation given by Pries at al. (1992).; (c) Ratio of the discharge hematocrits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 (a) Analytical solution of normalized cell flow rate as a function of hematocrit for channels with different widths; (b) Simulation results for cell flow rate in straight channel vs channel hematocrit; (c) Simulation results for cell flow rate in wavy channel vs channel hematocrit; (d) Ratio between RBC flux in wavy channel vs straight channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Suspension of RBCs flowing in a channel with wavy walls. The walls are shown in green, RBCs in red, and the edges of the CFL are shown in black. For this channel, W e f f = 57 µm and φ = 0.40. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Normalized CFL thickness, δ/W as a function of hematocrit for the straight-walled channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 RBC morphologies and CFLs (in yellow) for the smallest channel width W = 12 µm where the cell area fraction varies from φ = 0.05 -0.25 . . . . . . . . . . . . 4.11 Top: Initial RBC distribution for W = 22 µm, φ = 0.30 in a 22 µm × 120 µm channel. Bottom: RBC distribution at a later time. . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Dependence of the CFL width on the channel width and hematocrit for a straight channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 (a) CFL ratio as a function of channel width and hematocrit for the wavy-walled case; (b) Ratio of CFLs: wavy vs straight . . . . . . . . . . . . . . . . . . . . . . . . .

x 4 .

 4 14 (a) Evolution of the center of mass of N RBC = 4 in a straight channel of width W = 12µm and φ = 0.05; (b) Evolution of the center of mass of N RBC = 4 in a wavy channel of effective width W = 12µm and φ = 0.05 . . . . . . . . . . . . . . . 4.15 (a) Evolution of the center of mass of N RBC = 29 in a straight channel of width W = 12µm and φ = 0.40; (b) Evolution of the center of mass of N RBC = 29 in a wavy channel of width W = 12µm and φ = 0.40 . . . . . . . . . . . . . . . . . . . . 4.16 (a) Evolution of the center of mass of N RBC = 17 in a straight channel of width W = 57µm and φ = 0.05; (b) Evolution of the center of mass of N RBC = 17 in a wavy channel of width W = 57µm and φ = 0.05 . . . . . . . . . . . . . . . . . . . . 4.17 (a) Evolution of the center of mass of N RBC = 138 in a straight channel of width W = 57µm and φ = 0.40; (b) Evolution of the center of mass of N RBC = 138 in a wavy channel of width W = 57µm and φ = 0.40 . . . . . . . . . . . . . . . . . . . . 4.18 (a) Evolution of the elongation index of N RBC = 4 in a straight channel of width W = 12µm and φ = 0.05; (b) Evolution of the elongation index of N RBC = 4 in a wavy channel of width W = 12µm and φ = 0.05 . . . . . . . . . . . . . . . . . . . . 4.19 (a) Evolution of the elongation index of N RBC = 29 in a straight channel of width W = 12µm and φ = 0.40; (b) Evolution of the elongation index of N RBC = 29 in a wavy channel of width W = 12µm and φ = 0.40 . . . . . . . . . . . . . . . . . . . . 4.20 (a) Evolution of the elongation index of N RBC = 17 in a straight channel of width W = 57µm and φ = 0.05; (b) Evolution of the elongation index of N RBC = 17 in a wavy channel of width W = 57µm and φ = 0.05 . . . . . . . . . . . . . . . . . . . . 4.21 (a) Evolution of the elongation index of N RBC = 138 in a straight channel of width W = 57µm and φ = 0.40; (b) Evolution of the elongation index of N RBC = 138 in a wavy channel of width W = 57µm and φ = 0.40 . . . . . . . . . . . . . . . . . . . 4.22 WSS profile along the top and bottom walls of the channel (black). The CFL is shown in blue and the RBC membranes in red. The simulation is for a flat-walled channel with W = 12 µm and φ = 0.40. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Ratios of CFL and WSS for a hematocrit of φ = 0.40 to those for a RBC-free channel, τ 0 . The results are shown separately for the top and bottom walls of a 12 µmwide straight channel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 (a) Velocity field of subsection of simulation domain. Black dots represent RBC membrane coordinates. Arrows represent local velocity field and are colored according to streamwise velocity component; (b) Streamwise velocity profiles across the cell-wall gap at locations: x = 29 µm, x = 30 µm, x = 31 µm, and x = 32 µm. Dashed lines represent the linear velocity distribution according to Eq. 4.15. . . . . 4.25 Average perturbation WSS ratio as a function of channel width and hematocrit for a straight channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi 4 .

 4 26 (a) Snapshot of RBCs (red) flowing in a W = 57 µm channel with φ = 0.40 and instantaneous CFL (blue); (b) WSS profile (black) on the top and bottom walls (dashed and solid, respectively) and instantaneous CFL along the top and bottom walls for comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.27 (a) Velocity field of subsection of simulation domain. Black dots represent RBC membrane coordinates. Arrows represent local velocity vectors and are colored according to streamwise velocity component; (b) Streamwise velocity profiles across the cell-wall gap at locations: x = 72 µm, x = 73 µm, x = 74 µm, and x = 75 µm. Dashed lines represent the linear velocity distribution according to Eq. 4.15. . . . . 4.28 (a) Average perturbation WSS ratio as a function of channel width and hematocrit for a wavy-walled channel; (b) Ratio of average perturbation WSS wavy vs. straight108 5.1 (a) Concentration contour plot of 2D Gaussian hill at t = 200∆t with initial peak at x 0 = (200, 200) T ∆x (represented by blue star) for Pe = 0, D = 1.5∆x 2 /∆t; (b) Concentration profile at y = 200∆x. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 (a) Concentration contour plot of 2D Gaussian hill at t = 200∆t with initial peak at x 0 = (200, 200) T ∆x (represented by blue star) for Pe = 0, D = 0.0043∆x 2 /∆t; (b) Concentration profile at y = 220∆x. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 (a) Solute concentration initialized in the channel domain in the absence of an RBC (i.e. pure Poiseuille case) such that (xx c ) 2 + (yy C ) 2 ≤ R 2 where (x c , y c ) are the coordinate positions of the center and R = 2.6 µm; (b) Solute concentration initialized within RBC where α ≈ 1; (c) Solute concentration initialized within RBC where α = 0.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 (a) Flowfield relative to the laboratory reference frame; (b) Flowfield relative to the moving RBC with α ≈ 1; (c) Flowfield relative to the laboratory reference frame; (d) Flowfield relative to the moving RBC with α = 0.7 . . . . . . . . . . . . . . . . . 5.5 (a) Residence time probability distributions for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7. t = t/(W/2) 2 /D; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of higher Pe in the presence of RBC is to increase peak height. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Mean residence time. Without RBC, t m is independent of Pe. Increasing Pe results in a significant reduction in t m in the presence of RBC. . . . . . . . . . . . . . . . . . 5.7 Total remaining concentration of solute in the channel as a function of time. . . . . xii 5.8 (a) Probability distribution of the reaction position for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of higher Pe is to increase the peak height. . . . . . . 5.9 (a) Mean reaction position as a function of Pe; (b) Width of reaction position as a function of Pe. Higher Pe exhibit much more dispersion in their distributions than lower Pe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 (a) Residence time probability for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7. t = t/(W/2) 2 /D; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of higher Pe in the presence of RBC is to increase peak height. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Mean residence time as a function of Pe. . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 (a) Total remaining concentration in the channel domain as a function of time; (b) Total remaining concentration as a function of time for Pe = 1000. . . . . . . . . . . 5.13 (a) Total remaining concentration in the channel domain as a function of time; (b) Total remaining concentration as a function of time for Pe = 1000. . . . . . . . . . . 5.14 (a) Probability distribution of the reaction position for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of higher Pe is to increase the peak height. . . . . . . 5.15 (a) Mean reaction position as a function of Pe; (b) Width of reaction position as a function of Pe. Higher Pe exhibit much more dispersion in their distributions than lower Pe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Release positions for different angles of release, θ r , on the membrane for (a) α ≈ 1 and (b) α = 0.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.17 Mean residence time as a function of release angle . . . . . . . . . . . . . . . . . . . 5.18 Width of reaction position as a function of angle of release for (a) Pe = 0 (b) Pe = 100 (c) Pe = 1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii List of Tables 3.1 Minimum and maximum wall shear stress in the channel, influencing range R, and time duration T for the three different cases represented in Fig. 3.8. Lowercase and italicized t and b represent top and bottom wall, respectively. . . . . . . . . . . 50

  Fig 1.1 shows a schematic of the ramified cardiovascular system.

Figure 1

 1 Figure 1.1: Schematic of the ramified cardiovascular system showing artery branching into arterioles and capillaries followed by venules and veins which bring the blood back to the heart

Figure 1 . 2 :

 12 Figure 1.2: Scanning electron microscope showing the main components of blood: RBCs, WBCs, and platelets.

Figure 1

 1 Figure 1.3: (a) Biconcave RBC morphology; (b) Staining electron microscopy image of spectrin network; (c) Schematic of RBC membrane and cytoskeleton network

Figure 1

 1 Figure 1.4: Schematic showing the development of a spherocytotic RBC due to the deficiencies in the structural connections in the RBC. The deficiency in these structural connections results in a spheric morphology with reduced surface area to volume ratio. (Perrotta et al. 2008)

Figure 1

 1 Figure 1.5: A phase diagram taken from Kaoui et al. (2011) showing five distinct regions depending on flow (represented by the dimensionless number measuring the flow strength to the cell deformability) and the degree of confinement (represented by the vesicle diameter over the channel width)

  that induce active changes in the contractile state of vascular smooth muscle cells include neural, hormonal, metabolic, secreted, and mechanical factors. The endothelium releases nitric oxide (NO) into the neighboring smooth muscle cells which causes them to relax when a rise in shear stress is detected by the endothelium. These processes whereby mechanical cues are translated into biochemical responses by the cell are referred to as mechanotransduction. It is also hypothesized that RBCs are also involved in activating vasodilation, responding to mechanical cues or local hypoxic conditions by producing ATP which upon reaching the endothelium reacts to produce the vasodilator NO. It has been shown that when submitted to a large enough shear stress inducing a sufficient shape deformation, RBCs release ATP.[START_REF] Forsyth | Multiscale approach to link red blood cell dynamics, shear viscosity, and ATP release[END_REF] find ATP release constant below a shear stress threshold, whereas above the threshold, ATP release is increased and accompaned by large cellular deformations. They propose that Pannexin 1 hemichannel is the main avenue for ATP release independent of the shear stress level, and that the cystic fibrosis transmembrane conductance regulator (CFTR) only contributes to deformation-dependent release. Fig 1.6 shows a schematic for ATP release pathways from a RBC.

Figure 1 . 6 :

 16 Figure 1.6: ATP release pathways from a RBC. Pannexin 1 (Px 1) and Piezo1 can be activated by local shear stress; CFTR is activated by freed actin and upregulates Px1 (Zhang et al. 2018)

Fig 2 .

 2 1 graphically illustrates how post-collisional distributions at a given node fk are streamed to their new locations along the lattice, f k ( x + c k ∆t, t + ∆t).

Figure 2 . 1 :

 21 Figure 2.1: Graphical illustration of streaming step in LBM (Bao and Meskas 2014).

Fig

Figure 2 . 2 :

 22 Figure 2.2: The two-dimensional, nine-velocity (D2Q9) lattice model. There are 9 discrete velocities used to stream the distribution functions at each time step. c 0 is stationary; c 1-4 point in the east, north, west, and south directions, respectively; c 5-8 point in the north-east, north-west, south-west, and south-east directions, respectively.

Figure 2 .

 2 Figure 2.3: Schematic of the bounce-back boundary condition used in order to recover the no-slip boundary condition. The wall is located halfway between the solid nodes (depicted as open circles) and the adjacent fluid nodes (filled circles). The incoming post-collisional populations f 3 , f 6 , f 7 are reflected back equal and opposite into the fluid domain. (Bao and Meskas 2014).

Figure 2 . 4 :

 24 Figure 2.4: The discrete delta functions for the 3-point and the 4-point schemes (Shen 2016)

Figure 2 . 5 :

 25 Figure 2.5: Depiction of velocity interpolation on membrane node X n (t) (left) and force spreading to fluid node x (right). The hollow points denote the Eulerian fluid nodes, while the red dots depict the membrane nodes. The velocity is interpolated on node n from the velocity of the fluid nodes falling in the shaded region, while the force at x is determined from the forces of the membrane nodes falling in the shaded region.

Figure 2 . 6 :

 26 Figure 2.6: The membrane of the RBC is modelled as a collection of springs.[START_REF] Tsubota | Particle method for computer simulation of red blood cell motion in blood flow[END_REF]

  in the absence of flow are shown in Fig 2.7.

Figure 2 .

 2 Figure 2.7: Equilibrium shapes obtained for various reduced areas where ν represents the value of reduced area of the RBC (Shen 2016)

FigFigure 2

 2 Figure 2.8: Schematic of conditions for planar Poiseuille flow

  Fig 2.9 shows that the computed velocity profile closely agrees with the analytical solution.

Figure 2 .

 2 Figure 2.9: Parabolic velocity profile of planar Poiseuille flow showing close agreement of analytical and computed solutions

Figure 2 .

 2 Figure 2.10: (a) Unconfined slipper morphology (Ca = 10; Cn = 0.2); (b) Parachute morphology (Ca = 10; Cn = 0.45); (c) Confined slipper morphology (Ca = 10; Cn = 0.7)

  showed endothelial height profiles of around 3 µm for unsheared human umbilical vein endothelial cells (HUVECs); the resulting AFM images and cell profiles are shown in Fig 3.1.

Figure 3

 3 Figure 3.1: AFM images and cell profiles of unsheared human endothelial cells (HU-VECs) taken from Oberleithner, 2004. The image insets show lines corresponding to the height profile displayed at the bottom of each image. The study looked at the effect of the hormone aldosterone on HUVECs.

Figure 3 . 2 :

 32 Figure 3.2: 1D projection of regular lattice and curved boundary illustrating the Bouzidi interpolation bounceback scheme for curved boundaries

Figure 3

 3 Figure 3.3: An example of a domain with wavy boundaries

Figure 3

 3 Figure 3.4: A schematic of the domain showing part of the curved boundary in green. Nodes falling underneath the boundary are defined as solid nodes and are shown in green, and fluid nodes are shown in blue. The fluid node of interest is shown with a blue star. The line associated to the link 7 from this fluid node is shown as a solid red line, which passes through the bordering solid node denoted by a green star. The coordinates of the point of intersection must be known in order to calculate the distance necessary for applying the Bouzidi boundary condition.

Figure 3 . 5 :

 35 Figure 3.5: Bilinear interpolation. The value at the black spot is the sum of the value at each colored node multiplied by the area of the rectangle of the same color, divided by the total area of all four rectangles.

Chapter 3 .

 3 Figure 3.6: (a) Streamlines and velocity flow field corresponding to induced flow by a single RBC with a parachute morphology (C n = 0.45, Ca = 10, α = 0.7). Colorbar corresponds to relative magnitude of streamwise velocity to maximal velocity, U, in the channel; (b) Velocity vectors in the co-moving frame with the RBC showing zero velocity inside the RBC. Colorbar corresponds to relative magnitude of streamwise velocity in the co-moving frame to maximal velocity.

  ; Aouane 2015). The phase diagram showing different morphologies adopted as a function of C n and Ca is shown in Fig 3.7.

Figure 3 . 7 :

 37 Figure 3.7: Phase diagram showing different morphologies adopted as a function of the confinement ratio, C n , and the flow strength, C k W/R 0 , where C k is the capillary number, Ca, defined previously where the shear rate is defined as γ = R 0 U/(W/2) 2 . Figure taken from Aouane thesis (2015).

Chapter 3 .

 3 Figure 3.8: (a) Large channel (small confinement) and cell exhibiting unconfined slipper morphology and resulting wall shear stress on the top and bottom channel walls (C n = 0.2, CaW/R 0 = 15); (b) Increased confinement and cell exhibiting parachute morphology. Since the parachute is an axisymmetric shape, the wall shear stress curves for the top and bottom walls lie on top of each other (C n = 0.45, CaW/R 0 = 15); (c) Cell very confined and exhibiting confined slipper morphology (C n = 0.7, CaW/R 0 = 15). Note the different scales for shear stress for optimal visibility of the effect. x cm denotes the x-position of the center of mass of the RBC.

  Figure 3.9: Pressure gradient along the top and bottom walls of the channel. (a) Large channel (small confinement) with RBC exhibiting the unconfined slipper morphology (C n = 0.2, CaW/R 0 = 15); (b) Intermediate confinement with RBC exhibiting the parachute morphology. Since the parachute is an axisymmetric shape, the curves for the top and bottom walls lie on top of each other (C n = 0.45, CaW/R 0 = 15); (c) Small channel (large confinement) with RBC exhibiting the confined slipper morphology (C n = 0.7, CaW/R 0 = 15)

  Figure 3.10: Effect of reduced area on RBC morphology and wall shear stress for C n = 0.42 (W = 12 µm) and Ca = 0.1. (a) Stationary shapes obtained for different values of reduced area; (b) Maximum, τ max , and minimum, τ min , shear stress experienced in the channel with a single RBC normalized by the Poiseuille solution. The values for both the top and bottom of the channel are shown; (c) Plot of the normalized inverse minimum distance d -1 min W of the cell to the respective closest wall

  Figure 3.11: Effect of reduced area on RBC morphology and wall shear stress for C n = 0.42 (W = 12 µm) and Ca = 10. (a) Stationary shapes obtained for different values of reduced area; (b) Maximum and minimum shear stress experienced in the channel with a single RBC normalized by the Poiseuille solution. The values for both the top and bottom of the channel are shown; (c) Plot of the inverse minimum distance d -1 min W of the cell to the respective closest wall

  Figure 3.12: Effect of reduced area on RBC morphology and wall shear stress for C n = 0.72 (W = 7 µm) and Ca = 10. (a) Stationary shapes obtained for different values of reduced area; (b) Maximum and minimum shear stress experienced in channel with one RBC normalized by the Poiseuille solution. The values for both the top and bottom of the channel are shown; (c) Plot of the inverse of the normalized minimum distance, d -1 min W , of the cell to the respective closest wall

Figure 3

 3 Figure 3.13: Recasting results from Figs. 3.10c, 3.11c, and 3.12c such that τ max /τ Pois vs d -1 min W

Figure 3 .

 3 Figure 3.14: Normalized wall shear stress along a channel with a boundary given by a sinusoidal function where λ = 30 µm and A = 2 µm

  Figure 3.15: (a) One period of the wall geometry; (b) Average, maximum, and minimum WSS averaged in time a for channel containing a flowing RBC with α = 0.6; (c) RBC with α ≈ 1 and Ca = 0.1 and Ca = 10. The corresponding morphologies as the RBC transverses the undulating wall are shown in Fig. 10a and 11a and Fig. 10e.

  Figure 3.16: (a) Maximum and minimum wall shear stress along the top and bottom walls of a wavy-walled channel as a function of reduced area with C n = 0.42 (W e f f = 12 µm ) and Ca = 0.1 and Ca = 10; (b) The inverse of the normalized minimum distance as a function of reduced area

  Figure3.16a shows that the maximum wall shear stress in the wavy wall case is less sensitive to changes in reduced area than in the flat wall case. In extreme, the maximum wall shear stress only varies by ∼ 10%, whereas in the flat wall case it can vary by ∼ 30% among reduced areas. Waviness reduces the variation in maximum shear stresses observed for different reduced areas; however, the values are in general higher than the equivalent in the straight case. The inverse of the minimum distance normalized to the effective width of the channel is not predictive of the maximum wall shear stress as is suggested by Fig.3.16b.

  Figure3.17 shows snapshots of one RBC as it flows over one period of wall undulation. We see that what would be a stationary morphology in the case of a flat-walled channel is no longer so due to the variation of the flow field along the undulation. Because the RBC is effectively "rigid," it does not undergo a noticeable deformation in either the streamwise or perpendicular streamwise directions as it passes through the constriction and expansion of one wall undulation. Instead, we see that the tilt angle of the RBC changes over one period. As the cell approaches the maximal constriction of the two walls, the RBC tilts towards the horizontal axis, and as the cell exits the constriction, the tilt moves toward the vertical axis.The results above were for a rigid RBC (Ca=0.1). For Ca=10, representing a softer cell, all morphologies obtained are axisymmetric, so the maximum shear stress along the top and bottom walls are identical, as shown in Fig.3.16a.

  Figure 3.19: Five snapshots of the morphology of a single RBC as it passes over one wall undulation for different dephasing angles between the top and bottom walls (a) n = 0; (b) n = 1/5; (c) n = 2/5; (d) n = 3/5

Figure 3 .

 3 Figure 3.20: Maximum wall shear stress normalized to the Poiseuille solution (for a flat-walled channel having the same effective width) along the channel walls as a function of the degree of dephasing, n

  Figure 3.21: (a) Channel where amplitude is divided in half (i.e., A/W = 0.2). We no longer see the recirculation zones for Re = [1, 10]; (b) Re = 0.5; (c) Re = 5; (d) Re = 50; (e) Re = 125; (f) Height of the center of the recirculation zone, h, normalized to the channel amplitude, A, as a function of Re

Chapter 3 .

 3 Figure 3.22: Trajectories of the center of mass of a "rigid" RBC (Ca = 0.25) initially placed at (a) y0 = 25, y0 = 30, y0 = 50, y0 = 100 in the valley; (b) y0 = 100, y0 = 145; (c) y0 = 150 with the RBC positioned with a 0 • , 45 • , or 90 • angle to the horizontal; (d) Snapshot of trapped RBC y0 = 25

  Figure 3.23: Trajectories of the center of mass of a "soft" RBC (Ca = 5) initially placed at (a) y0 = 25, y0 = 30, y0 = 50, y0 = 100; (b) y0 = 100, y0 = 145; (c) y0 = 150 with the RBC positioned with a 0 • , 45 • , or 90 • angle to the horizontal

Chapter 3 .

 3 Figure 3.24: (a) Morphology exhibited in valley region; (b) Morphology exhibited in peak region

Figure 4 . 1 :

 41 Figure 4.1: Relative apparent viscosity as a function of channel width , W, and channel hematocrit, φ, for a channel with straight walls. Dashed lines represent the empirical relation of Pries et al.(1992)

Figure 4 . 2 :

 42 Figure 4.2: Top: Placement of cells in the channel for a hematocrit of φ = 0.15; Bottom: Placement of cells in the channel for a larger hematocrit of φ = 0.20. The RBCs are less extended in the perpendicular streamwise direction for φ = 0.20.

Figure 4 .

 4 Figure 4.3 illustrates the streamwise averaged velocity profiles in the channel for the two hematocrits compared to the Poiseuille solution.

Figure 4 . 3 :

 43 Figure 4.3: Streamwise averaged velocity profiles

Figure 4 . 4 :

 44 Figure 4.4: Relative apparent viscosity as a function of hematocrit in capillaries.

Figure 4 . 6 :

 46 Figure 4.6: (a) Ratio of channel hematocrit, φ, to discharge hematocrit, φ d , versus channel hematocrit for the straight case; (b) Ratio of tube hematocrit, φ, to discharge hematocrit, φ d , versus tube hematocrit for the wavy-walled case. Dashed lines in (a) and (b) represent empirical relation given by Pries at al. (1992).; (c) Ratio of the discharge hematocrits

Chapter 4 .

 4 Figure 4.7: (a) Analytical solution of normalized cell flow rate as a function of hematocrit for channels with different widths; (b) Simulation results for cell flow rate in straight channel vs channel hematocrit; (c) Simulation results for cell flow rate in wavy channel vs channel hematocrit; (d) Ratio between RBC flux in wavy channel vs straight channel

  Fig 4.8.

Figure 4 . 8 :

 48 Figure 4.8: Suspension of RBCs flowing in a channel with wavy walls. The walls are shown in green, RBCs in red, and the edges of the CFL are shown in black. For this channel, W e f f = 57 µm and φ = 0.40.

Figure 4

 4 Figure 4.9: Normalized CFL thickness, δ/W as a function of hematocrit for the straight-walled channels

  initial RBC distribution contained 6 columns of RBCs, while ours, where the channel length is nearly double that of Zhang et al. contains 12 as shown in the top panel of Fig. 4.11.

Figure 4 .

 4 Figure 4.11: Top: Initial RBC distribution for W = 22 µm, φ = 0.30 in a 22 µm × 120 µm channel. Bottom: RBC distribution at a later time.

  Figure 4.13: (a) CFL ratio as a function of channel width and hematocrit for the wavywalled case; (b) Ratio of CFLs: wavy vs straight

  Figure 4.14: (a) Evolution of the center of mass of N RBC = 4 in a straight channel of width W = 12µm and φ = 0.05; (b) Evolution of the center of mass of N RBC = 4 in a wavy channel of effective width W = 12µm and φ = 0.05

Figure 4 .

 4 Figure 4.15: (a) Evolution of the center of mass of N RBC = 29 in a straight channel of width W = 12µm and φ = 0.40; (b) Evolution of the center of mass of N RBC = 29 in a wavy channel of width W = 12µm and φ = 0.40

  Figure 4.16: (a) Evolution of the center of mass of N RBC = 17 in a straight channel of width W = 57µm and φ = 0.05; (b) Evolution of the center of mass of N RBC = 17 in a wavy channel of width W = 57µm and φ = 0.05

  Figure 4.17: (a) Evolution of the center of mass of N RBC = 138 in a straight channel of width W = 57µm and φ = 0.40; (b) Evolution of the center of mass of N RBC = 138 in a wavy channel of width W = 57µm and φ = 0.40

  Figure 4.18: (a) Evolution of the elongation index of N RBC = 4 in a straight channel of width W = 12µm and φ = 0.05; (b) Evolution of the elongation index of N RBC = 4 in a wavy channel of width W = 12µm and φ = 0.05

  Figure 4.19: (a) Evolution of the elongation index of N RBC = 29 in a straight channel of width W = 12µm and φ = 0.40; (b) Evolution of the elongation index of N RBC = 29 in a wavy channel of width W = 12µm and φ = 0.40

  Figure 4.20: (a) Evolution of the elongation index of N RBC = 17 in a straight channel of width W = 57µm and φ = 0.05; (b) Evolution of the elongation index of N RBC = 17 in a wavy channel of width W = 57µm and φ = 0.05

  Figure 4.21: (a) Evolution of the elongation index of N RBC = 138 in a straight channel of width W = 57µm and φ = 0.40; (b) Evolution of the elongation index of N RBC = 138 in a wavy channel of width W = 57µm and φ = 0.40

Figure 4 .

 4 Figure 4.22: WSS profile along the top and bottom walls of the channel (black). The CFL is shown in blue and the RBC membranes in red. The simulation is for a flatwalled channel with W = 12 µm and φ = 0.40.

Figure 4 .

 4 Figure 4.25: Average perturbation WSS ratio as a function of channel width and hematocrit for a straight channel

  Figure 4.27: (a) Velocity field of subsection of simulation domain. Black dots represent RBC membrane coordinates. Arrows represent local velocity vectors and are colored according to streamwise velocity component; (b) Streamwise velocity profiles across the cell-wall gap at locations: x = 72 µm, x = 73 µm, x = 74 µm, and

  Figure 4.28: (a) Average perturbation WSS ratio as a function of channel width and hematocrit for a wavy-walled channel; (b) Ratio of average perturbation WSS wavy vs. straight

  WSS ratio in Fig 4.28a largely follows the inverse trend of the average CFL thickness, shown in Fig 4.13a.

  Figure 5.1: (a) Concentration contour plot of 2D Gaussian hill at t = 200∆t with initial peak at x 0 = (200, 200) T ∆x (represented by blue star) for Pe = 0, D = 1.5∆x 2 /∆t; (b) Concentration profile at y = 200∆x.

  Fig. 5.2 shows results for this case when t = 200∆t.

  Figure 5.2: (a) Concentration contour plot of 2D Gaussian hill at t = 200∆t with initial peak at x 0 = (200, 200) T ∆x (represented by blue star) for Pe = 0, D = 0.0043∆x 2 /∆t; (b) Concentration profile at y = 220∆x.

  Figure 5.3: (a) Solute concentration initialized in the channel domain in the absence of an RBC (i.e. pure Poiseuille case) such that (xx c ) 2 + (yy C ) 2 ≤ R 2 where (x c , y c ) are the coordinate positions of the center and R = 2.6 µm; (b) Solute concentration initialized within RBC where α ≈ 1; (c) Solute concentration initialized within RBC where α = 0.7

Figure 5

 5 Figure 5.4: (a) Flowfield relative to the laboratory reference frame; (b) Flowfield relative to the moving RBC with α ≈ 1; (c) Flowfield relative to the laboratory reference frame; (d) Flowfield relative to the moving RBC with α = 0.7

  Figure 5.5: (a) Residence time probability distributions for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7. t = t/(W/2) 2 /D; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of higher Pe in the presence of RBC is to increase peak height.

Figure 5 . 6 :

 56 Figure 5.6: Mean residence time. Without RBC, t m is independent of Pe. IncreasingPe results in a significant reduction in t m in the presence of RBC.

Figure 5 . 7 :

 57 Figure 5.7: Total remaining concentration of solute in the channel as a function of time.

  Figure 5.8: (a) Probability distribution of the reaction position for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, andwith RBC α = 0.7. The effect of higher Pe is to increase the peak height.

  Figure 5.9: (a) Mean reaction position as a function of Pe; (b) Width of reaction position as a function of Pe. Higher Pe exhibit much more dispersion in their distributions than lower Pe.

  Figure 5.10: (a) Residence time probability for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7. t = t/(W/2) 2 /D; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of higher Pe in the presence of RBC is to increase peak height.

Figure 5 .

 5 Figure 5.11: Mean residence time as a function of Pe.

  Figure 5.12: (a) Total remaining concentration in the channel domain as a function of time; (b) Total remaining concentration as a function of time for Pe = 1000.

Figure 5 .

 5 Figure 5.13: (a) Total remaining concentration in the channel domain as a function of time; (b) Total remaining concentration as a function of time for Pe = 1000.

Chapter 5 .

 5 Figure 5.14: (a) Probability distribution of the reaction position for all cases: without RBC, with RBC α ≈ 1, and with RBC α = 0.7; (b) Pe = 10 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7. The effect of RBC is negligible; (c) Pe = 100 w/o RBC, with RBC α ≈ 1, and with RBC α = 0.7; (d) Pe = 1000 w/o RBC, with RBC α ≈ 1, andwith RBC α = 0.7. The effect of higher Pe is to increase the peak height.

  Figure 5.15: (a) Mean reaction position as a function of Pe; (b) Width of reaction position as a function of Pe. Higher Pe exhibit much more dispersion in their distributions than lower Pe.

Figure 5

 5 Figure 5.16: Release positions for different angles of release, θ r , on the membrane for (a) α ≈ 1 and (b) α = 0.7

Figure 5 .

 5 Figure 5.17: Mean residence time as a function of release angle

  Figure 5.18: Width of reaction position as a function of angle of release for (a) Pe = 0 (b) Pe = 100 (c) Pe = 1000

  .8, where min(τ/τ Pois ) corresponds to the minimum Chapter 3. Results: Red Blood Cells under Channel Flow: Effect of Endothelium and Red Blood Cell Morphology

Table 3 .

 3 1: Minimum and maximum wall shear stress in the channel, influencing range R, and time duration T for the three different cases represented in Fig.3.8.Lowercase and italicized t and b represent top and bottom wall, respectively.

	Case min(τ/τ Pois ) max(τ/τ Pois ) IR[µm]	T[ms]
	a	0.98 (t)	1.01 (t)	-	-
	b	0.91	1.18	15	5.80e-2
	c	0.83 (b)	1.26 (t)	12.2 (t); 14.2 (b) 4.70e-2 (t); 5.50e-2 (b)

(a) (b)

(a) (b)

Chapter 4

Suspensions in wavy-walled channels

We model blood as a suspension of RBCs flowing in microchannels with undulating walls in order to examine the cell-free layer, shear stress, and hemorheological properties. We model the RBCs as two-dimensional deformable biconcave vesicles using experimentally-measured properties. Through our simulations, we are able to observe several important characteristics of microscopic blood flows observed experimentally including the cell-free layer, Fahraeus-Lindqvist, and the Fahraeus effects. Channel widths range from 10 -60 µm and hematocrit values from 5 -40%, corresponding to conditions which might be present in the microvasculature. Section 4.1 will introduce the important flow anamolies found in blood flow in the microvasculature. Section 4.2 presents the results of our simulations with subsections dedicated to the Fahraeus-Lindqvist Effect (aparent viscosity), Fahraeus Effect (discharge hematocrit), cell-free layer, and wall shear stress. The conclusions of this chapter are summarized in Section 4.3.

Introduction

Blood is a suspension consisting of RBCs, white blood cells, platelets, and plasma containing various molecules and ions. RBCs make up the largest volume fraction among the various components at approximately 45% of the total volume, while white blood cells constitute ∼ 0.7%, and platelets make up ∼ 1/800 of total blood volume. The non-Newtonian dynamics of blood arise from the presence of discrete RBCs. An accurate description of blood flow must take them into consideration, especially in the microvasculature, for their presence leads to important consequences for hemodynamic and hemorheological properties.

When flowing in microvessels, it was observed that RBCs migrate toward the center region of the vessel, resulting in a core of high RBC density and leaving a layer of plasma adjacent to the vessel wall, known as the cell free layer (CFL). The CFL is an important property of blood flow in the microcirculation, since it has significant effects on flow resistance as well as transport of agents at the vessel wall [START_REF] Popel | Microcirculation and Hemorheology[END_REF]. Directly related to the CFL, another important phenomenon when considering the particulate nature of blood flow is the observation that the apparent viscosity decreases with decreasing diameter of vessel, a phenomenon known as the Fahraeus-Lindqvist Effect. The Fahreus-Lindqvist Effect was discovered in the seminal We recall similar conclusions made about the instantaneous CFL in Chapter 3 of this thesis and in previous works: we observe a peak in the wall shear stress when the cell is further from the wall and a valley when the cell is close to the wall, features which we can observe in Fig. 4.23; for instance, at x ≈ 100 on the top wall, we observe a local minimum in the WSS and a dip in the CFL, while at x ≈ 60 on the top wall, we observe that a large cell-wall gap corresponds to a local peak in WSS. [START_REF] Namgung | Effect of cell-free layer variation on arteriolar wall shear stress[END_REF] studying the relationship between WSS and flow rate and RBC aggregation, proposed the following Couette flow-type formula to obtain the local instantaneous WSS, τ, from their CFL measurements:

where µ ext is the plasma viscosity, V edge is the mean edge velocity of the RBC core, and δ is the local, instantaneous CFL thickness. However, this relationship would predict the inverse of what we have observed: that is, it predicts a low WSS at a location with a wide gap and vice-versa. To further understand, we examined the simulated velocity field in the gap between RBCs and the wall. A part of the simulation domain from Chapter 5

Advection-Diffusion of Passive Scalar Released by Red Blood Cells in Flow

In this chapter, we describe a 2D advection-diffusion solver with moving particles. It represents a general numerical model for the study of of advection-diffusion problems coupled with flow and moving particles. We use the model to study the problem of solute release from particles in the flow and the subsequent advection and diffusion in the flow to the channel walls. We simulate realistic RBC morphologies and conditions in the microvasculature. Two modes of release are investigated: one where the solute is located internally within the RBC and flux across the membrane is generated by the gradient of concentration and another where the solute concentration is initialized on the membrane of the RBC and is released to the plasma domain at a constant rate. We investigated the effect of flow and RBC shape on the residence time of the solute in the flow before reaction and its diluteness upon reaction with the endothelium. Section 5.1 briefly discusses potential applications of our numerical model. Section 5.2 presents the LBM framework, and Section 5.3 presents the results of brenchmarking of our advection-diffusion code. Section 5.4 introduces the fundamentals of residence time theory from which we can understand the two main deliverables we seek: residence time and diluteness at the wall. Section 5.5 presents the results for the diffusion of an internal solute, while Section 5.6 presents those for the solute release at the membrane. Section 5.7 summarizes our conclusions.

Chapter 5. Advection-Diffusion of Passive Scalar Released by Red Blood Cells in Flow

Residence Time Theory

At time t = 0, the solute species is introduced inside the RBC. The membrane is perfectly permeable and the species will be advected and diffused in the flow before reaching the channel walls where it is immediately and completely absorbed (i.e., perfectly absorbing walls). Depending on the initial spatial distribution and flow conditions, the solute particles will reach and react with the channel walls at different times. We can define a distribution function representing the residence time distribution, E(t), such that it has the properties of a probability distribution: it is always non-negative and ∞ 0 E(t)dt = 1. E(t) describes the amount of time the solute has spent within the channel domain such that:

(5.12)

where the numerator, C exit (t), is the total solute concentration leaving the system at any time t, which can be found from the integral of the flux across the channel walls, and the denominator is the cumulative sum of the solute concentration leaving the system for infinite time, which in the case of perfectly absorbing walls is equal to the initial amount of solute injected into the system, C 0 . Rewriting Eq. 5.12, we find

where j walls is the flux through the top and bottom channel walls. This can be calculated from our simulations where the solute is introduced into the system at time t = 0 by measuring the flux of the solute concentration across the channel walls as a function of time. The mean residence time is given by the first moment of the residence time distribution:

Another quantity of interest is the reaction position on the channel walls for which we define a probability distribution of reaction position similarly to the residence time distribution such that:

which is non-negative and ∞ -∞ P(x)dx = 1. We look at the mean reaction position as well as the second moment of the reaction position distribution (i.e., standard deviation of the distribution), which is a measure of the diluteness of the solute upon reaction with the wall and is calculated as the second moment of the reaction position distribution such that:

(5.16)