Numerical Simulations of Red-Blood Cells in Fluid Flow: A Discrete Multiphysics Study

In this paper, we present a methodological study of modelling red blood cells (RBCs) in shear-induced flows based on the discrete multiphysics (DMP) approach. The DMP is an alternative approach from traditional multiphysics based on meshless particle-based methods. The proposed technique has been successful in modelling multiphysics and multi-phase problems with large interfacial deformations such as those in biological systems. In this study, we present the proposed method and introduce an accurate geometrical representation of the RBC. The results were validated against available data in the literature. We further illustrate that the proposed method is capable of modelling the rupture of the RBC membrane with minimum computational difficulty.


Introduction
The human red blood cell (RBC; erythrocyte) is a homogeneous biconcave disk-like microparticle that is surrounded by a fluidic incompressible lipid bilayer membrane that is underpinned by a thin elastic cytoskeleton [1]. The membrane and its underlying cytoskeleton have been naturally selected to be sufficiently flexible to allow passage of the RBC through capillaries of 4 µm diameter in all tissues; this is approximately half of the cell's main diameter. The main physiological function of the RBC is the transfer of oxygen from the lungs to all tissues and the return of CO 2 to the lungs where it is exhaled. The RBC actively metabolises glucose that provides free energy from its bond cleavage to phosphorylate ADP to make ATP, the 'energy currency' of the cell. In the RBC, an ionic disequilibrium exists across the membrane with a much higher Na + concentration outside in the blood plasma, while inside the concentration of K + is relatively high. The steadystate of these concentrations is maintained by the membrane protein Na, K-ATPase that catalyses the hydrolysis of one molecule of ATP for three Na + ion ejected and two K + taken up. The detailed kinetics of RBC metabolism and ATP turnover are encapsulated in a comprehensive computer model [2][3][4].
Other functions of the RBC are mooted to involve the transfer of ATP and ADP agents inward and outward in the cell to release various molecules in blood flow through cell deformation in shear-induced strain fields [5]. These purines are proposed by some to be key effectors of arteriolar smooth muscle tone and hence affect peripheral resistance to blood flow and thence blood pressure. The dynamics and deformation of the RBC in shear flows are crucial in numerous biomedical applications including the detection of diseases associated with tissue perfusion, the experimental separation of healthy cells from infected parasitised ones, and studying human metabolism in vivo, and using nuclear magnetic resonance spectroscopy (MRS) [6]. The cell might also undergo membrane rupture (haemolysis) under conditions of extreme deformation. To have a deep understanding of these cellular responses, the dynamics of the RBC require investigation.
Furthermore, mechanobiology is concerned with the relationship between cell shape and energy consumption. Experimental methods to study this interrelationship are few, but NMR spectroscopy of RBCs distorted in stretched or compressed hydrogels is one recent example in this field. Specifically, stretching human RBCs by a factor of two enhances their rate of glucose consumption by approximately the same factor [7]. Therefore, being able to model the extent of altered membrane geometry will be key to understanding the role of mechanical or flow induced shape changes of RBCs (and other cells) in whole body energy metabolism.
There are some notable examples of numerical modelling for investigating the complex fluid dynamics of the RBC in shear-induced flows. RBCs in this context have been modelled by many methods including the lattice Boltzmann method (LBM ) [8][9][10] and the finite element (FE) method [11][12][13] in which the fluid-solid interface is captured by the immersed boundary method (IBM). Particle-based methods such as the smoothed particle hydrodynamics (SPH) [14][15][16], discrete particle dynamics (DPD) [17], and moving-particle semi-implicit (MPS) [18] methods have also been used for modelling RBC shape in shear induced fluid flow. The Lagrangian feature of particle-based methods is one of its advantages by which fluid-solid interfaces such as the RBC membrane can be explicitly captured. In all these methods, it is challenging to represent the RBC shape accurately. It is this that undergoes large deformations while leading to no structural damage. The RBC displays in-plane area incompressibility as observed in experimental studies [19,20].
In general, there are two main approaches to representing the RBC membrane, namely shell-based, and spring-based models. In the shell-based model [21], the membrane is represented by a zero-thickness quasi-two-dimensional structure. This view provides a continuum representation of the membrane that can be coupled with several constitutive laws including the neo-Hookean [22][23][24] and Skalak [13,25,26] laws for evaluating in-plane tension and transverse shear. A detailed comparison and investigation of the effect of the requisite constitutive laws has been conducted by Lac et al. [27] for spherical capsules under shear-induced deformations. On the other hand, spring-based models e.g., the latticespring model (LSM) and mass-spring model (MSM) represent the membrane through a (structured) geometrical network of points and particles that are inter-connected by means of virtual springs. In the latter approach, the mechanical response of the membrane depends on the geometry of the structured network (i.e., triangular, tetrahedral, etc.) and it is specified through mathematical correlations that relate the spring constant of the connecting bonds with the overall elastic properties of the membrane. For specifying the in-plane area incompressibility of the membrane in mathematical correlations, we can use linear models that employ Hookean spring characteristics [28], or nonlinear correlations that include virtual dampers [29], or angular springs [30,31]. In addition to these two main approaches (LSM and MSM), other models have also been used for capturing the deformation of the RBC membrane; these include an area-difference-elasticity-theory [32,33]; linear finite element (LFE) method [34], and high fidelity simulations of RBCs [35,36].
As mentioned above, the modelling of RBC behaviour under conditions of flow might also include the sporadic occurrence of haemolysis (membrane rupture) under specific conditions. Taking into account that haemolysis might occur in extreme conditions, a comprehensive model should be able to accommodate all these possible scenarios. On the other hand, and to the best of our knowledge, all models thus far in the literature are not fully comprehensive in this regard. Therefore, this motivated our use of hybrid multiphysics numerical modelling, which is fully capable of capturing all possible scenarios of the RBCs under a shear field. Therefore, we present a methodological study of modelling RBCs in shear-induced flows based on using what we call the discrete multiphysics (DMP) approach [37][38][39]. DMP is an alternative approach from traditional multiphysics by combining particle-based methods, e.g., the SPH for solving fluid flow and the MSM for capturing the elastic deformation of flexible solid materials, and it has been particularly successful in modelling biological systems [40][41][42]. The SPH method is a Lagrangian particle-based technique, which was developed for astrophysical studies [43,44]. The method was later extended for modelling fluid flow [45,46]. In SPH modelling the fluid domain is discretised into Lagrangian computational particles that possess fluid properties viz. density, viscosity, and internal energy, while interacting with their neighbouring particles in a way that is described mathematically by a kernel function [47]. The Lagrangian nature of the SPH method makes it suitable for simulating systems that undergo large deformations, such as occur in free surface flows [46,48], droplet and bubble dynamics [49][50][51], and fluid-solid interactions [52,53]. Within the DMP framework, the MSM was used to generate the RBC membrane through a linear mathematical correlation. In order to accommodate haemolysis in our model, we extended the DMP by setting a limited strain threshold beyond which the virtual bonds, that connect the membrane particles are dispelled, thus representing membrane rupture.

