SPH-DEM approach to numerically simulate the deformation of three-dimensional RBCs in non-uniform capillaries
© The Author(s) 2016
Published: 28 December 2016
Blood continuously flows through the blood vessels in the human body. When blood flows through the smallest blood vessels, red blood cells (RBCs) in the blood exhibit various types of motion and deformed shapes. Computational modelling techniques can be used to successfully predict the behaviour of the RBCs in capillaries. In this study, we report the application of a meshfree particle approach to model and predict the motion and deformation of three-dimensional RBCs in capillaries.
An elastic spring network based on the discrete element method (DEM) is employed to model the three-dimensional RBC membrane. The haemoglobin in the RBC and the plasma in the blood are modelled as smoothed particle hydrodynamics (SPH) particles. For validation purposes, the behaviour of a single RBC in a simple shear flow is examined and compared against experimental results. Then simulations are carried out to predict the behaviour of RBCs in a capillary; (i) the motion of five identical RBCs in a uniform capillary, (ii) the motion of five identical RBCs with different bending stiffness (K b ) values in a stenosed capillary, (iii) the motion of three RBCs in a narrow capillary. Finally five identical RBCs are employed to determine the critical diameter of a stenosed capillary.
Validation results showed a good agreement with less than 10% difference. From the above simulations, the following results are obtained; (i) RBCs exhibit different deformation behaviours due to the hydrodynamic interaction between them. (ii) Asymmetrical deformation behaviours of the RBCs are clearly observed when the bending stiffness (K b ) of the RBCs is changed. (iii) The model predicts the ability of the RBCs to squeeze through smaller blood vessels. Finally, from the simulations, the critical diameter of the stenosed section to stop the motion of blood flow is predicted.
A three-dimensional spring network model based on DEM in combination with the SPH method is successfully used to model the motion and deformation of RBCs in capillaries. Simulation results reveal that the condition of blood flow stopping depends on the pressure gradient of the capillary and the severity of stenosis of the capillary. In addition, this model is capable of predicting the critical diameter which prevents motion of RBCs for different blood pressures.
RBCs develop within the bone marrow . When a RBC is initially produced, it contains a nucleus inside the cell. However, RBCs eject their nuclei in the early stages of maturity before entering into the main blood stream [2, 3]. Healthy RBCs exhibit a biconcave discoidal shape with a mean diameter of about 8 µm and a mean thickness of about 2 µm at rest . In the cardiovascular network, blood continuously flows through millions of blood vessels, including smallest blood vessels (capillaries), which are even smaller than the mean diameter of an average healthy RBC. Due to the complex three-dimensional geometric structure of RBCs, they exhibit various types of motion and deformations when they flow in the capillaries . Studying the motion and deformation of RBCs is somewhat difficult due to the micro-dimensions of RBCs and the complexity of blood vessels. In this context, numerical modelling techniques have high potential for explaining and predicting the behaviour of RBCs in the capillaries. Among numerical modelling techniques, recently developed meshfree particle methods (MPM) are most effective for analysing problems with large deformations, such as those associated with RBCs . In particular, smoothed particle hydrodynamics (SPH), one of the popular and well-established meshfree particle approaches, has been used widely by researchers for analysing micro-scale hydrodynamics problems [7, 8]. In addition, spring network models are extensively used to model the elastic membrane of RBCs. However, the combination of these two techniques to model RBC flow has not been properly considered, which is one of the innovations in this study.
Recent developments in computational and numerical techniques have made solving the behaviour of RBCs in three-dimensional domains possible. Imani et al.  developed a three-dimensional numerical model to simulate the malaria infected red blood cell (IRBC). In their model, all blood components are modelled by discrete particles, while the malaria parasites inside the RBCs are represented by cluster of rigid particles. For this simulation, the biconcave shaped healthy RBC and the spherical shaped IRBC are used to qualitatively examine the behaviour of RBCs in a narrow 6 µm square channel. Results revealed that IRBC cannot flow through the narrow channels since the IRBCs are stiffer and less deformable . The membrane of the IRBC is modelled by a two-dimensional spring network and the RBCs are considered as two-dimensional disks in a three-dimensional channel. However, the motion and deformation of the RBCs is highly three-dimensional, as the cells exhibit three-dimensional deformations in microvessels . This model was not able to capture the three-dimensional behaviour of the RBCs. Tsubota and Wada  proposed a three-dimensional spring network model to estimate the elastic membrane force of a RBC membrane during its tank treading motion. In their model, the RBC membrane is discretised into triangular elements. Assuming a simple shear flow, a small external force was introduced on each node to reproduce the tank treading motion. This model was further improved by Nakamura et al.  to simulate the mesoscopic blood flow. However, they assumed that RBCs do not disturb the surrounding flow and a one-way coupling was implemented for the flow-RBC by pre-defining the macroscopic flow field.
Nagayama and Honda  developed a three-dimensional model to simulate the behaviour of the RBCs in blood vessels like capillaries the using moving particle semi-implicit (MPS) method. They developed a momentum equation to define the motion of RBCs, considering the inter-particle force, viscous diffusion and external force without solving the Navier–Stokes equations. They studied the motion and deformation of multiple RBCs in bent capillaries. However, this model was not employed to carry out a comprehensive study on motion and deformation of RBCs in capillaries. Pozrikidis [3, 14] developed a three-dimensional model to explain the flow-induced deformation of RBCs. The numerical instabilities of the model and the assumption of axisymmetric behaviour made the model impractical . Recently, the dissipative particle dynamics (DPD) method was employed by Ye et al.  to develop a three-dimensional RBC model to predict the flow through a tube containing interacting RBCs. However, they employed only two RBCs with different properties in order to simulate and investigate the effect of the infected RBC on the motion and deformation of the other RBC.
In this paper, a set of (up to five) RBCs is employed to predict the motion and deformation of RBCs more accurately. An advanced numerical modelling technique combining smoothed particle hydrodynamics (SPH) and the discrete element method (DEM) is used to model the motion and deformation of the set of three-dimensional RBCs in microvessels. The RBC membrane is modelled by a three-dimensional spring network using fundamentals of DEM and the RBC membrane is discretised into a finite number of particles. Each particle represents a finite mass with associated with density and pressure. The application of DEM allows to model the larger deformation of the moving RBCs. The forces acting on the RBC membrane are determined based on the minimum energy concept [5, 16–18]. First, we investigate the motion and deformation of five identical RBCs through a stenosed capillary. Then, the effect of RBC membrane bending stiffness (K b ) on the motion and deformation of the RBCs in a stenosed capillary is explored. Specifically, this study aims to predict the asymmetric motion and three-dimensional deformation of RBCs, when they have an uncharacteristic membrane bending stiffness due to infection by a disease like malaria. Furthermore, we explore how the motion and deformation behaviour of a set of RBCs in blood vessels with the diameters smaller than those of the blood vessels. Finally, the critical diameter for a stenosed capillary to prevent the motion of RBCs is investigated. With the aid of this model the behaviour of three-dimensional RBCs is predicted, with particular focus on blood flow rate under pathological conditions.
Numerical model and solution methodologies
Three dimensional RBC model
Smoothed particle hydrodynamics approach
Simulation results and discussion
Key simulation parameters for the model
Spring constant for stretching/compression
1 × 10−6 N/m
Spring constant for bending
1 × 10−10 N
Area expansion modulus for local area
3 × 10−3 N/m
Area expansion modulus whole membrane area
2 × 10−3 N/m
Density of the RBC membrane particles
Density of the cytoplasm particles
Density of the plasma particles
RBC membrane viscosity
1 × 10−3 Pa s
5 × 10−3 Pa s
1 × 10−3 Pa s
where µ and h are the dynamic viscosity of the plasma and height of the channel in y-direction respectively.
Deformation behaviour of multiple RBCs through a stenosed capillary
Due to the hydrodynamic interaction between RBCs, a lower DI for the 3rd RBC it is expected compared to the 2nd RBC during motion through the stenosed section of the capillary. However, the maximum DI of the 3rd RBC does increase compared to the value of the 2nd RBC. Similarly, the 3rd, 4th and 5th RBCs show higher DIs compared with their preceding RBC in the stenosed section (see Fig. 8). It is not possible to explain this phenomenon with the aid of the hydrodynamic interaction between RBCs and further studies have to be done to describe this behaviour (this phenomenon will be discussed and explained in next section “Effect of the initial set up of the problem domain on the deformation behaviour of RBCs through a stenosed capillary”. Moreover, the DI of the 1st RBC reduces significantly when it exits the stenosed section (see Fig. 8; after t = 0.018 ms). However, the DI gradually increases again with time (see Fig. 8) and the 1st RBC shows higher DI compared with the other RBCs due to the hydrodynamic interaction between RBCs.
Effect of the initial set up of the problem domain on the deformation behaviour of RBCs through a stenosed capillary
In reality, the blood continuously flows and RBCs exhibit deformation (deformed shapes) at all times. However, for the numerical simulations, an initial condition (t = 0) is assumed. For the simplicity of this study, it is assumed that the RBC begins with its typical biconcave shape and there is no deformation at t = 0. For the above three cases, if the horizontal distance from the RBC’s centre to the stenosed section (l 2 ) increases considerably, the maximum DI of all the RBC would have been same, when it passes through the stenosed section. However, it is computationally very expensive to increase the horizontal distance from the RBC’s mass centre to the stenosed section (l 2 ), since it significantly increases the number of particles in the problem domain. The growth in the number of particles in the problem domain would take a longer time to solve the problem and it would be very inefficient. Therefore, the problem domain is controlled for the conditions as explained in earlier sections.
The variation of the maximum DI of five RBCs, when they pass though the stenosed capillary, can be explained (see Fig. 8) using the above argument. The 1st RBC of five RBCs shows the highest maximum DI when it passes through the stenosed section due to the hydrodynamic interaction between the RBCs. The maximum DI of the 2nd RBC is lower than that value of the 1st RBC as expected again due to the hydrodynamic interaction between RBCs. However, the maximum DI of the 3rd RBC shows a greater value, compared with the 2nd RBC. This phenomenon happens due to the difference in the horizontal distance from the RBC’s mass centre to the stenosed section (l 2 ). Initially (at t = 0) the 3rd RBC is set farther away from the stenosed section compared to the 2nd RBC and the 3rd RBC takes longer time to enter the stenosed section of the capillary. During that time RBC experiences a significant deformation.
The increase in the DI of the 3rd RBC just before entering into the stenosed section of the capillary is thus higher than that value of the 2nd RBC. The DI of both RBCs increases significantly, while they pass though the stenosed region. However, due to that change in DI of the 2nd and 3rd RBCs just before entering the stenosed section of the capillary, the 3rd RBC exhibits higher maximum DI when it passes through the stenosed section compared with that value of the 2nd RBC. Similarly, trailing RBCs exhibit higher maximum DI compared with their preceding RBCs (except the 1st RBC). Therefore, it can be concluded that the initial horizontal distance from the RBC’s centre to the stenosed section (l 2 ) is a crucial parameter in numerical simulations, when the simulations are carried out to capture the behaviour of RBCs in stenosed capillary. This parameter has to be chosen properly, to obtain reliable enough results without affecting the computation cost too much. However, this study was limited to the conditions discussed for Fig. 7.
Deformation behaviour of the RBCs with different bending stiffness values in a stenosed capillary
In this study, the effect of RBC membrane bending stiffness on motion and deformation is studied. Five identical RBCs are used to simulate the motion and deformation behaviour of the RBCs in a stenosed capillary. A capillary with a total length of (L) of 57.2 µm is used for this study. The inlet (d i ) and outlet diameters (d 0 ) of the capillary are set to 10.0 µm. In this study the severity of the stenosed section is further increased in order to clearly compare the effect of the membrane bending stiffness of the RBC. The diameter of the stenosed area (d c ) is set to 5.2 µm. The inlet pressure is set to 500 Pa, while the outlet pressure is set to zero. The membrane bending stiffness of the RBCs changes when they are infected by diseases like malaria and cancers. In order to investigate the behaviour of diseased RBCs, the membrane bending stiffness of all the RBCs is changed from the typical value, K b (1 × 10−10 N) to 0.1 K b , 10 K b , 20 K b , 30 K b and 40 K b . The five identical RBCs and all the other simulation parameters are set as described earlier.
Deformation behaviour in narrow capillaries
Critical diameter of the stenosed section to stop the motion of blood flow
The blood flow rate in a capillary reduces when the capillary has a stenosed section and it causes to reduce the overall blood rate through the cardiovascular network. Depending on the severity of the stenosis, there is a high risk of microvascular blockage, which may lead to completely stoppage of the blood flow in that capillary . In this section, the critical stenosed diameter of a capillary is investigated to halt the blood flow in the capillary with the aid of five identical RBCs. In this study the diameter of the stenosed section is changed to 8.4, 7.6, 6.8, 6.0, 5.2, 4.4, 3.6 and 2.8 µm in the capillary and the radius of each round corner, R is adjusted accordingly. The inlet (d i ) and outlet diameters (d 0 ) of the capillary are set to 10.0 µm, while the total length of the capillary (L) is set to 57.2 µm. The inlet and outlet pressures are set to 1000 Pa and zero respectively.
However, in reality, when the blood flow stops, the blood pressure builds up before the stenosed section. In order to investigate the effect on the pressure, the simulations are carried out for different inlet pressure values. Similar to the previous case, the critical diameter of the stenosed section of the capillary, which stops the blood flow, is found for different inlet pressure values. In this study the inlet pressure of the stenosed capillary is changed to 500, 1000, 1500, 2000 and 2500 Pa while keeping the outlet pressure of the capillary to zero. All the other parameters are kept constant.
The critical diameter of the stenosed section of the capillary for different inlet pressures
Inlet pressure of the stenosed capillary (Pa)
Critical diameter of the stenosed section (µm)
A three-dimensional spring network model based on DEM is used in combination with the SPH method to model the motion and deformation of RBCs in capillaries. From this study, the below conclusions are drawn.
In numerical simulation, initial setting of the RBCs directly affects the deformation behaviour of the RBCs. The lengths of the capillaries should be long enough to obtain reliable enough results without affecting the computation cost too much.
When the membrane bending stiffness of the RBCs increase like in malaria infected RBCs they show highly asymmetrical deformed shapes and rolling motions. On the other hand, the RBCs with lower membrane bending stiffness values exhibit wrinkles on the membrane when they are deforming.
Irrespective of the membrane bending stiffness of the RBCs, RBCs deform a certain amount in order to pass through the stenosed section of the capillary.
The RBCs exhibit bullet-like shapes when they flow through the capillaries with narrower sections, which are narrower than the diameter of the RBCs at rest. However, they show parachute shapes when the diameter of the section they are moving through, increases.
There is a certain critical diameter for a given stenosed capillary and for a given pressure gradient which completely stops the motion of blood with RBCs, which leads to microvascular blockages.
HNPG coded the model, performed the numerical simulation and the analysis of results. SCS and ES contributed to conception of fluid mechanics and CFD, carried out the literature review, and helped the analysis of results. RF contributed the biomedical conception. WS contributed to the SPH model. YTG contributed to concept of models, carried out the literature review and performed the analysis of results. All authors have been involved in drafting the manuscript or revising it critically for important intellectual content and have given final approval of the version to be published. Each author has participated sufficiently in the work to take public responsibility for appropriate portions of the content. All authors read and approved the final manuscript.
Support provided by the ARC grants (FT100100172, DP150100828 and LP150100737), the High Performance Computer (HPC) resources in Queensland University of Technology (QUT) and the Queensland Cyber Infrastructure Foundation (QCIF) are gratefully acknowledged. The fourth author acknowledges that the Australian Governments fully fund the Australian Red Cross Blood Service for the provision of blood products and services to the Australian community. Ms. Sarah Barns is acknowledged for proofreading the manuscript.
The authors declare that they have no competing interests.
About this supplement
This article has been published as part of BioMedical Engineering OnLine Volume 15 Supplement 2, 2016. Computational and experimental methods for biological research: cardiovascular diseases and beyond. The full contents of the supplement are available online http://biomedical-engineering-online.biomedcentral.com/articles/supplements/volume-15-supplement-2.
Availability of data and materials
The data related to this research is available to all interested researchers upon reasonable requests.
Publication charges for this article have been funded by Australian Red Cross Blood Service.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Perrod DBC. Epigenetic PU. 1 silencing in myeloid leukemia by mimicrying a T cell specific chromatin loop. Berlin: Humboldt-Universität zu Berlin; 2013.Google Scholar
- Freund JB, Zhao H. A high-resolution fast boundary-integral method for multiple interacting blood cells. In: Pozrikidis C, editor. Computational hydrodynamics of capsules and biological cells. Boca Raton: CRC Press; 2010. p. 71.View ArticleGoogle Scholar
- Pozrikidis C. Numerical simulation of the flow-induced deformation of red blood cells. Ann Biomed Eng. 2003;31(10):1194–205.View ArticleGoogle Scholar
- Dupire J, Socol M, Viallat A. Full dynamics of a red blood cell in shear flow. Proc Natl Acad Sci. 2012;109(51):20808–13.View ArticleGoogle Scholar
- Polwaththe Gallage HN, Saha SC, Gu Y. Formation of the three-dimensional geometry of the red blood cell membrane. ANZIAM J. 2014;55:C80–95.MathSciNetView ArticleGoogle Scholar
- Liu Y, Liu WK. Rheology of red blood cell aggregation by computer simulation. J Comput Phys. 2006;220(1):139–54.MathSciNetView ArticleMATHGoogle Scholar
- Tanaka N, Takano T. Microscopic-scale simulation of blood flow using SPH method. Int J Comput Methods. 2005;2(04):555–68.View ArticleMATHGoogle Scholar
- Polwaththe-Gallage HN, et al. Numerical investigation of motion and deformation of a single red blood cell in a stenosed capillary. Int J Comput Methods. 2015;12:1540003.MathSciNetView ArticleGoogle Scholar
- Imai Y, et al. Three-dimensional simulation of blood flow in malaria infection. In: 13th international conference on biomedical engineering. Berlin: Springer; 2009.Google Scholar
- Secomb TW, Styp-Rekowska B, Pries AR. Two-dimensional simulation of red blood cell deformation and lateral migration in microvessels. Ann Biomed Eng. 2007;35(5):755–65.View ArticleGoogle Scholar
- Tsubota K-I, Wada S. Elastic force of red blood cell membrane during tank-treading motion: consideration of the membrane’s natural state. Int J Mech Sci. 2010;52(2):356–64.View ArticleGoogle Scholar
- Nakamura M, Bessho S, Wada S. Spring-network-based model of a red blood cell for simulating mesoscopic blood flow. Int J Numer Methods Biomed Eng. 2013;29(1):114–28.MathSciNetView ArticleGoogle Scholar
- Nagayama K, Honda K. 3D particle simulations of deformation of red blood cells in micro-capillary vessel. Rijeka: INTECH Open Access Publisher; 2012.View ArticleGoogle Scholar
- Pozrikidis C. Numerical simulation of cell motion in tube flow. Ann Biomed Eng. 2005;33(2):165–78.View ArticleGoogle Scholar
- Ye T, et al. Dissipative particle dynamics simulations of deformation and aggregation of healthy and diseased red blood cells in a tube flow. Phys Fluids. 2014;26(11):111902.View ArticleGoogle Scholar
- Tsubota K-I, Wada S, Yamaguchi T. Simulation study on effects of hematocrit on blood flow properties using particle method. J Biomech Sci Eng. 2006;1(1):159–70.View ArticleGoogle Scholar
- Pan TW, Wang T. Dynamical simulation of red blood cell rheology in microvessels. Int J Numer Anal Model. 2009;6:455–73.MathSciNetMATHGoogle Scholar
- Polwaththe-Gallage HN, et al. Numerical simulation of red blood cells’ motion: a review. In: 4th international conference on computational methods (ICCM 2012). Crowne Plaza, Gold Coast, QLD. 2012.Google Scholar
- Schauf B, et al. The laser diffractoscope—a new and fast system to analyse red blood cell flexibility with high accuracy. Lasers Med Sci. 2003;18(1):45–50.View ArticleGoogle Scholar
- Shi X, et al. A lattice Boltzmann fictitious domain method for modeling red blood cell deformation and multiple-cell hydrodynamic interactions in flow. Int J Numer Methods Fluids. 2013;72(8):895–911.MathSciNetView ArticleGoogle Scholar
- Pozrikidis C. Axisymmetric motion of a file of red blood cells through capillaries. Phy Fluids. 2005;17(3):031503.MathSciNetView ArticleMATHGoogle Scholar
- Omori T, et al. Tension of red blood cell membrane in simple shear flow. Phys Rev E. 2012;86(5):056321.View ArticleGoogle Scholar
- Boryczko K, Dzwinel W, Yuen DA. Dynamical clustering of red blood cells in capillary vessels. J Mol Model. 2003;9(1):16–33.View ArticleGoogle Scholar
- Cooke BM, Mohandas N, Coppel RL. The malaria-infected red blood cell: structural and functional changes. Adv Parasitol. 2001;50:1–86.View ArticleGoogle Scholar
- Sun C, Munn LL. Particulate nature of blood determines macroscopic rheology: a 2-D lattice Boltzmann analysis. Biophys J. 2005;88(3):1635–45.View ArticleGoogle Scholar
- Le D-V, et al. An implicit immersed boundary method for three-dimensional fluid-membrane interactions. J Comput Phys. 2009;228(22):8427–45.MathSciNetView ArticleMATHGoogle Scholar
- Frcitas RA. Exploratory design in medical nanotechnology: a mechanical artificial red cell. Artif Cells Blood Substit Biotechnol. 1998;26(4):411–30.View ArticleGoogle Scholar
- Fedosov DA, Caswell B, Karniadakis GE. A multiscale red blood cell model with accurate mechanics, rheology, and dynamics. Biophys J. 2010;98(10):2215–25.View ArticleGoogle Scholar