- Research
- Open Access

# SPH-DEM approach to numerically simulate the deformation of three-dimensional RBCs in non-uniform capillaries

- Hasitha-Nayanajith Polwaththe-Gallage
^{1}, - Suvash C. Saha
^{1}, - Emilie Sauret
^{1}, - Robert Flower
^{2}, - Wijitha Senadeera
^{1}and - YuanTong Gu
^{1}Email author

**15 (Suppl 2)**:161

https://doi.org/10.1186/s12938-016-0256-0

© The Author(s) 2016

**Published:**28 December 2016

## Abstract

### Background

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.

### Methods

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.

### Results

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.

### Conclusions

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.

## Keywords

- Blood flow
- Computational biomechanics
- Critical diameter
- Discrete element method
- Meshfree method
- Multiple red blood cells
- Smoothed particle hydrodynamics

## Background

RBCs develop within the bone marrow [1]. 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 [4]. 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 [5]. 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 [6]. 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. [9] 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 [9]. 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 [10]. This model was not able to capture the three-dimensional behaviour of the RBCs. Tsubota and Wada [11] 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. [12] 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 [13] 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 [12]. Recently, the dissipative particle dynamics (DPD) method was employed by Ye et al. [15] 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

*E*

_{ S }is calculated by:

*NS*is the number of springs, equal to 2856,

*K*

_{ S }is the spring constant for stretching/compression, and

*L*

_{ n }and

*L*

_{ n0 }are the present length and original length of the nth spring respectively. Energy associated with the bending deformation,

*E*

_{ B }is calculated by:

*NB*is the number of neighbouring triangles, equal to 2856,

*K*

_{ B }is the spring constant for bending, and

*θ*

_{ n }is the angle between the nth neighbouring triangles. Since the number of lipids per area of the RBC membrane is constant, the membrane area should be conserved locally and as a whole. The energy generated due to the local area changes (

*E*

_{ a }) and total area changes (

*E*

_{ A }) are calculated by:

*A*

_{ n }and

*A*

_{ n0 }are the current and original areas of the considered triangular element respectively, while

*A*and

*A*

_{ 0 }are the current and the reference value of the whole RBC membrane area. Here,

*A*

_{ 0 }is equal to the average area of a healthy matured RBC.

*K*

_{ a }and

*K*

_{ A }in Eqs. (3) and (4) are the area expansion moduli for local area and whole membrane area respectively. The total enclosed volume by the RBC membrane is conserved and energy generated in the membrane, due to the change in total enclosed volume (

*E*

_{ V }) is calculated by:

*V*and

*V*

_{ 0 }are the current and reference volume of the RBC respectively. The reference volume is equal to the average volume of a healthy matured RBC. Finally,

*K*

_{ V }is the penalty coefficient to maintain the

*V*as

*V*

_{ 0 }.

*E*) of the RBC is then calculated by taking the sum of the above energies;

*E*=

*E*

_{ S }+

*E*

_{ B }+

*E*

_{ A }+

*E*

_{ a }+

*E*

_{ V }. Finally, the forces acting on each particle are determined based on the principle of virtual work as:

**F**

_{ i }is the vectorial force acting on the ith particle and

**r**

_{ i }is the position vector of the ith particle.

*V*

_{ 0 }) is set to 60% of the initial volume of the sphere and the reference area of the RBC (

*A*

_{ 0 }) is assumed to be the area of the sphere with the radius of 3.1 µm. When the total energy of the RBC is minimised the typical discoidal biconcave shape of a healthy matured RBC is obtained (see Fig. 2).

## Smoothed particle hydrodynamics approach

**v**,

*m*,

*p*,

*ρ*and

*µ*are velocity, mass, pressure, density and dynamic viscosity of the SPH particles respectively. In the SPH method, any field function of the

*i*th particle is approximated by the same field function values of neighbouring

*j*th particles. Finally,

**F**

_{ i }and

*W*are the external force acting on the particles and the smoothing function.

## Simulation results and discussion

**v**but in opposite directions. Periodic boundary conditions are applied to the inlet and the outlet of the channel, such that a particle leaving the outlet will re-enter the fluid domain through the inlet. However, the properties of that particle is recalculated, using the properties of the neighbouring particles at inlet. Due to the motion of the top and bottom plates of the flow channel, plasma particles start to move and generate a pressure acing on the RBC. As a result the RBC elongates and shows a deformed shape (see Fig. 4). The energy constants of the RBC membrane and other parameters are set as in Table 1.

Key simulation parameters for the model

Parameter | Definition | Value | Reference |
---|---|---|---|

| Spring constant for stretching/compression | 1 × 10 | [11] |

| Spring constant for bending | 1 × 10 | [12] |

| Area expansion modulus for local area | 3 × 10 | [12] |

| Area expansion modulus whole membrane area | 2 × 10 | [12] |