General Equations
An incompressible Newtonian fluid in laminar flow is governed by: where u is the velocity vector, while ρ, µ, and p are the fluid properties density, kinematic viscosity, and pressure, respectively. D Dt is the material time derivative where t denotes time, and F denotes denotes the volumetric body force.
We note once again that the present numerical approach was based on the DMP model [38], which combines the SPH method for modelling fluid flow, and the MSM for representing the elastic properties of membrane particles. It is relevant to outline these approaches along with how their coupling has been achieved mathematically.

The SPH Method
The governing equations of motion i.e., Equations (1) and (2) are solved by the SPH method over computational particles using a kernel function, W(r ij , h). The kernel function, which can be concisely represented by W ij , relates particle i with its neighbouring particle j [46,54], based on their distance r ij = |r i − r j |, and the smoothing length, h as the horizon limit for the interactions with neighbouring particles. Figure 1 represents the interactions between two-dimensional SPH particles based on their relative proximity, through a smoothing kernel function. Among several kernel functions in the literature, the Lucy kernel function [44] was used in this study: where q = (r ij /κh) and β = (105/16 πh 3 ) for three-dimensional simulations. In the SPH method, any field variable f can be approximated using the summation over discrete particles as: where m j and J i denote the mass and the number of computational neighboring particles for particle i. The SPH method discretizes the continuity and momentum equations, Equations (1) and (2) into the following form as: In Equation (6), the first term is the pressure gradient and the second is the dissipation term that conforms to a laminar viscosity model [55]. To evaluate the pressure term in the SPH method, one may solve a Poisson equation to obtain an incompressible solution [50,56,57] or might consider a weakly compressible SPH (WCSPH) approach [48,58,59], in which the pressure is obtained by an equation of state (EOS) using density variations. We used this approach by using the so-called Tait EOS [47,48]: where, c 0 is the reference speed of sound. In order to keep density variations below 1%, it is recommended that the speed of sound should be set at least one order of magnitude larger than the maximum velocity in the domain [45]. ρ 0 is the reference density (1000 kg m −3 ), and γ is empirically set to be 7 in this power law equation [45,46].

