We start by summarizing the state-of-the-art physical models employed in the description of blood flow in chapter 2, because the choices made here affect all higher levels. The models that are analyzed in-depth in this thesis are the Stokes equation (modeling the hydrodynamics) and the Helfrich model (for the bending resistance of red blood cells). Our contribution to the
1.3 Contributions and outline
Physical phenomena (real world) Physical models (ch.2, [Pub2])
Mathematical considerations (ch.3, [Pub3]) Algorithms (ch.4and5, [Pub1] and [Pub2])
Numerical implementation (ch.5.3)
Validation (ch.5.3, [Pub1] and [Pub5]) Simulations
→Fundamental questions (ch.6.1, [Pub4] and [Pub5])
→Application-oriented questions (ch.6.2, [Pub3]) Explanation and/or new phenomena
Fig. 1.1: The various layers when doing research via numerical simulations. The present thesis makes contributions to each one, located in the chapters listed within the diagram and further detailed in the dedicated publications.
model layer is a rigorous comparison between the common Canham-Helfrich bending model and so-called linear bending models [Pub2].1
We then proceed in chapter3with presenting the mathematics of the hydrodynamic formula-tion, which is based on the Stokes equation. We extend the standard boundary integral equation to include volume-changing objects [Pub3]. The development of the extension was necessary in order to include oscillating and deformable microbubbles in blood flow. Existing simulation methods are not capable of handling red blood cells together with such bubbles.
On the algorithmic part, we shed for the first time light on the performance of various algorithms for the bending forces in chapter 4.3 ([Pub1] and [Pub2]1). Bending forces are required here for the proper simulation of red blood cells, but are also very important for related fields of research (capsules, vesicles, etc.) and even unrelated ones (computer graphics etc.). An intermediate component (the mean curvature) is also necessary for the microbubbles.
Going one layer up to the actual implementation, the author wrote the base of the code for the boundary integral method already during his master thesis [12]. Creating a new code was necessary as no suitable and publicly available BIM code existed. During his doctorate, he imple-mented necessary extensions (such as for the bubbles) and several performance optimizations.
He also implemented an MPI parallelization in the course of a KONWIHR research project, as outlined in chapter5.2.
Continuing, the validation layer happened to some degree also during the author’s master thesis.
Nevertheless, our discussion on the bending algorithms in chapter 4.3.3 provides additional evidence for the correctness of the code [Pub1]. Most importantly, the research on the behavior of single RBCs and the comparison with new experiments by our collaborators (chapter6.1) puts all layers to the test – successfully [Pub5].
1Despite being a review to large extents, publication [Pub2] contains two novel contributions: First, the rigorous comparison between the linear bending models and the Canham-Helfrich model and second, the computation of the bending errors when using a spherical harmonics expansion.
1 Introduction
Finally, we contribute to the highest layer and thus to the knowledge on blood flow by means of two studies. The first one is the just mentioned research, where we consider individual red blood cells flowing through a microchannel. Hence, it is important for the fundamental understanding of blood but also for certain applications. Our comparison between numerical and experimental measurements constitute the first of its kind with such details. We recover the two well-known common shapes assumed by the cells, croissants and slippers. Our most important finding is bistability and its systematic analysis, i.e. the coexistence of these shapes at the same set of parameters. It is presented in chapter 6.1 [Pub5]. The second study in chapter 6.2 [Pub3]
considers lipid-coated microbubbles in blood flow. As such, it is more about a specific application rather than fundamental insights. We show that application of an ultrasound (causing the bubbles to oscillate) together with the special properties of the lipid coating and the interaction with the red blood cells causes the bubbles to migrate to the vessel wall (margination). Thus we can conclude that such bubbles constitute an efficient agent in targeted drug delivery protocols. The results also highlight that e.g. in-vitro experiments must take red blood cells into account in order to arrive at conclusions that can be translated to the in-vivo case.
Since the present thesis is in the format of a cumulative dissertation, the following chapters will serve as a guide along the path from the lowest to the highest layer, while providing context and summaries of our individual contributions. The full details can be found in the attached publications (partII).
2 Biology and models of blood flow
In order to capture the physical behavior of blood flow, the two most important components that need to be modeled are the red blood cells and the fluid inside and outside of the cells. The following sections will shortly outline the biology and the corresponding models.
2.1 Structure and modeling of red blood cells
Overview Red blood cells (erythrocytes) are created in the bone marrow and lack a nucleus, i.e.
consist of only a membrane and an internal fluid (cytoplasm) [13]. The cytoplasm contains a high concentration of hemoglobin, a protein capable of binding oxygen molecules. It is therefore imperative for the cells in order to fulfill their main purpose, the delivery of oxygen from the lungs to the rest of the body. This inner fluid is Newtonian [14] and has a dynamic viscosity of µRBC = 5 – 15 mPa s for healthy human red blood cells at 37◦C [13]. The precise value depends to a large extent on the hemoglobin concentration, with older cells exhibiting higher concentrations and, therefore, also higher viscosities [13].
The only structural part of the cell, its membrane, is formed by a lipid bilayer and an underlying cytoskeleton, both tethered together by several types of transmembrane proteins [13]. More precisely, the cytoskeleton consists of a complex and possibly dynamic [15] network of spectrin proteins which forms a triangular mesh when expanded [16]. This network is responsible for providing the membrane with its shear elasticity. The lipid bilayer, on the other hand, is composed of two layers of phospholipids and cholesterol. It endows the cell with some bending rigidity [17].
Furthermore, the area of the composite membrane stays nearly constant under deformations (at
≈140µm2[18,19]), i.e. it has a very high surface area dilatation modulus. This is to a large extent due to the lipid bilayer, while the cytoskeleton’s dilatation resistance is notably smaller [17,20]. Longer-lasting changes of the area above≈4 % lead to cell rupturing [13,21], while short lasting (≈50µs) changes of up to approximately 40 % can be coped with [22].
Modeling the in-plane elasticity Since the membrane has a thickness of typically less than 100 nm [16] while the RBC in equilibrium has a typical diameter of 8µm [23], erythrocytes in flow simulations are usually modeled as inert objects having a homogeneous thin shell with zero thickness. Skalak et al. [24] introduced one of the most often used models for RBCs, the so-called Skalak model, which we will also employ in this thesis. It consists of a strain-hardening [25] component to model the shear elasticity with the corresponding shear modulus2 κS≈5×10−6N/m [27,28], and an area dilatation component with a modulusκA. As explained above,κAshould be high enough to get an almost constant surface area, with appropriate values of aroundκA = 103κS to105κS [21,24]. Such high moduli can impede a large performance penalty in numerical simulations, and hence smaller values are often chosen. Additional surface
2Note that different conventions exist for the definition ofκS. Here we use the one by Krüger [26].
2 Biology and models of blood flow
forces or constraints can be introduced to limit the surface deviations to large degrees [26,29].
Naturally, the Skalak model is not the only possible choice [30,31], with spring-network models having become another popular alternative [32–34].
Modeling the bending rigidity The bending rigidity of the RBC membrane is usually taken into account via the famous Helfrich (or Canham-Helfrich) model [35–37], which is prescribed in the form of a surface energy. A typical bending modulus isκB≈3×10−19N m [28,38]. The Helfrich model can be amended by an “area-difference elasticity” (ADE) term [21,39], but it was apparently not used so far for RBCs in flow except in reference [40]. Bending forces are important to correctly capture equilibrium shapes when the shear stresses are nearly zero [21,41]
or to describe the wrinkles [42,43] and shapes [11] that can appear during larger deformations.
They can also improve the stability of numerical simulations [25]. Although the Helfrich model is used in a large part of the literature, a proper algorithm is still not established and different results are obtained in practice depending on the algorithm as we show in [Pub1].
Several publications [21,42,44–52] do not employ the Helfrich model but rather a so-called
“linear bending model”. Different flavors exist. All of them have in common that the constitutive equation is not an energy, but rather the bending moments Mαβ are prescribed directly. The reference state aside, this is eitherMαβ ∼Haαβ orMαβ ∼bαβ whereaαβ is the metric tensor, bαβ the curvature tensor andHthe mean curvature. Their relationship to the Helfrich model was so far not completely clear. Only Pozrikidis [46, p. 279] [48, ch. 2.8.2] made some rather cryptic statements, suggesting that some of them should be equivalent in the 2D case.
In publication [Pub2] we show mathematically thatMαβ ∼Haαβ andMαβ ∼bαβ lead to the same traction jump (force per unit area) in 3D, but that this traction jump matches with the Helfrich result only to leading order, i.e. small deformations. Moreover, the higher orders lead to an additional qualitative difference in the tangential component of the traction jump: It is zero for the Helfrich model but non-zero for the linear bending models. Figure2.1exemplifies the traction jump as obtained from the Helfrich and the linear bending models for the biconcave discocyte shape, showing deviations of up to 35 %. Linear bending models are only equivalent if they prescribe in-plane tensions correctly, something which is usually not the case.
0 50 100 150 200 250
0 π/8 π/4 3π/8 π/2
(a)
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
0 π/8 π/4 3π/8 π/2
(b)
Tractionjump|4f|
Θ Helfrich model
Lin. bend.
Relativedeviation
Θ
Fig. 2.1: (a) The traction jump4fas obtained from the full Helfrich model in comparison with the linear bending models that neglect in-plane tensions for the discocyte shape (compare figure2.2). Units measured inκBand the large radius of the discocyte. (b) Relative deviation between the data from (a).Θ := arccos x3/p
x21+x22+x23 is the polar angle as measured from the center of the discocyte. Reprinted from publication [Pub2] with permission from IOP Publishing.
2.1 Structure and modeling of red blood cells
Stress-free shape An open question in modeling the cell elasticity is that of the stress-free (or reference) shape. Due to the twofold nature of the RBC membrane, the question applies to both the lipid bilayer as well as the cytoskeleton. For the first, the reference shape is typically encoded in the Helfrich model via the so-called spontaneous (or reference) curvature c0 (or H0). The asymmetric distribution of the phospholipids between the two layers [13], which are themselves homogeneous on the µm-scale, indicates a non-zero value for c0 that is spatially constant.
Unfortunately, experimental values are non-existent so far. The most common choice is therefore c0 = 0.
For the second constituent of the membrane, the cytoskeleton, a spatially inhomogeneous reference shape is required to explain results by Fischer [53]. He showed that the membrane is endowed with a “shape memory”: Membrane elements before a deformation and after equili-brating the deformed cell again are found at the same location. Two suggestions are currently discussed in the literature as possible candidates for the stress-free shape [41]: Either the disco-cyte equilibrium shape or an oblate spheroid that is nearly (but not completely) a sphere. So far the debate has not been settled, but more and more recent studies claim that the oblate shape is the correct one [54–61], although converse arguments also exist [62]. Analyses are made more complicated by the observation that both stress-free shapes result in distinct dynamics only for certain parameter ranges. If the ratio λ of the inner and outer viscosities is in the physiological range (λ≈5), the differences can be negligible [57]. In our studies in chapter6we use the discocyte reference state. We make a more in-depth comment on this issue in the outlook (chapter7.2.1).
Further details The volume of the RBC is set by osmotic balance. Hence, it stays constant as long as the environment does not change [21,63]. Simulations of cells in flow therefore assume a constant volume of around≈100µm3[18,19,64]. Still, although the underlying equations already describe an impenetrable membrane, an artificial volume drift can occur in practice as a numerical artifact. It can be countered by introducing an additional ad-hoc volume force [26], exploiting the no-flux condition [29] or by a simple volume rescaling procedure [29].
In the absence of external forces, a healthy RBC equilibrates to a biconcave shape called discocyte with a diameter of approximately 8µm and a thickness of 2.6µm [23] as depicted in figure2.2. A convenient mathematical formula to describe this geometry is given by Evans and Fung [23]. This shape is determined by the available excess area compared to a sphere with the same volume and the membrane forces [21]. In non-physiological environments [21, 54] or in case of diseases [65–67], other shapes such as stomatocytes or echinocytes are also possible. Application of external stresses (for example by placing the cells into a flow) leads to large reversible deformations as discussed in publication [Pub5] (chapter6.1).
Fig. 2.2: The biconcave discocyte is the typical equilib-rium shape of a healthy red blood cell. Half of it was made transparent for illustration purposes.
Red blood cells show additional features that were not mentioned so far and are not taken into account in the present work. The first one is membrane viscosity [7,28,38,68–71], which is
2 Biology and models of blood flow
often neglected in simulations [28] if the chosen method does not incorporate it as an intrinsic feature (such as in dissipative particle methods [33]). At least in some cases, the effect of membrane viscosity can also be modeled by a higher viscosity for the inner fluid [40,60,72].
Second, it was recently shown that the small fluctuations of the membrane cannot be explained only by the finite temperature of the environment, but rather that some part of it comes from active processes in the membrane [73]. The existence of fluctuations, active or passive, plays an important role in the determination of the elastic moduli [38]. For cells in flow they can, however, usually be neglected [28] except at very low flow rates where diffusion effects become important [11, 74] and near transition boundaries [75,76]. As for the membrane viscosity, they are mostly (but not always [77,78]) taken into account only if the method provides them naturally, although an actual comparison with and without fluctuations is missing so far. The third one is the≈6 nm [79], i.e. very thin brush-like glycocalyx covering the membrane surface. Due to its thinness, it was not yet taken into account and is not expected to influence the dynamics of fast flowing cells, although it might play a role in the formation of rouleaux [80,81].
We finally remark that the above given values are onlytypicalvalues forhuman healthyred blood cells in physiological conditions. The values obtained from different experiments can scatter a lot (compare the overviews given in [19, 21, 26,28]) and can also depend on the environment [21]. Pathological [3,67,82,83], aged [13] and cells from other mammals [62,84, 85] often have different properties.