| Penalty coefficient | 50 N/m | [12] |

| Density of the RBC membrane particles | 1098 kg/m | [25] |

| Density of the cytoplasm particles | 1050 kg/m | [26] |

| Density of the plasma particles | 1025 kg/m | [27] |

| RBC membrane viscosity | 1 × 10 | Set |

| Cytoplasm viscosity | 5 × 10 | [28] |

| Plasma viscosity | 1 × 10 | [28] |

*DI)*of the RBC is defined as the ratio between the lengths of the RBC in

*z*-direction to

*y*-direction and is calculated for different shear stress values. Simulation results reveal that the

*DI*increases with the shear stress as shown in Fig. 5. It agrees with the previous experimental results [19] with less than 10% difference. In this study Eq. (8) is used to calculate the shear stress (

*τ*),

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

*d*

_{ i }and

*d*

_{ 0 }respectively) are set to 10.0 µm, while the minimum diameter of the stenosed area (

*d*

_{ c }) is set to 6.8 µm. The total length of the capillary (

*L*) is 57.2 µm and the horizontal distance from the 5th RBC’s centre of mass to the inlet boundary (

*l*

_{ 1 }) and from the 1st RBC’s centre of mass to the narrowest part of the stenosed (

*l*

_{ 6 }) are 3.4 and 10.6 µm respectively (see Fig. 6). The radius of each round corner,

*R*is set to 3.2 µm (see Fig. 6). The distances between two consecutive RBCs (

*l*

_{ 2 },

*l*

_{ 3 },

*l*

_{ 4 }, and

*l*

_{ 5 }) are set to 4 µm. The inlet pressure is set to 1000 Pa, while the outlet pressure is set to zero.

*t*= 0.01 ms). However, at this stage, the 1st RBC shows more deformation compared with the following RBCs due to the hydrodynamic interaction between RBCs [20, 21]. The 1st RBC moves though the stenosed section at about

*t*= 0.02 ms and it experiences severe deformation during that time (see Fig. 7c). When

*t*= 0.025 ms the 1st RBC exits from the stenosed section and it recovers typical deformed parachute shape. When

*t*= 0.025 ms, the 2nd RBC moves through the stenosed section and similar to the 1st RBC, the 2nd RBC also undergoes a large deformation. Similar to the 1st and 2nd RBCs, all the RBCs experience severe deformation, when they pass through the narrowest section of the capillary and they recover their typical deformed parachute shape after the stenosed section. Therefore, similar behaviour would be expected for the all the RBCs if they flow further after the stenosed section of the capillary.

*DI*of five RBCs with time. As expected, due to the hydrodynamic interaction between RBCs, the 1st RBC shows the maximum

*DI*(when

*t*= 0.018 ms). The

*DI*s of the following RBCs are lesser than that of the 1st RBC when they pass through the stenosed section.

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

*DI*of the RBCs compared with preceding RBCs when it passing through the stenosed section, further simulations are conducted. In this study, a single RBC is employed with the different horizontal distances from the RBC’s mass centre to the inlet boundary (

*l*

_{ 1 }) (see Fig. 10). Consequently the horizontal distances from the RBC’s mass centre to the stenosed section (

*l*

_{ 2 }) is changed accordingly. Three cases are studied with

*l*

_{ 1 }= 3.4, 11.4 and 19.4 µm. A capillary with the stenosed diameter of 6.8 µm is used and all the other simulation conditions are kept same as described earlier.

*DI*of the RBCs with time for the three different cases is analysed. The simulation results show that the RBC, initially positioned closer to the stenosed section (

*l*

_{1}= 19.4 µm) experiences less deformation compared to the other RBCs, initially set farther away from the stenosed section (see Fig. 11). As can be seen in Fig. 11, the RBC relevant to

*l*

_{1}= 3.4 µm (initially set farther from the stenosed section) undergoes a considerable deformation before it enters to the stenosed section. Moreover, this RBC is subjected to a further deformation while it passes through the stenosed section and it shows the maximum

*DI*among the three RBCs. On the other hand, the RBC set closer to the stenosed section (

*l*

_{1}= 19.4 µm) does not have enough time to deform before entering to the stenosed section. The deformation of this RBC mainly occurs when it passes through the stenosed region. Therefore, it can be concluded that the RBCs set closer to the stenosed section experience less deformation compared with the deformation of the RBCs set farther away from the stenosed section.

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.

*K*

_{ b }). However, they show local uneven deformation with wrinkles on the membrane (see Figs. 12a, 13) when they pass through the stenosed region and just after the stenosed section. On the other hand, the five RBCs with higher bending stiffness values do not show any wrinkles on the deformed membrane. However, the 1st RBC of five RBCs with bending stiffness of 10

*K*