The MSM
As mentioned, the MSM [60], which is sometimes referred to as coarse grained molecular dynamics (CGMD) [38], was employed to implement elasticity of the RBC membrane. In the MSM [61], a network of harmonic bonds connects computational particles, allowing them to deform and stretch according to the Newtonian equations of motion using linear spring bonds. Figure 2 shows the structure of the MSM for a quasi-two-dimensional RBC membrane. To account for Hookean elasticity, harmonic bonds are used as follows: where k b and r 0 are the coefficient and the equilibrium distance of Hookean springs, respectively. In this methodological approach, our emphasis is on the discretization and implementation of the model, so we adopted the linear Hookean model for the sake of simplicity and lower computational costs, while other nonlinear models can be easily replaced without further difficulties. Another advantage of the MSM is inherited in the modelling of membrane rupture, for which we nullified the springs once their length exceed the ultimate threshold set equal to 1.05 of their initial value.

Coupling of SPH and MSM
To model the physical properties at the fluid-solid interface, we considered three types of boundary conditions to be taken into consideration [62], i.e., the no-penetration, the no-slip, and continuity of stresses at the interface. In continuum mechanics these conditions are often represented, respectively, as: and where n represents the unit normal vector. u s and u f denote the solid displacement and fluid velocity, respectively, while the solid and fluid stresses are represented by σ s and σ f . Within the DMP framework, we considered ghost SPH particles, which are assigned to MSM particles at the fluid-solid interface. To satisfy all three boundary conditions, these particles have dual identity; they interact with the SPH fluid particles for solving fluid properties but they interact with other MSM particles as membrane particles.

Numerical Algorithm
The velocity verlet (VV) algorithm was used to integrate over time by using a firstorder Euler approach with a variable time step, specified by the following stability condition, ∆t = ζh 2 /ν. Here, ν represents the dynamic viscosity equal to ν = ρ/µ and ζ = 0.125 provided satisfactory results [38]. Considering that ( * ) and (n) superscripts denote variables at the intermediate and nth time step, the particle velocities were calculated using the VV algorithm as: Then, particle densities were calculated according to: In the above equation,ρ is the rate of density variations according to Equation (5), in which the density should be updated based on the velocity difference between the particles u ij . However, using u ij causes numerical instabilities arising from the lag of the velocity in the VV algorithm. So, we use an extrapolated velocity in the calculation of the rate of density variations, which can be represented as: where u ij = (u i − u j ). Once particle velocities are updated, computational particles are moved to their new positions using: Then,ρ (n+1) i , and f (n+1) i were calculated at the (n + 1) time step using Equations (5) and (6), respectively. Finally, the velocity and density were calculated at the (n + 1) time step, respectively, as: and, We have extensively used the above-mentioned numerical algorithm in our previous studies [28,63,64] where we can refer interested readers to for more detailed discussions.

