A particle-based model for computing fluid flows and cell dynamics by Seyed Majid Hosseini Amin B.Sc., Sharif University of Technology, 2002 M.Sc., Sharif University of Technology, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Chemical and Biological Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2012 c© Seyed Majid Hosseini Amin 2012Abstract The connection between certain human diseases and changes in the me- chanical properties of living cells is well established, e.g. in the cases of malaria and cancer. However, the mechanism for the mechanical modifica- tions, which tend to facilitate the pathogenesis of such diseases, is not always clear. For instance, the overall loss of deformability of malaria-infected red blood cells (RBCs) corresponds to a 10-fold increase in the rigidity of the cell membrane. On the other hand, micropipette aspiration has only measured a 3-fold increase in the elastic modulus. In this thesis, a particle-based model is developed to explore the inter- play between the underlying microstructures and the behavior of the cell as a whole. The research consists of three related projects. The first project deals with the long-standing problem of Smoothed Particle Hydrodynamics (SPH) method with open boundaries and solid walls. We propose a “rotational pressure-correction scheme” with a consistent pressure boundary condition that leads to a large improvement in accuracy of calculated pressure and the drag coefficient on solid obstacles. The second and third projects concern developing a 2D and then 3D particle-based model for RBCs to explore the parasite-driven changes in malaria-infected RBCs. In our models the cell membrane is replaced by a set of discrete particles connected by linear or nonlinear springs. In addition, a linear bending elasticity is implemented using the deviation of the local curvature from the innate curvature of the biconcave shape of a resting RBC. The cytoplasm and the external liquid are modelled as homogeneous Newtonian fluids, and discretized by particles as in standard SPH solution of the Navier-Stokes equations. The malaria parasite is modelled as an aggregate of particles constrained to rigid-body motion. We argue that the discrepancy in the estimated elastic modulus of the membrane is caused by the presence of the rigid parasite particles inside infected cells, and have carried out numerical simulations to demon- strate this mechanism. Our three-dimensional simulation of RBC stretching tests by optical tweezers accurately demonstrates the compensating effects between the existence of malaria parasites and the elevated stiffness of the membrane on the overall deformability of infected RBC. iiPreface This thesis entitled “A particle-based model for computing fluid flows and cell dynamics” presents the research that I led and performed during my PhD study. The research conducted in this thesis was supervised, identified, and initiated by Professor J. J. Feng. In this preface, the contributions and collaborations to the papers published or submitted for publication from current thesis are briefly explained. • A version of chapter 5 has been published. S. M. Hosseini and J. J. Feng (2009), A particle-based model for simulating erythrocyte de- formation and transport through capillaries. Chemical Engineering Science, 64: 4488-4497. In close collaboration with J. J. Feng, I de- veloped the numerical model and perform simulations and drafted the paper. J. J. Feng helped me to prepare the final version of paper. • A version of chapter 4 has been published. S. M. Hosseini and J. J. Feng (2011), Pressure boundary conditions for computing incompress- ible flows with SPH. Journal of Computational Physics. 230: 7473- 7487. In close collaboration with J. J. Feng, I developed the numerical algorithm and carried out the simulations and drafted the paper. J. J. Feng helped me to prepare the final version of paper. • A version of chapter 6 has been submitted for publication. S. M. Hosseini and J. J. Feng (2012), How malaria merozoites reduce the deformability of infected RBC. I developed the particle based model for RBC to three-dimensional space and perform simulation for different stages of malaria infection under the supervision of J. J. Feng. The final version of the paper was prepared with help of J. J. Feng. iiiTable of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Thesis objectives . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Contribution of current study and its importance . . 3 1.3 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Human red blood cells . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Structure of RBC membrane . . . . . . . . . . . . . . 7 2.1.2 Lifespan and function of RBC . . . . . . . . . . . . . 8 2.1.3 Deformability of RBC . . . . . . . . . . . . . . . . . . 8 2.1.4 Measurement of cell deformability . . . . . . . . . . . 9 2.2 Human malaria . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1 Public health impact of malaria . . . . . . . . . . . . 13 2.2.2 The life cycle of malaria parasite . . . . . . . . . . . . 13 2.3 Malaria infected red blood cells (iRBCs) . . . . . . . . . . . 14 2.3.1 Molecular and structural changes of iRBCs . . . . . . 16 2.3.2 Morphology and mechanical properties of iRBCs . . . 17 2.3.3 Deformability loss in malaria iRBC . . . . . . . . . . 18 2.3.4 Rheological changes in blood containing iRBCs . . . 20 ivTable of Contents 3 Research methodology . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 Modelling cell mechanics . . . . . . . . . . . . . . . . . . . . 22 3.2 Continuum models . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3 Discrete models . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.1 Spectrin network models . . . . . . . . . . . . . . . . 24 3.3.2 Mesoscopic scale models . . . . . . . . . . . . . . . . 25 3.3.3 DPD for RBC modeling . . . . . . . . . . . . . . . . . 26 3.3.4 SPH for RBC modeling . . . . . . . . . . . . . . . . . 27 3.4 Smoothed particle hydrodynamics (SPH) . . . . . . . . . . . 28 3.4.1 Fundamentals . . . . . . . . . . . . . . . . . . . . . . 29 3.4.2 Kernel approximation of functions and their gradients 29 3.4.3 Particle approximation . . . . . . . . . . . . . . . . . 32 3.4.4 SPH for simulating fluid flow . . . . . . . . . . . . . . 33 3.4.5 Weakly compressible SPH (WCSPH) method . . . . . 34 3.4.6 Truly incompressible SPH algorithm . . . . . . . . . . 35 3.4.7 SPH approximation of the flow equations . . . . . . . 36 3.4.8 Limitations of SPH for open-boundary flows . . . . . 38 4 Correcting the pressure boundary conditions of SPH . . . 41 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 Projection schemes . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.1 Homogeneous projection scheme . . . . . . . . . . . 42 4.2.2 Rotational projection scheme . . . . . . . . . . . . . . 43 4.2.3 Natural boundary condition . . . . . . . . . . . . . . 44 4.3 SPH implementation . . . . . . . . . . . . . . . . . . . . . . 45 4.3.1 SPH interpolation . . . . . . . . . . . . . . . . . . . . 45 4.3.2 Solution algorithm . . . . . . . . . . . . . . . . . . . . 47 4.3.3 Boundary conditions . . . . . . . . . . . . . . . . . . 48 4.3.4 Tensile instability . . . . . . . . . . . . . . . . . . . . 51 4.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . 53 4.4.1 Poiseuille flow driven by body force in periodic domain 53 4.4.2 Fully developed Poiseuille flow with open boundaries 53 4.4.3 Channel flow past a square cylinder . . . . . . . . . . 56 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5 A 2D particle-based model for RBC . . . . . . . . . . . . . . 70 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.2 Physical model and numerical scheme . . . . . . . . . . . . . 71 5.2.1 Constitutive models for the membrane . . . . . . . . 72 5.2.2 Governing equations for fluid flow . . . . . . . . . . . 74 vTable of Contents 5.2.3 Discretization using SPH . . . . . . . . . . . . . . . . 75 5.2.4 Solution algorithm . . . . . . . . . . . . . . . . . . . . 76 5.2.5 Bending moment . . . . . . . . . . . . . . . . . . . . 78 5.3 Results and analysis . . . . . . . . . . . . . . . . . . . . . . . 79 5.3.1 Cell in shear flows . . . . . . . . . . . . . . . . . . . . 80 5.3.2 Cell in Poiseuille flows . . . . . . . . . . . . . . . . . 84 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6 How malaria merozoites reduce the deformability of iRBC 89 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.2 Physical model and numerical scheme . . . . . . . . . . . . . 90 6.2.1 Elasticity of the RBC membrane . . . . . . . . . . . . 92 6.2.2 Numerical method . . . . . . . . . . . . . . . . . . . . 95 6.3 Results and analysis . . . . . . . . . . . . . . . . . . . . . . . 98 6.3.1 Stretching of healthy RBC . . . . . . . . . . . . . . . 99 6.3.2 Stretching of RBC in the ring stage . . . . . . . . . . 100 6.3.3 Stretching of RBC in the trophozoite stage . . . . . . 103 6.3.4 Stretching of RBC in the schizont stage . . . . . . . . 106 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7 Conclusions and recommendations . . . . . . . . . . . . . . . 112 7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.1.1 Improving SPH numerical algorithm . . . . . . . . . . 112 7.1.2 Presenting a particles based model for RBC . . . . . 113 7.1.3 Studying the malaria iRBC by our numerical model . 113 7.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 viList of Figures 2.1 Schematic structure of cell membrane. Reprinted from Ref. [23] with permission from Elsevier. . . . . . . . . . . . . . . . . . . . 7 2.2 Schematic model of the red cell cytoskeleton. Reprinted from Ref. [154] with permission from John Wiley and Sons. . . . . . . . . . 8 2.3 Snapshots of the healthy RBC [94] (a) in the rest state and (b) passing through splenic cords to the splenic sinus. . . . . . . . . . 9 2.4 Micropipette aspiration experiment for red blood cell. Reprinted from Ref. [62] with permission from Elsevier. . . . . . . . . . . . 10 2.5 The images of deformed RBC by different stretching forces applied by two silica beads in diametrical position. Reprinted from Ref. [59] with permission from Elsevier. . . . . . . . . . . . . . . . . . 12 2.6 A schematic view of the life cycle of the malaria parasite. Reprinted from Ref. [51] with permission from Nature Publishing Group. . . 15 2.7 A schematic view of merozoite containing major organelles and cel- lular structures. Reprinted from Ref. [18] with permission from Elsevier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.8 Transmission, confocal, and 3D images of a trophozoite-stage IRBC (a), and schizont-stage IRBC (b and c). Reprinted from Ref. [35] with permission from Elsevier. . . . . . . . . . . . . . . . . . . . 18 2.9 The micropipette aspiration experiment shows progressive loss of deformability according to disease progression. Reprinted from Ref. [86] with permission from Elsevier. . . . . . . . . . . . . . . . . . 19 2.10 Optical images of healthy (H-RBC), ring stage (Pf-R-pRBC), tropho- zoite stage (Pf-T-pRBC) and schizont stage (Pf-S-pRBC) in rest state (left column) after tensile stretching by optical tweezers at a constant force of 68 ± 12 pN (middle column) and at a constant force of 151 ± 20 pN (right column). Reprinted from Ref. [142] with permission from Elsevier. . . . . . . . . . . . . . . . . . . . 20 3.1 Schematic representation of the compact support of a kernel func- tion that determines the finite interaction length between particles. 31 viiList of Figures 4.1 Computational domain for flow past a square cylinder of side D in a channel of width H , showing the three types of boundaries to be encountered in the simulation. . . . . . . . . . . . . . . . . . . . 49 4.2 (a) Tensile instability caused by the particles forming anisotropic strings in flow around a square cylinder. (b) The particle shift- ing approach [160] largely removes the undesirable clustering and maintains a more or less uniform particle distribution. . . . . . . 52 4.3 Development of the Poiseuille velocity profile in time: comparison between SPH simulation and the analytical solution. Two numer- ical solutions are shown, corresponding to initial particle spacing L0 = 0.05H and 0.1H . Time is made dimensionless by H/U0. . . . 54 4.4 Deviation of the computed pressure field p(r) from the analyti- cal solution pA(r) (i.e., a linear distribution along the flow direc- tion). (a) Homogeneous projection scheme. (b) Rotational pro- jection scheme. The pressure error has been scaled by 12ρU20 and zeroed at the center of the domain. . . . . . . . . . . . . . . . . . 55 4.5 Temporal evolution of the drag coefficient CD and its components due to pressure and viscous friction. The spatial resolution is at L0 = H/80 with 21300 particles. The drag coefficient is scaled by the finite-element benchmark value CFED and time is made dimen- sionless by D/U0. . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.6 Drag coefficient CD computed by the rotational and homogeneous projection schemes at five particle numbers. The data are shown as fractional deviations from the benchmark FE drag coefficient CFED . 60 4.7 (a) The centerline velocity profile u(x) computed by finite elements. (b) Deviation of the SPH velocity profiles from the FE profile com- puted using increasing numbers of particles. . . . . . . . . . . . . 61 4.8 (a) The centerline pressure profile p(x) computed by finite ele- ments. (b) Deviation of the SPH pressure profiles from the FE profile computed using increasing number of particles. Pressure has been scaled by 12ρU 2 0 and zeroed at the exit of the channel (x = 12D). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.9 Comparison between (a) velocity profiles and (b) pressure pro- files along the centerline computed by homogenous and rotational schemes. We plot the deviation of the SPH results from the FE benchmarks in Figs. 4.7(a) and 4.8(a). . . . . . . . . . . . . . . . 64 4.10 Pressure profile along the front (a) and the back (b) of the square cylinder computed by finite elements and the homogenous and ro- tational schemes. In all three computations p is zeroed at the exit. 66 viiiList of Figures 4.11 Drag coefficient CD on the square cylinder as a function of the Reynolds number Re. The SPH computation with the rotational scheme uses 68000 particles. The lattice-Boltzmann results of Breuer et al. [11] are shown for comparison. . . . . . . . . . . . . . . . . 68 5.1 (a) Particle-based model for the red blood cell, with fluid parti- cles representing the cytoplasm and the suspending plasm, and spring-connected solid particles representing the cell membrane. (b) Schematic of extensional and bending elasticity between par- ticles in the membrane. . . . . . . . . . . . . . . . . . . . . . . . 71 5.2 Element bending groups (EBGs) used to convert the bending mo- ment to pairs of forces acting on particles in the membrane. The line and arrow styles distinguish pairs of forces used to replace the bending moment on each line segment. . . . . . . . . . . . . . . . 78 5.3 Deformation of an RBC in shear flow. C = 2 × 104, G = 0.234 and ˆEB = 2.63× 10−3. The snapshots are for dimensionless times kt = 0, 2, 3, 4 and 10, and the last frame depicts the flow field surrounding the cell in the steady state. . . . . . . . . . . . . . . 81 5.4 Effect of increasing the bending elasticity on cell deformation. ˆEB = 1.31× 10−2 (first row), 2.62 × 10−2 (second row) and 5.26× 10−2 (third row) correspond to a bending elasticity 5, 10 and 20 times that of the real RBC. The three columns are for kt = 3, 5 and 7. G and C have the same value as in Fig. 5.3. . . . . . . . . . . . . 82 5.5 The steady shape of the RBC gets more elongated with increasing G values. All other dimensionless parameters have the same value as in Fig. 5.3. Note the quantitative agreement with the numerical result of [167] at G = 0.234, who used the immersed boundary method and a continuum model for the cell. . . . . . . . . . . . . 83 5.6 Deformation of a red blood cell in a Poiseuille flow. The flow is from left to right, and the five snapshots are taken at dimensionless time Umt/a = 0, 1.28, 3.85, 6.40 and 8.95, the last approaching the steady state. Despite the considerable change in shape, the circumference of the cell has increased a mere 1.05% relative to the undeformed state. . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.7 Comparison of the steady-state cell shape in a Poiseuille flow be- tween our simulation (black dots) and the result of [127] (continuous curve). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.8 Effects of (a) G and (b) ˆEB on the steady state deformation of the RBC in a Poiseuille flow. . . . . . . . . . . . . . . . . . . . . . . 86 ixList of Figures 5.9 Deformation of a red blood cell as it goes through a contraction in the channel. The flow is from left to right. The width ratio be- tween the narrow and wide channels is 0.4. Defined using the mean velocity in the wide channel, all flow and membrane parameters are the same as in Fig. 5.6. The snapshots are at dimensionless time Umt/a = 0, 0.641, 1.60, 1.92, 2.24 and 2.88. . . . . . . . . . . . . 87 6.1 Discretization of a healthy RBC by SPH particles. The membrane is represented by elastic springs connecting neighboring particles into a more or less uniform network. The total number of particles in this case is 13500. The RBC is centered at (0, 0, 0) and the coor- dinates are in unit of meters. The arrows represent the stretching force applied by optical tweezers. . . . . . . . . . . . . . . . . . . 91 6.2 Convergence test: elongation of the major axis of a model RBC as a function of the total number of SPH particles. The RBC contains no merozoites inside, and has a membrane shear modulus Gs = 15 µN/m and a bending modulus kb = 6 × 10−19 J. The stretching force F = 100 pN. . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.3 Variation of the longitudinal and transverse diameters of the healthy RBC as a function of the stretching force. The results of our nu- merical simulation are compared with experimental data [142] and prior computations of Imai et al. [71] and Fedosov et al. [41]. . . . 100 6.4 Three of the configurations of the merozoite used in our simulations. The top row shows the top view, and the bottom the corresponding side view. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.5 Variation of the axial and transverse diameters of the RBC as a function of the stretching force in ring stage. The solid curves are computation at Gs = 15µN/m without any merozoites. The error bar on the computed result at 150 pN indicates variations due to inclusion of the merozoite in the different configurations of Fig. 6.4 The experimental data of [142] are also plotted for comparison. . . 102 6.6 (a) Micrograph of an iRBC in trophozoite stage, with the mero- zoite outlined in translucent brown and the hemozoin in purple. Reprinted from Ref. [57] with permission from Elsevier. (b) Ax- isymmetric model for the iRBC used in the computations, with a spherical merozoite located in the thicker part of the cell. The lengths are in µm, and the parasite is centered at Y = −1.65µm. . 103 xList of Figures 6.7 (a) Deformation surface of trophozoite-stage iRBC as a function of membrane elastic modulus Gs and the merozoite radius Rm. The horizontal plane is at 39% deformation, the experimentally observed amount. (b) The intersection of the 39%-deformation plane and the deformation surface, plotted on the (Gs, Rm) plane. . . . . . . . . 105 6.8 (a) Scanning electron micrograph of an IRBC in the schizont stage. K indicates a knob. Reprinted from Ref. [150] with permission from Elsevier. (b) An oblate spheroid model for the schizont-infected iRBC, with length in µm. The parasite is represented by an oblate spheroid of the same aspect ratio as the whole cell. . . . . . . . . 106 6.9 Closely packed daughter cells of P. falciparum in the schizont stage. Reprinted from Ref. [58] with permission from Elsevier. . . . . . 107 6.10 (a) Deformation surface of schizont-stage iRBC as a function of membrane elastic modulus Gs and the and merozoite-to-iRBC vol- ume ratio v. The horizontal plane is at 22% deformation, the ex- perimentally observed amount. (b) The intersection of the 22%- deformation plane and the deformation surface, plotted on the (Gs, v) plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 xiAcknowledgements First and foremost, I would like to express my sincere gratitude to my out- standing supervisor, Professor James J. Feng for his invaluable support dur- ing my PhD studies. Working with Professor Feng was a priceless oppor- tunity for me that allowed me to learn a lot about the scientific approach in doing a solid research. Most importantly I learnt that it is impossible to build a solid academic career without academic integrity and ethics. His encouragements when I got disappointed were a magnificent help to tackle the difficulties. Being so serious and criticizing during the scientific discus- sions while being so friendly after the meeting defined a perfect supervisor in my mind. Last but not least, he showed me what it means to be a truly scientific character. I want to thank Dr. Peng Gao and Dr. Giovanni Ghigliotti who helped me a lot during the progress of my research. They patiently spent lots of time to discuss the mathematical modeling in my projects. I am indebted to my lovely wife, Sara Saberi, for her exclusive support and continuous encouragement and love which undoubtedly was one of the main reasons of every success that I had in my grad life in UBC. xiiDedication This thesis is dedicated to my father, Hashem Hosseini, who taught me being patient with difficulties and enjoy from each single moment of my life. It is also dedicated to my mother, Zahra Etedali, who supports me with her invaluable love in each single step of my life. I also would like to express my sincere gratitude to my lovely wife, Sara Saberi, who has always been a great source of love to me. She was a truly wise and supportive friend throughout the tough days of my grad life in UBC. xiiiChapter 1 Introduction 1.1 Motivation The pathogenesis of some severe human diseases such as malaria and cancer is crucially dependent on changing the mechanical properties of infected cells [79, 142]. In cancer, some genetic mutations modify the underlying micro structures of the cancerous cells and drastically change their mechanical properties. The experimental studies measured the Young’s modulus of cancerous cells to be about one-tenth of healthy cells [79]. As a result, the cancerous cell are very flexible and capable of migrating through size- limited pores in the basal membrane and endothelium during metastasis. In malaria, the parasite derived proteins initiate some biochemical reactions that decrease the deformability of red blood cell (RBC) by one order of magnitude. Some researchers estimated an increase in the elastic modulus of the infected RBC (iRBC) by more than a factor of 10 [141]. The infected RBCs are too stiff to squeeze through tiny capillaries to deliver oxygen to different tissues of body as their normal task. Instead they will block the capillaries, and may cause coma and death. Studying the exact mechanical changes of diseased cells may lead to understanding the dominant mechanism of disease progression, which is extremely important for designing new drugs. In addition, knowledge of the mechanical properties of infected cells helps in designing new diagnostic devices and techniques to detect diseases at an early stage. Here we specifically focused on studying the pathogenicity of Plasmod- ium falciparum within the human RBC. Malaria is a global health issue that kills more than 781,000 people and affecting more than 216 million people per year according to the World Health Organization’s 2011 Report. Malaria is usually caused by being bitten by an infective female Anopheles mosquito. Infecting sporozoites migrate to the liver right after entering the bloodstream from the mosquito’s saliva. After multiplying to thousands of merozoites in the infected liver they will escape to the bloodstream to target the red blood cells. The erythrocytic phase of malaria parasite starts by the merozoite’s penetrating an RBC. The malaria parasite exports proteins to 11.2. Thesis objectives the host cell membrane to modify its properties. In addition, the parasite undergoes growth in size and multiplies to great numbers, swelling the host cell. Eventually the RBC membrane ruptures some 48 hours after invasion to release the merozoites into the blood to infect new RBCs. There are many unanswered question about the disease progression dur- ing the whole infection process, starting from invasion of RBC to its final destruction. The mechanism for the progressive loss of deformability is one of the major question due to its importance in disease progression. For ex- ample, the current understanding attributes the loss of RBC deformability to a 10-fold increase in elastic stiffness of membrane, which is in turn due to extra cross-linking among proteins in the cytoskeleton. But certain mi- cropipette aspiration tests reported less than 3-fold increase in shear elastic modulus of membrane. An explanation of this discrepancy will shed light on the dominant factor causing rigidification of infected RBCs, which may ultimately be useful for designing new anti-malaria drugs. This example highlights the need to develop a rational and quantitative understanding of evolving cellular mechanics during disease progression, which has motivated the computational efforts undertaken in this thesis. 1.2 Thesis objectives The objective of this thesis is to develop a particle-based model for cell me- chanics and use it to explore the interplay between microstructural remodel- ing inside the cell and its mechanical behavior as a whole. Malaria infected RBC is one of the most prominent examples of such structural-property- disease coupling, a dynamic relationship that tends to evolve during disease progression. The pathogenicity of Plasmodium falciparum is largely owing to its altering the mechanical properties of iRBCs. The biochemical reac- tions initiated by parasite-derived proteins change the architectures of some microstructures inside the RBC, which can drastically change the overall be- havior of iRBC. Representing the RBC with Lagrangian particles may give us the unique capability to explore the contribution of different underlying microstructures on the overall behavior of the cell. In such models a clus- ters of particles having common properties that differ from the surroundings may represent various organelles, and the interaction among different parti- cles may be dynamically evolved to reflect the dominant factor in changing the macro-mechanical properties of RBC in disease state. In this PhD dissertation, I focus on malaria iRBC and develop a particle based model to simulate the mechanical changes in different disease stages. 21.2. Thesis objectives It is well known that the severe symptoms of malaria disease are mainly due to stiffening the infected RBCs. The current understanding ascribes the loss of RBC deformability to a 10-fold increase in elastic stiffness of membrane due to extra cross-linking among proteins in the cytoskeleton. Meanwhile micropipette aspiration tests reported less than 3-fold increase in shear elastic modulus of membrane. We believe the discrepancy is caused by the presence of the rigid parasite particles inside infected cells, and have carried out numerical simulations to demonstrate this mechanism. 1.2.1 Contribution of current study and its importance The major contributions of the current PhD thesis can be summarized as follows: • Introducing the “rotational pressure-correction scheme” to the stan- dard incompressible SPH method. We proposed to use a consistent pressure boundary condition that relates the normal pressure gradient to the local vorticity. We demonstrated that the previous method of enforcing incompressibility by a projection method with an artificial homogeneous Neumann boundary condition for the pressure Poisson equation produces a numerical boundary layer that compromises the solution near boundaries. Difficulties related to inconsistent pressure boundary condition at solid walls and inflow and outflow boundaries have previously kept a large class of flow problems inaccessible to SPH simulations. We showed that our proposed scheme computes the pres- sure and velocity accurately near open boundaries and solid objects, and extends the scope of SPH simulation beyond the usual periodic boundary conditions. • Developing a two-dimensional (2D) particle-based model for the red blood cell, and using it to compute cell deformation in simple shear and pressure-driven flows. In this model the cell membrane was replaced by a set of discrete particles connected by nonlinear springs; the spring law enforces conservation of the membrane area to a high accuracy. In addition, a linear bending elasticity was implemented using the devia- tion of the local curvature from the innate curvature of the biconcave shape of a resting red blood cell. The cytoplasm and the external liq- uid were modeled as homogeneous Newtonian fluids, and discretized by particles as in standard smoothed-particle-hydrodynamics (SPH) solution of the Navier–Stokes equations. We employed the discrete 31.3. Thesis outline particles not only as a numerical device for solving the partial differ- ential equations, but also as a means for incorporating microscopic physics into the model. In other words, the particles are a numerical device for solving the partial differential equations as well as a vehicle for incorporating microscopic physics. Numerically, the fluid flow and membrane deformation were computed, via the particle motion, by a two-step explicit scheme. The model parameters were determined from experimental measurements of cell viscosity and elastic moduli for shear, areal dilatation and bending. • Extending the particle-based model for the red blood cell to three- dimensional (3D) space to explore the parasite-driven changes in in- fected RBCs. The particle-based model was able to account for the contributions of the cell membrane elasticity, cytosol viscosity and par- asite shape and volume on the dynamic behavior of iRBC. Our discrete model sheds light on the discrepancy in literature data on the increased shear elastic modulus of RBC membrane in mature stages of infection among different studies. We showed that the discrepancy is caused by the presence of the rigid parasite particles inside infected cells, and have carried out numerical simulations to demonstrate this mecha- nism. The cell membrane is represented by a set of discrete particles connected by linear springs that represent shear and bending elasticity. The cytosol is modeled as a homogeneous Newtonian fluid, and dis- cretized by particles as in standard smoothed particle hydrodynamics models. The malaria parasite is modeled as an aggregate of particles constrained to rigid-body motion. Our 3D simulation of RBC stretch- ing tests by optical tweezers demonstrates a “compensating effects” between the existence of malaria parasites and the elevated stiffness of the membrane on the overall deformability of infected RBC. The particle-based model can be easily extended to other types of living cells. 1.3 Thesis outline Chapter 2 presents a general background on the RBC structure and its functionality, the life cycle of P. falciparum, and existing knowledge on the gradual loss of deformability in malaria-infected red cells. Chapter 3 presents details of the research methodology, including the rationale behind choosing the discrete model over continuum models, and a comparison among popular discrete models. Finally we give a detailed 41.3. Thesis outline discussion of the SPH method as a key element in both our physical model and our numerical scheme. Chapter 4 presents the details of our proposed pressure boundary condi- tions for computing incompressible flows with SPH. We introduce the “ro- tational pressure-correction scheme” into the standard incompressible SPH method. The proposed scheme was tested by simulating the channel flow past a square cylinder and the results showed markedly improvement in the accuracy of the solution. Our numerical experimentation has shown the ro- tational projection scheme as robust, accurate and efficient in dealing with open boundaries and flow around solid obstacles. Chapter 5 presents the details of our 2D particle-based model for the red blood cell. The model predicted an elongated shape for RBC in a sim- ple shear flow with the membrane and cytoplasm undergoing tank-treading motion. In a Poiseuille flow, we predicted a parachute shape for the RBC. The model predictions of tank-treading in shear and the parachute shape in Poiseuille flows are in excellent agreement with experimental data and prior continuum-based computations. Chapter 6 presents the details of our 3D particle-based model for ex- ploring the changes in the malaria infected RBC. The infection process is usually divided into three stages: ring, trophozoite and schizont. For each stage we adopted a different morphology for the iRBC and the merozoite according to experimental observations. The 3D model was used to sim- ulate the RBC stretching tests by optical tweezers in different stages of disease. Our numerical simulations demonstrated that the presence of the rigid parasite particles inside infected cells has a great effect in reducing the overall deformability of the iRBC. Previous schemes that attributed the observed cell rigidification entirely to membrane stiffening were flawed, and they overestimated the membrane modulus. Finally, Chapter 7 summarizes the thesis, outlines the limitations of the current work, and makes recommendations on future work. 5Chapter 2 Background 2.1 Human red blood cells Human red blood cell (RBC) is a bright red, biconcave, disk-shaped cell approximately 8 µm in diameter. RBCs account for about 40–45 percent of human blood by volume and are responsible for the blood’s red color. Their principal task is to take up oxygen (O2) in the lungs and deliver it to the body tissues by squeezing through the body’s capillaries. The circulatory system of blood flow transfers oxygen to different parts of body. The total number of red blood cells at any given time is roughly 2 ∼ 3×1013 (20-30 trillion), which is approximately one quarter of the total human body cell number for adult humans [116]. The red blood cells (RBCs) are also called erythroid cells or erythrocytes (from Greek erythros for “red” and kytos for “hollow”, with cyte translated as “cell” in modern usage). Human erythrocyte has a very simple structure, without a cell nucleus and organelles. That is the reason for its being so flexible and able to squeeze and pass through the microscopic capillaries. Around 2.4 million new erythrocytes are produced per second in the bone marrow [124]. The lifespan of normal RBC is around 120 days. Approximately a third of the total RBC volume comprises hemoglobin. The spectral properties of the hemic iron ions in hemoglobin produces the RBC’s red color. The role of hemoglobin is very crucial for transporting the oxygen by RBC. The healthy RBC in an undeformed state has a biconcave shape, 6 – 8 µm in diameter and 2 µm in thickness. It has a surface area around 136 µm2 and a volume around 90 µm3. If swollen up to a sphere shape while conserving the membrane surface area, the sphere’s volume would be 150 µm3. The fact that the RBC’s surface area is much larger than a sphere of equal volume helps it to easily deform under flow condition and external forces. It is well known that if RBCs are forced to go through capillaries with diameters in the range of 3 µm to 13 µm its shape will deform to a parachute or slipper-like one [109, 135]. 62.1. Human red blood cells Figure 2.1: Schematic structure of cell membrane. Reprinted from Ref. [23] with permission from Elsevier. 2.1.1 Structure of RBC membrane The RBC membrane comprises of a viscous and nearly area-incompressible lipid bilayer with an attached cytoskeleton. The elastic properties of the membrane can be entirely attributed to the cytoskeleton, which consists of a spectrin protein network linked by short actin filaments. Because of the lipid bilayer, the RBC membrane is highly conservative against areal expansion or contraction. In the meantime, it is very flexible with respect to shear defor- mation and bending, such that the RBC can deform readily to pass through narrow capillaries. Such properties are thanks to the multi-layered structure of the membrane. Figure 2.1, taken from [23], schematically depicts the main components of the RBC membrane, including the following components: the phospholipid bilayer, cholesterol molecules, transmembrane proteins and the underlying spectrin network. The cytoskeleton comprises a multi-protein complex formed by struc- tural proteins including α and β spectrin, ankyrin, protein 4.1 and actin. The mechanical properties and integrity of RBC are achieved by interaction between the membrane skeleton proteins and lipid bilayer and transmem- brane proteins. The side-to-side attachment of α and β spectrin forms flex- ible rod-like heterodimers. The tetramers will be formed by head-to-head self-association of the mentioned rod-like heterodimers. Ankyrin connects the tetramers to the cytoplasmic domain of the integral membrane protein band 3. Junctional complexes are the connection points between multiple spectrin tetramers with actin protofilaments, tropomyosin, tropomodulin and adducin. A schematic model showing the spatial relationship among these components is shown in Fig. 2.2 [154]. 72.1. Human red blood cells Figure 2.2: Schematic model of the red cell cytoskeleton. Reprinted from Ref. [154] with permission from John Wiley and Sons. 2.1.2 Lifespan and function of RBC The life cycle of RBC starts in the bone marrow. It takes about 7 days to develop a mature erythrocyte from a committed stem cell in a process which is called erythropoiesis. The mature RBC stays in blood circulation for about 100–120 days, before finally being trapped in the spleen, with its components recycled by macrophages. RBCs lose nucleus during their mat- uration which make them unable to perform protein synthesis and mitosis. In addition, RBCs can not repair any structural damage or metabolic fail- ures from synthesis of proteins or lipids as a result of lacking the nucleus and cellular organelles. But a cell without nucleus is extremely flexible and highly specialized for their main role of delivering the oxygen from the lungs to the tissues and then transporting carbon dioxide back from the tissues to the lungs. This task is accomplished by hemoglobin as an iron-containing biomolecule that can bind oxygen and carry it to the destination. 2.1.3 Deformability of RBC The deformability is a crucial property of RBC in transporting oxygen to various parts of body through micro-vessels. A healthy RBC can squeeze through microvasculature with cross sections one-third its own diameter by 82.1. Human red blood cells Figure 2.3: Snapshots of the healthy RBC [94] (a) in the rest state and (b) passing through splenic cords to the splenic sinus. remarkable deformation from its initial biconcave shape. The most drastic shape change of RBC happens during its transit in vivo from the splenic cords to the splenic sinus. As depicted in Fig. 2.3 [94], the RBC has to deform from its initial biconcave shape to a totally different one to pass through the spleen in each circulation. Such level of deformability is an important and distinguishing feature of RBC. With aging, the RBC gradually loses its ability to deform. After around 120 days the aged and less flexible RBC will be trapped in the spleen and removed from blood recirculation.Various diseases, including diabetes, sickle cell anemia and malaria, compromise the mechanical properties of the RBC, and cause it to lose it deformability. This gives great importance to the study of the mechanics of healthy and diseased cells, especially on the relationship between intracellular structures and deformability. 2.1.4 Measurement of cell deformability A number of experimental techniques have been developed to measure the deformability of RBC. Four of the most successful are briefly described be- low. Micropipette aspiration Evans et al. [40] proposed to use micropipette aspiration test for probing the deformability of RBC and extracting the elastic properties of membrane in 1976. This method is still being used for experimental studies of RBC 92.1. Human red blood cells Figure 2.4: Micropipette aspiration experiment for red blood cell. Reprinted from Ref. [62] with permission from Elsevier. [73]. In micropipette aspiration, the leading edge of cell surface is being recorded while its surface is aspirated into a small glass tube by a suction pressure. The displacement of the tongue is recorded as a function of the as- piration pressure. Using some elastic constitutive equation, one can match this function by theoretical or computational results, thus extracting the mechanical properties of the membrane. The displacement of the tongue can be tracked with a very high accuracy by a light microscope and the aspiration pressure can be as small as 0.1-0.2 pN/µm2 [62]. An image from micropipette aspiration of RBC is depicted in Fig. 2.4. A tongue with- drawal time can be obtained by pausing the suction pressure. The tongue withdrawal time (te), membrane viscosity (ηm), membrane elasticity (µ), viscosity of cytoplasm(ηi) and the erythrocyte thickness (d) are related to each other by following equation [105]: te = (ηm + ηid) µ (2.1) If the cell undergoes an elastic deformation, the tongue will be completely reintegrated into the cell body after release of the suction pressure. But a plastic deformation caused by a very large suction pressure will result in partial reintegration. Atomic force microscopy The deformability of RBC can also been probed by atomic force microscopy (AFM) which is a technique for quantitative assessment of microdeforma- tion. AFM can successfully combine high-resolution scanning with nano- indentation of AFM probe. In this method the AFM probe is vibrating over a very small piece of RBC membrane by a piezoelectric ceramic posi- tioner that deforms by applying a prescribed voltage. The AFM probe is 102.1. Human red blood cells a flat cantilever with a tip which come into contact with the sample while the cantilever behaves as a spring. The deflection of the probe can be ac- curately recorded by an optical system consisting of a laser diode and a position-sensitive detector. The tip movement can be tracked with a very high accuracy [105] and the local value of the Young’s modulus can be calcu- lated by measuring the cantilever deflection during the indenting the cell sur- face. Recently, the AFM method has been widely utilized for morphological study of living cells; it produces very high-resolution data of topographical structures of living cells. Optical tweezers Stretching cells by optical tweezers is a popular method for probing the deformability of the cells and extracting their elastic properties. Henon et al. [59] were among the first to use this method. In their experiments, two silica beads attached to opposite ends of the cell are trapped by laser beams and moved in opposite directions. The process is illustrated by Fig. 2.5. More recent experiments [24, 85, 93, 142] were done by adhering one bead to the glass surface and trapping another one with a laser beam. The cell is directly stretched by moving the trapped bead with the laser beam. A very small stretching force will impose a finite deformation on the axial and transverse diameter of RBC. This process will be repeated for several times with progressively increasing stretching force to obtain a correlation between deformed axial and transverse diameters and the corresponding stretching force. Then a mathematical model and computer simulation are needed to match the experimental data and extract the mechanical properties of the cell. Of course, the choice of constitutive laws and the accuracy of numerical solutions may affect the final calculated mechanical parameters. Such protocols have been applied to stretch cancer cells, healthy RBCs and iRBCs in the ring, trophozoite and schizont stages of malaria infection. The measured deformation is fitted to computations to yield an effective membrane modulus that varies with the progression of the disease. Microfluidics analysis of red blood cell The deformability of RBC has also been probed by microfluidic systems. One attractive feature of such devices is that they can represent to a degree the morphology and connectivity of the microvascular network. Another is their capability to test a large number of RBCs in flow fields under physi- ologically relevant flow conditions. The micro-channels are usually made of 112.1. Human red blood cells Figure 2.5: The images of deformed RBC by different stretching forces applied by two silica beads in diametrical position. Reprinted from Ref. [59] with permission from Elsevier. transparent substrata, such as polydimethylsiloxane (PDMS) and glass, so the RBCs behavior could be recorded by video microscopy and examined later to extract the mechanical parameters. Based on the size of the micro- channel the cell velocity and shape can be quantitatively described as a func- tion of the applied pressure drop. Tomaiuolo et al. [152] proposed a method to measure cell membrane viscoelastic properties in converging/diverging flow which are in good agreement with previous measured values in liter- ature. Recently, microchannels-based measurements are becoming popular for accurate estimating of areas and volumes of individual cells [47]. 122.2. Human malaria 2.2 Human malaria 2.2.1 Public health impact of malaria According to the World Health Organization’s 2011 World Malaria Re- port, more than 216 million cases of malaria were reported worldwide, with 781,000 deaths from malaria in 2010. Approximately 91% of reported malaria-related deaths occur in sub-Saharan Africa. More than 86% of malaria deaths globally happened to children under 5 years of age. Al- though the incidents of malaria globally has declined by 17% since 2000, it is still affecting a very large number of people and thus constitutes a global health issue. Malaria happens commonly in countries with a high percent- age of population living in poverty, and unfortunately it can itself be a cause of poverty and a major hindrance to economic development in such coun- tries. Many species of malaria parasites have been discovered so far, which include Plasmodium falciparum, Plasmodium malariae, Plasmodium ovale, Plasmodium vivax and Plasmodium knowlesi [104]. Among these, Plasmod- ium falciparum causes the most severe form of malaria and is responsible for about 90% of reported malaria-related deaths [139]. Thus, P. falciparum is also among the most thoroughly studied of this genus of parasites. 2.2.2 The life cycle of malaria parasite Malaria is usually caused by being bitten by an infective female Anophe- les mosquito. The female mosquitoes of the Anopheles genus are acting as transmission vectors and get infected by feeding on an infected human car- rier in the first place. Getting a blood meal from an infected human will transmit the malaria infection from human being to the mosquito. Further differentiation of parasite gametocytes into male or female gametes will hap- pen in the mosquito’s body and then fuse in the mosquito’s gut. Eventually they will turn into oocysts in the gut wall which releases sporozoites upon ruptures. The sporozoites migrate to the salivary glands of the mosquito, and become ready to transmit to humans by skin biting [145]. Development of malaria parasite in the human body comprises an exoerythrocytic phase and an erythrocytic phase (Fig. 2.6). In the exoerythrocytic phase, the in- fecting sporozoites migrate to the liver right after entering the bloodstream from the mosquito’s saliva. They will multiply to thousands of merozoites in the infected liver during a period of 8 ∼ 30 days. They then escape to the bloodstream to target the red blood cells. Malaria merozoites pen- etrate into RBC to start the erythrocytic phase, which involves infection of the erythrocytes [7]. An invading merozoite undergoes 4 or 5 rounds of 132.3. Malaria infected red blood cells (iRBCs) mitotic DNA replication which results in 16-20 daughter merozoites by the end of the erythrocytic phase. Approximately 48 hours after invasion, the RBC membrane will become destabilized and eventually rupture. Then the merozoites break out of the cell and begin to invade other healthy RBCs. It is very hard for the immune system to detect the malaria parasites in the body. Because in the exoerythrocytic phase and the erythrocytic phase malaria parasites reside in the liver and inside RBCs, respectively, they are protected from being destroyed by white blood cells. However, as will be explained shortly, the infected RBCs are far less deformable than healthy RBCs, and this enables the spleen to detect and filter them from bloodstream. To avoid the risk of being trapped in the spleen, merozoites export adhesive proteins onto the surface of the infected RBCs, which cause the infected RBCs to stick to the walls of small blood vessels [14]. 2.3 Malaria infected red blood cells (iRBCs) Merozoites released from the liver look like lemons that are approximately 1.5 µm long and 1 µm wide as schematically shown in Fig. 2.7. These in- vading merozoites are equipped with complex organelles including a nucleus, an endoplasmic reticulum, a rudimentary Golgi, and a single mitochondrion [58]. There are many unanswered questions about the process of merozoite attachment and invasion of the RBC. It is clear that merozoites use the api- complexan invasion organelles (apical complex, pellicle and surface coat) to recognize and enter the host erythrocyte. Merozoite hits the RBC surface in a random orientation which is followed by reorientation of merozoite such that the apical complex is in proximity to the erythrocyte membrane. They tightly attached to each other in the junction point and the merozoite starts to penetrate inside the RBC. The parasite forms a parasitophorous vesicle or parasitophorous vacuole (PV) while entering the host RBC to create a controlled environment for its development inside the erythrocyte [18]. It is well known that the RBC does not have phagocytic capability. In fact, the 2D submembrane cytoskeleton in the RBC prevents endocytosis. Therefore, the formation of the parasitophorous vacuole is solely coming from the parasite. Many proteins on the surface of the merozoite have un- known functionality, but some of them are well-known motor proteins that are believed to be responsible for driving the invasion process. For instance, discovery of an actin capping and a depolymerising agent (cytochalasin B) on the parasite surface suggested an actomyosin motor [117]. The erythro- cytic development of malaria parasite happens within the PV formed during 142.3. Malaria infected red blood cells (iRBCs) the invasion step. The parasite exports various proteins to the host RBC to optimize the chemistry of the environment for its development. So far more than 400 proteins have been detected including kinases, lipases, adhesins, proteases and chaperone-like proteins. These proteins produce profound Figure 2.6: A schematic view of the life cycle of the malaria parasite. Reprinted from Ref. [51] with permission from Nature Publishing Group. 152.3. Malaria infected red blood cells (iRBCs) Figure 2.7: A schematic view of merozoite containing major organelles and cellular structures. Reprinted from Ref. [18] with permission from Elsevier. remodeling of the infected erythrocyte, and they are believed to play ma- jor roles in decreasing the RBC deformability while greatly elevating the permeability of the RBC membrane to low-molecular-weight solutes. This increased transport across the membrane is responsible for altering the ionic composition of the host cell cytoplasm [74]. 2.3.1 Molecular and structural changes of iRBCs Malaria is caused by mosquito-borne parasites of the genus Plasmodium, of which Plasmodium falciparum is a commonly studied species. Malaria merozoites roughly 1 µm in size penetrate healthy red blood cells (RBCs) and initiate a series of biochemical and morphological changes. A schematic view of erythrocytic phase of malaria parasite is shown in Fig. 2.6. The whole infection process of RBC can be divided into ring, trophozoite and schizont stages [79]. • Ring stage starts about 30 min after merozoite invasion of the RBC and lasts about 20 hours. The entry of malaria merozoite is facilitated by its ovoid shape and pointed tip, but shortly after invasion it flattens into a thin disc. The parasite has a diameter of 2–3 µm and a thickness 162.3. Malaria infected red blood cells (iRBCs) of 0.5 µm by the mid-ring stage, which is 10–15 hours after invasion [58]. • Trophozoite stage starts about 20 hours after invasion. In the trophozoite stage, the merozoites export parasite proteins to the sur- face membrane of the RBC that cause crosslinking of the spectrin net- work and stiffening of the RBC membrane. As the parasite matures it continues to modify the host cell to both facilitate access to nutrients and to protect the infected RBC from being removed by spleen. To achieve the second goal, the plasma membrane of iRBCs becomes dis- torted by structures known as knobs, which serve to promote adhesion of the infected RBC to the vascular endothelium [58]. Gruenberg et al. [54] showed that the knob size and density on the RBC surface could be used to characterize different stages of parasite development in the infected RBC. • Schizont stage starts about 40 hours after invasion. At this late stage the parasite undergoes nuclear division to produce 12 to 20 mero- zoites, which pack themselves closely within the RBC, swell the host to a spherical shape, and further export parasite proteins to the RBC membrane. The overall deformability of the RBC, subject to stretch- ing by optical tweezers, decreases by more than a factor of 10 [142]. The rigidified RBCs can block capillaries and cause severe anemia, coma and even death. At the end of the schizont stage, about 48 hours after invasion, the RBC membrane bursts and the merozoites break out of the cell and begin to invade other healthy RBCs. 2.3.2 Morphology and mechanical properties of iRBCs It is well known that pathological changes in the mechanical properties of cells is the main cause of several human diseases [79, 142]. In malaria, the biochemical reactions initiated by invading parasites change the internal structure and mechanical behavior of living cells. For instance, the malaria infected RBCs will be so stiffened in the mature stages that they can no longer deform sufficiently to traverse tiny capillaries to deliver oxygen to various parts of the body. Instead they block the capillaries and disrupt the blood flow, which potentially could lead to coma and even death. Some researchers reported more than 10-fold increase in elastic modulus of infected RBCs [141]. Clearly in such cases the disease progression is facilitated by microstructural changes which results in changing the mechanical behavior of living cells as a whole. 172.3. Malaria infected red blood cells (iRBCs) Figure 2.8: Transmission, confocal, and 3D images of a trophozoite-stage IRBC (a), and schizont-stage IRBC (b and c). Reprinted from Ref. [35] with permission from Elsevier. As mentioned before, the exported proteins from parasites change the permeability of RBC membrane and thus allow modification of the cyto- plasmic content and volume. Meanwhile the mechanical parameters of the membrane, notably the various elastic moduli, change as well. These result in drastic morphological changes of the iRBC, especially in late trophozoite and schizont stage as shown in Fig. 2.8. By the end of the trophozoite stage, the parasite occupies more than one third of cell volume. Enlarging the parasite in schizont stage and multiplication to 16-20 daughter mero- zoites, along with influx through the now porous membrane, will swell the biconcave shape of RBC to a near spherical shape by the late schizont stage. Variation in the surface area to volume ratio is another factor in modifying the deformability of infected RBC. 2.3.3 Deformability loss in malaria iRBC As mentioned above, the progress of malaria disease is closely related to the the loss of deformability of iRBC. Several researchers have carried out experimental studies of malaria iRBC by different techniques. Shelby et al. 182.3. Malaria infected red blood cells (iRBCs) Figure 2.9: The micropipette aspiration experiment shows progressive loss of de- formability according to disease progression. Reprinted from Ref. [86] with per- mission from Elsevier. [132] performed a microfluidic study on the deformability of malaria iRBC under flow condition. They showed that a healthy RBC with diameter of 8 µm can squeeze into microchannels with a diameter as small as 2 µm. Meanwhile, the majority of trophozoites iRBCs cannot pass through 4 µm channels and the schizont iRBCs will easily block 6 µm channels. The loss of deformability of iRBC has also been probed by micropipette aspiration and cell stretching by optical tweezers as illustrated in Fig. 2.9 and Fig. 2.10, respectively. Nash et al. [107] performed micropipette as- piration test on the malaria iRBC in ring, trophozoite and schizont stage. They observed that ring stage iRBCs took about 50% longer to enter 3 microns pipettes in comparison with non-parasitised cells. Similar increase was observed in the critical pressure required to induce cell entry. The much higher loss of deformability in trophozoite and schizont stage result in four to six fold increase in entry time and critical required pressure. The shear modulus of the membrane was measured by pipette aspiration of small mem- brane tongues. They reported a 13% increase in the ring stage, with a slight stiffening of the iRBC membrane, and a much larger increase in the shear modulus for late schizont stage, from 4 – 7 µN/m (for healthy RBC) to 9 – 14 µN/m. 192.3. Malaria infected red blood cells (iRBCs) Figure 2.10: Optical images of healthy (H-RBC), ring stage (Pf-R-pRBC), tropho- zoite stage (Pf-T-pRBC) and schizont stage (Pf-S-pRBC) in rest state (left column) after tensile stretching by optical tweezers at a constant force of 68±12 pN (middle column) and at a constant force of 151 ± 20 pN (right column). Reprinted from Ref. [142] with permission from Elsevier. Suresh et al. [142] probed the altered deformability of malaria iRBC by optical tweezers test. The deformability curve was generated by record- ing the deformation of the axial and transverse diameters for a range of stretching force values. They built a finite-element mechanical model to extract the mechanical parameters from their experimental results. They estimated the elastic shear modulus of the membrane to be 5.3, 16, 21.3 and 53.3 µN/m for healthy RBC, iRBC in ring, trophozoite and schizont stages, respectively. Separate studies have given even higher estimations of the schizont-stage shear modulus, ranging from 89 µN/m [72] to 119 µN/m [170]. These shear-modulus values estimated from whole-cell deformation are far above the locally measured value of Nash et al. [107]. This discrep- ancy partly motivated our study. 2.3.4 Rheological changes in blood containing iRBCs Blood flows like a non-Newtonian fluid and increasing the shear stresses will result in decreasing of its effective viscosity. In severe malaria the infected cells tend to attach to the blood vessels due to knob-like structures on the RBC membrane produced by parasite-exported proteins. The anaemia case by attachment of RBCs to blood vessels turns into a drastic drop of hema- tocrit which can exponentially decrease the whole-blood viscosity in larger 202.3. Malaria infected red blood cells (iRBCs) vessels. At the same time the plasma fibrinogen levels will increase in se- vere malaria cases, which actually counterbalances the effects of anaemia. The two competing mechanisms may lead to only a moderate change in the measured whole-blood viscosity in severe malaria cases [27]. Blood flow properties in capillaries with diameters less than 5 µm are determined by the cellular level properties. In particular, the loss of RBC deformability in late intraerythrocytic development stages (trophozoites and schizonts) will cause a sharp increase in resistance to flow in small capillar- ies. The small branches of the circulatory system, such as high endothelial venules, can be completely blocked by the attachment of masses of rigid, sticky RBCs. Severe symptoms such as cerebral malaria and coma are caused by the blockage of these vessels. 21Chapter 3 Research methodology 3.1 Modelling cell mechanics Certain human diseases are closely associated with changes in the mechan- ical behavior of living cells [79, 142]. Cancer and malaria are two examples of such diseases, which are brought on, respectively, by internal factors (ge- netic mutation) and external factors (malaria parasite). In both diseases the biochemical reactions initiated by internal or external factors change the me- chanical properties of affected cells, and this modification greatly facilitates the disease progression. For instance, the ability of the cancerous cells to migrate through size-limited pores in the basal membrane and endothelium during metastasis is greatly increased by being much more deformable than the healthy cells. Experimental studies measured 10-fold decrease in the Young’s modulus of cancerous cells in comparison with healthy cells [79]. On the other hand, the effective elastic modulus of the red blood cell will increase by more than a factor of 10 after infection by malaria parasite [141]. The stiffened infected RBCs can no longer deform sufficiently to traverse narrow capillaries to deliver oxygen to various parts of the body. Instead they increase the blood flow resistance and block the micro capillaries and leading to coma and even death. Therefor studying the alteration in the mechanical behavior of infected cells can help us to understand the exact mechanism behind the disease progression. Modelling of cell mechanics is very useful for exploring the effects of underlying structure and organelles on the overall behavior of cell as a whole. After validation by in vitro experiments, mechanical models of living cells could be used for prediction of in vivo behavior of infected cells. Such important information provided by mechanical modeling will help us to design new diagnostic devices and techniques to detect diseases in an early stage. The mechanical models could be very useful in the fight against the diseases by simulating the effects of new medications on the behavior of infected cells in the body. Theoretical modeling has followed either the continuum approach or the microstructural approach to simulate the mechanical behavior of living cells 223.2. Continuum models [84]. The advantages and disadvantages of each approach are discussed in the following subsections. 3.2 Continuum models In the continuum approach the components of cell or the whole cell are treated as homogeneous materials which behave based on some constitutive equations. The constitutive equations and the relevant model parameters are basically the ones that best fit experimental data. The Newtonian liquid drop model is one of the simplest continuum models; it treats the whole cell as a liquid drop with a constant cytoplasmic viscosity and a constant surface tension [29, 163]. Later on the cell membrane elasticity was added to the model. More complex constitutive equations have also been proposed to account for shear-thinning and viscoelasticity. The recent continuum models simulate the cell mechanics by an elastic or viscoelastic membrane encapsulating the liquid cytoplasm. Still more sophisticated continuum models seek to reflect the multiphase nature of the cellular components, either through mixture-type “biphasic models” [60] or an explicitly two-phase treatment of the cell as a deformable capsule, with liquid cytoplasm enclosed by an elastic or viscoelastic membrane. The latter approach is related to fluid-structure interactions in mechanical systems, and thus is especially appealing to fluid dynamicists [32, 82, 113, 114, 122]. Being amenable to computational methods which are already developed for fluid and solid mechanics is one the most attractive aspects of contin- uum approach. However, They are not able to make any connection between the microstructural changes inside the cell and the behavior of the cell as a whole. One example is the activation of white blood cells, which may dras- tically alter the cells’ mechanical properties [161, 162]. Such changes cannot be faithfully represented by, say, the relaxation time of a viscoelastic model of leukocytes [28]. As mentioned before, understanding of such connection is valuable in understanding severe disease like cancer and malaria. 3.3 Discrete models Discrete models can be built based on the microstructures and micro-organelles inside the cells. The interaction between microstructures will be postulated so as to produce the overall behavior as seen in experimental studies. Dis- crete models have been developed to approximate the living cells at the spectrin molecular level [8, 80] or at the mesoscopic scale [30, 109, 118]. 233.3. Discrete models The most attractive feature of discrete models is their capability to ac- count for the interplay between microstructural remodeling inside the cell and its mechanical behavior as a whole. To explore the effects of some dis- eases, it is possible to make some changes in microstructures of the living cells and study the overall behavior of the infected cell as a whole. The models to be used in this thesis fall in this category. In the following we will first briefly review the mostly widely used discrete models, and then focus on one based on smoothed particle hydrodynamics (SPH), which is central to the present study. 3.3.1 Spectrin network models So far, the network models have been the most popular and successful among discrete cell-mechanical models. Stamenovic and Coughlin [140] modeled the cytoskeleton by a network of struts and cables that represent the mi- crotubules and actin filaments. Boey et al. [8] proposed a three-dimensional network model inspired by the spectrin network in erythrocytes. An elas- tically based network free energy may then be minimized to determine the equilibrium shape of the red cell, as well as its deformation subject to ex- ternal stretching [80]. Allowing the cytoskeletal network to undergo active polymerization, Herant et al. [60] incorporated the kinetics and network- membrane interaction through convection-diffusion equations. Here, the network is described as a pseudo-continuum having statistical properties. One obvious advantage of network models is the knowledge, to a certain degree of approximation, of the configuration of the intracellular compo- nents. If the network is discrete with a large but finite number of links, that physical discretization can be used directly as a numerical discretization for computation [10, 129]. In view of the mechanical behavior of the spectrin polymer, a wormlike chain model is adopted to describe its nonlinear elas- ticity. The network configuration can be more or less regular or randomized to mimic experimentally observed patterns. Bending elasticity is introduced based on a spontaneous curvature and the difference in orientation between neighboring elements of the network. On this level, Boey et al. [8] were able to simulation the micropipette aspiration of red cells. Li et al. [80] deter- mined the equilibrium shape of the red cell by minimizing the free energy, and also computed its deformation subject to external forces. Some network models tried to incorporate more intracellular details by allowing the cytoskeletal network to undergo active polymerization [60]. So the interaction between the kinetics and network-membrane is governed by convection-diffusion equations. Obviously the capability of network models 243.3. Discrete models in making connection between the physical discretization and the numerical discretization is a great advantage in modeling pathological changes on the cellular level [10, 129]. However, high computational cost is a major limi- tation of network models that include very detailed information of spectrin molecules. 3.3.2 Mesoscopic scale models Mesoscopic models coarse-grain the structural details of the intracellular components while striving to retain key elements or effects important for cellular-level mechanics. Several such models have been proposed for simulating the living cell and in particular RBCs [30, 31, 109, 118]. An earlier model in this category is that of Dzwinel et al. [31], which comprises a volume of elastic material with an inner skeleton. Their model did not include the RBC membrane and the internal cytosol which made it impossible to simulate the dynamics behavior of RBC in shear flow such as the tumbling and tank-treading. The next generation of mesoscopic methods treat the RBC membrane as a network of elastic springs which constraints the surface area. By adding the bending rigidity they were able to produce very promising results. The more sophisticated models account for coupling between the RBC membrane and the inside and outside fluids, as well as interaction among multiple cells in a suspension [30, 41]. Recently, particle-based representations of cells are being actively pur- sued as a new class of discrete microstructural models. Boryczko et al. [9] and Dzwinel et al. [31] used a rectangular lattice of particles connected by elastic springs to model the elastic behavior of red blood cell. The interaction between particles including conservative and dissipative forces among them can be postulated in different methods such as dissipative particle dynamics (DPD) method or the smoothed particle hydrodynamics (SPH) method or moving-particle semi-implicit (MPS) method. Pivkin and Karniadakis [118] successfully simulate the dynamic behavior of RBC by a coarse-grained RBC model which used Dissipative Particle Dynamics (DPD) [64]. Thanks to their meshless nature, Lagrangian particle-based methods such as SPH and MPS are particularly suitable for simulating the large RBC deformations in micro-capillaries [146, 155, 157]. So far the DPD and SPH methods have been widely used with promising results. In the following, their advantages and disadvantages will be discussed in more details. 253.3. Discrete models 3.3.3 DPD for RBC modeling The Dissipative Particle Dynamics (DPD) method was first introduced by Hoogerbrugge and Koelman [64] two decades ago for simulating microscopic hydrodynamic phenomena. They combined features from molecular dynam- ics and lattice-gas automata to present a novel particle-based method. They showed that DPD can solve the isothermal Navier-Stokes flow with relatively few particles, which makes it much faster than molecular dynamics and much more flexible than lattice-gas automata schemes. Using DPD for simulating the dynamic behavior of a coarse-grained RBC model was proposed by Pivkin and Karniadakis [118]. By coarse-graining the spectrin network model, they represented the RBC membrane by a 2D triangulated network of particles. Following the spectrin network models [8, 26, 81], they wrote the potential energy of the system as a sum of the in-plane elastic energy Vin−plane (stretching springs), the bending energy Vbending, the area conservation Varea and volume conservation constraints Vvolume: V(xi) = Vin−plane + Vbending + Varea + Vvolume. (3.1) The above expression of free energy will result in the following nodal forces: fi = −∂V(xi)/∂xi, i ∈ 1...N. (3.2) The in-plane free energy includes the spring energy Us and the elastic energy stored in the lipid membrane and other protein materials: Vin−plane = ∑ Us(lj) + ∑ Cq Aqk , (3.3) where lj is the length of the spring j, Ak is the area of the kth triangle, and the constant Cq and exponent q are model parameters to be properly selected to produce acceptable results. The bending energy was defined as Vbending = ∑ j∈1..Ns kb[1− cos(θj − θ0)], (3.4) where kb is the bending constant, θj is the instantaneous angle between two adjacent triangles having the common edge j, and θ0 is the spontaneous angle. It is worth mentioning that θ0 is not the spontaneous angle between two adjacent triangles that lie on the undeformed biconcave RBC. Rather it is the angle consistent with the curvature of the sphere to which the mem- brane would return if free of stress. Thus, it is a model parameter common 263.3. Discrete models to all vertices. However, since in a discrete representation the angle between adjacent triangles depends on the resolution of the tessellation, θ0 was set according to the total number of vertices for any specific triangulation of the RBC surface [43]: θ0 = cos−1 (√ 3(N − 2)− 5pi√ 3(N − 2)− 3pi ) . (3.5) The potential terms regarding the area and volume conservation con- straints were defined as follows: Varea = ka(A−Atot0 )2 2Atot0 + ∑ j∈1..Nt kd(Aj −A0)2 2A0 , (3.6) Vvolume = kv(V − V tot0 )2 2V tot0 , (3.7) where the model parameters ka, kd and kv were described as global area, local area, and volume constraint coefficients, respectively [42, 43]. Using the above models, Karniadakis and coworkers have simulated the stretching of red cells in different stages of malaria infection by optical tweez- ers [41, 118], single RBC dynamics in Poiseuille flows [118], cell adhesion to walls of blood vessels [41, 44], cell microrheology [43], cell-cell aggregation, structure formation and overall rheology in suspension of cells [41]. 3.3.4 SPH for RBC modeling SPH is a fully Lagrangian particle-based method well suited for comput- ing fluid-structure interaction with very large deformations. Inspired by the SPH method, Tanaka and Takano [146] and Tsubota et al. [155, 157] presented particle-based models for RBC including membrane particles and fluid particles. The elastic behavior of the membrane was represented by two sets of springs that produce resistance to both stretching and bending. The motion of fluid particles was governed by the Navier-Stocks equations following the classical SPH discretization. More recently, Imai et al. [71] and our group [68] have extended the SPH models to compute red cell motion and deformation in complex flow fields. Because the SPH models represent a coarse-grained version of the spec- trin network, the model parameters are usually seen as phenomenological. They are chosen such that the membrane as a whole exhibits desirable elas- tic behavior with prescribed moduli. This differs from the spectrin network 273.4. Smoothed particle hydrodynamics (SPH) models [8, 26, 81] that seek to account for actual spectrin properties and network configurations. In this regard, it is interesting to note that non- linear springs can be chosen in an SPH model so that the overall behavior approximates the Skalak model for a continuum membrane [4, 68]. The elas- tic reaction of a regular hexagonal network can be calculated analytically in the limit of small deformations [43, 112]. This way, the macroscopic elas- tic properties of the membrane, such as the shear, dilatational and Young’s moduli, can be related to parameters in the SPH model. Through work carried out in this thesis, we have come to appreciate the unique advantages of discrete particle models for cell mechanics. The particles offer far more freedom in modeling the internal structure of cells than continuum. For instance, various organelles can be represented by a group of particles having common properties that differ from the surrounding particles and their interaction may be dynamically evolved to reflect the remodeling of the intracellular structures triggered by internal or external factors. The work to be discussed in the rest of the thesis consists of the develop- ment of a SPH computational scheme, and applications of SPH in simulating 2D and 3D physical models of the red cell. Since SPH is at the foundation of all these, we devote the rest of this chapter to a review of the generalities of the method. 3.4 Smoothed particle hydrodynamics (SPH) Smoothed particle hydrodynamics (SPH) is a mesh-free Lagrangian method that discretizes the materials with moving particles. From a mathematical point of view, the particles are just interpolation points to approximate so- lutions of partial differential equations (PDEs). From a physical point of view, the particles are material particles which can have properties that dif- fer from the surrounding particles. Before representing the materials with moving particles, one need to define the interactions between them so as to faithfully reproduce the equations of fluid dynamics or continuum mechanics [100]. Gingold and Monaghan [48] and Lucy [88] proposed a kernel estima- tion technique that can represent approximately an arbitrary function at any point by the values of the function at the neighbor particles. The first applications were for a class of astrophysical problems. Later an improved algorithm was developed that conserves the linear and angular momentum exactly for a compressible nondissipative fluid [49]. The improvement made a clearer connection between SPH and molecular dynamics [65, 66]. 283.4. Smoothed particle hydrodynamics (SPH) An effective technique for approximating continuum partial differential equations, SPH is an attractive computational method for a very wide range of applications, such as shock waves simulation [98], incompressible flow [21], waves breaking on arbitrary structures [17], multiphase flows [101], liquid metal moulding [17] and elasticity and fracture [6]. Recently, SPH has been used as an effective method for biological application such as blood flow simulation in arteries [134] and mechanical modeling of RBC in healthy and infected state [146, 155, 157]. 3.4.1 Fundamentals SPH provides an approximation for the values of the field functions such as the density, velocity, energy and their derivatives at any point in the problem domain. The particle-wise discretized approximation of PDEs will turn into a set of ordinary differential equations (ODEs) with respect only to time due to the Lagrangian nature of SPH. Any standard integration routines can be utilized to solve this set of discretized ODEs. Solving any problem with SPH technique requires the following steps: • Filling the computational domain with a set of arbitrarily distributed particles. No connectivity information is needed for the particles thanks to meshless nature of SPH method. • Formulating the field variables in the integral representation. This is usually called the SPH kernel approximation. • Replacing the integral approximation of field variables with summa- tion over neighboring particles. This is known as the SPH particle approximation. The interaction length is defined based on the com- pact support of a kernel function. • Replacing PDEs from governing equations with ODEs with respect only to time by virtue of Lagrangian particle approximation. • Time marching by solving the ODEs with an explicit integration al- gorithm borrowed from conventional finite difference. 3.4.2 Kernel approximation of functions and their gradients Two steps of approximation are needed to derive any SPH formulation. First, the original PDEs are written in the integral representation form. Second, one needs to derive the particle approximation of these integral 293.4. Smoothed particle hydrodynamics (SPH) equations. The integral representation of an arbitrary function f can be expressed as f(r) = ∫ f(r′)δ(r − r′)dr′, (3.8) where δ is the Dirac delta function. The first approximation step starts with replacing the Dirac delta function by a smooth kernel function W [87]: f(r) ' ∫ f(r′)W (r − r′, h)dr′. (3.9) This integral representation of f(r) is called the kernel approximation. The kernel function is chosen to be an even function having the following desir- able properties: • Compactness: W (r − r′, h) = 0 when |r − r′|/h > κ • Approximation of delta function: limh−→0 W (r − r′, h) = δ(r − r′) • Normalization: ∫ W (r − r′, h)dr′ = 1 • Smoothness. A kernel function is schematically shown in Fig. 3.1. The compact sup- port h defines the finite interaction length between particles or the influence area of each particle. κ is a constant that depends on the type of ker- nel functions and determines its the effective (non-zero) area. Because of the compact support, the integration is limited to the neighboring particles within the compact support of the kernel function. The kernel approxima- tion operator is commonly indicated by the angle bracket <> in the SPH literatures [87]. Therefore, the equation (5.2.3) maybe rewritten as < f(r) >= ∫ f(r′)W (r − r′, h)dr′. (3.10) Taylor expansion may be used to estimate the errors in the SPH kernel representation: < f(r) >= ∫ [f(r) + f ′(r)(r′ − r) +O((r′ − r)2)]W (r − r′, h)dr′ = f(r) ∫ W (r − r′, h)dr′ + f ′(r) ∫ (r′ − r)W (r − r′, h)dr′ +O(h2) = f(r) +O(h2), (3.11) 303.4. Smoothed particle hydrodynamics (SPH) Figure 3.1: Schematic representation of the compact support of a kernel function that determines the finite interaction length between particles. where O() stands for the residual. The second integral vanishes because W is an even function of r:∫ (r′ − r)W (r − r′, h)dr′ = 0. (3.12) We have also used the unity condition of the kernel function. Therefore, the SPH kernel approximation of a function is of second order accuracy. The approximation of the spatial derivatives of a function can be achieved by replacing f(r) with the derivatives in equation (3.10), and then perform an integration by parts. For example: < ∇ · f(r) >= ∫ [∇ · f(r′)]W (r − r′, h)dr′ = − ∫ f(r′) · ∇W (r − r′, h)dr′. (3.13) The surface integrals have been put to zero in the above since the kernel has compact support, and is zero on the domain boundary. An exception occurs for particles near the domain boundary, and special treatments have to be incorporated to address the boundary deficiency. 313.4. Smoothed particle hydrodynamics (SPH) 3.4.3 Particle approximation In SPH model the continuous integral representations in kernel approxi- mation equation is computed by summation over neighbor particles within the compact support. This discretized process is the second major approx- imation step in SPH method. After discretizing the whole computational domian with finite number of particles, the infinitesimal volume dr ′ at the location of particle j may be replaced by the finite volume of the particle δVj = mjρj . wheremj and ρj are mass and density of the particle, respectively. Thus, the discretized particle approximation of any continuous function can be expressed as f(r) ' ∫ f(r′)W (r − r′, h)dr′ ∼= N∑ j=1 f(rj)W (r − rj, h)δVj = N∑ j=1 mj ρj fjW (r − r′j, h). (3.14) For function values on any point ri, we will have f(ri) ≈ N∑ j=1 mj ρj fjWij , (3.15) where Wij = W (r− r′j, h). The same procedure can be used to compute the particle approximation of spatial derivatives of an arbitrary function: ∇f(ri) ≈ N∑ i=1 mj ρj (fj)∇W (r − r′j , h), (3.16) ∇ · f(ri) ≈ N∑ j=1 mj ρj f(rj) · ∇iWij , (3.17) where ∇iWij is taken with respect to particle i: ∇iWij = xi − xj rij ∂W (ri − rj, h) ∂rij = xij rij ∂Wij ∂rij . (3.18) However, the above expression of derivatives will not vanish exactly if f has a constant function. Using an identity, we can recast the particle approx- imation so as to ensure a zero value for derivative of a constant function. For any differentiable function Φ, ∇f = 1Φ[∇(fΦ)− f∇(Φ)]. (3.19) 323.4. Smoothed particle hydrodynamics (SPH) The SPH particle discretized form of this equation is ∇f(ri) = 1Φi N∑ i=1 mj ρj Φj(fj − fi)∇iWij , (3.20) which does vanish for a constant f . Taking Φ = 1 and Φ = ρ will give us two formulas for the particle approximations of gradients: ∇f(ri) ≈ N∑ i=1 mj ρj (fj − fi)∇iWij, (3.21) and ∇f(ri) ≈ 1ρi N∑ i=1 mj(fj − fi)∇iWij. (3.22) These are identical to the order of approximation of the SPH scheme. 3.4.4 SPH for simulating fluid flow In a Lagrangian framework, the governing equations for a viscous flow prob- lem include the continuity and the momentum equations: Dρ Dt = −ρ∇ · v, (3.23) Dv Dt = g + 1 ρ ∇ · τ − 1 ρ ∇p, (3.24) where t is time, g is the body force, p is pressure, v is the velocity vector, τ is the viscous stress tensor, D/Dt refers to the material derivative computed by following the motion of particles. These equations need to be supplemented by an incompressibility condition ∇ · v = 0. (3.25) However, a weakly compressible model is often used, in which the above is replaced by a postulated equation of state. This makes possible a simple explicit time-marching scheme. In the following we give more details about the two different numerical approaches. 333.4. Smoothed particle hydrodynamics (SPH) 3.4.5 Weakly compressible SPH (WCSPH) method Monaghan [96] and Morris et al. [102] showed that the incompressibility in low Reynolds number flows can be approximated by assuming a compressible fluid with a large sound speed. Their results demonstrated an acceptable accuracy by using Mach number of M < 0.1. This approach is usually called weakly compressible SPH (WCSPH). The slight compressibility causes some fluctuation in the density of fluid particles. Also to satisfy the CFL condition a very small time step needed to be used in WCSPH. In this scheme the system of governing equations including the continuity equation and the momentum equation are supplemented by following equation of state p p0 = ( ρ ρ0 )γ − 1, (3.26) where γ = 7, c is the speed of sound and ρ0 is the initial reference density, and p0 = c2ρ0/γ. In theory, the compressibility effect is on the order of O(M2). So, having the M ≤ 0.1 will limit the density variation to the order of 1% [21]. The benefit of relaxing the incompressibility condition is to enable an explicit scheme for solving the Navier-Stokes system. This can be implemented in a two-step prediction-correction scheme. The two-step algorithm strives to restore the initial density of each par- ticle and simulate an incompressible fluid flow. In the prediction step, an intermediate velocity v˜ is calculated solely based on body force g and the viscous force in the momentum equation, and the particle positions are up- dated accordingly. We have not imposed any incompressibility constraint so far and we are expecting to see deviation from the initial density for some particles. The intermediate density of each particle can be computed di- rectly from the updated particle position based on equation (3.24). In the correction step of the algorithm, a pressure is calculated from the perturbed density field using the equation of state. Thus the particles in denser areas will have a higher density and as a result a higher pressure, and vice versa. This pressure gradient is used to compute a velocity correction: vˆi = −∆t ∑ j mj ( pi ρ2i + pj ρ2j ) ∇iWij, (3.27) which tends to disperse or aggregate particles from areas of higher or lower densities, respectively. Thus, the new velocity v = v˜ + vˆ at the end of the time step should be approximately divergence-free. Finally, the position of 343.4. Smoothed particle hydrodynamics (SPH) the particles is updated by a central differencing scheme: ri(t+∆t) = ri(t) + ∆t2 [vi(t) + vi(t+∆t)] . (3.28) To have a stable solution, the time-step ∆t should be constrained by the CFL condition [102]: 4t ≤ 0.25h c (3.29) and another condition based on viscous diffusion: 4t ≤ 0.125h 2 ν . (3.30) 3.4.6 Truly incompressible SPH algorithm Simulating the fully confined moderate and high Reynolds number flows with WCSPH usually result in large pressure fluctuations. This is because the small density fluctuations tend to be amplified by the power index γ. Following the classical projection scheme of Chorin [15] and Temam [149] for incompressible Navier-Stokes equations, an SPH projection method has been proposed to deal with high Reynolds number flows [21]. This truly incom- pressible SPH (ISPH) method treats the pressure and the velocity as primi- tive variables and can present a very good approximation of a divergence-free velocity field. ISPH employs a projection scheme. The momentum equation will be spilt into two parts: a prediction step based on the viscosity term and the body force terms v˜ − vn 4t = g + 1 ρ ∇ · τ , (3.31) and a correction step based on the pressure gradient to enforce the divergence- free condition on the velocity: vn+1 − v˜ 4t = − 1 ρ ∇pn+1. (3.32) The intermediate velocity v˜ usually is not divergence free. So, applying the divergence operator on the two sides of equation (3.32) will result in ∇ · ( 1 ρ ∇pn+1 ) = ∇ · v˜ 4t . (3.33) 353.4. Smoothed particle hydrodynamics (SPH) With a constant density, the above equation can be rewritten as the pressure Poisson equation: ∇2pn+1 = ρ4t∇ · v˜, (3.34) where ∇2 is the Laplacian operator. The pressure field will be calculated implicitly and the velocity field will be updated by the pressure gradient as vn+1 = v˜ − ( 1 ρ ∇pn+1 ) 4t. (3.35) Thus, the incompressibility is implemented at the cost of solving the pressure- Poisson equation. 3.4.7 SPH approximation of the flow equations Using either WCSPH or ISPH, the particle approximation of the various terms need to be evaluated, and these formulas are summarized below. In ISPH, one needs to compute the divergence of velocity: ∇ · vi ≈ 1 ρi N∑ i=1 mj(vj − vi) · ∇iWij (3.36) as well as a formula for the Laplacian of pressure. Lee et al. [78] proposed a very stable discritized form of pressure Poisson equation: ∇2pi ≈ 2 ρi N∑ i=1 mj pijrij · ∇iWij (rij2 + 0.01h2) , (3.37) where pij = pi − pj. In WCSPH, we need to evaluate the variation in particle density subject to disturbance to the flow field or particle position. Two equivalent formulas have been developed. The interpolant form can be obtained with replacing the function f by density in the equation (3.15), thus giving us [95] ρi = N∑ j=1 mjWij . (3.38) This is convenient in computing ρ from updated particle positions. The second formula is derived by discretizing the continuity equation (3.24)[95]: dρi dt = N∑ j=1 mjvij · ∇iWij, (3.39) 363.4. Smoothed particle hydrodynamics (SPH) where vij denotes vi − vj. The momentum equation contains a pressure gradient term and a viscous term. By using the equation (3.22), the particle approximation of pressure gradient term may be written as ∇pi = 1 ρi N∑ j=1 mj(pj − pi)∇iWij (3.40) This form of pressure gradient will vanish for a constant pressure field but it is not able to conserve the linear and angular momentum [95]. To overcome this problem, we first transform the pressure gradient term as follows: ∇p ρ = ∇ ( p ρ ) + ( p ρ2 ) .∇ρ (3.41) Then the pressure gradient term in equation (3.24) can be rewritten in a symmetric form: ∇p ρ = N∑ j=1 mj ( pj ρ2j + pi ρ2i ) ∇iWij. (3.42) The linear and angular momenta will be conserved by equation (3.42) due to a symmetric central force between pairs of particles [95]. The effects of constant external pressure can be cancelled by replacing the value of pressure p by p− pex in all of mentioned equations. The approximation of the viscosity term has a somewhat tortuous his- tory. At the beginning, SPH was widely used to model compressible fluids at high Reynolds number (Re > 103) [2]. An artificial viscosity has been widely used to handle high-Mach number shock waves. The artificial viscosity was designed for the same purpose as the von Neumann-Richtmyer viscosity in finite-difference methods. The most widely used form was introduced by Monaghan [95]: Πij = { −αcµ˜ij+βµ˜ij2 ρij if vij · rij < 0 0 otherwise (3.43) µ˜ij = hvij · rij rij2 + 0.01h2 , (3.44) where ρij = (ρi + ρj)/2 and the term 0.01h2 guarantees a nonzero value for the denominator. Some researchers have taken this form of artificial viscosity to model the real viscosity [2]. Morris et al. [102] showed the inaccuracy of velocity 373.4. Smoothed particle hydrodynamics (SPH) profiles produced by using such a viscosity in low Reynolds number flows. In its place, they proposed a more realistic form of viscosity for approximating the viscous diffusion, which is modeled after expressions used in [97] for heat conduction. Since then, the viscous term of Morris et al. [102] has been widely adopted for viscous flow simulations, with good accuracy: {(1 ρ ∇ · µ∇)v}i = N∑ j=1 mj(µi + µj)rij · ∇iWij ρiρj(rij2 + 0.01h2) vij . (3.45) This form of viscosity diffusion conserves the linear momentum exactly but the angular momentum only approximately [102]. Finally the SPH momen- tum equation (3.24) including the simplified form of viscous equation (3.45) can be rewritten as dvi dt = gi − N∑ j=1 mj ( pj ρ2j + pi ρ2i ) ∇iWij + N∑ j=1 mj(µi + µj)vij ρiρj ( 1 rij ∂Wij ∂ri ) . (3.46) 3.4.8 Limitations of SPH for open-boundary flow simulations As mentioned before, the SPH method originated from astrophysical prob- lems some three decades ago [48, 88] and in recent years has evolved into a unique tool for computational fluid dynamics. Its meshless formalism makes it a very unique tool to simulate large strains and morphological changes [67, 110, 148]. Nevertheless, SPH is still a relatively new tool in compu- tational fluid dynamics, and is not as well developed as finite difference or finite element [99, 100]. Some fundamental issues, such as the treatment of incompressibility and boundary conditions, have not been fully resolved. Traditionally, SPH methods employ a weakly compressible formalism with an artificial equation of state that specifies pressure as an algebraic function of density. A large speed of sound is used to maintain an ac- ceptable density variation [95, 102, e.g.]. This weakly compressible SPH (WCSPH) algorithm is easy to program but has several disadvantages. The high speed of sound imposes a severely restrictive CFL criterion on the time step [21]. The pressure is subject to large and non-physical fluctuations, es- pecially for relatively coarse spatial resolution. Such fluctuations may cause numerical instability. More important, they corrupt the result in problems where pressure is of physical interest, as in computing the drag on a solid 383.4. Smoothed particle hydrodynamics (SPH) object [78]. To bypass these difficulties, Cummins and Rudman [21] intro- duced truly incompressible SPH (ISPH) algorithms based on the projection scheme. Numerical results have demonstrated that ISPH largely eliminates the density and pressure fluctuations, produces more accurate predictions of velocity and forces on solids, and is in general more efficient than WCSPH [21, 78]. However, these ISPH algorithms suffer from an inconsistency in the pres- sure boundary condition. In the pressure-correction projection scheme, pres- sure is computed from a Poisson equation, and is then used to correct the velocity field to make it divergence-free [15, 149]. Solving the Poisson equa- tion requires boundary conditions for the pressure, which are not required by the original Navier-Stokes equation. Following the pioneers of the origi- nal projection scheme [15, 149], Cummins and Rudman [21] and Lee et al. [78] both employed a homogeneous Neumann condition: n · ∇p|Γ = 0, (3.47) where n is the normal vector to the boundary Γ. In many flow situations, of course, the normal pressure gradient does not vanish on the boundary. Important examples are open-boundary problems where the fluid enters and exits the computational domain driven by a pressure gradient, and flow around bluff bodies or through channels with variable cross sections. In such cases, the simple homogeneous condition of Eq. (3.47) produces a numerical boundary layer at the boundary and corrupts the accuracy of the solution. Partly for this reason, open-boundary flows are rarely computed using SPH. The standard practice is to impose periodic conditions on the inlet and the exit, and replace the pressure gradient driving the flow by a body force [78]. In fact, Cummins and Rudman [21] and Lee et al. [78] both concluded their paper by highlighting the need to develop a more suitable pressure boundary condition. This is a well-known problem in the projection scheme, and is by no means specific to SPH. A physically consistent boundary condition on pres- sure has long been sought after [53]. Recently Guermond et al. [55] pub- lished a careful examination of the various projection methods, and recom- mended the “rotational pressure-correction scheme” of Timmermans et al. [151], with a nonhomogeneous pressure condition for the Poisson equation. Their finite-element numerical tests show that this largely removes the arti- ficial boundary layer caused by Eq. (3.47) and preserves the accuracy of the solution throughout the domain. We implement in SPH the rotational pressure-correction scheme, with its nonhomogeneous pressure boundary condition. The details of the method 393.4. Smoothed particle hydrodynamics (SPH) and the validations were explained in chapter 4. Our work is inspired on the one hand by the need for a suitable pressure boundary condition in ISPH [21, 78], and on the other by the successful finite-element computations using the rotational scheme [55]. We have validated our proposed algorithm by computing pressure-driven Poiseuille flow with open boundaries and flow past a square cylinder. In each case, the SPH result is in very good agreement with analytical or finite- element benchmarks. Thus, the new pressure boundary condition enhances the capability of SPH in computing incompressible flows; previously inac- cessible problems with open and solid boundaries can now be computed routinely. 40Chapter 4 Pressure boundary conditions for computing incompressible flows with SPH 4.1 Introduction As discussed in section 3.4.8, previous incompressible SPH formalisms suf- fer from an inconsistent boundary condition when simulating open-boundary flows and flow around obstacles. This is rooted in the traditional projec- tion scheme in the solution algorithm. As a result, open-boundary flows are rarely computed using SPH, and they are usually replaced by a problem with periodic conditions on the inlet and exit, with the flow being driven by a body force. The whole problem stems from the lack of a suitable pressure boundary condition; the simple homogeneous condition of Eq. (3.47) pro- duces a numerical boundary layer at the boundary and corrupts the accuracy of the solution. We propose here to use a rotational pressure-correction scheme in SPH method, with a nonhomogeneous pressure boundary condition for the pressure- Poisson equation. Our results show the superiority of this scheme over the traditional scheme. In this chapter the details of implementation of the new algorithm are explained. The new ISPH algorithm is validated by simu- lating pressure-driven Poiseuille flow with open boundaries and flow past a square cylinder. In each case, the SPH result is in very good agreement with analytical or finite-element benchmarks. Thus, the new pressure boundary condition enhances the capability of SPH in computing incompressible flows; previously inaccessible problems with open and solid boundaries can now be computed routinely. 414.2. Projection schemes 4.2 Projection schemes 4.2.1 Homogeneous projection scheme The classical projection scheme of Chorin [15] and Temam [149] features the homogeneous pressure boundary condition (Eq. 3.47), and will be called hereafter the homogeneous projection scheme. Consider the flow of an in- compressible fluid governed by the Navier-Stokes equations: ∇ · u = 0, (4.1) du dt + 1 ρ ∇p− ν∇2u = g, (4.2) where g is a body force. For the time being, let us assume a Dirichlet boundary condition on the velocity: u|Γ = uΓ. (4.3) Neumann conditions on n·∇u, as appropriate for the exit in open-boundary flows, will be considered in subsection 4.2.3. The projection scheme is a prediction-correction technique. In the pre- diction step, we compute an intermediate velocity field u˜ from the momen- tum equation (Eq. 4.2) without the pressure term. This velocity satisfies the boundary condition (Eq. 4.3), but in general not the continuity equa- tion (Eq. 4.1). In the correction step, we compute a pressure field from the intermediate velocity, and then use the pressure gradient to correct u˜ into a solenoidal velocity u. Using second-order time stepping, the scheme can be written as [55]: 1 2∆t(3u˜ k+1 − 4uk + uk−1)− ν∇2uk = g, u˜k+1|Γ = uΓ, (4.4) 1 2∆t(3u k+1 − 3u˜k+1) + 1 ρ ∇pk+1 = 0, ∇ · uk+1 = 0, n · uk+1|Γ = n · uΓ, (4.5) where k and k + 1 indicate the old and new time steps, respectively. In the SPH implementation, the temporal differencing gives the material derivative and thus no advection term appears. For our purpose here, the most inter- esting aspect is that pk+1 is computed from a Poisson equation by taking the divergence of the first equation of Eq. (4.5): ∇2pk+1 = 3ρ 2∆t∇ · u˜ k+1 . (4.6) 424.2. Projection schemes Equation (4.5) implies an homogeneous Neumann boundary condition for pk+1: n · ∇pk+1|Γ = 0, (4.7) which is evidently inconsistent with the momentum equation (Eq. 4.2). Un- less the physical problem is such that the normal pressure gradient indeed vanishes at the boundary, this homogeneous Neumann condition produces a numerical boundary layer in the solution and corrupts its accuracy [55]. 4.2.2 Rotational projection scheme In the prediction step, one may retain a pressure gradient based on the prior time step when computing the intermediate velocity u˜k+1 from the momentum equation. Then what appears in the correction step will be a pressure difference between the current and prior steps. The essence of the rotational scheme is to modify this pressure difference so as to make the pressure Neumann condition consistent with the momentum equation [151]. Following Guermond et al. [55], we write the scheme in this form: 1 2∆t(3u˜ k+1 − 4uk + uk−1)− ν∇2u˜k+1 + 1 ρ ∇pk = g, u˜k+1|Γ = uΓ(4.8) 1 2∆t(3u k+1 − 3u˜k+1) + 1 ρ ∇φk+1 = 0, ∇ · uk+1 = 0, n · uk+1|Γ = n · uΓ, (4.9) where φk+1 is the modified pressure defined as: φk+1 = pk+1 − pk + µ∇ · u˜k+1, (4.10) µ = ρν being the dynamic viscosity. The following observations can be made of this projection scheme. 1. Retaining the old pressure gradient ∇pk in the prediction step in- creases the accuracy of the scheme [158]. 2. Substituting Eq. (4.10) into Eq. (4.9) and then adding it to Eq. (4.8), we arrive at the proper discretization of the momentum equation: 1 2∆t(3u k+1 − 4uk + uk−1)− ν∇2uk+1 + 1 ρ ∇pk+1 = g, (4.11) where we have used the identity ∇(∇· u˜k+1)−∇2u˜k+1 = ∇×∇×u˜k+1 and the fact that ∇ × ∇ × u˜k+1 = ∇ × ∇ × uk+1 = −∇2uk+1 since ∇× u˜k+1 = ∇× uk+1 and ∇ · uk+1 = 0 from Eq. (4.9). 434.2. Projection schemes 3. The use of φk+1 for correcting the velocity, instead of pk+1, is the key difference from the homogeneous scheme. φk+1 can be solved from a Poisson equation with its own homogeneous Neumann condition implied by Eq. (4.9): ∇2φk+1 = 3ρ 2∆t∇ · u˜ k+1 , n · ∇φk+1|Γ = 0. (4.12) Translated into pressure, these become ∇2pk+1 = ∇2pk − µ∇2(∇ · u˜k+1) + 3ρ 2∆t∇ · u˜ k+1 , n · ∇pk+1|Γ = n · ∇(pk − µ∇ · u˜k+1)|Γ. (4.13) Using Eq. (4.8) and the identities following Eq. (4.11), the nonhomo- geneous Neumann boundary condition for pk+1 can be rewritten as n · ∇pk+1|Γ = n · (ρg + µ∇2uk+1)|Γ, (4.14) which is consistent with the momentum equation (Eq. 4.11). It was called the “rotational” pressure boundary condition because the vis- cous stress term µ∇2uk+1 = −µ∇×∇× uk+1 was viewed as the curl of the vorticity vector [55]. 4.2.3 Natural boundary condition So far, our discussion is limited to Dirichlet boundary conditions for the velocity. With open-boundary flows, it is more appropriate to impose a natural boundary condition at the exit of the flow domain. Denoting by Γ1 the portion of the boundary with Dirichlet velocity condition and Γ2 the exit, the boundary conditions are u|Γ1 = uΓ, (4.15)[ −pn+ µn · (∇u+∇uT)]|Γ2 = −pan, (4.16) 444.3. SPH implementation where pa is the ambient pressure. Following Guermond et al. [55], we modify the rotational scheme as follows: 1 2∆t(3u˜ k+1 − 4uk + uk−1)− ν∇2u˜k+1 + 1 ρ ∇pk = g u˜k+1|Γ1 = uΓ, [−pkn+ µn · (∇u˜k+1 + (∇u˜k+1)T)]|Γ2 = −pan, (4.17) 1 2∆t(3u k+1 − 3u˜k+1) + 1 ρ ∇φk+1 = 0, ∇ · uk+1 = 0, n · uk+1|Γ1 = n · uΓ, φk+1|Γ2 = 0 (4.18) φk+1 = pk+1 − pk + µ∇ · u˜k+1. (4.19) Note that Eq. (4.18) implies that on Γ1 the homogeneous Neumann con- dition is imposed on φk+1 as before, but on Γ2 a Dirichlet condition is used. Equivalently, pk+1 satisfies the non-homogeneous Neumann condition Eq. (4.14) on Γ1 and a non-homogeneous Dirichlet condition on Γ2. The latter removes the indeterminacy in the pressure-Poisson equation. The ho- mogeneous scheme (Eqs. 4.4, 4.5) can be similarly modified, with a Dirichlet pressure condition at the exit: pk+1|Γ2 = 0. 4.3 SPH implementation To some degree, SPH particles may be viewed as representing blobs of ma- terial. But they are for the most part a numerical device for discretizing the Navier-Stokes equations in Lagrangian coordinates. The particles are interpolation points from which the fluid properties may be calculated by smoothing over neighboring particles. Although the particles have mass and move according to Newton’s law of motion, the forcing terms stem directly from discretizing the continuum governing equations. Solid boundaries are discretized by solid particles, and the fluid-solid interaction is defined so as to satisfy the no-slip boundary condition. 4.3.1 SPH interpolation The SPH method allows any property to be interpolated from its values at a set of discrete points—the SPH particles—using a kernel or weighting function W (r − r′, h), which specifies the contribution to any field variable at position r by a particle at r′ that lies within the compact support of the kernel function. The range of the compact support, indicated by the length scale h, determines the maximum interaction length between particles. The 454.3. SPH implementation weighting function is normalized such that∫ V W (r − r′, h)dr′ = 1, lim h→0 W (r − r′, h) = δ(r − r′), (4.20) V being the entire space. If a field variable A(r) is known only at a dis- crete set of particles r′j , then its value at any spatial location r can be approximated by: 〈Ah (r)〉 = ∫ V A(r′)W (r − r′, h) dr′ = ∑ j A(r′j)W (r − r′j , h)Vj = ∑ j mj ρj A(r′j)W (r − r′j, h),(4.21) where Vj is the volume element at r′j , and has been replaced by the ratio between the mass and density of the jth particle: Vj = mj/ρj . The sum- mation is over all particles that lie within a circle or sphere with the radius of compact support of the kernel function. In our ISPH formalism, mj and ρj are constant for all particles. The subscript j is retained to be general and consistent with convention. The same is true for the representation of viscosity µ in formulae below. We have tested the cubic and quintic spline kernels [102] as well as the truncated Gaussian kernel [17]. The Gaussian kernel provides a good com- promise between computational cost and accuracy, and has been used in all the results to be presented. In terms of the distance s = |r − r ′|/h, it can be written as W (s) = { e−s 2 −e−9 pih2(1−10e−9) 0 ≤ s ≤ 3 0 s > 3. (4.22) We have used h = 1.2L0, L0 being the initial spacing between neighboring particles in a regular square lattice. In classical SPH methods, the gradient approximations is based on con- volving the variables with the kernel function [100]: 〈∇Ah(r)〉 = ∫ V A(r′)∇W (r − r′, h) dr′ =∑ j mj ρj A(r′j)∇W (r − r′j, h). (4.23) The interpolation among neighboring particles is equivalent to a Taylor ex- pansion that theoretically has second-order accuracy [111]. However, as the particles move and their spatial distribution becomes irregular, the accu- racy in approximating the gradient vector can be much compromised. To 464.3. SPH implementation improve the accuracy of the kernel-based approximation, a number of cor- rective procedures have been proposed [13, 111, e.g.]. In our simulations, the normalization technique proposed by Oger et al. [111] works well. They constructed a correction matrix that accounts for irregular distribution of neighbors around each particle, which in two dimensions is written as L(r) = [ ∑ j Vj(xj − x)∂Wij∂x ∑ j Vj(xj − x)∂Wij∂y∑ j Vj(yj − y)∂Wij∂x ∑ j Vj(yj − y)∂Wij∂y ] −1 , (4.24) with Vj = mj/ρj and Wij is a shorthand for W (ri − rj , h). Multiplying L(r) into the gradient of the kernel produces high accuracy even with highly irregular particle distributions [111]. For example, the pressure gradient and velocity divergence are computed as ∇pi = ∑ j mj ρj (pj − pi)L(ri) · ∇iWij , (4.25) ∇ · ui = 1 ρi ∑ j mj(uj − ui) ·L(ri) · ∇iWij, (4.26) where pi = p(ri), ui = u(ri), and ∇i indicates differentiation with respect to ri. Note that the right-hand-sides of Eqs. (4.25) and (4.26) use the “symmetric” form of the interpolation formula for the gradient [100], an oft-used alternative to Eq. (4.23). In simulating the flow around a square cylinder, in particular, the particle distribution becomes highly disordered at the front and back of the obstacle. The classic formulae would produce large pressure fluctuations, which are mostly suppressed by the corrected forms of the gradient and divergence formulae. 4.3.2 Solution algorithm For both the homogeneous and rotational schemes, the algorithm reflects the same two-step structure. We will outline the solution procedure by referring to the equations of the rotational scheme, and a similar one is followed for the homogeneous scheme. The prediction step consists in solving a Helmholtz equation for the intermediate velocity u˜k+1 (Eq. 4.8 or 4.17), while the correction step a Poisson equation for the modified pressure φk+1 (Eq. 4.12 or its counterpart for natural boundary conditions at the exit). In a time-dependent setting, we have tested implicit and explicit schemes in both steps. Balancing stability and efficiency, we have settled on a scheme similar to that of Lee et al. [78] of solving for the intermediate velocity 474.3. SPH implementation explicitly in the prediction step and for the modified pressure implicitly in the correction step. Of course, in the end we correct the velocity to make it divergence-free via Eq. (4.9) or Eq. (4.18). The Laplacian of velocity and modified pressure are discretized in SPH using the following formulae [102, 160]: ν∇2ui = ∑ j mj(µi + µj)rij ·L(ri) · ∇iWij ρiρjr2ij (ui − uj), (4.27) ∇2φi = ∑ j 2mj ρj rij · L(ri) · ∇iWij r2ij (φi − φj), (4.28) where rij = ri−rj and rij = |rij |. In solving the pressure-Poisson equation, we combine Eq. (4.28) for all the particles to construct a linear system of equations for φk+1. If all boundary segments of the domain have Dirichlet boundary condition on u and Neumann boundary condition on φ, there will be an indeterminacy in φ and the resulting linear system will be degenerate. This is to be expected since in the incompressible Navier-Stokes system, the absolute value of pressure is immaterial. The indeterminacy is removed by replacing one of the diagonal elements of the coefficient matrix by a very large number, thereby setting that φ value to nil [147]. The linear system is solved by using the stabilized version of the bi-conjugate gradient method (Bi-CGSTAB). 4.3.3 Boundary conditions Because of the Lagrangian nature of SPH particles, implementing boundary conditions at fixed spatial boundaries is less straightforward than in mesh- based methods. To be specific, we discuss three types of boundaries that appear in our simulations: solid walls, inflow boundary (entry) and outflow boundary (exit). These come from the geometry of flow around an obstacle in a channel (Fig. 4.1), and are common to a broad class of flow simulations. Solid walls For velocity, the no-slip boundary condition on solid walls has been imple- mented in earlier work with image particles [21, 83], dummy particles with zero velocity [70, 78], dummy particles with extrapolated velocity [102], and a normalization procedure to compensate for the boundary deficiency of neighboring particles [123]. Here we adopt the method of Morris et al. [102] of extrapolating the fluid velocity onto layers of dummy particles inside the 484.3. SPH implementation Figure 4.1: Computational domain for flow past a square cylinder of side D in a channel of width H , showing the three types of boundaries to be encountered in the simulation. solid wall. Initially we deploy a square lattice of SPH particles with spacing L0 to cover the fluid domain as well as the solid surfaces. Furthermore, three layers of dummy particles are drawn inside the solid, also at spacing L0. Our Gaussian kernel has a cut-off range of 3h (Eq. 4.22) and we set h = 1.2L0. Thus three layers of dummy particles, in addition to the wall particles, cover the compact support of all fluid particles. At the beginning of the prediction step, we assign a velocity to each dummy particle that is the negative of the interpolated velocity of its image inside the fluid domain. These dummy particles will serve as neighbors in solving the momentum equation for the fluid particles, but they will not move from one time step to the next. The velocity of the wall particles remains zero, of course. This setup is used for the channel walls as well as the solid obstacle in the middle of the domain. For dummy particles along the diagonal inside the corners of the square cylinder (cf. Fig. 4.1), we follow the scheme of Lee et al. [78] by averaging between their neighbors. The extrapolation of velocity into the solid wall is essential for accurate evaluation of wall shear stresses and the overall drag [144]. On solid walls, the homogeneous Neumann condition for p (Eq. 4.7) or φ (Eq. 4.12) applies. Following Lee et al. [78], we simply propagate the p or φ value of the wall particle to the dummy particles along the normal direction. This is implemented numerically by manipulating the relevant entries in the linear system. 494.3. SPH implementation Entry Most SPH simulations up-to-date have replaced entry and exit conditions by periodic conditions [78, 102]. These are rather easy to implement; par- ticles exiting the computational domain are simply reinserted at the entry, carrying the same properties. Lastiwka et al. [76] have implemented a gen- uine entry condition, for compressible flows, by attaching an inflow zone upstream of the computational domain. Properties of particles inside the inflow zone are assigned to reflect the analytical boundary conditions to be realized at the entry. For our incompressible flow, the proper entry condition is a prescribed unidirectional velocity profile (uin(y), 0), say a parabolic profile representing fully developed Poiseuille flow. Our approach is similar in spirit to Lastiwka et al. [76]. We replicate the first line of fluid particles at the entry, called the entry layer, three times upstream at spacing L0. The entry layer and the three layers of image particles are all assigned the prescribed velocity profile uin. Thus, the fluid particles immediately downstream of the entry layer are preceded by four particles upstream, all bearing the prescribed velocity uin, and the boundary deficiency problem is solved. The Neumann condition n · ∇φ = 0 (or n · ∇p = 0 for the homogeneous scheme) is easily implemented by assigning the φ (or p) value of an entry-layer particle onto its images upstream. When a particle of the entry layer moves beyond L0 downstream of the inlet, it becomes an interior particle and its upstream image particle becomes part of the entry layer. Meanwhile, another image particle is injected far upstream to be the third image particle. Exit Under certain circumstances, it is possible to treat the exit the same way as the entry, with a prescribed velocity profile. For instance, this is reasonable for Poiseuille flow in a straight channel and for more complex channel flows with an exit far downstream of flow disturbances. We have tested these cases in numerical experiments. However, the natural boundary condition (Eq. 4.16) is more appropriate in most cases. We will deal with the simple case of equilibrated pressures at the exit (i.e., no suction or blowing), so that the normal stress balance of Eq. (4.16) reduces to a simple Neumann condition for the velocity n · ∇u = 0. To compensate for the boundary deficiency at the exit, we again copy the last layer of fluid particles (the exit layer) onto three layers of image particles 504.3. SPH implementation downstream of the exit with spacing L0, each bearing the same velocity as their progenitor in the exit layer. Now the velocity of the near-exit fluid particles can be updated as for the inner particles, and the viscous stress and velocity divergence can be evaluated in the standard way. For pressure, we impose the Dirichlet conditions p = 0 and φ = 0 at the exit for the homogeneous and rotational schemes, respectively. To implement this, the image particles downstream of the exit take on p or φ values extrapolated from the upstream particles, which are incorporated into the discretization of the Poisson equation by properly modifying the matrix. When a particle of the exit layer crosses the outlet boundary, it ceases to be a fluid particle and becomes an image particle, bearing the same velocity. Meanwhile the upstream fluid particle becomes part of the exit layer. The p or φ value of the new image particle is again extrapolated from the upstream fluid particles with the condition that p = 0 or φ = 0 at the exit. Finally, we mention a boundary-deficiency issue in computing ∇ · u˜ for the rotational scheme, which appears in the Poisson equation for φ and at the end of the correction step when extracting the physical pressure p from φ (Eq. 4.10). Computing ∇ · u˜ near the boundaries requires u˜ for the dummy or image particles outside the domain. At the entry and exit, we assign the u˜ value of the entry or exit layer to the image particles as described above, and thus ∇ · u˜ can be computed. Near solid boundaries, however, using the linearly extrapolated velocity for dummy particles, as is done in implementing the no-slip boundary condition, would produce large spurious pressure on the dummy particles and corrupt the solution near the solid. This is because the pressure correction scheme in Eq. (4.10) is designed for fluid particles and does not account for the peculiar property of the dummy particles in ostensibly carrying a velocity but not really moving. Instead, it is more reasonable to assign u˜ = 0 on the dummy particles inside solid walls in computing ∇ · u˜. Numerical experiments show that this more accurately predicts the pressure field near solid walls. Interestingly, the homogeneous scheme is very insensitive to u˜ values on the dummy particles. The pressure field differs only by some 2% between the two choices. This may be because in the homogeneous scheme, pressure pk+1 is computed directly from the Poisson equation with no further correction by ∇ · u˜. 4.3.4 Tensile instability The SPH representation of a fluid of constant and uniform density relies on a more or less uniform distribution of the particles. Because of how the particles interact via the kernel function, however, clustering of particles is 514.3. SPH implementation (a) X Y -1 0 1 -1 -0.5 0 0.5 1 (b) X Y -1 0 1 -1 -0.5 0 0.5 1 Figure 4.2: (a) Tensile instability caused by the particles forming anisotropic strings in flow around a square cylinder. (b) The particle shifting approach [160] largely removes the undesirable clustering and maintains a more or less uniform particle distribution. liable to arise, especially when the material is subject to tensile deformation. This gives rise to the well-documented tensile instability [99]. An example from our own computations is illustrated in Fig. 4.2(a) for flow around a square cylinder to be discussed in detail in section 4.4.3. Among the various remedies proposed in the literature [13, 52, 99, 108, 160], we have found the particle shifting strategy of Xu et al. [160] most ef- fective and efficient for our ISPH projection scheme. Originally proposed by Nestor et al. [108] for a Finite Volume Particle Method, the idea consists in shifting particles slightly away from their streamlines and correcting their velocity and pressure according to a first-order interpolation. The direction and amount of shifting are determined from the arrangement of nearby parti- cles, and we have adopted the formula of Xu et al. [160]. Figure 4.2(b) shows that this scheme effectively suppresses the formation of particle strings and averts the onset of tensile instability. 524.4. Numerical results 4.4 Numerical results 4.4.1 Poiseuille flow driven by body force in periodic domain As a baseline, we compute the transient development of a Poiseuille flow in a straight channel driven by a body force. The SPH simulation is set up in the usual way as in the literature [102], with the pressure gradient set to zero, and periodic boundary conditions imposed on both ends of the domain. Thus, n · ∇p does vanish on all boundary segments, and the homogeneous and rotational projection schemes become identical for this simple problem. We have two objectives. The first is to validate our code for both projection methods against this benchmark problem, which has an analytical solution [102]. The second is to examine the accuracy of the solution with increasing numbers of particles, and determine the necessary level of spatial resolution. The channel has a width H and a length 2H. Initially the particles are arranged in a uniform square pattern of spacing L0, with zero initial velocity. At t = 0, the particles start to move under a constant body force g. The temporal development of the velocity profiles is recorded and compared with the analytical solution. The final steady-state centerline velocity U0 = gH2/(8ν) corresponds to a Reynolds number Re = U0H/ν = 1. Figure 4.3 compares the computed velocity profiles at four times with the analytical solution, the last one t = 1 approaching steady state. In this unidirectional flow, the particles on the same streamline move with identical velocity, and each data point in the plot corresponds to one such row of particles. Length, velocity and time are made dimensionless by H, U0 and H/U0, respectively. Numerical solutions at two spatial resolutions are shown; initial particle spacing L0 = 0.05H and L0 = 0.1H correspond to a total of 1120 and 360 particles. Note first the excellent agreement between the SPH solutions and the analytical solution. This serves as a validation of our methodology and code. Second, the numerical solutions at both spatial resolutions are practically indistinguishable; this suggests that L0 = 0.1H is sufficient for such a simple problem. Finer resolutions are required for more complex flow geometries, as will be shown shortly. 4.4.2 Fully developed Poiseuille flow with open boundaries To investigate the effect of the pressure Neumann boundary condition in the homogenous and rotational schemes, we simulate the steady-state Poiseuille flow in a straight channel with entry and exit boundary conditions. The channel width and length are both H. The velocity and pressure boundary 534.4. Numerical results y/H u /U 0 -0.5 -0.25 0 0.25 0.5 0 0.2 0.4 0.6 0.8 1 Analytical solution SPH, L0 = 0.05H SPH, L0 = 0.1H t = 1.0 t = 0.15 t = 0.05 t = 0.01 Figure 4.3: Development of the Poiseuille velocity profile in time: comparison between SPH simulation and the analytical solution. Two numerical solutions are shown, corresponding to initial particle spacing L0 = 0.05H and 0.1H . Time is made dimensionless by H/U0. conditions are imposed as described in section 4.3.3. Specifically, the velocity assumes a parabolic profile at the inlet: uin(y) = U0H2 (H 2 − 4y2), (4.29) and vanishes on the side walls. The Reynolds number Re = U0H/ν = 1. On both the inlet and side walls, p or φ, in the two projection schemes respectively, satisfies the homogeneous Neumann condition. At the exit, the velocity has zero normal gradient, and the Dirichlet conditions p = 0 and φ = 0 are imposed for the homogeneous and rotational schemes. For this calculation we used a total of 2000 particles with L0 = 0.025H, which provides sufficient resolution. Once the initial transient dies out, both schemes accurately reproduce the parabolic velocity profile throughout the domain, and their difference 544.4. Numerical results (a) (b) Figure 4.4: Deviation of the computed pressure field p(r) from the analytical solution pA(r) (i.e., a linear distribution along the flow direction). (a) Homogeneous projection scheme. (b) Rotational projection scheme. The pressure error has been scaled by 12ρU 2 0 and zeroed at the center of the domain. lies in the pressure field. To illustrate the pressure boundary layers clearly, we subtract the analytical linear pressure distribution from the computed 554.4. Numerical results pressure field to produce a pressure error, with its value zeroed at the center of the domain for a clearer view. Figure 4.4 depicts the two-dimensional distribution of such pressure errors for the two schemes using a snapshot of the SPH particles. As expected, the homogeneous scheme produces a prominent boundary layer at the entry (Fig. 4.4a). However, there is also appreciable error at the exit. Unlike at the entry, the p = 0 condition here and the linear extrapolation of p for downstream image particles are con- sistent with the expected linear pressure profile. Thus the exit error arises from a different source than the entry boundary layer; it is due to errors in the velocity since a weaker Neumann condition n ·∇u = 0 is imposed at the exit. In Fig. 4.4(b), the rotational scheme effectively eliminates the entry boundary layer, and reduces the maximum error by an order of magnitude. Note that even in this case, the entry and exit tend to produce larger errors due to extrapolation for the image particles outside the domain. Given the analytical pressure gradient in a Poiseuille flow dp/dx = 8µU0/H2 and the domain length of H, the scaled pressure drop between the entry and exit would be 16/Re = 16. Relative to this, the maximum pressure error in Fig. 4.4(a) for the homogeneous scheme amounts to 0.3%, which is quite acceptable. The rotational scheme has a maximum error of about 0.04%. The results have clearly demonstrated the pressure boundary layer suf- fered by the homogeneous projection scheme and the capacity of the rota- tional scheme to eliminate it. But for this simple problem, the ill effect is small and limited to a small area of the domain, and does not seriously cor- rupt the quality of the solution on the whole. This is no longer the case for flow around solid objects, where the incorrect pressure boundary condition seriously compromises quantities of physical interest such as the drag force. 4.4.3 Channel flow past a square cylinder The last problem of our numerical experimentation is the flow around a square cylinder in a two-dimensional channel (Fig. 4.1). As a test of the rota- tional projection scheme, this problem is attractive for several reasons. First, flow around obstacles is known to cause difficulties for WCSPH schemes [78]. The artificial pressure is liable to large spurious oscillations near the solid surface, thus hampering an accurate computation of the drag. Such prob- lems are where the ISPH, with the rotational pressure boundary condition, is potentially most advantageous and useful. Second, this is an oft-used benchmark problem in the literature for testing the capability of various numerical methods to accurately calculate the pressure and velocity fields and the drag force [11, e.g.]. For the geometry and parameters of our SPH 564.4. Numerical results simulation, we have carried out a high-resolution finite-element solution as a benchmark. Finally and most interestingly, Lee et al. [78] have computed a similar problem using ISPH. In their setup, the fluid flow is driven by a body force in a channel containing the square cylinder, on which the homogeneous Neumann condition is imposed for the pressure. Periodic boundary condi- tions are used between the inlet and the outlet; this necessitates a very long section downstream of the obstacle to allow the flow to relax to the same condition at the inlet. Even with such precautions, the calculated drag coefficient is nearly 40% lower than their finite-volume benchmark. Thus, solving this problem using our schemes provides an opportunity to investi- gate the benefits of the non-homogeneous Neumann boundary condition on pressure. The geometry of the problem is schematically shown in Fig. 4.1. The channel width is H and its length is 3H, and a square cylinder of side D = H/8 is situated at the center of the channel. Thus, the entry is at x = −12D and the exit at x = 12D. Initially SPH particles are deployed on a regular square lattice of side L0. The square cylinder is represented by wall particles and interior dummy particles at the same spacing. The boundary conditions are imposed as explained in section 4.3.3, with a parabolic velocity at the entry and zero normal-velocity-gradient at the exit. For the homogeneous projection scheme, p = 0 at the exit and n · ∇p = 0 on the other three sides of the channel and on the surface of the square cylinder. For the rotational scheme, we require φ = 0 at the exit and n·∇φ = 0 on all the other boundary segments. The sharp corners of the square cylinder tend to cause spurious disturbances to the pressure of nearby particles. This is alleviated by local smoothing: in computing the pressure gradient ∇pk in the prediction step, the corner pressure is assumed equal to that of the nearest wall particle. In presenting results in the following, we scale length by the width of the square cylinder D, velocity by the centerline velocity U0, time by D/U0 and pressure by 12ρU20 . Most computations to be reported below correspond to a Reynolds number Re = U0D/ν = 1. Higher Re results will be compared with prior studies in subsection 4.4.3. As a benchmark, we have solved the same problem using our in-house finite-element software AMPHI [164]. The finest finite-element mesh size is approximately 0.06D and a total of 39973 elements are used. Mesh refinement shows that this spatial resolution is sufficient; doubling the total number of elements causes a 0.7% change in the drag. 574.4. Numerical results Convergence with number of SPH particles To examine how the SPH results depend on spatial resolution, we have car- ried out computations using the rotational projection scheme at five spatial resolutions, with L0 = H/40, H/58, H/80, H/120 and H/160 and the total number of particles being 5490, 11664, 21300, 46260 and 79860, respec- tively. The results are then compared with the finite-element (FE) solution to establish convergence with respect to increasing particle numbers, and to determine the level of resolution needed for an accurate solution that does not depend on the number of particles. The SPH solutions at the five spa- tial resolutions show the same qualitative trend. Starting with zero velocity throughout the domain except for the parabolic velocity profile imposed at the entry, the flow develops in time and approaches a steady state at t ≈ 1. This “steady state”, of course, is in an average sense that allows fluctuations due to the motion of individual particles. In the following, we will mostly discuss steady-state results. First, we examine how the drag on the square cylinder converges with increasing spatial resolution. The drag force fd, along the flow direction (x), is computed by integrating the pressure and viscous stresses around the surface of the square cylinder: fd = ∫ Sx ( −p+ 2µ ∂u ∂x ) nxdy+ ∫ Sy µ ∂u ∂y nydx = − ∫ Sx pnxdy+µ ∫ Sy ∂u ∂y nydx, (4.30) where Sx refers to the left and right faces of the cylinder with nx = ±1 being the x component of the outward normal vector, and Sy the top and bottom faces of the cylinder with ny = ±1 being the y component of the outward normal. The viscous normal stress vanishes at the solid surface Sx because ∂u/∂x = −∂v/∂y = 0. We use the pressure of the surface particles on the square cylinder for p. The viscous stress (or velocity gradient) is evaluated using fluid and wall particles as well as dummy particles inside the cylinder, the latter bearing velocity linearly extrapolated from that of the fluid and wall particles. This is known to produce the shear stress accurately on solid walls [144]. A drag coefficient is then defined as CD = fd 1 2ρU20D2 . (4.31) As an example, Fig. 4.5 shows the temporal variations of the drag force and its pressure and viscous components at a spatial resolution of L0 = H/80. Note that the forces have been made dimensionless into drag 584.4. Numerical results t C D /C DF E 0 2 4 6 8 100 0.5 1 1.5 2 CD Total / CD FE CD Pressure / CD FE CD Viscous / CD FE Figure 4.5: Temporal evolution of the drag coefficient CD and its components due to pressure and viscous friction. The spatial resolution is at L0 = H/80 with 21300 particles. The drag coefficient is scaled by the finite-element benchmark value CFED and time is made dimensionless by D/U0. coefficients, which are then scaled by the drag coefficient CFED computed by finite elements under the same conditions. The steady state is attained around t = 1. The subsequent fluctuations are small, with a magnitude on the order of 2% of the mean. The total drag coefficient falls below the finite- element benchmark by some 5% at this spatial resolution, and the pressure contribution is about twice that due to viscous friction. A time-averaged drag has been computed over a time period after attain- ing the steady state (1 ≤ t ≤ 10) to smooth out the fluctuations associated with the passage of individual particles around the cylinder. Figure 4.6 de- picts this averaged drag coefficient for the five particle numbers in terms of its deviation from the FE benchmark CFED . For the rotational scheme, CD oscillates around CFED and tends to converge to the latter as the number of 594.4. Numerical results Particle numbers (C D - C DF E )/ C DF E 0 20000 40000 60000 80000 100000 -0.3 -0.2 -0.1 0 0.1 Rotational Homogenous Figure 4.6: Drag coefficient CD computed by the rotational and homogeneous projection schemes at five particle numbers. The data are shown as fractional deviations from the benchmark FE drag coefficient CFED . particles increases. The maximum deviation is 5%. At the finest resolution, with 79860 particles, the drag is predicted within 2% of the benchmark. For comparison, CD computed using the homogeneous scheme is also plotted and will be discussed in the next subsection. Now we examine convergence of the velocity and pressure fields inside the computational domain. The time-averaged centerline profiles of the hor- izontal velocity u(x), computed using increasing numbers of particles, are compared with the FE profile in Fig. 4.7. In this and subsequent centerline profiles, the time-averaging is done as follows. We record the longitudinal position x, velocity u and pressure p of particles within 0.025L0 of the cen- terline for a time period. Then we divide the centerline into a number of bins of equal width and average the u and p for particles within each bin to 604.4. Numerical results (a) x / D u /U 0 -10 -5 0 5 10 0 0.2 0.4 0.6 0.8 1 FE (b) x/D (u - u FE )/ U 0 -10 -5 0 5 10 -0.08 -0.06 -0.04 -0.02 0 0.02 79860 particles 21300 particles 5490 particles Figure 4.7: (a) The centerline velocity profile u(x) computed by finite elements. (b) Deviation of the SPH velocity profiles from the FE profile computed using increasing numbers of particles. 614.4. Numerical results produce the centerline profiles. For clarity we omitted the results for 11664 and 46260 particles, as they conform to the same trend of convergence. The convergence with increasing number of particles is evident, with the max- imum error relative to the characteristic velocity U0 being around 7% for 5490 particles, 1% for 21300 particles and below 0.5% for 79860 particles. The convergence with increasing number of particles is also shown by the pressure profiles of Fig. 4.8. The greatest errors in pressure occur immedi- ately upstream and downstream of the square cylinder. Relative to the FE pressure drop across the obstacle, the maximum pressure error in SPH com- putations drops from about 13% for 5490 particles to 5% for 21300 particles and finally to 3% for 79860 particles. Note that this is several times greater than the velocity errors. But in comparison with the weakly compressible SPH scheme, the pressure is more accurately computed by a wide margin; WCSPH underpredicts the centerline pressure behind the square cylinder by up to 30% [78]. Based on the above discussions, it is clear that the SPH result converges with increasing number of particles, and L0 = H/80, with 21300 particles, provides sufficient spatial resolution. Results presented hereafter are all for 21300 particles. Superiority of the nonhomogeneous pressure boundary condition To ascertain the superiority of the rotational scheme over the homogeneous one for flows more complex than unidirectional channel flows, we have used both schemes to simulate the channel flow around the square cylinder at a spatial resolution of L0 = H/80, with 21300 particles. For both schemes, the overall behavior of the solution is the same. Starting with the parabolic ve- locity profile at the entry and zero velocity in the domain, the flow undergoes an initial transient before reaching a steady state. Although the transient is not of interest here, we note that the rotational scheme is more robust and allows twice as large a time step than the homogenous scheme without incurring numerical instability. This is a potential benefit for simulating transient flows. The time-averaged centerline velocity and pressure profiles are compared between the two schemes in Fig. 4.9, using the FE profiles of Figs. 4.7(a) and 4.8(a) as the baseline. The maximum velocity error is roughly 2% for the homogeneous scheme, and 0.8% for the rotational scheme. Similarly the pressure error relative to the FE pressure drop across the cylinder (Fig. 4.8a) decreases from roughly 9% for the homogeneous scheme to 5% for the rota- tional one. For the homogeneous scheme, the larger errors at the front and 624.4. Numerical results (a) (b) x / D (p - p F E )/ (0 . 5 ρ U 2 0 ) -10 -5 0 5 10 -1.5 -1 -0.5 0 0.5 1 1.5 79860 particles 21300 particles 5490 particles Figure 4.8: (a) The centerline pressure profile p(x) computed by finite elements. (b) Deviation of the SPH pressure profiles from the FE profile computed using increasing number of particles. Pressure has been scaled by 12ρU 2 0 and zeroed at the exit of the channel (x = 12D). 634.4. Numerical results (a) (b) Figure 4.9: Comparison between (a) velocity profiles and (b) pressure profiles along the centerline computed by homogenous and rotational schemes. We plot the deviation of the SPH results from the FE benchmarks in Figs. 4.7(a) and 4.8(a). 644.4. Numerical results the back of the square cylinder can be ascribed to the homogeneous Neu- mann boundary condition dp/dn = 0. At the entry of the channel, there is a pressure boundary layer similar to Fig. 4.4(a). But its magnitude is so much smaller than that at the solid cylinder as to be invisible in the plot. We further compare the rotational and homogeneous schemes by the pressure distribution along the front and back faces of the square cylin- der (Fig. 4.10). For the SPH results, each data point represents the time- averaged pressure on a wall particle. Within the central part of the solid surface, roughly −0.2 < y/D < 0.2, both SPH profiles stay fairly close to the FE profile. The velocity is low immediately upstream and downstream of the cylinder, and there is no recirculating eddy at Re = 1. Thus, p is relatively uniform in these “dead zones”. At the corner regions, however, the FE solution features sharp peaks in the front and valleys in the back, corresponding to the rapid turning of the streamlines. To a good degree, the rotational scheme captures these features while the homogeneous scheme misses out completely. But note the oscillations in the rotational profiles, indicative of the greater sensitivity of the scheme to errors in enforcing the solenoidality of the velocity. This difference in the surface pressure bears directly on the computed drag force on the cylinder. As the homogeneous scheme underpredicts p at the front and overpredicts it at the back, it produces a drag that is between 10% and 27% below the correct value in Figure 4.6, depending on spatial resolution. In contrast, the prediction of the rotational scheme falls within 5%. Based on the results discussed in this subsection, we may conclude that the rotational scheme, with the physically consistent nonhomogeneous Neumann condition for pressure, has accomplished the goals set out at the beginning of this work. Higher Reynolds numbers So far we have discussed numerical results at Re = 1 for flow around the square cylinder. This problem was chosen partly because Breuer et al. [11] and Lee et al. [78] have published similar simulations using finite volume, lattice-Boltzmann and SPH methods. Their computations are for Reynolds numbers up to a few hundred. In the following, we present computations using the rotational projection scheme at higher Reynolds numbers and compare them with prior results as another test of the nonhomogeneous Neumann condition for pressure. We have obtained accurate steady-state solutions at Re = 5, 10, 20, 30 and 50. The rotational projection scheme works equally well in all cases; 654.4. Numerical results (a) (b) Figure 4.10: Pressure profile along the front (a) and the back (b) of the square cylinder computed by finite elements and the homogenous and rotational schemes. In all three computations p is zeroed at the exit. the nonhomogeneous pressure boundary condition (Eq. 4.13) is apparently unaffected by the increasing inertia. However, higher Re does put a more 664.4. Numerical results stringent demand on the SPH scheme itself. For example, increasing the flow speed reduces the allowable time step via the CFL condition. At higher Re, fine wake structures may develop that demand higher spatial resolution, and a longer wake will require a longer domain. Both increase the magnitude of the computation and slow down convergence in our iterative solver for the pressure Poisson equation. In addition, we have also noticed that the simulation tends to become noisier at higher Re. For example, for Re = 1 there is no need to apply the XSPH velocity smoothing among neighboring particles [100], a remedy often used in the literature, especially in WCSPH scheme [102]. For Re ≥ 20, however, we notice small spatial fluctuations in the wake of the cylinder that can be removed by velocity-smoothing with an XSPH coefficient of 0.02. Finally, the natural boundary condition at the exit, Eq. (4.16), tends to generate small lateral velocity disturbances at higher Re. The Dirichlet boundary condition is generally more rigid than the Neumann condition, and this observation is not specific to SPH. For Re above 10, therefore, we have imposed the parabolic Poiseuille velocity profile at the exit. That the domain length is sufficient for this profile to develop at the exit has been confirmed separately by finite-element computations. The highest Reynolds number that we have tested is Re = 100, at which the SPH scheme fails to yield a converged solution. One expects an oscil- lating wake with vortex shedding in this case. Our scheme seems to incur a numerical instability with the oscillation growing in time, until the implicit solver for pressure fails. This may indicate that the wake now requires finer temporal and spatial resolutions and a longer domain length. Both the ho- mogeneous and rotational schemes fail at this Re, suggesting that the failure is due to factors other than the projection scheme and pressure boundary condition. Therefore we have not investigated this further. In the following, we will compare the steady-state solutions, obtained using the rotational projection scheme, with those in the literature. Figure 4.11 compares the drag coefficient on the square cylinder with previous predictions of Breuer et al. [11] using lattice-Boltzmann method, for Re up to 50. Our computational geometry is matched to that of [11], with H = 8D. Our domain length is 32D, shorter than their 50D, but the results are insensitive to this. In general, the SPH drag tends to fall slightly below that of the lattice-Boltzmann drag; the percentage error ranges from 2% for Re = 1 to under 7% for Re = 50. This is consistent with the accuracy noted in Figs. 4.5 and 4.6. Therefore, we consider the agreement satisfactory. In a slightly different geometry, with a channel width H = 5D and a length of 32D, Lee et al. [78] computed the flow around the square cylin- 674.4. Numerical results Re C D 0 10 20 30 40 50 10 20 30 40 50 Our SPH Breuer et al. Figure 4.11: Drag coefficient CD on the square cylinder as a function of the Reynolds number Re. The SPH computation with the rotational scheme uses 68000 particles. The lattice-Boltzmann results of Breuer et al. [11] are shown for compar- ison. der using WCSPH and ISPH, and compared the results with a finite-volume benchmark. At Re = 30, they obtained a pressure drag coefficient C pD = 2.46 by finite volume. The ISPH yielded CpD = 1.55, a 37% underprediction, which is consistent with the error of our homogeneous scheme in Fig. 4.6. The WCSPH scheme predicted a much larger CpD = 6.49, symptomatic of the spurious pressure in the weakly compressible treatment. Under identical conditions, our SPH method with the rotational scheme predicts C pD = 2.14, much closer to the benchmark than their ISPH using the homogeneous pro- jection scheme. Note that Lee et al. [78] defined Re and C pD using the aver- age velocity, and we have converted their values according to our definitions using the centerline velocity at the entry. 684.5. Conclusions 4.5 Conclusions In this work, we have implemented a rotational projection scheme to com- pute incompressible flows using SPH. The main difference from conventional projection schemes is a nonhomogeneous Neumann boundary condition for the pressure Poisson equation that is consistent with the momentum balance. This resolves a well-known difficulty in projection schemes wherein impos- ing an artificial boundary condition of vanishing normal pressure gradient produces pressure boundary layers at the inlet and outlet to flow domains and on solid surfaces. Applying the rotational projection scheme to two-dimensional channel flows and flow around a square cylinder, we have demonstrated that the new pressure boundary condition removes the spurious pressure boundary layers and markedly improves the accuracy of the solution. In particular, it predicts a drag coefficient on the cylinder within 7% of the correct value for Reynolds numbers up to 50, as compared with errors on the order of 20% incurred by the artificial pressure boundary condition. Our numerical experimentation has shown the rotational projection scheme as robust, accurate and efficient in dealing with open boundaries and flow around solid obstacles. Based on these findings, we recommend it as a supe- rior algorithm to the homogeneous projection scheme in computing truly in- compressible flows using SPH. This removes a shortcoming in the formalism that has at times limited SPH simulations to periodic boundary conditions and hampered accurate evaluation of the pressure and drag force on solid obstacles. Although the scheme needs to be tested in more complex flow situations, it has shown promise in extending the capacity of incompressible SPH to interesting and important problems previously inaccessible to SPH. 69Chapter 5 A 2D particle-based model for RBC 5.1 Introduction As explained in chapter 3, we see a unique potential for discrete particle models in simulating the biological phenomena on cellular level. In com- parison with the network, the particles offer far more freedom in modeling the internal structure of cells. For instance, clusters of particles having common properties that differ from the surroundings may represent various organelles, and the interaction among different particles may be dynami- cally evolved to reflect, say, the remodeling of the cytoskeleton triggered by internal or external factors as alluded to at the beginning. From a computa- tional standpoint, particle methods are meshless, and particularly suitable for simulating large deformation of soft matter. Thus, it seems worthwhile to develop the ideas of discrete particle models further, and potentially into a tool for studying the structural-property-disease coupling that sometimes comes under the name of nanobiomechanics [79]. This chapter represents a first step in that direction. We adapt the SPH picture of particles, which have overlapping regions of influence for smooth- ing and interpolation, by adding extensional and bending elasticity between the particles representing the RBC membrane. As indicated above, the dis- crete nature of the model allows one to go beyond the continuum framework to probe micro- and nanostructural responses to external stimulation. How- ever, we limit ourselves in this chapter to establishing a basic cell model without active reconfiguration of the microstructures. Thus, the following objectives will be pursued in this chapter: (a) to present a 2D particle-based model for the erythrocyte with a Newtonian-liquid cytoplasm and an elastic membrane having both extensional and bending moduli; (b) to validate this model by computing RBC motion and deformation in 2D shear and pressure- driven flows, and comparing the numerical predictions to experimental data and continuum-based computations in the literature. 705.2. Physical model and numerical scheme Figure 5.1: (a) Particle-based model for the red blood cell, with fluid particles representing the cytoplasm and the suspending plasm, and spring-connected solid particles representing the cell membrane. (b) Schematic of extensional and bending elasticity between particles in the membrane. 5.2 Physical model and numerical scheme Consider an RBC suspended in a Newtonian plasma. Similar to previous discrete-particle models, we deal only with 2D planar geometry here, and leave the extension to 3D to a future endeavor. We view the RBC as a cap- sule made of an elastic membrane enclosing a Newtonian cytoplasm. The plasma and the cytoplasm flow according to the Navier-Stokes equations. Thus each is discretized by “fluid particles” in the SPH sense, the two kinds having different properties tuned to reflect the density and viscosity of the two liquids. The membrane is replaced by a collection of “solid particles”, connected by elastic springs obeying a nonlinear spring law (Fig. 5.1). In addition, we introduce a separate bending elasticity that controls the varia- tion of the membrane curvature. Note that in 2D, shearing in the tangential plane of the membrane is excluded; membrane deformation is limited to stretching and compression tangential to the membrane and bending. It seems appropriate at this point to briefly review the ideas underlying the SPH method and clarify the relationship between these and our cell model. SPH originated in astrophysics some three decades ago [48, 88], and was later adapted for solving hydrodynamic problems [102]. In this form, SPH is mostly a numerical device for discretizing the Navier-Stokes equations in Lagrangian coordinates. The particles are interpolation points from which the fluid properties may be calculated by “smoothing” over neighboring particles. Although the particles have mass and move according 715.2. Physical model and numerical scheme to Newton’s law of motion, the forcing terms stem directly from discretizing the continuum governing equations. Solid boundaries are discretized by solid particles, and the fluid-solid interaction may be defined so as to satisfy the no-slip boundary condition. Being completely meshless, the method does not require connectivity data as do finite volume and finite element methods. Thus, SPH is convenient in dealing with complex flows exhibiting large deformations, and especially fluid-solid interactions with large interfacial deformation [67]. In our cell model of Fig. 5.1, the fluid components inside and outside the membrane are represented by classical SPH particles. For the mem- brane particles, however, we insert additional physics on the particle level in the form of elastic forces that on the cell level amount to the appropriate elastic properties of the membrane. Specifically, we introduce a set of non- linear springs between each pair of neighboring particles on the membrane to account for compression and stretching (Fig. 5.1b). The spring law, to be detailed in the following, reflects the areal conservation and strain-hardening nature of the RBC membrane [137]. In addition, we introduce a bending elasticity that penalizes deviations of the local curvature from that of the biconcave “resting shape” of RBCs. This way, the particles representing the membrane not only play their conventional role in SPH as interpolating points, but also carry new physics that define membrane elasticity. As men- tioned above, it is the latter role that makes the discrete-particle approach particularly suitable for modeling the coupling between cell mechanics and microstructural evolution inside the cell. 5.2.1 Constitutive models for the membrane The lipid-bilayer structure of the red blood cell membrane endows it with a very high modulus against areal expansion or contraction, so the surface area is essentially conserved [32, 137]. In the meantime, the membrane is very flexible with respect to shear deformation and bending, such that the RBC can deform readily to pass through narrow capillaries. The membrane is so thin compared with the cell size that it may be considered a 2D elastic shell. A constitutive equation for such a membrane can be obtained by adapting 3D elasticity or by postulating a 2D relationship directly. [4], among others, have compared the constitutive equations that have been proposed for RBC membrane. The 2D version of Hooke’s law is the simplest constitutive law; it assumes a linear dependence of tension on surface deformations, and is ap- plicable only for small deformations. The neo-Hookean and Mooney-Rivlin models are classical hyperelastic models for rubber-like materials. But their 725.2. Physical model and numerical scheme lack of membrane areal conservation and strain-softening behavior for large deformation make them inappropriate for RBC membrane modeling [4]. So far, the most successful constitutive model for RBC membranes ap- pears to be that of [137]. When a 2D elastic shell is subject to in-plane stretching, with extensional ratio λ1 and λ2 along the principal directions, [137] proposed the following strain-energy function: WS = ES 4 ( 1 2 I21 + I1 − I2 ) + ED 8 I 2 2 , (5.1) where ES and ED are the shear and dilatation modulus of the membrane, and the strain invariants I1 = λ21+λ22−2 and I2 = λ21λ22−1. The extensional stress components in the membrane follow from WS : T1 = λ1 λ2 [ ES 2 ( λ21 − 1 ) + ED 2 λ22 ( λ21λ22 − 1 )] , (5.2) T2 = λ2 λ1 [ ES 2 ( λ22 − 1 ) + ED 2 λ21 ( λ21λ22 − 1 )] . (5.3) Any changes in the area produces a deviation of λ1λ2 from unity, which will result in a large elastic tension due to the large magnitude of ED relative to ES . As a result, the membrane area is kept approximately constant during deformation. In addition, the Skalak constitutive model exhibits strain-hardening for large deformation [4]. Based on these features, we have decided to adopt the Skalak constitutive model for the extensional springs in our particle model. In our 2D simulation, the RBC membrane is a closed planar curve instead of a 2D curved surface. By putting the out-of-plane stress component T2 to zero, one obtains an expression of λ2 in terms of λ1. When plugged into Eq. (5.2), this yields the one-dimensional stress for stretching the membrane [4]: T = ES 2 λ1(λ21 − 1) √ 1 + Cλ41 1 + Cλ21 [ 1 + C(1 + Cλ21) (1 + Cλ41)2 ] , (5.4) where C = ED/ES , the ratio between the dilatation and shear moduli, is above 104 for healthy RBC [137]. For our extensional springs between the membrane particles, λ1 will be the ratio between the deformed and resting lengths. The bending elasticity is reflected by a resistance against deviation from the equilibrium membrane curvature corresponding to the biconcave natural 735.2. Physical model and numerical scheme shape of the RBC. Following [39], we describe this biconcave shape by the following equation: y¯ = 0.5 (1− x¯2)1/2(c0 + c1x¯2 + c2x¯4), −1 ≤ x¯ ≤ 1, (5.5) where c0 = 0.207, c1 = 2.002, c2 = 1.122, and the non-dimensional coordi- nates (x¯, y¯) are scaled by the radius of a human RBC a = 3.9 µm. Almost all prior continuum models for the RBC membrane have used a linear bending elasticity [3, 120, 167]: m = EB(κ− κ0), (5.6) where m is an external bending moment exerted on a infinitesimal segment of the membrane, EB is the bending modulus, and κ and κ0 are the local curvature in the deformed and resting states. The sign conventions for m and κ are that m is positive if it squeezes the outside of the membrane and stretches the inside, and κ is positive if the membrane is locally concave, with the center of the osculating circle lying outside the cell. This is the case at the center of the RBC, while near the edge κ < 0. To implement Eq. (6.3) for our particle-based membrane, κ will be computed from the position of neighboring particles and m must be converted to nodal forces. Details will be given below in subsection §5.2.5. 5.2.2 Governing equations for fluid flow In the SPH algorithm, incompressibility is approximated by a small arti- ficial compressibility. Thus, the governing equations may be written in a Lagrangian framework as: Dρ Dt = −ρ∇ · v, (5.7) Dv Dt = g + 1 ρ ∇ · τ − 1 ρ ∇p, (5.8) p = p0 [( ρ ρ0 )γ − 1 ] , (5.9) where t is time, g is the gravitational acceleration, p is pressure, v is the ve- locity vector, τ is the viscous stress tensor and D/Dt refers to the material derivative. Although body force will not be important for the computa- tions to be presented here, we will include g in discussing the algorithm for completeness. 745.2. Physical model and numerical scheme The artificial compressibility is a device for coupling the particle pressure to their motion [95]. The motion of particles, if not observing the constraint ∇ · v = 0, produces a variation in the particle density ρ through Eq. (5.7). In turn, this causes a pressure disturbance through the artificial equation of state (Eq. 5.9), which can then be used to correct the velocity field and make it solenoidal. In the equation of state [5], ρ0 and p0 are reference quantities, and the large exponent γ = 7 produces a strong pressure response to density variations and keeps the density variations negligibly small (below 1%), even at high Reynolds numbers [16, 102]. 5.2.3 Discretization using SPH The SPH method allows any function to be interpolated from its values at a set of discrete points—the SPH particles—using a kernel or weighting function W (r − r′, h), which specifies the contribution to any field variable at r by a particle at r′ that lies within 2h of r. The weighting function is normalized such that [95]∫ V W (r − r′, h)dr′ = 1, lim h→0 W (r − r′, h) = δ(r − r′), (5.10) V being the entire space. If a field variable A(r) is known only at a dis- crete set of particles r′j , then its value at any spatial location r can be approximated by: < Ah(r) > = ∫ V A(r′)W (r − r′, h) dr′ = N∑ j=1 A(r′j)W (r − r′j, h)∆Vj = N∑ j=1 mj ρj AjW (r − r′j , h), (5.11) where ∆Vj is the volume element at r′j , and has been replaced by the ratio between the mass and density of the jth particle: ∆Vj = mj/ρj . The summation is over all particles that lie within a circle of radius 2h centered at r, and Aj is a shorthand for A(r′j). The gradient ∇A(r) is evaluated through an integration by parts to transfer the gradient operator onto W [70]: ∇Ah(r) = N∑ j=1 mj ρj Aj∇W (r − r′j, h). (5.12) Note that calculating the spatial derivatives in SPH requires no mesh in- formation, and this gives a straightforward way to construct gradient of a 755.2. Physical model and numerical scheme function from its values at the SPH particles. In this chapter we adopt a popular spline-based kernel function W (r, h) = σhν × 1− 32s 2 + 34s 3 if 0 ≤ s < 1 1 4 (2− s)3 if 1 ≤ s < 2 0 if 2 ≤ s (5.13) where s = |r|/h, ν is the number of dimensions (ν = 2 here) and σ = 2/3, 10/7pi or 1/pi in one, two and three dimensions, respectively. This kernel has compact support so that its interactions are exactly zero for r > 2h. The second derivative of this kernel is continuous and the leading-order error in an interpolation is O(h2). Higher-order splines can be used, but they interact at larger distances and thus are computationally costlier. 5.2.4 Solution algorithm Through the interpolation operation outlined above, any partial differential equation can be discretized into ordinary differential equations governing the motion of SPH particles. In particular, the momentum equation (Eq. 5.8) may be discretized for fluid particles as [95]: Dvi Dt = g − ∑ j mj( piρ2i + pj ρ2j )∇iWij +Πij , (5.14) where Wij is a shorthand for W (ri−r′j , h) and ∇i designates the derivative with respect to ri. Πij represents the viscous stress term, and we employ the formula suggested by [102]: Πij = ∑ j mj (µi + µj)rij · ∇iWij ρiρj r2ij vij , (5.15) with vij = vi − vj , rij = ri − rj and rij = |rij| is the distance between particles i and j, whose viscosities µi and µj may differ if they represent different phases or fluid components. The momentum equation has three forcing terms on the right-hand-side: body force, the pressure gradient and the viscous force. These must be treated properly, along with the continuity equation and equation of state, to approximate incompressibility. At each time step, the governing equations are solved for each particle to update its position, velocity and pressure. The sequence in which the forcing terms are incorporated can differ from one algorithm to another. Here we use a fully explicit two-step algorithm 765.2. Physical model and numerical scheme [70]. In the first step, the momentum equation is solved with body force g and viscous force Π but not the pressure gradient. Thus, an intermediate velocity v˜ is generated and the particle positions are updated accordingly. Because this step is not subject to the incompressibility constraint, we expect it to perturb the density of some particles away from the reference value ρ0. The density variation can be computed directly from the updated particle position [95] ρ(r) =∑ j mj W (r − r′j , h). (5.16) This is equivalent, within the interpolation errors of the scheme, to the continuity equation Eq. (5.7) [95, 102]. In the second step of the algorithm, the pressure is calculated, using the equation of state (Eq. 5.9), from the perturbed ρ field. Thus, areas with denser particles have a larger ρ and a higher p that would tend to dis- perse them, and vice versa. The latter action is achieved through correcting the velocity field by solving the momentum equation again, with only the pressure gradient term on the right-hand-side: vˆi = ∑ j mj ( Pi ρ2i + Pj ρ2j ) ∇iWij. (5.17) For each particle, the velocity vi = v˜i + vˆi is taken to be the new velocity at the end of the time step, and it should be approximately divergence-free. Finally, the position of the particles are updated by a central differencing scheme: ri(t+∆t) = ri(t) + ∆t2 [vi(t) + vi(t+∆t)]. (5.18) The procedure is repeated for the next time step till a specified time is reached. Numerical stability of the explicit scheme puts a limit on the time step. The following criterion, due to [102], works well for our computations: ∆t ≤ 0.125ρh 2 µ . (5.19) In addition, SPH algorithms are susceptible to a well-known numerical insta- bility known as the tensile instability, whereby particles tend to form clumps and cause unrealistic fracture in the material when it is being stretched [100]. We suppress this instability by using the scheme of [96], which introduces a small repulsive force between nearby particles when they are in a state of tensile stress. 775.2. Physical model and numerical scheme Figure 5.2: Element bending groups (EBGs) used to convert the bending moment to pairs of forces acting on particles in the membrane. The line and arrow styles distinguish pairs of forces used to replace the bending moment on each line segment. 5.2.5 Bending moment In representing the bending elasticity in our particle-based model, the mo- ment m needs to be transformed to forces acting on the membrane particles. This is accomplished by using the element bending group (EBG) idea for elastic shells [168]. In our 2D membrane, the EBG for a membrane particle is made of two adjacent line segments connecting it to the two neighbor- ing particles. Thus, a membrane with n particles has n line segments and n overlapping EBGs. Figure 5.2 depicts the EBG centered at membrane particle P3 that involves 2 line segments connecting 3 membrane particles. At each time step the membrane curvature at each particle is calculated by passing a circle over three neighboring particles. For instance, the curvature at point P3 will be the inverse of the radius of the circle passing through points P2, P3 and P4. Equation (6.3) then gives the bending moment m3 at P3. In the EBG scheme, m3 acts on both line segments that meet at P3. For the segment between P2 and P3, m3 is replaced by a pair of equal and opposite forces: F32 = m3/r23 on P3 and −F32 on P2. For the segment between P3 and P4, similarly, m3 generates F34 on P3 and −F34 on P4. This amounts to two forces on P3. But P3 is also the end point of two other EBGs centered at P2 and P4. Thus, m2 and m4 will produce a force −F23 and −F43, respectively, on P3. In the end, the particle P3, and every other membrane particle, receives 4 nodal forces as a result of bending elasticity. 785.3. Results and analysis Finally, the equation of motion for the membrane particles is: Dvi Dt = g − ∑ j mj( piρ2i + pj ρ2j )∇iWij +Πij + 4∑ n=1 Fi,n + 2∑ n=1 Ti,n, (5.20) where the index n refers to neighboring particles on the membrane. This is similar to Eq. (5.14) except for the 4 bending-based nodal forces Fi,n and 2 extensional spring forces Ti,n along the line segments. In summary, the main idea underlying our model is the same as in prior models by [146] and [157], which is to produce realistic dynamics on the cell level by manipulating the physics on the particle level. However, our model represents advances beyond prior work in two key aspects: (i) Our model properly implements liquid-membrane interactions, in a way that is consistent with the SPH algorithm. The membrane particles interact with the liquid particles in the cytoplasm and the suspending plasma according to Eq. (5.20), with an inter-particle pressure that moderates the repulsion and attraction between particles. As such, the fluid motion and membrane deformation is fully coupled [155]. Moreover, since the membrane particles are held together by extensional elasticity, the particle pressure prevents the fluid particles from leaking through the membrane [106]. (ii) We use more realistic constitutive equations for the extensional and bending elasticity, as well as model parameter values that correspond to physiologically realistic values in the literature. As will be shown in the next section, this makes possible direct comparisons with experimental data as well as numerical re- sults based on continuum models. Such comparisons were not possible with the earlier particle models. For instance, [146] had to treat the membrane bending elasticity as an adjustable parameter when comparing with experi- ments. 5.3 Results and analysis In this section, we study two benchmark problems to demonstrate the ca- pability of our particle-based model: cell deformation in shear flows and Poiseuille flows. The cell deformation can be viewed as the outcome of the competition between viscous forces from the external flow and elastic resistance of the membrane, which is embodied by the dimensionless group G = µUm ES , (5.21) 795.3. Results and analysis where µ is the plasma viscosity, and the characteristic velocity Um is the mean velocity for a Poiseuille flow and ka for a simple shear at shear rate k, a being the radius of the RBC. G can be likened to the capillary number in drop dynamics [169, e.g.]. We can also define two ratios between the moduli: ˆEB = EB a2ES , C = ED ES . (5.22) Note that in our 2D simulation, ES and ED have the dimension of force over length, while EB has that of energy. The last independent dimensionless group is the Reynolds number Re = ρUma µ , (5.23) which is on the order of 10−4 for RBC motion in microcirculation as well as in the computations to be described below. Thus, inertia will be negligible. 5.3.1 Cell in shear flows Fischer et al. [46] observed the deformation of an RBC in simple shear in an experimental setup that keeps the center of the RBC stationary. For sufficiently high shear rates, the cell deforms from its biconcave rest shape to an ellipsoid and then to an elongated spindle-like shape oriented at an angle with the undisturbed flow direction. Meanwhile, the membrane and cytoplasm execute a rotating motion around the center of the cell, which is well known as “tank-treading”. To mimic that experimental setup, we place an RBC in the center of an 8a × 4a rectangular domain with top and bottom walls moving in opposite directions. Periodic boundary conditions are imposed at the left and right boundaries such that particles exiting from one end will emerge from the other. The simulations to be presented use some 12,000 SPH particles in the domain, and 96 particles on the membrane. In this and previous calculations [70], we have confirmed that the spatial resolution is adequate; doubling the number of particles causes a change in the result on the order of 1%. The computations have used model parameters from physiological data for real human RBC, with ES = 5.0 × 10−6 N/m, ED = 0.1 N/m and EB = 2.0 × 10−19 N·m [136]. The large ratio of C = 2 × 104 ensures that little surface areal dilation occurs and the membrane area is essentially conserved. We have examined differing viscosities for the cytoplasm and the suspending fluid, but the results presented will have equal viscosity for the two liquids: µ = 6.0 × 10−3 Pa·s. These are the baseline parameters for 805.3. Results and analysis Figure 5.3: Deformation of an RBC in shear flow. C = 2 × 104, G = 0.234 and ˆEB = 2.63 × 10−3. The snapshots are for dimensionless times kt = 0, 2, 3, 4 and 10, and the last frame depicts the flow field surrounding the cell in the steady state. the simulations. We have also systematically varied G and ˆEB to probe the effects of the membrane elasticity; in these exercises C = 2 × 104 is kept constant. Figure 5.3 depicts a typical simulation by a sequence of snapshots of the particles in the domain. One of the membrane particles is drawn in larger size to serve as a marker to illustrate the tank treading motion. In the fi- nal frame, which is essentially steady state, the cell has an aspect ratio of 6.22 and assumes an angle of θ = 13◦ with respect to the far-field flow, and exhibits the tank-treading motion. Note that the cell deforms on the flow time scale 1/k, rather than on a time scale defined by membrane elastic- ity, say µa3/EB . This is because at G/ ˆEB = 90, the bending elasticity is overwhelmed by the flow. Indeed, the final shape of the cell is an almost symmetric cigar shape; the native curvature of the membrane is barely man- ifested. The steady-state circumference of the cell has increased 3.6% from the rest state. The 2D version of the Skalak constitutive law (Eq. 5.4) as- sumes that in-plane extension is accompanied by out-of-plane contraction. Thus areal conservation does not imply constant circumference of the cell in the plane. The steady-state configuration of the cell is sensitive to ˆEB and G. In- creasing ˆEB amounts to stronger resistance to bending of the membrane. As a result, the deformed shape of the RBC bears a more distinct signature 815.3. Results and analysis Figure 5.4: Effect of increasing the bending elasticity on cell deformation. ˆEB = 1.31 × 10−2 (first row), 2.62 × 10−2 (second row) and 5.26 × 10−2 (third row) correspond to a bending elasticity 5, 10 and 20 times that of the real RBC. The three columns are for kt = 3, 5 and 7. G and C have the same value as in Fig. 5.3. of the biconcave resting shape (Fig. 5.4). In fact, for higher ˆEB values, the tank-treading amounts to a periodic rather than steady solution; the high curvature at the edge of the undeformed RBC convects around the mem- brane. At sufficiently large ˆEB (e.g., ˆEB = 0.262), tank-treading can no longer be achieved; instead the cell tumbles as if it were a rigid particle. In comparison, the effect of G is less spectacular. Increasing G while keeping all other parameters unchanged tends to increase the aspect of the elongated cell, and decrease its angle of inclination (Fig. 5.5). Both are rather mild effects as G varies by a factor of 25. [167] computed the behavior of an RBC in simple shear using a con- tinuum model. The membrane is treated as a neo-Hookean viscoelastic material, and the fluid-membrane interaction is accounted for by the im- mersed boundary method. Although the membrane constitutive equations differ between their study and ours, both used measured membrane prop- 825.3. Results and analysis Figure 5.5: The steady shape of the RBC gets more elongated with increasing G values. All other dimensionless parameters have the same value as in Fig. 5.3. Note the quantitative agreement with the numerical result of [167] at G = 0.234, who used the immersed boundary method and a continuum model for the cell. erties to determine the model parameters. Thus, there is close agreement between the two studies. For instance, the evolution of the cell shape and the tank-treading motion in Fig. 5.3 are essentially identical to the predictions of [167]. The steady-state cell shape, computed for identical dimensionless parameters, is in quantitative agreement between the two (Fig. 5.5). This serves as a validation of our particle-based model as well as our SPH algo- rithm. Another feature of the simulation that can be quantitatively compared with previous work is the tank-treading frequency f , defined as the inverse of the period of tank treading. The data of [153] and [45] show the frequency f , scaled by the shear rate k, to be in the range f/k = 0.02 – 0.038. Our results show a higher frequency; for example, f/k = 0.163 for the conditions in Fig. 5.3. An obvious cause of this discrepancy is the two-dimensionality in our simulation. In a Stokes flow, a 2D solid cylinder rotates with a period of 2/k while a sphere has a much longer period 4pi/k [19, 119]. The factor of 2pi is roughly the difference in f/k between our 2D computation and 3D measurements. In addition, our representation of the membrane is simplistic. Membrane viscosity may have also been a factor. 835.3. Results and analysis 5.3.2 Cell in Poiseuille flows Pressure-drive flow in a tube is a close analogy for blood flow in capillaries. Thus the deformation of an RBC in Poiseuille flow not only constitutes a benchmark problem for our model and numerical method, but has direct relevance to microcirculation. It is the latter connection that has motivated numerous previous studies. Perhaps the most profound discovery is that the red blood cell deforms into a characteristic parachute shape in order to traverse capillaries smaller than its undeformed diameter [135]. Numeri- cal computations have mostly employed the continuum representation of an elastic shell enclosing a viscous liquid. [166] published one of the earliest finite-element computations of RBC deformation. [127] modeled RBC de- formation in a capillary using lubrication theory, for both axisymmetric and fully 3D geometries. The amount of deformation is studied as a function of the tube diameter and the flow velocity. [121] computed the axisymmetric motion of a file of red blood cells in the Stokes regime using the boundary- integral method. The cell membrane was modeled by an elastic shell obeying the Skalak constitutive law [137]. Against this backdrop, we test our particle-based model by computing the deformation of an initially biconcave RBC in a pressure-driven flow in a 2D channel. The channel width is 2.15a, 7.7% wider than the diameter of the undeformed RBC, and its length is 10a. The RBC is initially placed in the middle of the channel with its broad side facing the flow direction (Fig. 5.8), and a parabolic velocity field is imposed initially throughout the domain. For the two ends of the domain, periodic conditions are used such that particle exiting one end re-enters the other. The number of SPH particles is comparable with that for the shear-flow computations, and we adopt similar physical parameters based on experimental measurements of cell membrane moduli. For simplicity, we again assume equal viscosity for the cytoplasm and the suspending liquid. Figure 5.6 depicts the deformation of an RBC in a capillary, with the dimensionless parameters being C = 2 × 104, G = 0.024 and ˆEB = 1.97 × 10−3. These values have been chosen to match those of [127], who assumed a somewhat smaller bending modulus EB = 1.5 × 10−19 N·m. The higher flow velocity in the middle of the channel causes the center of the cell to bulge forward. This counters and quickly overcomes the native curvature of the membrane in the front (Umt/a = 3.85) and eventually leads to the 2D version of the parachute shape in steady state. Note that most of the deformation occurs within the first few cell radii that the cell center travels. Similar to the simulation in Fig. 5.3, the bending elasticity is weak and the 845.3. Results and analysis Figure 5.6: Deformation of a red blood cell in a Poiseuille flow. The flow is from left to right, and the five snapshots are taken at dimensionless time Umt/a = 0, 1.28, 3.85, 6.40 and 8.95, the last approaching the steady state. Despite the considerable change in shape, the circumference of the cell has increased a mere 1.05% relative to the undeformed state. cell deforms with the fluid more or less affinely. Since the parameters G and ˆEB were chosen to be consistent with real RBC properties and typical flow conditions in the capillary, the shape evolution in Fig. 5.6 being in general agreement with in vivo observations [135, e.g.] is encouraging. One may also notice the slight top-bottom asymmetry in the cell shape. This is a numerical artifact due to updating the particle positions sequentially in an explicit scheme, and diminishes with reduced time step. A more quantitative comparison is done in Fig. 5.7 with the numerical result of [127]. Both studies used the same dimensionless parameters G, ˆEB , C and the same cell-to-tube diameter ratio. The steady-state RBC shape is in close agreement between the two. However, Secomb’s cell appears to enclose more area than ours. This is a geometric effect. His calculation was for an axisymmetric cell while ours is planar. Thus, the conservation of membrane area corresponds to differing contour lengths in the plane. Besides, the membrane curvature is computed differently; the axisymmetric membrane has two principal curvature while our planar contour has one. Thus, the bending elasticity is represented differently in the two geometries, even though the material parameters are matched. Finally, [127] used the viscoelastic Kelvin model for the membrane while we used the nonlinear elastic Skalak model. These factors have given rise to the small differences in Fig. 5.7. We have also varied G and ˆEB from physiologically based values to explore their effects on the deformed RBC shape. As expected, greater G 855.3. Results and analysis Figure 5.7: Comparison of the steady-state cell shape in a Poiseuille flow between our simulation (black dots) and the result of [127] (continuous curve). (a) (b) Figure 5.8: Effects of (a) G and (b) ˆEB on the steady state deformation of the RBC in a Poiseuille flow. produces more deformation, with the parachute taking on a deeper dome shape (Fig. 5.8a). Since inertia is negligible, changing G can be visualized 865.3. Results and analysis Figure 5.9: Deformation of a red blood cell as it goes through a contraction in the channel. The flow is from left to right. The width ratio between the narrow and wide channels is 0.4. Defined using the mean velocity in the wide channel, all flow and membrane parameters are the same as in Fig. 5.6. The snapshots are at dimensionless time Umt/a = 0, 0.641, 1.60, 1.92, 2.24 and 2.88. as changing the velocity and viscous shearing in the external flow, the direct cause of cell deformation. The effect of ˆEB is subtler. While smaller bending modulus leads to more pointed ends at the edge of the cell, the center also becomes thicker, and in a sense less deformed (Fig. 5.8b). In the undeformed resting shape, the cell has a “dimple” at the center. As the bending modulus is weakened, it becomes easier to override the innate concavity at the dimple. Thus, as the stronger shear near the channel walls produces narrow and pointed edges, the cytoplasm is squeezed into the central part and causes it to swell. Finally, Fig. 5.9 shows a simulation of an RBC entering a contraction in the capillary. This geometry is relevant to both the microcirculatory network and recent studies of cell mechanics in microfluidic channels [133, 162]. The channel width is 2.15a in the upstream portion and 0.86a in the narrower part downstream. The extensional flow at the contraction deforms the cell 875.4. Summary into a boomerang. As the RBC enters the narrow tube, the two wings fold toward each other, increasing the bending of the membrane at the crotch in between. This causes the crotch to move back, shortening the wings and lengthening the front part of the RBC as it settles into a steady shape in the narrower tube. This steady shape bears a close resemblance to in vivo observations [128, e.g.]. Note that from the third frame (Umt/a = 1.60) onward, the thin gap between the RBC and the corner or inner wall of the channel contains only one layer of particles. Thus the local details of the flow are not adequately resolved in the computation. 5.4 Summary We have proposed a particle-based model for the red blood cell as the basis for developing more general discrete-particle microstructural models for var- ious cell types. The main goals of this chapter are to describe the model and present solutions on benchmark problems as validations. The cell membrane is represented by particles connected by nonlinear springs. A linear bend- ing elasticity is implemented by using the local curvature. The inner and outer liquids are discretized by particles in the standard smoothed-particle- hydrodynamics (SPH) procedure. Thus, the model is an adaptation of SPH ideas; the particle-level physics is designed so as to produce the proper mem- brane elasticity. The particles are a numerical device for solving the par- tial differential equations as well as a vehicle for incorporating microscopic physics. The model predictions of tank-treading in shear and the parachute shape in Poiseuille flows are in excellent agreement with experimental data and prior continuum-based computations. The model parameters are determined according to physiological measurements of cell viscosity and membrane elas- ticity, and no curve-fitting is involved. This agreement provides support for both the particle-based RBC model and the SPH-based numerical algorithm. However, it is necessary to mention the limitations imposed by our 2D model. Approximating the 3D RBC with a 2D model is equivalent to as- suming an infinite depth in the perpendicular direction. The conservation of RBC surface area imposed by Skalak law for a 3D RBC is not equivalent to conservation of circumference of 2D RBC membrane. The 2D model cannot account for the stress and strain in the perpendicular direction that may arise in the real 3D RBC. Nevertheless, The model does capture the most important qualitative dynamics of the 3D cell. 88Chapter 6 How malaria merozoites reduce the deformability of infected RBC 6.1 Introduction The most important changes in mechanical properties of the iRBC are in- creased tendency to adhere to the vascular walls and to other cells, and progressive loss of deformability. The latter, in particular, leads to capillary occlusions and serious anemia, coma and even death [20, 132]. We did a series of simulation to correlate the rigidity of the infected red cell to its structural components, especially the merozoites. The details of this study have been described in this chapter. The prevailing thinking in the literature is that the higher rigidity of the iRBC is caused by remodeling of the RBC membrane by parasite-derived proteins [89, 90, 92, 103, 126, e.g.]. More recently, network [8, 81], continuum [142] and dissipative particle dynamics models [41] have been developed to compute the loss in overall deformability of iRBCs when stretched by optical tweezers. This is expressed in terms of an elevated effective shear modulus for the membrane Gs. In the late stages of infection, Es is estimated to be between 40 and 60 µN/m [25, 41, 72, 141, 142], up from roughly 6 µN/m for healthy RBCs [34, 107]. When the membrane elasticity is probed locally by micropipette aspi- ration of small membrane tongues, Nash et al. [107] recorded much less increase of the membrane modulus due to infection by P. falciparum. In the late schizont stage Gs is found to be in the range of 9 – 14 µN/m. This is too low to account for the reduced deformation of the whole cell measured by stretching. Aside from membrane rigidification, Rathod and coworkers [1, 61] noted three other factors that may influence whole-cell deformability: the shape of the cell (surface-volume ratio), the cytoplasmic viscosity and the membrane viscoelasticity. Conceptually, the last factor encompasses 896.2. Physical model and numerical scheme Gs if elasticity is seen as one ingredient of the membrane viscoelasticity. Moreover, Nash et al. [107] argued that “the main factor in the drastic loss of deformability of the trophozoites and schizonts was the presence of the large, very resistant parasite itself.” If proved to be true, this would explain the discrepancy between the relatively modest increase in membrane rigid- ity and the marked loss of whole-cell deformability. In light of the above speculations, it would be desirable to probe the effects of each factors sepa- rately. However, this task presents great experimental difficulties, and Antia et al. [1] noted that “most studies that measure changes in RBC deformabil- ity are unable to distinguish between the underlying structural factors that contribute to these changes”. This chapter describes our attempt to use numerical simulations to quan- tify the relative importance of the membrane elasticity and the merozoites in iRBCs. The other two parameter, cell shape and cytoplasmic viscosity, will be chosen based on measured data and remain fixed, even though these data often contain large variations. In a smoothed-particle-hydrodynamics (SPH) framework, we represent the elastic membrane, cytosol, merozoites and suspending plasma by discrete particles. Then we simulate the stretch- ing test of RBCs using optical tweezers [24, 85] over a range of Gs values and merozoite sizes and spatial arrangements. The main finding is that the presence of a sizeable merozoite or cluster of merozoites greatly reduces the ability of the iRBC to deform under stretching. This confirms the supposi- tion of Nash et al. [107]. A prescribed amount of loss of deformability can be realized by different combinations of membrane modulus and merozoite size. Neglecting the merozoite contribution will lead to gross overestimation of the membrane modulus. In this sense, the conventional way of interpreting cell-stretching data is flawed. 6.2 Physical model and numerical scheme Mechanical models for biological cells typically involve dynamic coupling between viscous fluids and soft solids. Thanks to the merozoites inside the red cell, the problem at hand also features multiple solid domains that may undergo large displacement and deformation. This makes a discrete physical model attractive. As examples of previous work on such discrete models, Discher et al. [8, 26] developed a spectrin network model for the erythrocyte cytoskeleton, and simulated micropipette aspirations. Li et al. [80] used a similar model for stretching of red cells by optical tweezers. More recently, Karniadakis and coworkers [43, 118] coarse-screened the network 906.2. Physical model and numerical scheme Figure 6.1: Discretization of a healthy RBC by SPH particles. The membrane is represented by elastic springs connecting neighboring particles into a more or less uniform network. The total number of particles in this case is 13500. The RBC is centered at (0, 0, 0) and the coordinates are in unit of meters. The arrows represent the stretching force applied by optical tweezers. model of the membrane and coupled it to fluid flow by the dissipative particle dynamics (DPD) method. This allowed them to simulate deformation and motion of red cells subject to complex flow fields, even with solid merozoites inside [41, 44]. Yamaguchi et al. [71, 156] and Hosseini and Feng [68] have adopted a similar approach with the fluid flow being solved by smoothed particle hydrodynamics (SPH) method. In this work, we adopt a 3D particle-based model for the P. falciparum- infected RBC, which is a generalization of our previous 2D SPH model [68]. The RBC is represented by an elastic membrane enclosing a Newtonian viscous cytosol and possibly merozoites. Each component is discretized by particles. The membrane consists of a more or less uniformly distributed set of particles connected by elastic springs (Fig. 6.1). The cytosol is dis- cretized by classical SPH particles that flow according to the Navier-Stokes equations. The merozoite is a stiff particle that resides inside a parasitopho- rus vacuole (PV) [58]. For simplicity we replace the PV-merozoite complex 916.2. Physical model and numerical scheme by a single solid object. Numerically this is realized by constraining a set of SPH particles to move as a rigid body [100]. When simulating the stretching experiments [24, 85, 93, 142], we ignore the viscous flow of the ambient fluid because it plays no role in the final stretching. A constant ambient pressure is applied on the RBC membrane. Overall, our model is similar in spirit to those of Yamaguchi et al. [71, 156] and Karniadakis et al. [44, 118]. But there are important differences in the treatment of the elasticity of the membrane, as will be explained shortly. In our computational scheme, the membrane and merozoite particles not only play their conventional role in SPH as interpolating points, but also carry new physics that define membrane elasticity and merozoite rigidity. Compared with continuum models for cell mechanics, e.g. [24, 84, 142], the discrete particle-based model offers far more freedom in modelling the internal structure of cells and in representing the interplay between such intracellular microstructures and the mechanical behavior of the whole cell [68]. 6.2.1 Elasticity of the RBC membrane In our physical model, the RBC membrane includes the bilipid layer as well as the underlying spectrin network. Therefore, a general viscoelastic model should be used [43]. In the current work, however, we focus on steady- state stretching under a fixed load. Besides, the viscosity of the membrane does not vary significantly with the progress of disease [107]. Therefore, we have neglected membrane viscosity and treated the membrane as purely elastic. Specifically, we endow the membrane with elasticity against in-plane deformation and bending. We start with an analytical formula for the surface of a biconcave healthy red blood cell [37]: z = ± R0 2 [ 1− x2 + y2 R20 ]1/2 [ c0 + c1 x2 + y2 R20 + c2 ( x2 + y2 R20 )2] , (6.1) where R0 = 3.91 µm is the cell radius on the central plane, c0 = 0.207161, c1 = 2.002558 and c2 = −1.122762. This surface is discretized into roughly equilateral triangles, with a prescribed mesh size matching that in discretiz- ing the bulk fluids by SPH particles. This conformation is inspired by the spectrin network and has been used in almost all prior network models [8, 41, 56, 71, 81, e.g.]. We deploy SPH particles on the nodes of the mem- brane, and connect each pair of neighbors by a linear spring. The resting 926.2. Physical model and numerical scheme length of the springs equals their natural length on the undeformed RBC surface. When the RBC deforms, the springs will develop a tension f ij = k(rij − r0,ij)rijrij , (6.2) where rij = |rj − ri| is the length of the spring connecting particles i and j, r0,ij is its resting length and k is the spring constant. Note that the same k is assigned to all springs though their r0,ij may differ slightly. For bending elasticity, we have adopted the Helfrich bending energy expressed in terms of the local Gaussian curvature K and mean curvature H [130]. In this model, the bending resistance of the RBC membrane results in a surface force density [159]: τb = −kb[(2H − C0)(2H2 − 2K +C0H) + 2∇2sH]n, (6.3) where kb is the bending modulus, C0 is the spontaneous curvature and n is the normal vector. ∇2s is the surface Laplacian. The surface quantities are computed on the discrete triangular mesh using the scheme of Meyer et al. [91]. For the nearly uniform triangular meshwork on the undeformed RBC, this yields highly accurate results for the surface curvatures and bending resistance. When the RBC is stretched, so do the surface triangles, and the accuracy tends to decline. This puts a limit on the amount of stretching that can be computed with confidence; the results to be presented are for stretching forces up to 100 or 150 pN. It is interesting to note the connections and differences between our model and previous ones. The correspondence between elastic spring net- works and continuum hyperelastic shells has been examined thoroughly before [23, 43, 112, e.g.]. In our model, the surface tesselation produces largely uniform equilateral triangles, although to cover the curved surface of Eq. (6.1), some links are slightly shorter or longer. The overall elasticity of such a network, therefore, is expected to approximate that of a perfectly uni- form triangular network. For small strains, the shear modulus Gs, Young’s modulus Es and area dilatation modulus Ks are all proportional to the spring constant k: Es = 2/√3k, Gs = √3/4k, Ks = √3/2k, with a Poisson ratio of 1/3 [112]. For large strains, we have calculated the elasticity of the network following the scheme of Dao et al. [23]. Both the Young’s modulus and shear modulus show moderate strain-hardening, similar to the Skalak model [112, 137]. Of prior discrete models for the RBC membrane, the spectrin network models [8, 26, 43, 81, 118] and particle-based models [68, 71] have been most 936.2. Physical model and numerical scheme successful. In terms of in-plane elasticity, the spectrin network models strive to account for the characteristics of the spectrin polymer and the morphology of the network. A nonlinear wormlike chain model is used, which is more sophisticated than our linear spring. The particle model of Imai et al. [71] uses linear springs that connect SPH particles computed using the moving particle semi-implicit (MPS) method. A common feature of all these is that the undeformed biconcave shape is chosen as the “reference state” in which the in-plane elasticity is entirely relaxed. The bending elasticity is a more subtle property to represent using a discrete model. The spectrin network models define a bending energy based on the angle θ between pairs of triangles [43, 81]: Vb = kb[1− cos(θ − θ0)], where kb is a bending modulus and θ0 is the spontaneous angle. Similarly, Imai et al. [71] represent the membrane’s resistance to bending in terms of a force: Fb = kb tan(θ/2), which acts along the normal of each triangle. In both approaches, the angle between discrete elements of the membrane depends on the surface tesselation and would become smaller with finer discretization. Thus, kb is not a physically meaningful quantity, and must be modified when the network is coarsened or refined [43, 118]. Imai et al. [71] used kb as an adjustable parameter to produce the best fitting to data for a particular triangulation. Ultimately, the stretching of RBCs is dominated by the in-plane moduli, and bending elasticity has little effect [81]. But it does affect the equilibrium shape of the RBC as an energy minimizer. Since the Fb of Imai et al. [71] vanishes only for a flat sheet, the membrane will not assume the correct biconcave resting shape [81]. Finally, these existing models have incorporated an energy to penalize variation of the surface area of the RBC membrane [26, 43, 71]. The spectrin network models also employ a volume-conservation energy. In our model, in- compressibility of the cytosol ensures conservation of the cell volume, and no additional constraint is needed. Areal incompressibility is an important fea- ture of the RBC membrane. In biological terms, the area-incompressibility of the lipid bilayer is a local property. To impose an elastic energy penalty against change in the total area seems unphysical. It would allow the pos- sibility of compensating for the expansion of one region by the contraction of another, potentially far-away region. Thus, we decided not to include an additional area-elastic energy. For the parameters we have used here, the maximum area increase is 12% for the largest stretching forces. To im- prove this, one approach is to use nonlinear springs that is stiff against area dilatation. This has been attempted in our 2D model [68] and in 3D cap- sule simulations [4] where the springs are designed to produce an effective response of the membrane described by the Skalak model [137]. 946.2. Physical model and numerical scheme 6.2.2 Numerical method SPH originated in astrophysics some three decades ago, and was later adapted for simulation of incompressible hydrodynamic problems [102]. In this con- text, the particles serve mostly as a numerical device for discretizing the Navier-Stokes equations in Lagrangian coordinates. One of the attractive features of SPH is convenience in dealing with complex fluid-solid interac- tions with large interfacial deformation [67]. This strength is a direct result of its being completely meshless; there is no need to generate connectivity data as in finite volume and finite element methods. Thus, SPH is well suited to the problem at hand. The incompressibility condition can be handled in one of two ways. The truly incompressible algorithm follows the traditional projection method and uses a Poisson equation to compute the pressure field needed to make the final velocity solenoidal [69, 78]. The weakly compressible algorithm, on the other hand, uses an artificial equation of state to compute a pressure from density fluctuations [100]. The density variation can be controlled to within 1% if a numerical speed of sound is chosen to be an order of magnitude higher than the maximum flow velocity [96, 102]. Our current problem involves the slow stretching of red blood cell with vanishing Reynolds number, and the weakly compressible scheme is used for its convenience and simplicity. In a Lagrangian framework, the governing equations are written as fol- lows: Dρ dt = −ρ∇ · v, (6.4) Dv Dt = g + 1 ρ ∇ · τ − 1 ρ ∇p+ fm, (6.5) p p0 = ( ρ ρ0 )γ − 1, (6.6) where t is time, g is the body force, p is pressure, v is the velocity vector, τ is the viscous stress tensor, D/Dt refers to the material derivative and fm is the membrane force per unit mass, which is the sum of shearing and bending forces for membrane particles and is zero for fluid particles. In addition, particles representing the parasite are constrained to move as a rigid body according to the total hydrodynamic force on them. These special particles do not follow Eq. (6.5). Generalities of the SPH method can be found in previous reviews [100, e.g.]. Essentially we deploy a large number of particles in the domain and have each carry the dynamic properties of the field and move in a Lagrangian 956.2. Physical model and numerical scheme framework, their motion being determined by ordinary differential equations that result from discretizing the governing partial differential equations. For the ith particle, its continuity and momentum equations are discretized as follows: Dρi dt = ∑ j mjvij∇iWij, (6.7) Dvi Dt = g − ∑ j mj ( pi ρ2i + pj ρ2j ) ∇iWij +Πij + fm,i, (6.8) where the index j refers to neighbors of i within a compact support, m is the particle’s mass, vij = vi − vj is the relative velocity, and Wij = W (ri − rj) is the weighting function [100].The operator ∇i designates the derivatives with respect to ri. Πij represents the viscous stress term, and we employ the formula of Morris et al. [102]: Πij = ∑ j mj (µi + µj)rij · ∇iWij ρiρjr2ij vij , (6.9) where rij = |ri − rj | is the distance between particles i and j, whose vis- cosities µi and µj may differ if they represent different phases or fluid com- ponents. We advance the above equations explicitly in time, and in each time step solve the continuity equation, the momentum equation and the equation of state by a two-step prediction-correction algorithm [95, 96]. The two-step algorithm strives to restore the initial density of each particle and approx- imate an incompressible fluid flow. In the prediction step, an intermediate velocity v˜ is calculated from the momentum equation with the body force, viscous force and membrane force on the right hand side. Then the particle positions are updated according to v˜, and an intermediate density is com- puted for each particle from Eq. (6.7) using the new particle position. Since the pressure gradient and the incompressibility constraint are disregarded in this step, we expect the particle density to deviate from the initial constant density. In the correction step, a pressure is calculated from the perturbed ρ field using the equation of state Eq (6.6). Thus particles in denser areas will have a higher ρ and hence a higher p, and vice versa. This pressure gradient is used to compute a velocity correction vˆi = −∆t ∑ j mj ( pi ρ2i + pj ρ2j ) ∇iWij, (6.10) 966.2. Physical model and numerical scheme Number of particles Fi n a ld efo rm a tio n o fm a jor a xi s (% ) 0 5000 10000 15000 2000038 38.5 39 39.5 40 40.5 41 41.5 42 42.5 Figure 6.2: Convergence test: elongation of the major axis of a model RBC as a function of the total number of SPH particles. The RBC contains no merozoites inside, and has a membrane shear modulus Gs = 15 µN/m and a bending modulus kb = 6× 10−19 J. The stretching force F = 100 pN. which tends to disperse or aggregate particles from areas of higher or lower densities, respectively. Thus, the new velocity v = v˜ + vˆ at the end of the time step will be approximately divergence-free. Finally, the position of the particles is updated by a central differencing scheme: ri(t+∆t) = ri(t) + ∆t2 [vi(t) + vi(t+∆t)] . (6.11) Our 3D hydrodynamic solver is based on the open-source software SPHysics, v2.2 [22], and we have added the treatment of the cell membrane and float- ing merozoites inside the iRBC. To determine the required spatial resolution for accurate results, we have run a series of stretching simulations with in- 976.3. Results and analysis creasing numbers of particles. Figure 6.2 plots the steady-state elongation of the major axis as a function of the total number of solid and fluid parti- cles. Somewhat surprisingly, even a very coarse resolution with 361 particles captures the final deformation within 7%. This may be because the final deformation does not depend on the dynamics of the flow and is thus an easy quantity to compute. All results to be presented below have been computed using 13540 particles, with 2930 on the membrane. This corresponds to an inter-particle separation that is roughly 0.051R0. 6.3 Results and analysis Stretching cells by optical tweezers is a popular method for probing the deformability of the cells and extracting their elastic properties. Henon et al. [59] was among the first to use this method. In their experiments, two silica beads attached to opposite ends of the cell are trapped by laser beams and moved in opposite directions. More recent experiments [24, 85, 93, 142] stretched cancer cells, healthy RBCs and iRBCs in the ring, trophozoite and schizont stages of malaria infection. The measured deformation is fitted to computations to yield an effective membrane modulus. We simulate the application of the stretching force by the beads in a similar way to that by Karniadakis and coworkers [43, 118]. We take 2% of the membrane particles with the smallest y-coordinates and 2% with the largest y and mark them as the “contact particles” between the RBC and the beads. The stretching force, evenly divided among the particles on each end, is then applied to each contact particle (cf. Fig. 6.1). The stretching of the cell is simulated until a steady state is achieved, when the stretching force is balanced by the elastic stresses developed among the membrane particles. As a healthy RBC gets infected and undergoes the ring, trophozoite and schizont stages, the following geometric and physical attributes change accordingly: the shape, surface area (S) and volume (V ) of the whole cell, the shape, size and number of merozoites inside the cell, and the membrane elasticity. In the following simulations, we have chosen a cell shape, S and V based on experimental observations of each stage, and systematically examined how the membrane shear modulus Gs and merozoite size influence the deformability of the cell. The bending elasticity kb may vary as well, but we have found little data for iRBC in the literature. For healthy human RBC, reported kb values range from 10−19 to 9× 10−19 J [38, 80, 138]. We have fixed kb = 6× 10−19 J in all the simulations. Similarly, we do not have 986.3. Results and analysis data on how the cytosol density and viscosity vary with disease progression, and fixed them at 1050 kg/m3 and 4.1 mPa·s for all the simulations presented below. These are values cited for healthy RBCs in the literature [12, 77]. 6.3.1 Stretching of healthy RBC As a baseline for later simulations, we have first computed the stretching of a healthy RBC from the equilibrium shape given by Eq. (6.1). Engelhardt and Sachmann [34] measured a membrane shear modulus of Gs = 6.1 µN/m for healthy RBC. Nash et al. [107] reported values between 4 and 8 µN/m for clinical and cultured samples. Therefore it seems reasonable to take Gs = 6 µN/m. The steady-state longitudinal and transverse diameters of the stretched cell, as a function of the stretching force, are plotted in Fig. 6.3. Prior experimental data and two computations are also shown for comparison. Our numerical prediction of the longitudinal diameter L agrees closely with experimental data. But the shortening of the transverse diameter (width) W is considerably underpredicted. This discrepancy may be due to simplifications in our physical model and improper parameter values. For example, the elastic springs are assumed linear, and there is no rigid constraint on area conservation. In our computation, the area increases by 12% for the largest stretching force F = 100 pN. Healthy RBCs exhibit considerable variations in the kb and Gs values, and the above computation has used the mean values. Moreover, it is illuminating to compare our simulations with that of Imai et al. [71] using SPH and that of Fedosov et al. [41] using DPD. Imai et al.’s model differs from ours in the bending elasticity and the inclusion of a total areal conservation term. Their membrane elastic constants for shearing and bending have been chosen for closest fitting of the experimental data. On the other hand, our parameters are from measurements; it is notable that the two computations agree so closely. The prediction of Fedosov et al. [41] agrees better with experiments than ours. Their model is more sophisticated in several aspects. The elastic spring uses a nonlinear wormlike-chain potential, and they employ additional con- straints on the conservation of cell volume and surface area. Their effective shear modulus Gs = 6.3 µN/m is close to ours, but their bending modulus kb = 2.4× 10−19 J is much smaller than ours. Conceivably weaker bending elasticity allows for larger curvature in the transverse plane. 996.3. Results and analysis X X X X X X X X X Force (pN) (D - D 0 )/D 0 (% ) 0 20 40 60 80 100 -40 -20 0 20 40 60 80 100 120 140 Experimental Fedosov et al. Imai et al. Our results X Figure 6.3: Variation of the longitudinal and transverse diameters of the healthy RBC as a function of the stretching force. The results of our numerical simulation are compared with experimental data [142] and prior computations of Imai et al. [71] and Fedosov et al. [41]. 6.3.2 Stretching of RBC in the ring stage Upon penetrating an RBC, the malaria merozoite flattens into a thin disc. By the mid-ring stage (10–15 hours after initial invasion), it has a diameter of 2–3 µm and a thickness of about 0.5 µm. The shape and size of the in- fected RBC are little changed from those of healthy RBCs in the ring stage [58]. Some 5–15 % of the host volume is occupied by the parasitophorus vacuole (PV) containing the parasite [131]. As mentioned before, we neglect the PV in the simulations and assign its volume to the solid disc represent- ing the merozoite. The initial condition for the stretching simulation is the 1006.3. Results and analysis (a) (b) (c) Figure 6.4: Three of the configurations of the merozoite used in our simulations. The top row shows the top view, and the bottom the corresponding side view. equilibrium biconcave profile of the healthy RBC (Fig. 6.1). On the shear modulus Gs in the ring stage, there is considerable scatter in existing data [107, 115, 143]. The latest measurement [115] by membrane refractive index and nanoscale fluctuations gives Gs = 15.3 µN/m, while recent computa- tions [41, 142] have used similar values. We have adopted Gs = 15 µN/m for the simulations in this subsection. To probe the effect of the merozoite, we first run the stretching simulation without any inside the RBC. Then we have tested two sizes of the solid disc, with a thickness of 0.5 µm and diameters of 2 and 3 µm, positioned inside the RBC in different locations and orientations (Fig. 6.4). Our main finding is that with its small size, the merozoite in the ring stage has negligible effects on the deformation of the cell when stretched. The final deformation as a function of stretching force is plotted in Fig. 6.5 for the different merozoite sizes and configurations. In all situations, the effect of merozoite on the RBC deformation is less than 5%. This is understandable since, given their small size, the merozoites are free to move and adjust to the deformation of the membrane. Conceivably, under much larger stretching forces, the disc-shaped merozoite may come into near-contact with the membrane and interfere with cell deformation directly. The numerical predictions in Fig. 6.5, using Gs = 15 µN/m, also agree reasonably well with experimental data. Similar agreements have been noted by earlier computations using similar parameter values and the same bi- concave resting shape [41, 142]. However, observations and measurements 1016.3. Results and analysis Force (pN) (D - D 0 )/ D 0 (% ) 0 50 100 150 -60 -40 -20 0 20 40 60 Experimental Our results Figure 6.5: Variation of the axial and transverse diameters of the RBC as a func- tion of the stretching force in ring stage. The solid curves are computation at Gs = 15µN/m without any merozoites. The error bar on the computed result at 150 pN indicates variations due to inclusion of the merozoite in the different config- urations of Fig. 6.4 The experimental data of [142] are also plotted for comparison. suggest that the area to volume ratio S/V decreases in the ring stage by some 15% from that of the healthy RBC [36]. Such a reduction would be a separate mechanism for the stiffening of the iRBC [107], which has not been explored systematically in the present study. Note that this would challenge the apparent agreement with experiment in Fig. 6.5, though not the main conclusion of this subsection on the negligible effect of the merozoite. 1026.3. Results and analysis (a) (b) Figure 6.6: (a) Micrograph of an iRBC in trophozoite stage, with the merozoite outlined in translucent brown and the hemozoin in purple. Reprinted from Ref. [57] with permission from Elsevier. (b) Axisymmetric model for the iRBC used in the computations, with a spherical merozoite located in the thicker part of the cell. The lengths are in µm, and the parasite is centered at Y = −1.65µm. 6.3.3 Stretching of RBC in the trophozoite stage During the trophozoite stage the parasite is at its metabolically most ac- tive state. It secretes and exports proteins such as the knob-associated histidine-rich protein (KAHRP) [50] onto the RBC membrane that changes its structure and permeability. As a result, the membrane’s elastic modulus increases [90, 142], and the shape, volume and surface area of the cell change as well [36, 131]. In fact, the trophozoite stage is a transition for the RBC shape, from the biconcave in the ring stage to the near spherical shape in the schizont stage. Meanwhile, the parasite itself increases in size and becomes more spherical, eventually reaching roughly 4 µm across [58]. There is a great deal of variation in the published data on the shape and size for both the iRBC and the merozoite in the trophozoite stage. For example, one study reported around 30% loss of volume in late trophozoite stage [165], another reported little change [125], and more recent studies documented 10%–25% increase in volume by the late trophozoite stage [36, 131]. Regarding the surface area, there are reports of increase [36], decrease [165] or no change [107, 131]. Roughly speaking, the iRBC is more or less 1036.3. Results and analysis axisymmetric, but with asymmetry between the two ends of the elongated axis. Besides, there are many bumps and cavities on the surface. A typical image is shown in Fig. 6.6(a). Based on these observations, we have adopted an elongated axisymmetric shape for the RBC (Fig. 6.6b), the narrower end reflecting the cavities at the right end of the real RBC. We select the length scales of the model iRBC to match the volume and length measured by two experiments [36, 131], which are the most recent and in agreement. But the surface area is around 15% less than reported by Esposito et al. [36] for lack of the surface bumps and cavities. We represent the merozoite as a rigid sphere initially located in the thicker part of the iRBC (Fig. 6.6b). We have tested a range of merozoite sizes, with radius up to Rm = 2.25 µm and volume up to 47% of total iRBC volume [131]. The initial center of the parasite is always 1.65 µm below the center of the iRBC. As a result of the aforementioned structural and conformational changes, iRBCs in the trophozoite stage become much less deformable than their healthy counterparts. Conventionally, this loss of deformability is ascribed to an elevated membrane stiffness. Several groups have used an effective shear modulus Gs ranging from 21.3 to 40 µN/m to explain the reduced de- formation of the iRBC when stretched [41, 71, 141, 142]. On the other hand, micropipette aspiration of small areas of the membrane shows a much lower local modulus, 9 – 14 µN/m at the late schizont stage, up from around 6 µN/m [34, 107] for healthy RBC membrane. In the following, we explore the effect of the parasite particle on the overall deformability of iRBC with our numerical model, and demonstrate that the current practice of attributing the loss of deformability entirely to increased membrane modulus is flawed. We have carried out a series of simulations by varying the radius of the merozoite Rm from 0 to 2.25 µm, and separately the membrane shear modulus Gs from 10 to 27.5 µN/m. The stretching force is fixed at F = 150 pN in these simulations. The final axial length of the cell, after attaining steady-state stretching, is plotted in Fig. 6.7(a) as a function of Rm and Gs. The relative deformation, as a percentage of the undeformed length of the iRBC, decreases with both Gs and Rm. The dependence on Rm is relatively mild for smaller merozoites, reminiscent of the lack of an effect in the ring stage (Fig. 6.5). But for larger Rm, the deformability plunges precipitously, as the stretching entails narrowing of the transverse dimension of the cell with the membrane pressing directly on the rigid merozoite inside. The most interesting observation of Fig. 6.7(a) is a kind of compensa- tion between Gs and Rm; a certain amount of deformability can be reached by having a larger merozoite paired with a softer membrane, or a smaller 1046.3. Results and analysis (a) (b) R m [ µ m] G S [µ N/ m ] 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 10 15 20 25 30 Figure 6.7: (a) Deformation surface of trophozoite-stage iRBC as a function of membrane elastic modulus Gs and the merozoite radius Rm. The horizontal plane is at 39% deformation, the experimentally observed amount. (b) The intersection of the 39%-deformation plane and the deformation surface, plotted on the (Gs, Rm) plane. 1056.3. Results and analysis (a) (b) Figure 6.8: (a) Scanning electron micrograph of an IRBC in the schizont stage. K indicates a knob. Reprinted from Ref. [150] with permission from Elsevier. (b) An oblate spheroid model for the schizont-infected iRBC, with length in µm. The parasite is represented by an oblate spheroid of the same aspect ratio as the whole cell. merozoite with a stiffer membrane. To illustrate this compensation more quantitatively, we note that experiments have reported a 39% axial length- ening of iRBC when stretched by F = 150 pN in the trophozoite stage [142]. By cutting the “deformation surface” of Fig. 6.7(a) by the horizontal plane of 39% deformation, we obtain the intersection curve by linear interpolation and plot it in Fig. 6.7(b). Clearly, the 39% overall deformation observed in the experiment could have been produced by an infinite number of (Rm, Gs) combinations. In particular, taking the merozoite size Rm = 2.15 µm in late trophozoite stage [131], the 39% deformation requires a membrane shear modulus of Gs = 18.2 µN/m, comparable with the local modulus of 9–14 µN/m measured by aspiration of small membrane tongues [107]. On the other hand, extrapolating the curve toward Rm = 0 will give Gs values on the order of 30 µN/m, consistent with those obtained by interpreting the loss of deformability entirely as due to an increasing Gs [41]. Thus, the discrepancy in Gs in the literature is explained by the existence of the solid merozoite. 6.3.4 Stretching of RBC in the schizont stage The last erythrocytic stage of malaria parasite starts around 40 hours after invasion. During schizont stage the merozoite undergoes 4 or 5 rounds of mitosis that results in 16–20 daughter merozoites. The daughter merozoites and a remnant body of hemozoin crystals at the center of the cell are tightly 1066.3. Results and analysis Figure 6.9: Closely packed daughter cells of P. falciparum in the schizont stage. Reprinted from Ref. [58] with permission from Elsevier. enclosed within a compact digestive vacuole in the RBC in the schizont stage. Meanwhile, the volume of the iRBC shrinks from the trophozoite stage, with the iRBC taking on a more rounded shape, swelling and eventually bursting at the end of the schizont stage [36, 58, 131]. The task of defining a suitable geometric model for the iRBC in schizont stage is as challenging as that for the trophozoite stage. First, the iRBC shape is geometrically complex, with numerous small-scale patterns. An example is shown in Fig. 6.8(a) taken from [150], with a knob indicated by K. Second, the iRBC and merozoite configurations differ among cells in any cohort, and continually change with disease progression. Finally and consequently, there is a greater amount of variation among the observations and measurements. For example, there is considerable disagreement among experiments on how the iRBC surface area changes from the trophozoite to the schizont stage. Esposito et al. [36] reported an area reduction of 18%, while other saw little change in the surface area [131]. All seem to agree, however, on a significant volume reduction on the order of 20%. Similarly, care must be taken in representing the conformation of the merozoites. The mitotic proliferation of new daughter merozoites starts from the beginning of the schizont stage, leading to a closely packed array as depicted in Fig. 6.9. Various estimations put the total volume of the parasites at the end of the schizont stage from 40 µm3 up to 63 µm3 [33, 75, 131]. The upper bound amounts to 83% of total cell volume [131]. In view of the above, it would be difficult for any physical model to cap- ture all the geometric and structural features of the iRBC and merozoites in 1076.3. Results and analysis the schizont stage. Inspired by the image of Fig 6.8(a) and prior computa- tions [41], we have used an oblate spheroid to represent the iRBC in schizont stage. The three axes, 6.9 × 6.9 × 3.2 µm, are chosen to match the volume V = 81 µm3 and surface area S = 98.2 µm2 reported by Esposito et al. [36]. The tightly packed parasites will be modeled as a single rigid body suspended in the cytosol. For convenience we assume an oblate spheroidal shape with the same aspect ratio as the iRBC, initially centered inside the red cell (Fig. 6.8b). Its size will be varied in the simulations to demonstrate the effect of the parasites on the deformation of the stretched iRBC. Figure 6.10 depicts a deformation surface similar to that of Fig. 6.7(a); the percentage elongation of the iRBC is plotted as a function of the mem- brane shear modulus Gs and the volume of the solid parasite v as a per- centage of the total volume of the iRBC. The stretching force is fixed at 100 pN. The deformation decreases with both Gs and v, as expected. In partic- ular, the presence of the parasite plays a very important role in stiffening the iRBC as a whole, as speculated by Nash et al. [107]. A closer inspec- tion shows that the merozoite exerts a stronger effect in suppressing total deformation for smaller Gs. This is probably because a softer membrane deforms easily to hug the contour of the merozoite, thus allowing its effect to be fully manifested. A stiff membrane, in contrast, holds its shape well under stretching and partially shields the effect of the solid inclusion. The experiment reported a 22% elongation of the iRBC in the schizont stage under a stretching force of F = 100 pN [142]. Intersecting the deforma- tion surface by a horizontal plane at 22%, we obtain the curve of Fig. 6.10(b) that plots the Gs and v combinations that would produce the observed amount of deformation. Again we note a compensating effect between the two factors; a larger merozoite may compensate for a softer membrane and vice versa. Three points are particularly interesting. First, for the actual parasite volume ratio v ≈ 80% in the late schizont stage [131], the 22% de- formation requires a membrane shear modulus of Gs ≈ 16.8 µN/m. This is reasonably close to the range (9–14 µN/m) estimated from micropipette as- piration [107]. Second, the estimation of Gs ≈ 16.8 µN/m is also close to the value of 18.2 µN/m obtained from Fig. 6.7(b) in the trophozoite stage. This seems to suggest that the membrane modulus stays relatively unchanged in the later stages of infection. It is unclear whether this is consistent with real- ity. Finally, given the significant contribution of the merozoite to the loss of iRBC deformability, neglecting it and attributing the reduced deformation entirely to a stiffening of the membrane would lead to an overestimation of the membrane modulus. Figure 6.10(b) would suggest Gs = 35 µN/m in the absence of the merozoite, much above the locally measured values. This 1086.3. Results and analysis (a) (b) ν [ %] G S [µ N/ m ] 0 20 40 60 80 100 5 10 15 20 25 30 35 40 Figure 6.10: (a) Deformation surface of schizont-stage iRBC as a function of membrane elastic modulus Gs and the and merozoite-to-iRBC volume ratio v. The horizontal plane is at 22% deformation, the experimentally observed amount. (b) The intersection of the 22%-deformation plane and the deformation surface, plotted on the (Gs, v) plane. is comparable to Gs = 40 µN/m used by Fedosov et al. [41] to predict the proper amount of deformation without including the parasite, and further 1096.4. Summary below those used in earlier studies [141, 142] ranging between 50 and 60 µN/m. 6.4 Summary This study is based on a simple mechanical argument that the presence of a solid parasite inside the red cell should affect its overall deformation. If the solid is of a substantial size, it will conceivably interfere with the deformation of the membrane when the cell is stretched. This was hypothesized by Nash et al. [107] more than 20 years ago, but overlooked in subsequent experi- mental and computational studies on stretching RBC by optical tweezers. Using a relatively simple particle-based model, we have carried out detailed computations that have quantified this effect for the first time. In the ring stage, we show that the parasite volume is too small to have much effect on the cell deformation. Thus, the previous assumption that the increased RBC rigidity is entirely due to elevated membrane modulus is correct in this case. In the trophozoite and schizont stages, on the other hand, the red cell membrane tends to press against the solid parasite particle when sufficiently stretched. This leads to a significant reduction in the cell deformability. In each case, we show that the measured size of the parasite coupled with the locally measured membrane shear modulus will produce the observed degree of iRBC rigidification. Ignoring the parasite and attributing the loss of deformation entirely to the membrane rigidity would overestimate the membrane modulus by more than a factor of 2 in the trophozoite stage and up to a factor of 4 in the schizont stage. This explains a discrepancy between the membrane shear modulus locally measured by micropipette aspiration and globally estimated from deformation under stretching. Although our study has established the mechanical effect of solid para- sites in hampering overall cell deformability, we must point out the assump- tions and limitations of the model that may have affected the fidelity and accuracy of the results. The actual RBC membrane has a complex struc- ture, and only the grossest features are retained in our particle-linear-spring network. The most severe shortcoming is the lack of area conservation. One way to address this issue is to use nonlinear springs designed to be strongly resistant to areal dilatation or compression [4], in the spirit of the contin- uum Skalak model [137]. The linearity of the spring leads to only a moderate amount of strain-hardening at large deformation; a stiffer nonlinear spring will improve area conservation by strong strain-hardening. In addition, we have neglected the viscosity of the membrane partly because this work fo- 1106.4. Summary cuses on the steady-state stretching instead of the dynamic transient. The membrane should be treated as viscoelastic in dynamic situations. Although we have focused on a specific example of iRBC being stretched by optical tweezers, the model and methodology developed may be relevant to the more general and fundamentally important question of cell responses to forcing. An example is the transit of red cells in microcirculation and in microfluidic devices designed for cell sorting or cytometric analysis. In a recent review, Hoffman and Crocker [63] noted the difficulty in understand- ing how external force propagates through different intracellular structures, and called for the development of “a detailed model for the (dynamic) me- chanical response of all the relevant load-bearing structures in the cell and their interconnections”. Our particle-based model for P. falciparum-infected RBC is a modest effort toward that goal. 111Chapter 7 Conclusions and recommendations 7.1 Summary The main outcomes of this thesis can be summarized along three research directions as follows. 7.1.1 Improving SPH numerical algorithm We successfully resolved the long-standing problem of SPH method with open boundaries and solid walls. The traditional SPH scheme for comput- ing the incompressible fluid flows suffers from a numerical pressure boundary layer at solid walls and inflow and outflow boundaries, which compromises the accuracy of solution near such boundaries. We demonstrated that the well-known difficulty in projection schemes is caused by the artificial homo- geneous Neumann boundary condition for the pressure Poisson equation. We resolved this problem by utilizing a “rotational pressure-correction scheme” with a consistent pressure boundary condition that relates the normal pres- sure gradient to the local vorticity. As shown in chapter 4, this scheme computes the pressure and velocity accurately near open boundaries and solid objects, and extends the scope of SPH simulation beyond the usual periodic boundary conditions. The proposed rotational projection scheme was applied to two-dimensional channel flows and flow around a square cylinder. we have demonstrated that the new pressure boundary condition removes the spurious pressure boundary layers and markedly improves the accuracy of the solution. In particular, it predicts a drag coefficient on the cylinder within 7% of the correct value for Reynolds numbers up to 50, as compared with errors on the order of 20% incurred by the artificial pressure boundary condition. Our numerical experimentation has shown the rota- tional projection scheme as robust, accurate and efficient in dealing with open boundaries and flow around solid obstacles. Based on these findings, we recommend it as a superior algorithm to the homogeneous projection 1127.1. Summary scheme in computing truly incompressible flows using SPH. This removes a shortcoming in the formalism that has at times limited SPH simulations to periodic boundary conditions and hampered accurate evaluation of the pressure and drag force on solid obstacles. 7.1.2 Presenting a particles based model for RBC We have proposed, in chapter 5, a two-dimensional particle-based model for the red blood cell as the basis for developing more general discrete- particle microstructural models for various cell types. In our model the cell membrane is replaced by a set of discrete particles connected by nonlinear springs; the spring law enforces conservation of the membrane area to a high accuracy. In addition, a linear bending elasticity is implemented using the deviation of the local curvature from the innate curvature of the biconcave shape of a resting red blood cell. The cytoplasm and the external liquid are modeled as homogeneous Newtonian fluids, and discretized by particles as in standard smoothed-particle-hydrodynamics (SPH) solution of the Navier- Stokes equations. Thus, the discrete particles serve not only as a numerical device for solving the partial differential equations, but also as a means for incorporating microscopic physics into the model. The model parameters are determined from experimental measurements of cell viscosity and elastic moduli for shear, areal dilatation and bending. Our proposed model suc- cessfully predict the characteristic parachute shape of RBC in simple shear and its tank treading motion in pressure-driven flows with a high accuracy. 7.1.3 Studying the malaria iRBC by our numerical model We extended the particle-based RBC model to three dimensions (3D) to explore the parasite-driven changes in malaria-infected RBCs. Malaria par- asite initiates a series of biochemical modifications of the cell membrane. Meantime, the parasite itself undergoes continual changes in morphology and increase in size. Both factors conspire to cause a drastic decrease in the deformability of iRBC. The loss of deformability plays a crucial role in disease progression by blocking micro-capillaries. So, it is very important to understand the mechanism behind the progressive stiffening of iRBC during the infection process. Conventionally the loss of deformability is attributed entirely to an increase in the rigidity of the membrane through extra cross- linking caused by parasite proteins. Numerical simulations of optical tweezer stretching test estimated a roughly 10-fold increase in the elastic modulus of the membrane by the last stage of infection. This contrasts the moderate 1137.2. Limitations 3-fold increase in membrane modulus measured locally by micropipet aspi- ration. We utilized our particle-based model to explore the discrepancy in the estimated elasticity of the membrane, and recognized the solid parasite itself as a major contributor to the loss of deformability of iRBC. In our simulations the malaria parasite is modeled as an aggregate of particles constrained to rigid-body motion. The morphology of the iRBC in different stages of disease is designed based on experimental observations. Our three-dimensional simulation of RBC stretching tests have revealed a compensating effect between the existence of malaria parasites and the el- evated stiffness of the membrane on the overall deformability of infected RBC. Based on our simulations, the locally measured membrane modulus coupled with a parasite of appropriate shape and size can predict the mea- sured amount of overall deformation. This resolves the discrepancy about the membrane modulus in the literature, and points to a more realistic in- terpretation of the cell-stretching tests. 7.2 Limitations A downside of the particle-based method is its high computational costs. This in turn puts a limit on the size of the physical system that can be simulated, and on the complexity of the physical model underlying the sim- ulation. In particular, we note the following factors that limit the scope and depth of our study. • Simplicity of the model : Of necessity, the model neglects many de- tails of the cellular microstructure and its interplay with the overall behavior of the cell as a whole. Replacing the complex structure of the membrane with a simple 2D network of particles and springs is a drastic simplification, which has made it impossible to simulate inter- esting phenomena such as the invasion process by parasites, interaction of invading parasite with the membrane, the knob-like structure on the surface of iRBC in the trophozoite stage, pitting and rupture of the membrane at the end of the schizont stage. • Chemical-mechanical coupling: We hinted at the profound biochemical changes caused by P. falciparum infection. In this and many other cellular dynamic processes, the chemo-mechanical coupling plays key roles, and this has been completed ignored in the current work. A more comprehensive model, therefore, should contain a module dealing with the chemical changes in the cell, which may affect the mechanical 1147.2. Limitations properties via, say, modifying the spectrin polymers and their cross- linking. • Cell surface area conservation: Healthy RBCs are known to conserve its surface area under deformation. Our current model suffers from area variations on the order of 10%. Since physiologically the conser- vation is owing to local membrane properties, we did not adopt the conventional approach of using an energetic penalty on variations of the whole-cell area. Instead, we envision improvement of this aspect by implementing suitable nonlinear springs in the model. • Dynamic responses: We have ignored the viscosity of the membrane and treated it as purely elastic. In the cell-stretching simulations, we further neglected the viscosity of the surrounding liquid. These simplifications do not affect the amount of final deformation, which is the focus of Chapter 6, but would affect the dynamic response of the cell in flow situations considered in Chapter 5. Ideally, the mem- brane should be modeled as viscoelastic, with the parameters tuned to accurate measurements. • Simple geometries and flow conditions: The cell-deformation simula- tions are all in relatively simple geometries and flow setups. These fall far short of the physiological conditions in micro-capillaries, tis- sues and organs such as the spleen. In vitro experiments have also employed increasing complex designs (e.g. microfluidic channels) to mimic in vivo conditions. Our current limitations to small, simple flow systems can be expanded and even lifted by developing parallel computation, a task that is underway in our group. To end the thesis, we may say that the explorations reported herein have illustrated the possibilities and promises of particle-based models for representing complex cell dynamics. An ambitious generalization would be to use discrete particles to represent intracellular microstructures, and then to probe their evolution in response to chemical and mechanical stimuli. This approach promises a unique route to understanding the microstructural remodelling inside the cell as a consequence of biochemical reactions on the one hand, and as the cause for pathological changes of cell property on the other. 115Bibliography [1] M. Antia, T. Herricks, and P. K. Rathod. Microfluidic approaches to malaria pathogenesis. Cell. Microbio., 10:19681974, 2008. [2] P. Artymowicz and S. H. Lubow. Dynamics of binary-disk interaction. 1: Resonances and disk gap sizes. Astrophysical Journal, Part 1 (ISSN 0004-637X), 421:2: 651–667, 1994. [3] P. Bagchi. Mesoscale simulation of blood flow in small vessels. Biophys. J., 92:1858–1877, 2007. [4] D. Barth`es-Biesel, A. Diaz, and E. Dhenin. Effect of constitutive laws for two-dimensional membranes on flow-induced capsule deformation. J. Fluid Mech., 460:211–222, 2002. [5] G. K. Batchelor. An Introduction to Fluid Dynamics. Cambridge University Press, 1999. [6] W. Benz and E. Asphaug. Simulations of brittle solids using smoothed particle hydrodynamics. Comput. Phys. Commun., 87:253–65, 1995. [7] G.H. Bledsoe. Malaria primer for clinicians in the united states. South Med J., 98:12:1197–1204, 2005. [8] S. Boey, D. Boal, and D. Discher. Simulations of the erythrocyte cytoskeleton at large deformation. i. microscopic models. Biophys. J., 75:1573–1583, 1998. [9] K. Boryczko, W. Dzwinel, and D. A. Yuen. Dynamical clustering of red blood cells in capillary vessels. J. Mol. Model., 9:16–33, 2003. [10] D. C. Bottino. Modeling viscoelastic networks and cell deformation in the context of the immersed boundary method. J. Comput. Phys., 147:86–113, 1998. [11] M. Breuer, J. Bernsdorf, T. Zeiser, and F. Durst. Accurate computa- tions of the laminar flow past a square cylinder based on two different 116Bibliography methods: lattice-Boltzmann and finite-volume. International Journal of Heat and Fluid Flow, 21:186–196, 2000. [12] C.Y. Chee, H.P. Lee, and C. Lu. Using 3d fluidstructure interaction model to analyse the biomechanical properties of erythrocyte. Physics Letters A, 372:1357–1362, 2008. [13] J. K. Chen, J. E. Beraun, and C. J. Jih. An improvement for ten- sile instability in smoothed particle hydrodynamics. Computational Mechanics, 23:279–287, 1999. [14] Q. Chen, M. Schlichtherle, and M. Wahlgren. Molecular aspects of severe malaria. Clin. Microbiol. Rev., 13:3:439–450, 2000. [15] A.J. Chorin. Numerical solution of the Navier-Stokes equations. Math. Comput., 22:745–762, 1968. [16] P. Cleary, J. Ha, V. Alguine, and T. Nguyen. Flow modelling in casting processes. Appl. Math. Model., 26:171–190, 2002. [17] A. Colagrossi and M. Landrini. Numerical simulation of interfa- cial flows by smoothed particle hydrodynamics. J. Comput. Phys., 191:448–75, 2003. [18] A. F. Cowman and B. S. Crabb. Invasion of red blood cells by malaria parasites. Cell, 124:4:755–766, 2006. [19] R. G. Cox, I. Y. Z. Zia, and S. G. Mason. Particle motions in sheared. suspensions. XXV. Streamlines around cylinders and spheres. J. Col- loid Interface Sci., 27:7–18, 1968. [20] H.A. Cranston, C.W. Boylan, G.L. Carroll, S.P. Sutera, J.R. Williamson, I.Y. Gluzman, and D.J. Krogstad. Plasmodium falci- parum maturation abolishes physiologic red cell deformability. Sci- ence, 223:400–403, 1984. [21] S. J. Cummins and M. Rudman. An SPH projection method. J. Comput. Phys., 152:584–607, 1999. [22] R.A. Dalrymple, M. Gmez-Gesteira, B.D. Rogers, A. Panizzo, S. Zou, A.J.C. Crespo, G. Cuomo, and M. Narayanaswamy. Smoothed Particle Hydrodynamics for Water Waves, in Advances in numerical simula- tion of nonlinear water waves, volume 1. Ma, Q. ed., World Scientific Publishing, 2009. 117Bibliography [23] M. Dao, J. Li, and S. Suresh. Molecularly based analysis of deforma- tion of spectrin network and human erythrocyte. Materials Science and Engineering C, 26:1232–1244, 2006. [24] M. Dao, C.T. Lim, and S. Suresh. Mechanics of the human red blood cell deformed by optical tweezers. Journal of the Mechanics and Physics of Solids, 51:2259–2280, 2003. [25] M. Diez-Silva, M. Dao, J. Han, C. Lim, and S. Suresh. Shape and biomechanical characteristics of human red blood cells in health and disease. MRS BULLETIN, 35:382–388, 2010. [26] DE Discher, DH Boal, and SK Boey. Simulations of the erythrocyte cytoskeleton at large deformation. II. Micropipette aspiration. Bio- phys. J., 75:1584–1597, 1998. [27] A.M. Dondorp, P.A. Kager, J. Vreeken, and N.J. White. Abnormal blood flow and red blood cell deformability in severe malaria. Para- sitology Today, 16:228–232, 2000. [28] C. Dong and R. Skalak. Leukocyte deformability: finite element mod- eling of large viscoelastic deformation. J. Theor. Biol., 158:173–193, 1992. [29] JL Drury and M Dembo. Hydrodynamics of micropipette aspiration. Biophys. J., 76:110–128, 1999. [30] M. M. Dupin, I. Halliday, C. M. Care, L. Alboul, and L. L. Munn. Modeling the flow of dense suspensions of deformable particles in three dimensions. Phys. Rev. E, 75:066707–066724, 2007. [31] W. Dzwinel, K. Boryczko, and D. A. Yuen. A discrete-particle model of blood dynamics in capillary vessels. J. Coll. Interface Science, 258:163–173, 2003. [32] C. D. Eggleton and A. S. Popel. Large deformation of red blood cell ghosts in a simple shear flow. Phys. Fluids, 10:1834–1845, 1998. [33] D.A. Elliott, M.T. McIntosh, H.D. Hosgood, S. Chen, G. Zhang, P. Baevova, and K.A. Joiner. Four distinct pathways of hemoglobin uptake in the malaria parasite Plasmodium falciparum. Proc. Natl. Acad. Sci. USA, 105:2463–2468, 2008. 118Bibliography [34] H. Engelhardt and E. Sackmann. On the measurement of shear elas- tic moduli and viscosities of erythrocyte plasma membranes by tran- sient deformation in high frequency electric fields. Biophys Journal, 54:3:495–508, 1988. [35] A. Esposito, J. Choimet, J. Skepper, J. Mauritz, V. Lew, C. Kaminski, and T. Tiffert. Quantitative imaging of human red blood cells infected with Plasmodium falciparum. Biophysical Journal, 99:953–960, 2010. [36] A. Esposito, J. Choimet, J. Skepper, J. Mauritz, V. Lew, C. Kaminski, and T. Tiffert. Quantitative imaging of human red blood cells infected with Plasmodium falciparum. Biophysical Journal, 99:953–960, 2010. [37] E. Evans and Y. Fung. Improved measurements of the erythrocyte geometry. Microvascular Research, 4:4:335–347, 1972. [38] E. Evans and R. Skalak. Mechanics and thermodynamics of biomem- branes. CRC Press, Inc., 1980. [39] E. A. Evans and Y. C. Fung. Improved measurements of the erythro- cyte geometry. Microvasc. Res., 4:335–347, 1972. [40] E.A. Evans, R. Waugh, and L. Melnik. Elastic area compressibility modulus of red cell membrane. Biophysical Journal, 16:585–595, 1976. [41] D. Fedosov, B. Caswell, S. Suresh, and G. Karniadakis. Quantifying the biophysical characteristics of Plasmodium-falciparum-parasitized red blood cells in microcirculation. PNAS, 108:35–39, 2011. [42] D. A. Fedosov. Multiscale modeling of blood flow and soft matter. PhD thesis, Brown University, 2010. [43] D. A. Fedosov, B. Caswell, and G. E. Karniadakis. Systematic coarse- graining of spectrin-level red blood cell models. Computer Methods in Applied Mechanics and Engineering, 199:1937–1948, 2010. [44] D. A. Fedosov, B. Caswell, and G. E. Karniadakis. Wall shear stress- based model for adhesive dynamics of red blood cells in malaria. Bio- physical Journal, 100:2084–2093, 2011. [45] T. M. Fischer. Tank-tread frequency of the red cell membrane: De- pendence on the viscosity of the suspending medium. Biophys. J., 93:2553–2561, 2007. 119Bibliography [46] T. M. Fischer, M. Stohr-Lissen, and H. Schmid-Schonbein. The red cell as a fluid droplet: tank tread-like motion of the human erythrocyte membrane in shear flow. Science, 202:894–896, 1978. [47] S. C. Gifford, M. G. Frank, J. Derganc, C. Gabel, R. H. Austin, T. Yoshida, and M. W. Bitensky. Parallel microchannel-based mea- surements of individual erythrocyte areas and volumes. Biophysical Journal, 84:623–633, 2003. [48] R. A. Gingold and J. J. Monaghan. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon. Not. Roy. Astron. Soc., 181:375–389, 1977. [49] R. A. Gingold and J. J. Monaghan. Kernel estimates as a basis for general particle methods in hydrodynamics. J. Comput. Phys., 46:429– 453, 1982. [50] F.K. Glenister, R.L. Coppel, A.F. Cowman, N. Mohandas, and B.M. Cooke. Contribution of parasite proteins to altered mechanical prop- erties of malaria-infected red blood cells. Blood, 99:1060–1063, 2002. [51] D. E. Goldberg and A. F. Cowman. Moving in and renovating: export- ing proteins from plasmodium into host erythrocytes. Nature Reviews Microbiology, 8:9:617–621, 2010. [52] J. P. Graya, J. J. Monaghan, and R. P. Swift. Sph elastic dynamics. Computer Methods in Applied Mechanics and Engineering, 190:6641– 6662, 2001. [53] P. M. Gresho and R. L. Sani. On pressure boundary conditions for the incompressible navier-stokes equations. Int. J. Num. Methods Fluids, 7:1111–1145, 1987. [54] J. Gruenberg, D. R. Allred, and I. W. Sherman. Scanning electron microscope-analysis of the protrusions (knobs) present on the surface Plasmodium falciparum-infected erythrocytes. The Journal of Cell Biology, 97:795–802, 1983. [55] J. L. Guermond, P. Minev, and J. Shen. An overview of projection methods for incompressible flows. Comput. Methods Appl. Mech. En- grg, 195:6011–6045, 2006. 120Bibliography [56] J.C. Hansen, R. Skalak, S. Chien, and A. Hoger. An elastic network model based on the structure of the red blood cell membrane skeleton. Biophysical Journal, 70:146–166, 1996. [57] E. Hanssen, C. Knoechel, N. Klonis, N. Abu-Bakar, S. Deed, M. LeGros, C. Larabell, and L. Tilley. Cryo transmission X-ray imag- ing of the malaria parasite, P. falciparum. Journal of Structural Biol- ogy, 173:161–168, 2011. [58] E. Hanssen, P. J. McMillan, and L. Tilley. Cellular architecture of Plasmodium falciparum-infected erythrocytes. International Journal for Parasitology, 40:1127–1135, 2010. [59] S. Henon, G. Lenormand, A. Richert, and F. Gallet. A new deter- mination of the shear modulus of the human erythrocyte membrane using optical tweezers. Biophys. J., 76:1145–1151, 1999. [60] M Herant, WA Marganski, and M Dembo. The mechanics of neu- trophils: Synthetic modeling of three experiments. Biophys. J., 84:3389–3413, 2003. [61] T. Herricks, M. Antia, and P. K. Rathod. Deformability limits of Plasmodium falciparum-infected red blood cells. Cell. Microbio., 11:13401353, 2009. [62] R. M. Hochmuth. Micropipette aspiration of living cells. Journal of Biomechanics, 33:1522, 2000. [63] B. D. Hoffman1 and J. C. Crocker. Cell mechanics: Dissecting the physical responses of cells to force. Annual Review of Biomedical En- gineering, 11:259–288, 2009. [64] P. J. Hoogerbrugge and J. M. V. A. Koelman. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Euro- phys. Lett., 19:155–160, 1992. [65] W. G. Hoover. Isomorphism linking smooth particles and embedded atoms. Physica A, 260:244–54, 1998. [66] W. G. Hoover, C. G. Hoover, and E. C. Merritt. Smooth-particle applied mechanics: Conservation of angular momentum with tensile stability and velocity averaging. Phys. Rev. E, 69:016702: 1–10, 2004. 121Bibliography [67] S. M. Hosseini and N. Amanifard. Presenting a modified SPH al- gorithm for numerical studies of fluid-structure interaction problems. Int. J. Eng. Trans. B, 20:167–178, 2007. [68] S. M. Hosseini and J. J. Feng. A particle-based model for the transport of erythrocytes in capillaries. Chem. Eng. Sci., 64:4488–4497, 2009. [69] S. M. Hosseini and J. J. Feng. Pressure boundary conditions for computing incompressible flows with SPH. J. Comput. Phys., 230:74737487, 2011. [70] S. M. Hosseini, M. T. Manzari, and S. K. Hannani. A fully explicit three-step SPH algorithm for simulation of non-Newtonian fluid flow. Int. J. Num. Methods Heat Fluid Flow, 17:715–735, 2007. [71] Y. Imai, H. Kondo, T. Ishikawa, C. T. Lim, and T. Yamaguchi. Model- ing of hemodynamics arising from malaria infection. Journal ofBiome- chanics, 43:1386–1393, 2010. [72] G.Y. Jiao, K.S.W. Tan, C.H. Sow, M. Dao, S. Suresh, and C.T. Lim. Computational modeling of the micropipette aspiration of malaria in- fected erythrocytes. IFMBE Proceedings, 23:5: 1788–1791, 2009. [73] C. Kelemen, S. Chien, and G. M. Artmann. Temperature transition of human hemoglobin at body temperature: Effects of calcium. Bio- physical Journal, 80:26222630, 2001. [74] K. Kirk and K. J. Saliba. Targeting nutrient uptake mechanisms in plasmodium. Current Drug Targets, 8:75–88, 2007. [75] M. Krugliak, J. Zhang, and H. Ginsburg. Intraerythrocytic plasmod- ium falciparum utilizes only a fraction of the amino acids derived from the digestion of host cell cytosol for the biosynthesis of its proteins. Molecular and Biochemical Parasitology, 119:249–256, 2002. [76] M. Lastiwka, M. Basa, and N. J. Quinlan. Permeable and nonreflecting boundary conditions in SPH. Int. J. Num. Methods Fluids, 61:709– 724, 2008. [77] D.V. Le, J. White, J. Peraire, K.M. Lim, and B.C. Khoo. An implicit immersed boundary method for three-dimensional fluidmembrane in- teractions. Journal of Computational Physics, 228:8427–8445, 2009. 122Bibliography [78] E.-S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, and P. Stansby. Comparisons of weakly compressible and truly incompress- ible algorithms for the SPH mesh free particle method. J. Comput. Phys., 227:8417–8436, 2008. [79] G. Y. H. Lee and C. T. Lim. Bimechanics approaches to studying human diseases. TRENDS in Biotechnology, 25:111–118, 2007. [80] J. Li, M. Dao, C. Lim, and S. Suresh. Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte. Biophys. J., 88:3707–3719, 2005. [81] J. Li, M. Dao, C. T. Lim, and S. Suresh. Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte. Biophysical Journal, 88:3707–3719, 2005. [82] X. Li and K. Sarkar. Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic mem- brane. J. Comput. Phys., 227:4998–5018, 2008. [83] L. D. Libersky, A. G. Petschek, T. C. Carney, J. R. Hipp, and F. A. Allahdadi. High strain lagrangian hydrodynamics. J. Comput. Phys., 109:67–75, 1993. [84] C. T. Lim, E. H. Zhou, and S. T. Quek. Mechanical models for living cells - a review. Journal of Biomechanics, 39:195–216, 2006. [85] C.T. Lim, M. Dao, S. Suresh, C.H. Sow, and K.T. Chew. Large de- formation of living cells using laser traps. Acta Mater., 52:1837–1845, 2004. [86] C.T. Lim, E.H. Zhou, A. Li, S.R.K. Vedula, and H.X. Fu. Exper- imental techniques for single cell and single molecule biomechanics. Materials Science and Engineering C, 26:1278–1288, 2006. [87] G. R. Liu and M. B. Liu. SmoothPd Pdrticlp Hydrodyndmics; a mesh- free particle method, volume 1. World Scientific Publishing Co. Pte. Ltd., 2003. [88] L. B. Lucy. A numerical approach to the testing of the fission hypoth- esis. Astron. J., 82:1013–1024, 1977. [89] A. G. Maier, M. Rug, M. T. O’Neill, M. Brown, S. Chakravorty, T. Szestak, J. Chesson, Y. Wu, K. Hughes, R. L. Coppel, C. Newbold, 123Bibliography J. G. Beeson, A. Craig, B. S. Crabb, and A. F. Cowman. Exported proteins required for virulence and rigidity of plasmodium falciparum- infected human erythrocytes. Cell, 134:48–61, 2008. [90] A.G. Maier, B.M. Cooke, A.F. Cowman, and L. Tilley. Malaria par- asite proteins that remodel the host erythrocyte. Nature Reviews. Microbiology, 7:341354, 2009. [91] M. Meyer, M. Desbrun, P. Schroder, and A. H. Barr. Discrete differential-geometry operators for triangulated 2-manifolds. techni- cal note, 2002. [92] J. P. Mills, M. Diez-Silva, D. J. Quinn, M. Dao, M. J. Lang, K. S. W. Tan, C. T. Lim, G. Milon, P. H. David, O. Mercereau-Puijalon, S. Bon- nefoy, and S. Suresh. Effect of plasmodial RESA protein on deforma- bility of human red blood cells harboring Plasmodium falciparum. Proceedings of the National Academy of Sciences, 104(22):9213–9217, 2007. [93] J.P. Mills, L. Qie, M. Dao, C.T. Lim, and S. Suresh. Nonlinear elastic and viscoelastic deformation of the human red blood cell with optical tweezers. Mech. Chem. Biosyst., 3:169–180, 2004. [94] N. Mohandas and P. G. Gallagher. Red cell membrane: past, present, and future. BLOOD, 112:3939–3948, 2008. [95] J. J. Monaghan. Smoothed particle hydrodynamics. Ann. Rev. Astron. Astrophys., 30:543–574, 1992. [96] J. J. Monaghan. Simulating free surface flows with SPH. J. Comput. Phys., 110:399–406, 1994. [97] J. J. Monaghan. Heat conduction with discontinuous conductivity. Applied Mathematics Reports and Preprints 95/18,Monash University, Melbourne, Australia, 1995. [98] J. J. Monaghan. Sph and riemann solvers. J. Comput. Phys., 136:298– 307, 1997. [99] J. J. Monaghan. SPH without a tensile instability. J. Comput. Phys., 159:290–311, 2000. [100] J. J. Monaghan. Smoothed particle hydrodynamics. Reports on Progress in Physics, 68:8: 1703–1759, 2005. 124Bibliography [101] J.J. Monaghan and A. Kocharyan. Sph simulation of multi-phase flow. Computer Physics Communications, 87:225–235, 1995. [102] J. P. Morris, P. J. Fox, and Y. Zhu. Modeling low Reynolds number incompressible flows using SPH. J. Comput. Phys., 136:214–226, 1997. [103] C. A. Moxon, G. E. Grau, and A. G. Craig. Malaria: modification of the red blood cell and consequences in the human host. British Journal of Haematology, 154:670–679, 2011. [104] I. Mueller, P.A. Zimmerman, and J.C. Reeder. Plasmodium malariae and plasmodium ovalethe ”bashful” malaria parasites. Trends Para- sitol, 23:6:278–283, 2007. [105] M. Musielak. Red blood cell-deformability measurement:review of techniques. Clinical Hemorheology and Microcirculation, 42:47–64, 2009. [106] N. Nakamura, N. Tanaka, and T. Masuzawa. Improvement of model of red blood cell in analysis of blood flow using SPH method (in japanese). Nihon Kikai Gakkai Kanto Shibu. Seimitsu Kogakkai Ibaragi Koenkai Koen Ronbunshu, pages 45–46, 2006. [107] G. Nash, E. O’Brien, E. Gordon-Smith, and J. Dormandy. Abnormal- ities in the mechanical properties of red blood cells caused by Plas- modium falciparum. Blood, 74:855–861, 1989. [108] Nestor, R., Basa, M., and N. Quinlan. Moving boundary problems in the finite volume particle method. ERCOFTAC SIG SPHERIC IIIrd International workshop, Lausanne, Switzland, pages 109–114, 2008. [109] H. Noguchi and G. Gompper. Shape transitions of fluid vesicles and red blood cells in capillary flows. PNAS, 102:40:14159–14164, 2005. [110] G. Oger, M. Doring, B. Alessandrini, and P. Ferrant. Two-dimensional sph simulations of wedge water entries. Journal of Computational Physics, 213:803–822, 2006. [111] G. Oger, M. Doring, B. Alessandrini, and P. Ferrant. An improved sph method: Towards higher order convergence. Journal of Computational Physics, 225:1472–1492, 2007. [112] T. Omori, T. Ishikawa, D. Barth`es-Biesel, A.-V. Salsac, J. Walter, Y. Imai, and T. Yamaguchi. Comparison between spring network 125Bibliography models and continuum constitutive laws: Application to the large de- formation of a capsule in shear flow. Phys. Rev. E, 41918:1–11, 2011. [113] T. Pan and T. Wang. Dynamical simulation of red blood cell rheology in microvessels. [114] Vijay Pappu and Prosenjit Bagchi. 3D computational modeling and simulation of leukocyte rolling adhesion and deformation. Comput. Biol. Med., 38:738–753, 2008. [115] Y. Park, M. Diez-Silva, G. Popescu, G. Lykotrafitis, W. Choi, M. Feld, and S. Suresh. Refractive index maps and membrane dynamics of human red blood cells parasitized by Plasmodium falciparum. PNAS, 105:13730 – 13735, 2008. [116] F. Pierig, S. Serafini, L. Rossi, and M. Magnani. Cell-based drug delivery. Advanced Drug Delivery Reviews, 60:2: 286–295, 2008. [117] J. C. Pinder, R. E. Fowler, A. R. Dluzewski, L. H. Bannister, F. M. Lavin, G. H. Mitchell, R. J. M. Wilson, and W. B. Gratzer. Acto- myosin motor in the merozoite of the malaria parasite, plasmodium falciparum: implications for red cell invasion. Journal of Cell Science, 111:1831–1839, 1998. [118] I. V. Pivkin and G. E. Karniadakis. Accurate coarse-grained modeling of red blood cells. Phys. Rev. Lett., 101:118105–118109, 2008. [119] G. G. Poe and A. Acrivos. Closed-streamline flows past rotating single cylinders and spheres: inertia effects. J. Fluid Mech., 72:605–623, 1975. [120] C. Pozrikidis. Numerical simulation of the flow-induced deformation of red blood cells. Ann. Biomed. Eng., 31:1194–1205, 2003. [121] C. Pozrikidis. Axisymmetric motion of a file of red blood cells through capillaries. Phys. Fluids, 17:031503, 2005. [122] C. Pozrikidis. Numerical simulation of cell motion in tube flow. Ann. Biomed. Eng., 33:165–178, 2005. [123] R. W. Randles and L. D. Libersky. Smoothed particle hydrodynamics: some recent improvements and applications. Comput. Methods Appl. Mech. Eng., 139:375–408, 1996. 126Bibliography [124] E. Sackmann. Biological Membranes Architecture and Function, volume 1 of Handbook of Biological Physics, (ed. R.Lipowsky and E.Sackmann). Elsevier, 1995. [125] K.J. Saliba, H.A. Horner, and K. Kirk. Transport and metabolism of the essential vitamin pantothenic acid in human erythrocytes in- fected with the malaria parasite Plasmodium falciparum. Journal of Biological Chemistry, 273:10190–10195, 1998. [126] S. Sanyal, S. Ege, G. Bouyer, S. Perrot, I. Safeukui, E. Bischoff, P. Buf- fet, K. W. Deitsch, O. Mercereau-Puijalon, P. H. David, T. J. Temple- ton, and C. Lavazec. Plasmodium falciparum stevor proteins impact erythrocyte mechanical properties. Blood, 8, 2011. [127] T. W. Secomb. Mechanics of red blood cells and blood flow in narrow tubes. In C. Pozrikidis, editor, Modeling and Simulation of Capsules and Biological Cells, chapter 4, pages 163–191. Chapman and Hall, London, 2003. [128] T. W. Secomb, R. Hsu, and A. R. Pries. Tribology of capillary blood flow. J. Eng. Tribol., 220:767–774, 2006. [129] T. W. Secomb, B. Styp-Rekowska, and A. R. Pries. Two-dimensional simulation of red blood cell deformation and lateral migration in mi- crovessels. Ann. Biomed. Eng., 35:755–765, 2007. [130] U. Seifert. Configurations of fluid membranes and vesicles. Advances in Physics, 46:1:13–137, 1997. [131] Y. M. Serebrennikova, J. Patel, W. K. Milhous, and L. H. Garc´ıa- Rubio. Quantitative analysis of morphological alterations in Plas- modium falciparum infected red blood cells through theoretical inter- pretation of spectral measurements. Journal of Theoretical Biology, 265:493–500, 2010. [132] J. P. Shelby, J. White, K. Ganesan, P.K. Rathod, and D.T. Chiu. A microfluidic model for single-cell capillary obstruction by Plas- modium falciparum-infected erythrocytes. Proc Nat’l Acad Sci USA, 100:14618–14622, 2003. [133] J. Patrick Shelby, John White, Karthikeyan Ganesan, Pradipsinh K. Rathod, and Daniel T. Chiu. A microfluidic model for single-cell capillary obstruction by plasmodium falciparum-infected erythrocytes. PNAS, 100:14618–14622, 2003. 127Bibliography [134] M. Sinnott, P. W. Cleary, and M. Prakash. An investigation of pul- satile blood flow in a bifurcation artery using a grid-free method. Proc. Fifth International Conference on CFD in the Process Industries, Mel- bourne, 2006. [135] R. Skalak and P. I. Branemark. Deformation of red blood cells in capillaries. Science, 164:717–719, 1969. [136] R. Skalak, N. ¨Ozkaya, and T. C. Skalak. Biofluid mechanics. Ann. Rev. Fluid Mech., 21:167–204, 1989. [137] R. Skalak, A. Tozeren, R. P. Zarda, and S. Chien. Strain energy function of red blood cell membrane. Biophys. J., 13:245–264, 1973. [138] J. Sleep, D. Wilson, R. Simmons, and W. Gratzer. Elasticity of the red cell membrane and its relation to hemolytic disorders: an optical tweezers study. Biophys. J., 77:3085–3095, 1999. [139] R.W. Snow, C.A. Guerra, A.M. Noor, H.Y. Myint, and S.I. Hay. The global distribution of clinical episodes of plasmodium falciparum malaria. BLOOD, 434:214–217, 2005. [140] D. Stamenovic and M. Coughlin. A quantitative model of cellular elasticity based on tensegrity. J. Biomech. Eng., 122:39–43, 2000. [141] S. Suresh. Mechanical response of human red blood cells in health and disease: Some structure-property-function relationships. Journal of materials research, 21:1871–1877, 2006. [142] S. Suresh, J. Spatz, J.P. Mills, A. Micoulet, M. Dao, C.T. Lim, M. Beil, and T. Seufferlein. Connections between single-cell biomechanics and human disease states: gastrointestinal cancer and malaria. Acta Bio- materialia, 1:15–30, 2005. [143] R. Suwanarusk, B. M. Cooke, A. M. Dondorp, K. Silamut, J. Sat- tabongkot, N. J. White, and R. Udomsangpetch. The deformability of red blood cells parasitized by plasmodium falciparum and p. vivax. The Journal of Infectious Diseases, 189:190–194, 2004. [144] H. Takeda, S. M. Miyama, and M. Sekiya. Numerical simulation of viscous flow by smoothed particle hydrodynamics. Prog. Theor. Phys., 92:939–960, 1994. 128Bibliography [145] A. Talman, O. Domarle, F. McKenzie, F. Ariey, and V. Robert. Game- tocytogenesis: the puberty of plasmodium falciparum. Malaria Jour- nal, 1:3–24, 2004. [146] N. Tanaka and T. Takano. Microscopic-scale simulation of blood flow using sph method. Int. J. Comput. Methods, 2:558–568, 2005. [147] J. M. Tang and C. Vuik. On deflation and singular symmetric positive semi-definite matrices. J. Comput. Applied Math., 206:603–614, 2007. [148] A. M. Tartakovsky, A. L. Ward, and P. Meakin. Pore-scale simulations of drainage of heterogeneous and anisotropic porous media. Phys. Fluids, 19:103301, 2007. [149] R. Temam. Sur l’approximation de la solution des ´equations de Navier- Stokes par la me´thode des pas fractionnaires II. Arch. Ration. Mech. Anal, 33:377–385, 1969. [150] L. Tilley, M. W. A. Dixon, and K. Kirk. The Plasmodium falciparum- infected red blood cell. The International Journal of Biochemistry and Cell Biology, 43:6: 839–842, 2011. [151] L.J.P. Timmermans, P.D. Minev, and F.N. Van De Vosse. An ap- proximate projection scheme for incompressible flow using spectral elements. Int. J. Numer. Methods Fluids, 22:673–688, 1996. [152] G. Tomaiuolo, M. Barra, V. Preziosi, A. Cassinese, B. Rotolid, and S. Guido. Microfluidics analysis of red blood cell membrane viscoelas- ticity. Lab on a chip, 11:449–454, 2011. [153] R. Tran-Son-Tay, S. P. Sutera, and P. R. Rao. Determination of red blood cell membrane viscosity from rheoscopic observations of tank- treading motion. Biophys. J., 46:65–72, 1984. [154] W. T. Tse and S. E. Lux. Red blood cell membrane disorders. British Journal of Haematology, 104:1:2–13, 1999. [155] K. Tsubota, S. Wada, and T. Yamaguchi. Particle method for com- puter simulation of red blood cell motion in blood flow. Comput. Methods Progr. Biomed., 83:139–146, 2006. [156] K. Tsubota, S. Wada, and T Yamaguchi. Particle method for computer simulation of red blood cell motion in blood flow. Computer Method and Program in Biomedicine, 83:139–146, 2006. 129Bibliography [157] K. Tsubota, S. Wada, and T. Yamaguchi. Simulation study on ef- fects of hematocrit on blood flow properties using particle method. J. Biomech. Sci. Eng., 1:159–170, 2006. [158] J. van Kan. A second-order accurate pressure-correction scheme for viscous incompressible flow. SIAM J. Sci. Comput., 7:870–891, 1986. [159] P. M. Vlahovska, T. Podgorski, and C. Misbah. Vesicles and red blood cells in flow: From individual dynamics to rheology. Computer Methods in Applied Mechanics and Engineering, 10:775–789, 2009. [160] R. Xu, P. Stansby, and D. Laurence. Accuracy and stability in in- compressible sph (isph) based on the projection method and a new approach. Journal of Computational Physics, 228:6703–6725, 2009. [161] Belinda Yap and Roger D. Kamm. Cytoskeletal remodeling and cellu- lar activation during deformation of neutrophils into narrow channels. J. Appl. Physiol., 99:2323–2330, 2005. [162] Belinda Yap and Roger D. Kamm. Mechanical deformation of neutrophils into narrow channels induces pseudopod projection and changes in biomechanical properties. J. Appl. Physiol., 98:1930–1939, 2005. [163] A. Yeung and E. Evans. Cortical shell-liquid core model for passive flow of liquid-like spherical cells into micropipets. Biophys. J., 56:139– 149, 1989. [164] P. Yue, C. Zhou, J. J. Feng, C. F. Ollivier-Gooch, and H. H. Hu. Phase- field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing. J. Comput. Phys., 219:47–67, 2006. [165] M. A. Zanner, W. R. Galey, J. V. Scaletti, J. Brahm, and D. L. Vander Jagt. Water and urea transport in human erythrocytes infected with the malaria parasite plasmodium falciparum. Molecular and Biochem- ical Parasitology, 40:2:269–278, 1990. [166] P. R. Zarda, S. Chien, and R. Skalak. Interaction of viscous incom- pressible fluid with an elastic body. In T. L. Geers, editor, Computa- tional Methods for Fluid-Solid Interaction Problems, New York, 1977. American Society of Mechanical Engineers. 130Bibliography [167] J. Zhang, P. Johnson, and A. Popel. An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows. Phys. Biol., 4:285–295, 2007. [168] D. Zhou and R. H. Wagoner. Development and application of sheet- forming simulation. J. Mater. Process. Tech., 50:1–16, 1995. [169] D. Zhou, P. Yue, and J. J. Feng. Viscoelastic effects on drop deforma- tion in a converging pipe flow. J. Rheol., 52:469–487, 2008. [170] E.H. Zhou, C.T. Lim, K.S.W. Tan, and S.T. Quek. Finite element modeling of the micropipette aspiration of malaria-infected red blood cells. Proc. SPIE 5852, 763, 5852:763–767, 2005. 131
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A particle-based model for computing fluid flows and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A particle-based model for computing fluid flows and cell dynamics Hosseini Amin, Seyed Majid 2012
pdf
Page Metadata
Item Metadata
Title | A particle-based model for computing fluid flows and cell dynamics |
Creator |
Hosseini Amin, Seyed Majid |
Publisher | University of British Columbia |
Date | 2012 |
Date Issued | 2012-04-23 |
Description | The connection between certain human diseases and changes in the mechanical properties of living cells is well established, e.g. in the cases of malaria and cancer. However, the mechanism for the mechanical modifications, which tend to facilitate the pathogenesis of such diseases, is not always clear. For instance, the overall loss of deformability of malaria-infected red blood cells (RBCs) corresponds to a 10-fold increase in the rigidity of the cell membrane. On the other hand, micropipette aspiration has only measured a 3-fold increase in the elastic modulus. In this thesis, a particle-based model is developed to explore the interplay between the underlying microstructures and the behavior of the cell as a whole. The research consists of three related projects. The first project deals with the long-standing problem of Smoothed Particle Hydrodynamics (SPH) method with open boundaries and solid walls. We propose a "rotational pressure-correction scheme" with a consistent pressure boundary condition that leads to a large improvement in accuracy of calculated pressure and the drag coefficient on solid obstacles. The second and third projects concern developing a 2D and then 3D particle-based model for RBCs to explore the parasite-driven changes in malaria-infected RBCs. In our models the cell membrane is replaced by a set of discrete particles connected by linear or nonlinear springs. In addition, a linear bending elasticity is implemented using the deviation of the local curvature from the innate curvature of the biconcave shape of a resting RBC. The cytoplasm and the external liquid are modelled as homogeneous Newtonian fluids, and discretized by particles as in standard SPH solution of the Navier-Stokes equations. The malaria parasite is modelled as an aggregate of particles constrained to rigid-body motion. We argue that the discrepancy in the estimated elastic modulus of the membrane is caused by the presence of the rigid parasite particles inside infected cells, and have carried out numerical simulations to demonstrate this mechanism. Our three-dimensional simulation of RBC stretching tests by optical tweezers accurately demonstrates the compensating effects between the existence of malaria parasites and the elevated stiffness of the membrane on the overall deformability of infected RBC. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Collection |
Electronic Theses and Dissertations (ETDs) 2008+ |
Date Available | 2012-04-23 |
Provider | Vancouver : University of British Columbia Library |
DOI | 10.14288/1.0059285 |
URI | http://hdl.handle.net/2429/42186 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of Chemical and Biological Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2012-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_2012_spring_hosseini-amin_seyed-majid.pdf [ 5.67MB ]
- Metadata
- JSON: 1.0059285.json
- JSON-LD: 1.0059285+ld.json
- RDF/XML (Pretty): 1.0059285.xml
- RDF/JSON: 1.0059285+rdf.json
- Turtle: 1.0059285+rdf-turtle.txt
- N-Triples: 1.0059285+rdf-ntriples.txt
- Original Record: 1.0059285 +original-record.json
- Full Text
- 1.0059285.txt
- Citation
- 1.0059285.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
Canada | 18 | 0 |
United States | 14 | 8 |
Germany | 13 | 58 |
Russia | 7 | 0 |
China | 4 | 0 |
France | 2 | 0 |
Hong Kong | 2 | 0 |
Ukraine | 2 | 0 |
Mexico | 1 | 0 |
Romania | 1 | 0 |
Malaysia | 1 | 0 |
Iran | 1 | 1 |
United Kingdom | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 24 | 58 |
Burnaby | 7 | 0 |
Dartmouth | 6 | 0 |
Ashburn | 4 | 0 |
Beijing | 3 | 0 |
Berlin | 3 | 2 |
Saint-Quentin | 2 | 0 |
Washington | 2 | 0 |
Central District | 2 | 0 |
Winnetka | 2 | 0 |
Kuala Lumpur | 1 | 0 |
San Jose | 1 | 0 |
Bellefonte | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0059285/manifest