_{ b }shows some asymmetric behaviour compared to the other four RBCs when it reaches the downstream of the capillary (see Fig. 12c at

*t*= 0.152 ms). The asymmetric behaviour of that RBC is clearly evidenced when the RBCs have higher bending stiffness values; 20, 30 and 40

*K*

_{ b }(see Fig. 12d–f at

*t*= 0.152 ms). Interestingly, the RBCs with the highest bending stiffness values (30 and 40

*K*

_{ b }) do not show observable deformation during their motion before the stenosed section (see Fig. 12e, f at

*t*= 0.048 ms). However, they deform into bullet-like shapes when they pass through the stenosed region and the deformed RBCs are not symmetrical in shape (see Fig. 12e, f at

*t*= 0.096 ms). Furthermore, the RBCs show rolling motions after the stenosed region of the capillary (see Fig. 12e, f at

*t*= 0.152 ms). Additionally, the RBCs with the highest bending stiffness value (40

*K*

_{ b }) exhibit complicated asymmetric deformed shapes when passing through the narrowest section of the capillary.

*K*

_{ b }) value of the RBC as shown in Fig. 14. It can be seen from Fig. 14 that the average total membrane energy of one RBC increases significantly with the membrane bending stiffness (

*K*

_{ b }) of the RBCs, since

*K*

_{ b }is dominant on the total energy of the RBCs. Furthermore, the RBCs with higher membrane energies are very unstable in terms of energy and they are vulnerable to rupture [22].

*DI*to compare the amount of deformation, when RBCs show highly asymmetrical three-dimensional deformed shape. Furthermore, the average total membrane energy does not reflect the amount of bending or the

*DI*of the RBCs, since the higher membrane bending stiffness values always contribute for higher total membrane energy of the RBCs. Therefore, the average membrane energy (

*E*

_{ b }) is normalized using membrane bending stiffness (

*K*

_{ b }). Recalling Eq. (2), this

*E*

_{ b }/

*K*

_{ b }parameter expresses the amount of deformation or the deformability of the RBC. Therefore this parameter can be used to compare the deformability of RBCs, when they have different membrane bending stiffness values.

*E*

_{ b }/

*K*

_{ b }is employed to compare the amount of deformation of the RBCs. It can be seen in Fig. 15, that the RBC with the lowest membrane stiffness (0.1

*K*

_{ b }) shows a very high

*E*

_{ b }/

*K*

_{ b }value compared to RBCs with the typical membrane bending stiffness (

*K*

_{ b }). The

*E*

_{ b }/

*K*

_{ b }ratio decreases substantially when the membrane bending stiffness of the RBCs increases from

*K*

_{ b }to 10

*K*

_{ b }(see Fig. 15). Moreover,

*E*

_{ b }/

*K*

_{ b }further decreases gradually with the membrane bending stiffness of the RBC when the membrane bending stiffness of the RBCs increases from 10 to 50

*K*

_{ b }(see Fig. 15). However, that value reaches a more or less steady value, when the membrane bending stiffness increases from 50 to 100

*K*

_{ b }. Therefore, 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. It can be concluded that the

*DI*or the amount of deformation of the RBCs is governed by the diameter of the stenosed section. Moreover, in order to flow the blood through the cardiovascular network, all the RBCs deform and squeeze through stenosed sections, depending on the diameter of the stenosed section.

## Deformation behaviour in narrow capillaries

*L*) of 60.0 µm and the length of the narrow section (

*l*) of the capillary is set to 21.2 µm (see Fig. 16). The diameter of the narrow section (

*d*

_{ c }) is set to 6.0 µm, in order to represent a lung capillary. The radius of each round corner,

*R*is set to 3.5 µm (see Fig. 16). The inlet (

*d*

_{ i }) and outlet (

*d*

_{ i }) diameters of the capillary is set to 10.0 µm as can be seen in Fig. 16, while the inlet and outlet pressures are set to 1000 Pa and zero respectively.

*DI*s of the three RBCs are similar to each other (see Figs. 17, 18). However, after the narrow section, the 1st RBC exhibits a more deformed shape with a higher

*DI,*while the 3rd RBC shows a rounder shape with a lower

*DI*due to the hydrodynamic interaction between the RBCs. On the other hand, the 2nd RBC takes an intermediate

*DI*compared with the 1st and the 3rd RBCs. This difference in the

*DI*s of the RBCs occurs due to the hydrodynamic interaction between the RBCs.

## 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 [24]. 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) |
---|---|

500 | 4.472 |

1000 | 4.004 |

1500 | 3.985 |

2000 | 3.992 |

2500 | 3.989 |

## Conclusions

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.

## Declarations

### Declarations

**Authors’ contributions**

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.

### Acknowledgements

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.

### Competing interests

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.

### Funding

Publication charges for this article have been funded by Australian Red Cross Blood Service.

**Open Access**This 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.

## Authors’ Affiliations

## References

- 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