Problem Set-Up
The literature contains several models of the geometry of the RBC, including the frequently cited one by Evans et al. [25]. In this study, we used a model which we consider represents the dimples of the biconcave disc more accurately. The model consists of a set of three parametric equations that involve Jacobi elliptic functions [65] as: and where, SN , CN , and DN are the Jacobi elliptic functions. The coefficient A was defined by using three geometrical distances obtained from experimental electron microscope measurements.
This was used to define other intermediate parameters. Also, g = 1 − (DN (µ, m) 2 SN (ν 1 , m p ) 2 ) and, The key specific lengths of the RBC are the cell diameter at the plane splitting the biconcave into two cup-shaped halves d = 8 × 10 −6 m, the thickness of the cell at the thickest position (like the width of the tyre on a car wheel) h = 2.12 × 10 −6 m, and the thickness at the centre of the dimple, b = 1 × 10 −6 m. Considering these parameters, the surface area and volume of the RBC are found to be equal to S = 85.8 fL and V = 128 µm 2 respectively [65]. Figure 3 shows the shape of the RBC in a two dimensional projection of the three-dimensional body. In the present analysis the RBC was represented by a hybrid model, as described in Section 2 by which the SPH was used for modelling the fluid flow, and the RBC membrane was simulated by the MSM. In order to apply the surface area incompressibility of the membrane (conservation of area), the spring constants were taken to be sufficiently large so that a change in the net surface area was minimized. The spring constant could not be taken to be very large since large spring constants resulted in numerical instabilities. Here, the spring constant was set to k b = 0.0005 N/m. The shear modulus of the RBC was tuned by adjusting the angular spring constant. It has been seen that the angular constant of k a = 1 × 10 −20 Nm/rad led to convincing simulated behaviours in shear fields. Figure 4 represents the 3D simulation domain with height, width, and depth of H, W, and D, respectively. The top and bottom walls abide the no-slip boundary conditions, while periodic boundary conditions were implemented on the rest of boundaries. To maintain a constant shear rateγ, the top and bottom boundaries were given with velocity U and −U, respectively. The RBC was placed at the centre, equidistant from all domain boundaries. A uniformly-spaced Cartesian grid was used to populate the domain with fluid particles. All in all, there were four different types of particles in this model: (1) particles of the top boundary (white); (2) particles of the bottom boundary (red); (3) the RBC particles (blue); (4) the fluid particles (green), as illustrated in Figure 4. In total, the simulation domain contained more than one million computational particles comprising all four particle types. It should be noted that a particle-resolution study was carried out for four different particle resolutions in a shear flow [28] in our previous study, and it was observed that the results and hydrodynamics did not improve significantly beyond the resolution that is equivalent to the one we used in this study. The test-case was simulated on the UK's national super computing facilities (ARCHER with 2.7 GHz, 12-core Intel Xeon E5-2697). The computational cost was 23 h of wall-time on approximately 1500 CPU cores (MPI) for each test-case.
The domain undergoes shear in the x-direction by moving the top and bottom boundary particles in opposite directions. Subsequently, the shear imposed RBC deformation are described and analysed in the next section. . Computational domain for cell under shear rate; the RBC is illustrated by red particles, the top and bottom boundaries are blue, and the fluid particles are green (different particle size, particularly for green and red particles, is a post-processing trick for better illustration and has no further modelling or physical interpretations).

Results
In this section, we show that the presented model was capable of capturing the dynamics of the RBC with superior features, relative to other models reported in the literature. Thus, we show how the current model could predict membrane rupture that was induced by extreme flow conditions.
In order to characterise the deformation of the RBC, a quantitative measure was introduced. Specifically, this was the deformation index, D 12 = |L − H|/(L + H), where L and H were the length of the RBC in the plane passing through its origin in the direction of the main axis, and its transverse direction, respectively. Figure 5 represents the variation of the deformation index for three capillary numbers with respect to dimensionless simulation time and compares them with the numerical data in the literature [13]. It was shown that the present model produced accurate results with respect to the models in the literature and it captured the frequency of deformation with the sinusoidal variations of flow rate. Any discrepancy in the results was due to the different geometries used to model a RBC; these produced a larger amplitude in the sinusoidal variation of the deformation index [13]. Considering the variations of the deformation index shown in Figure 6, it was interesting to note that the RBC underwent a tumbling motion, and it rotated around its centre of mass when the shear rates were symmetrical. This is the reason why there was a sinusoidal trend in the deformation index. At small capillary numbers, the RBC did not deform very much, and accordingly the deformation index exhibited minimal variation over time. The variation of the deformation index was restricted between 0.56 ≤ D 12 ≤ 0.58 for the smallest capillary number. At moderate capillary numbers i.e., Ca = 0.02, the deformation index indicated sinusoidal behaviour of large amplitude. This was interpreted to be due to the resistance of the RBC against the applied shear. When the applied shear met the RBC when it was oriented such that its main diameter was along the flow direction, the RBC geometry was similar to its undisturbed shape. This was the point in time when the deformation index was at its peak value. On the other hand, when the main diameter of the RBC was transverse to the flow direction, the deformation index decreased considerably. For large capillary numbers, the shear rate was sufficiently high to stretch the RBC along the direction of shear. The variation of the deformation index with such high capillary numbers was also due to the re-orientation of the RBC during its tumbling motion. Figure 6 represents the RBC shape at different simulation times for three capillary numbers corresponding to Figure 6. Our numerical model was further compared with data in the literature [13] and it showed that the shape of the RBC was accurately captured. The shape of the RBC remained almost undisturbed during its tumbling motion. For moderate values of capillary number, the RBC was similar in shape to its undisturbed form whereby its major diameter was along the direction of the applied shear. For large capillary numbers, the RBC was stretched along it direction of shear. In many applications, once RBCs are exposed to high shear rates, the membranes of unhealthy cells might rupture. To simulate this biomedical eventuality is difficult but our method is one that can model the responses under shear of healthy cells and unhealthy or parasite-infected ones. On the other hand, current numerical techniques generally have systemic problems for modelling the membrane rupture of RBCs under shear fields. In summary, we present in Figure 7 membrane rupture of a RBC under flow-induced shear. It was clear that rupture could occur at relatively low shear rates in contrast to when the RBC was exposed to extreme shear such as those in Figure 5.
The modelling is a numerical depiction of RBCs under shear and serves to illustrate its capabilities for applications to particles of other geometries and elastic properties. In other words, in order to simulate realistic membrane rupture, careful experimental analysis and membrane characterisation were required to determine the threshold strain that must be in the numerical model. These characterisations lie in the realm of cellular biophysics of RBCs (and other cell types) that will be developed as we refine the current numerical modelling to help interpret experimental data.

Conclusions
In this study, we employed a novel numerical technique using the DMP model for simulating the deformation of RBCs in shear-induced flows. This methodological study describes a simple but accurate and flexible numerical technique for simulating a myriad of problems involving RBCs in in vivo and in vitro applications. In this study, we discussed the methodology in details from the mathematical, numerical, and technical perspectives, describing the governing equations, numerical algorithm, and implementation techniques.
The model has a Lagrangian particle-based nature and utilised a complex mathematical formulation derived from the Jacobi elliptic functions, which provided an accurate representation of the RBC geometry. Then, the proposed mathematical formulation was linked to the MSM, which modelled the RBC through a complex network of interconnected particles on its interface. The proposed model is employed for modelling RBC deformation under shear at various capillary numbers. It was observed that RBCs represented a sinusoidal behaviour in their deformation index factor which is accurately comparable with available data in the literature.
We have further extended the model to simulate the rupture of the RBC's membrane under extremely large shear rates, which could be used in many biological and pharmaceutical applications. The membrane rupture was implemented via breaking the spring bonds between particles of which their deformation exceeded a certain strain threshold. Another advantage of the proposed method is that it can be easily extended to other complex physics, and applied to other real-world problems.