Bridging the computational gap between mesoscopic and continuum modeling of red blood cells for fully resolved blood flow

TL;DR

This paper introduces a modular computational framework for simulating blood flow with fully resolved red blood cells, combining a lattice Boltzmann solver, a novel finite element method, and an immersed boundary approach, offering improved speed, stability, and flexibility.

Contribution

The paper presents a new nodal projective FEM solver for RBCs that is faster, more stable, and requires fewer parameters than existing methods, with enhanced membrane modeling capabilities.

Findings

The solver achieves stability and efficiency comparable to mesoscopic models.

It accurately captures RBC viscoelastic behavior and tank-treading frequencies.

The framework is validated against experimental data and is scalable for multiple cells.

Abstract

We present a computational framework for the simulation of blood flow with fully resolved red blood cells (RBCs) using a modular approach that consists of a lattice Boltzmann solver for the blood plasma, a novel finite element based solver for the deformable bodies and an immersed boundary method for the fluid-solid interaction. For the RBCs, we propose a nodal projective FEM (npFEM) solver which has theoretical advantages over the more commonly used mass-spring systems (mesoscopic modeling), such as an unconditional stability, versatile material expressivity, and one set of parameters to fully describe the behavior of the body at any mesh resolution. At the same time, the method is substantially faster than other FEM solvers proposed in this field, and has an efficiency that is comparable to the one of mesoscopic models. At its core, the solver uses specially defined potential energies, and builds upon them a fast iterative procedure based on quasi-Newton techniques. For a known material, our solver has only one free parameter that demands tuning, related to the body viscoelasticity. In contrast, state-of-the-art solvers for deformable bodies have more free parameters, and the calibration of the models demands special assumptions regarding the mesh topology, which restrict their generality and mesh independence. We propose as well a modification to the potential energy proposed by Skalak et al. 1973 for the red blood cell membrane, which enhances the strain hardening behavior at higher deformations. Our viscoelastic model for the red blood cell, while simple enough and applicable to any kind of solver as a post-convergence step, can capture accurately the characteristic recovery time and tank-treading frequencies. The framework is validated using experimental data, and it proves to be scalable for multiple deformable bodies.

Figures

Tables

Table 1. Table 1 : Overview of used parameters

Parameter	Value
$\Delta x$	$0.5\mu m$
relaxation time $\tau$	$[1,2]$
$\Delta t$	equation (54)
Skalak B,C,D	$5,5000,35pN/\mu m$
$k_{bending}$	$1pN\mu m$
$k_{volume}$	$250$
Rayleigh ${\alpha }_D,{\beta }_D$	$0,0.01$
${\kappa }_{\text{damping}}$	$0.6$
IBM cycles	$[1,5]$
RBC mesh resolution (surface vertices)	$\{66,130,258,514,719,1026\}$

Table 2. Table 2 : Comparison of npFEM with Hemocell [ 21 , 22 ] . The timing refers to the deformable bodies solvers.

SV: Surface Vertices	Hemocell SV:642	npFEM SV:258	npFEM SV:667
Average execution time per time step [ms/step]	0.58	1.35	5.19

Equations

$${\mathbf{v}}_{n+1}=\frac{{\mathbf{x}}_{n+1}-{\mathbf{x}}_n}{h},$$

$${\mathbf{F}}_{int}({\mathbf{x}}_{n+1})+{\mathbf{F}}_{ext}({\mathbf{x}}_n)-\mathbf{C}{\mathbf{v}}_{n+1}=\mathbf{M}\frac{{\mathbf{v}}_{n+1}-{\mathbf{v}}_n}{h},$$

$$\tilde{\mathbf{M}}({\mathbf{x}}_{n+1}-{\mathbf{y}}_n)=h^2{\mathbf{F}}_{int}({\mathbf{x}}_{n+1}),$$

$$\underset{{\mathbf{x}}_{n+1}}{\min }\frac{1}{2h^2}{\Vert {\tilde{\mathbf{M}}}^{\frac{1}{2}}({\mathbf{x}}_{n+1}-{\mathbf{y}}_n)\Vert }_F^2+\sum _i{E}_i({\mathbf{x}}_{n+1}),$$

$$E_i(\mathbf{x})=\underset{{\mathbf{S}}_i\mathbf{p}\in {\mathcal{M}}_i}{\min }{\tilde{E}}_i(\mathbf{x},{\mathbf{S}}_i\mathbf{p}),{\tilde{E}}_i(\mathbf{x},\mathbf{z})={\Vert {\mathbf{G}}_i\mathbf{x}-\mathbf{z}\Vert }_F^2,$$

$$E(\mathbf{x},\mathbf{p})=\underset{\mathbf{z}\in \mathcal{M}}{\min }{d}^{{}^{\prime }}(\mathbf{G}\mathbf{x},\mathbf{z})={\Vert \mathbf{G}\mathbf{x}-\mathbf{p}\Vert }_F^2.$$

$$g({\mathbf{x}}_{n+1})=\frac{1}{2h^2}tr(({\mathbf{x}}_{n+1}-{\mathbf{y}}_n{)}^T\tilde{\mathbf{M}}({\mathbf{x}}_{n+1}-{\mathbf{y}}_n))+\sum _i{w}_i\left(\underset{{\mathbf{S}}_i\mathbf{p}\in {\mathcal{M}}_i}{\min }{\tilde{E}}_i(\mathbf{x},{\mathbf{S}}_i\mathbf{p})\right),$$

$$\begin{array}{rl}g({\mathbf{x}}_{n+1})=& \frac{1}{2h^2}tr\left(({\mathbf{x}}_{n+1}-{\mathbf{y}}_n{)}^T\tilde{\mathbf{M}}({\mathbf{x}}_{n+1}-{\mathbf{y}}_n)\right)+\\ & \underset{\text{projective dynamics potential energies: quadratic formulation}}{\underbrace{\frac{1}{2}tr\left({\mathbf{x}}_{n+1}^T\mathbf{L}{\mathbf{x}}_{n+1}\right)-tr\left({\mathbf{x}}_{n+1}^T\mathbf{J}\mathbf{p}({\mathbf{x}}_{n+1})\right)+\frac{1}{2}tr\left({\mathbf{p}}^T({\mathbf{x}}_{n+1})\mathbf{S}\mathbf{p}({\mathbf{x}}_{n+1})\right)}},\end{array}$$

$$g(\mathbf{x})=\frac{1}{2h^2}tr((\mathbf{x}-\mathbf{y}{)}^T\tilde{\mathbf{M}}(\mathbf{x}-\mathbf{y}))+\sum _i{w}_i{E}_i^{PD}(\mathbf{x})+\sum _i{E}_i^{non-PD}(\mathbf{x}),$$

$$\nabla g(\mathbf{x})=\frac{1}{h^2}\tilde{\mathbf{M}}\left(\mathbf{x}-\mathbf{y}\right)+\underset{\text{−FintPD}}{\underbrace{\mathbf{L}\mathbf{x}-\mathbf{J}\mathbf{p}(\mathbf{x})}}+\underset{\text{−Fintnon−PD}}{\underbrace{\sum _i{\nabla }_{\mathbf{x}}{E}_i^{non-PD}(\mathbf{x})}}.$$

$$E(\mathbf{f},\mathbf{T})=\frac{k}{2}{\int }_S{\Vert {\nabla }_S\mathbf{f}-\mathbf{T}{\nabla }_S\mathbf{g}\Vert }_F^{2}dA,$$

$$E_i(\mathbf{x},\mathbf{T})=\frac{k_i}{2}{A}_i||{\mathbf{D}}_{\mathbf{s}}{{\mathbf{D}}_{\mathbf{m}}}^{-1}-\mathbf{T}|{|}_F^2,$$

$$E_i(\mathbf{x},\mathbf{T})=\frac{k_i}{2}{A}_i||\mathbf{F}-\mathbf{T}|{|}_F^2,$$

$$\mathbf{F}=\mathbf{U\Sigma }{\mathbf{V}}^T,$$

$$det(\mathbf{F})=det(\mathbf{\Sigma })=1.$$

$$\frac{k_i}{2}{A}_i||\mathbf{\Sigma }-{\mathbf{\Sigma }}^{\ast }|{|}_F^2\text{s.t.}\text{det}({\mathbf{\Sigma }}^{\ast })=1,$$

$$E(\mathbf{f},\mathbf{R})=\frac{k}{2}{\int }_S||{\mathbf{\Delta }}_S\mathbf{f}-\mathbf{R}{\mathbf{\Delta }}_S\mathbf{g}|{|}_2^{2}dA,$$

$$E_i\left(\mathbf{x},\mathbf{R}\right)=\frac{k_i}{2}{A}_i{\Vert {\mathbf{X}}_f\mathbf{c}-\mathbf{R}{\mathbf{X}}_g\mathbf{c}\Vert }_2^2,$$

$$V=\frac{1}{3}\underset{\mathcal{V}}{\iiint }\nabla \cdot \mathbf{x}d\mathbf{x}=\frac{1}{3}\underset{\partial \mathcal{V}}{\iint }{\mathbf{x}}^T\mathbf{n}d\mathbf{x},$$

$$V(\mathbf{x})=\frac{1}{3}\underset{\partial \mathcal{V}}{\iint }{\mathbf{x}}^T\mathbf{n}d\mathbf{x}=\frac{1}{9}\sum _i{A}_i({\mathbf{x}}_{i_1}+{\mathbf{x}}_{i_2}+{\mathbf{x}}_{i_3}{)}^T{\mathbf{n}}_i,$$

$$E_i(\mathbf{x},\mathbf{p})=\frac{w_i}{2}{\Vert \mathbf{p}-\mathbf{x}\Vert }_F^2\text{s.t.}\mathbf{p}\in C.$$

$$\begin{array}{rl}& \underset{\mathbf{p}}{\min }\frac{w_i}{2}{\Vert \mathbf{p}-\mathbf{x}\Vert }_F^2,\\ & \text{s.t.}C(\mathbf{p})=C(\mathbf{x}+\Delta \mathbf{x})=0\\ & \Delta \mathbf{x}=\lambda \nabla C(\mathbf{x}),\end{array}$$

$$\Delta \mathbf{x}=-\frac{C(\mathbf{x})}{{\Vert \nabla C(\mathbf{x})\Vert }_F^2}\nabla C(\mathbf{x}),$$

$$\nabla C(\mathbf{x})\approx \frac{1}{3}[{\bar{\mathbf{n}}}_1^T,...,{\bar{\mathbf{n}}}_n^T{]}^T,$$

$$E[\varphi ]={\int }_V\Psi [\mathbf{\varphi };\mathbf{X}]d\mathbf{X}={\int }_V\Psi [\mathbf{F}]d\mathbf{X}.$$

$$\mathbf{P}(\mathbf{F})=\partial \Psi (\mathbf{F})/\partial \mathbf{F}.$$

$$\Psi ({\lambda }_1,{\lambda }_2)=f({\lambda }_1)+f({\lambda }_2)+g({\lambda }_1{\lambda }_2),$$

$$\Psi =\frac{B}{4}\left(\frac{1}{2}{I}_1^2+I_1-I_2\right)+\frac{C}{8}{I}_2^2,$$

$$f({\lambda }_i)$$

$$g({\lambda }_i{\lambda }_j)$$

$$E_i={\int }_{{\mathcal{T}}_i}\Psi (\mathbf{F})d\mathbf{X}=W_i\cdot \Psi ({\mathbf{F}}_i),$$

Full text

1 Introduction

Blood is a complex suspension of red blood cells (RBCs), white blood cells and platelets, submerged in a Newtonian fluid, the plasma. Especially, the RBCs (else, erythrocytes) cover 40-45% of the whole blood volume, and their behavior has a direct impact on blood rheology. The accurate modeling of the collective transport of the cells in the plasma is of paramount importance since it can help us decipher in vivo phenomena, e.g., blood clotting, margination of platelets, or the cell-free layer (Mountrakis 2015 [1]). In addition, the deformability of RBCs is strongly linked to some pathological conditions, e.g., hereditary disorders (like spherocytosis, elliptocytosis, and stomatocytosis), metabolic disorders (like diabetes, hypercholesterolemia, and obesity), malaria, or sickle anemia as described by Tomaiuolo 2014 [2]. In more details, a human RBC has a biconcave discocyte shape with a liquid interior (cytoplasm, hemoglobin solution) of dynamic viscosity $6-10~{}cP$. A healthy red blood cell has an average surface area of $134~{}\mu m^{2}$, a volume of $94~{}\mu m^{3}$, a diameter of $7.82~{}\mu m$ and a varying thickness of $0.81~{}\mu m$ at the dimple to $2.57~{}\mu m$ at the periphery. An analytical formula for the RBC shape is proposed in Evans and Skalak 1980 [3]. The RBC membrane is a complex multi-layered structure consisting of an external lipid bilayer attached to an underlying cytoskeleton (spectrin-link network). These complicated anucleate structures transfer oxygen between blood and tissues thanks to their ability to undergo substantial deformations while maintaining their area and volume.

The absence of a universal numerical model capable of describing the viscoelastic behavior of a RBC is the motivation of this study. A universal model should fulfill a number of criteria, such as generality, robustness, accuracy and performance. By generality and accuracy, we mean the ability of the model to simulate any material and all the expected behavior as close as possible to the observed data (experiments). By robustness, we refer to the capability of the model to adapt to extreme physical circumstances such as large deformations, confined flows, or high flow shear rates, and to numerically unfavorable settings like large time steps. Last but not least, performance is a key requirement in the field of computational biophysics, since the model should be extended to cases including millions or even billions of deformable bodies. To give an example, a volume of $1~{}mm^{3}$ of blood at a hematocrit 40% contains approximately five million RBCs. For the modeling of RBCs, two main approaches are proposed in the literature. The first is a continuum-level approach, that views the membrane as a continuous medium which obeys specific partial differential equations and constitutive laws. The discretization of these equations is commonly done through a finite element method (FEM). Many researchers have successfully used this approach for the modeling of membranes, including Krüger et al. 2011 [4], MacMeccan et al. 2009 [5], Shrivastava and Tang 1993 [6]. Klöppel et al. 2011 [7] employed a two-layer model for an even more accurate description of the lipid bilayer and the underlying cytoskeleton. The aforementioned models use constitutive laws for the description of the membrane material, such as the material introduced by Skalak et al. 1973 [8], Neo-Hookean materials proposed in [9] or the material described by Yeoh 1990 [10] (for a thorough investigation of material models the reader is referred to Dimitrakopoulos et al. 2012 [11] and Siguenza et al. 2017 [12]). Among the many advantages of this approach is the guaranteed mesh-independence, which is a product of the way the equations are discretized. This modeling approach satisfies many of the above-mentioned criteria, namely generality, robustness and accuracy. However, they do not perform well in simulations of multiple deformable bodies due to their high computational cost. The second modeling approach operates at a mesoscopic level, and represents surface properties through a mass-spring system. More precisely, mesoscopic modeling of the RBCs spans from the spectrin-link to the coarse-grained spectrin-link levels. The spectrin-level approach tries to imitate the physics of the cytoskeleton of a RBC as a fine-grained mass-spring model. The spectrin level was introduced by Discher et al. 1998 [13] and extended further by Dao et al. 2003 & 2006 [14, 15] and Li et al. 2005 [16]. Though, this approach is limited by the high computational cost, given that it requires $23~{}867$ surface vertices [17] for the simulation of the whole spectrin network. For this reason, the coarse-grained models of Dupin et al. 2007 [18], Pivkin and Karniadakis 2008 [17], Fedosov et al. 2010 [19], and Reasor et al. 2011 [20] have gained substantial attention during the last years, where the surface discretization demands about $500$ vertices. Currently, state-of-the-art solvers like Hemocell [21, 22] and the solver of Fedosov [19] are based on this modeling approach. These coarse-grained methods outperform all the others in terms of computational time. Nevertheless, since they are mass-spring systems they inherit limitations on generality, robustness and accuracy.

As far as the simulation of the blood plasma is concerned, there exists a plethora of mature CFD approaches. For our study, we make use of the lattice Boltzmann method (LBM) [23] which indirectly solves the Navier-Stokes equations. Another candidate method widely used for the simulation of blood flows, is dissipative particle dynamics (DPD) [17, 19]. Both are particle-based methods which deploy either collision and streaming operations or Newtonian laws, respectively, for the time advancement of the fluid.

The coupling of the fluid with the solid phase is a critical point for an accurate and stable simulation. Peskin 1972 [24] developed the immersed boundary method (IBM) to model blood flow in combination with moving heart valves. The strength of the IBM is that the fluid solver does not have to be modified except from the addition of a forcing term, and the fluid mesh does not need to be adjusted dynamically. Moreover, the deformable body and its discrete representation do not need to conform to the discrete fluid mesh, which leads to a very efficient fluid-solid coupling. In the context of lattice Boltzmann, an alternative solution would be the deployment of moving bounce back nodes [23]. With this approach, the computational cost however steeply increases, given the need to track the lattice nodes that transition from solid to fluid state.

In a typical numerical framework for blood flow, the computational time is dominated by the structural solver for the deformable RBCs. Even with the relatively cheap mass-spring systems, Závodsky et al. 2017 [22] report for example that the deformable bodies solver constitutes over 95% of the total computational time. Since the numerical models for Newtonian fluid flow and for the fluid-solid coupling are more mature than the ones for the physics of RBCs, our main focus is to develop a novel solution which strictly fulfills the criteria of generality, robustness, accuracy and performance without any compromise. We propose a novel approach for deformable viscoelastic bodies based on nodal projective finite elements method (npFEM). Inspired by the research of Liu et al. 2017 [25] and Bouaziz et al. 2014 [26], our solver extends these theories to the field of computational biophysics. The expression “nodal” refers to the mass lumping technique, in which both the masses and the forces are lumped on the vertices of the discretized body, and therefore the finite elements act like generalized viscoelastic springs. The term “projective” stands for the use of specially designed potential energies that help us build a fast solver based on quasi-Newton optimization techniques [27, 28]. Lastly, we discretize the continuous potential energies using linear basis functions, leading to the standard Galerkin method [29] in the field of FEM. Our solver inherits the versatility and robustness of FEM and is almost as fast as plain mass-spring systems. To our knowledge, this is the first solver bridging the gap between the two approaches, continuum and mesoscopic, in the field of RBC modeling.

The present work is organized as follows: In section 2 we thoroughly introduce the npFEM method focusing on the modeling of RBCs, the lattice Boltzmann method and the immersed boundary method. Following, in section 3 we put all the different pieces together and show how the complete computational framework operates. Finally, in section 4 we present extensive validations of our computational framework and give insights about its performance and scalability to systems of multiple blood cells.

2 Methods

The focus of this study is to build a computational framework for fully resolved 3D blood flow. The constituents of this framework are the solid body solver (section 2.1), the fluid solver (section 2.2), and the fluid-structure interaction (FSI) (section 2.3). The solid body solver based on the nodal projective finite elements method (npFEM) is our main contribution to the field. Given the volume fraction of red blood cells compared to the other cells (white, platelets), and consequently their importance on blood rheology, we develop the npFEM solver around RBC simulations. The fluid solver, which is used for the simulation of the blood plasma, is based on the lattice Boltzmann method (LBM) as implemented in the open source software Palabos [30]. The coupling between the two solvers is done via an immersed boundary method known as multidirect forcing scheme [31].

2.1 Nodal projective FEM (npFEM)

This section provides a complete review of the methods introduced in the research of Liu et al. 2017 [25] and Bouaziz et al. 2014 [26], and some novel extensions proposed to adapt them to the field of RBC modeling. A validation of these adaptations is provided in later sections, see results section 4. A red blood cell is represented by its triangulated membrane, and its total mass is lumped on the membrane vertices (including the interior fluid). This section presents how we realize the dynamics of deformable bodies, as well as all the potential energies that describe the behavior of a red blood cell, namely local area and global volume conservation, bending and material.

Let us assume a surface mesh consisting of $n$ vertices with positions $\bm{x}\in\mathbb{R}^{n\times 3}$ and velocities $\bm{v}\in\mathbb{R}^{n\times 3}$. The evolution of the system in time follows Newton’s law of motion. The external forces are defined as $\bm{F}_{ext}\in\mathbb{R}^{n\times 3}$ and come from the fluid-solid interaction, while the internal forces are $\bm{F}_{int}\in\mathbb{R}^{n\times 3}$. The internal forces are position dependent, i.e., $\bm{F}_{int}(\bm{x})=-\sum_{i}\nabla E_{i}(\bm{x})$, where $E_{i}(\bm{x})$ is a scalar discrete elemental potential energy. Moreover, the summation of all the elemental potential energies of any type results in the total elastic potential energy of the body. Conservative internal forces, i.e, derived from potentials, are closely related to the assumption of hyperelastic materials [32]. The implicit Euler time integration leads to the following advancement rules (subscripts $n$ and $n+1$ refer to time $t$ and $t+h$, respectively):

[TABLE]

where $\mathbf{M}\in\mathbb{R}^{n\times n}$ is the mass matrix, $h$ is the time step and $\mathbf{C}\in\mathbb{R}^{n\times n}$ is the damping matrix acting like a proxy for viscoelasticity (Rayleigh damping, see section 2.1.4). The mass matrix is built by lumping the total mass of the body on the mesh vertices, resulting in a diagonal structure. The lumping can be done either by equally distributing the mass or by weighting, based on the corresponding area per vertex. Combining the above equations we derive

[TABLE]

where $\widetilde{\mathbf{M}}=\mathbf{M}+h\mathbf{C}$ and $\bm{y}_{n}=\bm{x}_{n}+h\widetilde{\mathbf{M}}^{-1}\mathbf{M}\bm{v}_{n}+h^{2}\widetilde{\mathbf{M}}^{-1}\bm{F}_{ext}$. Equation (3) can be turned into an optimization problem [26] as

[TABLE]

where $||\cdot||_{F}$ is the Frobenius norm. Indeed, setting the derivative of equation (4) to zero (thereby minimizing the objective function), we recover Newton’s second law of motion. The solution to the minimization problem gives $\bm{x}_{n+1}$. The choice of implicit Euler integration serves a dual purpose. Firstly, it provides unconditional stability for arbitrary time step. Secondly, this integration scheme is characterized by some numerical dissipation linked to the time step (the smaller the time step the smaller the numerical dissipation). This dissipation enhances the stability of the coupling (fluid & solid) by reducing energy oscillations that could potentially lead to instabilities. This dissipation term acts like the viscous dissipation provided by Rayleigh damping (matrix $\mathbf{C}$).

A remark at this point is that the variational form of implicit Euler (4) is just a reformulation of Newtonian dynamics, which in their turn can be viewed from the Hamiltonian dynamics perspective. Nevertheless, to accurately capture the trajectory and deformation of the bodies (dynamics) we need to guarantee the conservation of linear & angular momenta. For this, one has to be cautious with the selection of the potential energies and the viscoelastic terms. According to Noether’s theorem, momenta are conserved when the elastic energies are rigid motion invariant. For the viscoelastic terms, we give an extended discussion on the conservation of momenta in section 2.1.4. Taking into consideration all the above, we guarantee that we resolve as accurately as possible the dynamics of the simulated bodies, reproducing the correct physics as dictated by Newton’s laws.

2.1.1 Extended projective dynamics method and quasi-Newton optimization

Projective dynamics method [26] requires a special form of potential energies, based on the notion of constraint projection:

[TABLE]

where $\mathcal{M}_{i}$ is a constraint manifold defined by a desired undeformed state, $\bm{p}_{i}=\mathbf{S}_{i}\bm{p}\in\mathbb{R}^{m\times 3}$ is an auxiliary projection variable on $\mathcal{M}_{i}$ and $\mathbf{G}_{i}$ is a mapping from $\bm{x}$ to an element-wise deformation representation which depends on the mesh topology and the initial configuration. For an elastic body, there are many elemental energies of various types, therefore various projections. For compactness, we stack all the projections into $\bm{p}\in\mathbb{R}^{m\times 3}$ and define binary selector matrices $\mathbf{S}_{i}\in\mathbb{R}^{m\times m}$ such that $\bm{p}_{i}=\mathbf{S}_{i}\bm{p}$. The dimension $m$ is equal to $\sum_{i}c_{i}$ (summation over all quadratic energies), where $c_{i}$ is one for the bending energy, two for the local area energy and $n$ for the global volume conservation, i.e., it depends on the corresponding energy (see section 2.1.2). More in details, the matrix $\mathbf{G}_{i}\in\mathbb{R}^{m\times n}$ formulates the deformation descriptor and places it in the involved rows and columns, where the columns correspond to the involved vertices of the $i$-th energy and the rows refer to the corresponding projection $\bm{p}_{i}$. The term “constraint” is used because the energy penalizes the deviation from the undeformed state, and thus constraints the body from deformation, and the term “projection” corresponds to finding the minimum distance from the current state onto the manifold.

Generalizing the above, all the possible undeformed configurations define the so-called constraint manifold $\mathcal{M}$ (Cartesian product of the individual constraint manifolds $\mathcal{M}_{i}$), which mathematically is formulated as the zero level-set of the total potential energy. From a geometric point of view, the potential energy is a distance $d$ which quantifies how far the deformed state $\mathbf{G}\bm{x}$ is from the manifold $\mathcal{M}$ (rest state). The projection of the current state onto $\mathcal{M}$ is denoted by $\bm{p}$. There can be many user-defined manifolds. For example, a manifold for volume preservation does not guarantee area conservation and therefore, we have to define another manifold to guarantee the latter constraint. Recapitulating the above, the total potential energy can be formalized as

[TABLE]

The minimum is achieved when $\bm{p}$ is the projection of $\mathbf{G}\bm{x}$ onto $\mathcal{M}$. Figure 1 summarizes the definition of the distance and its correspondence to a potential energy. The generic potential energy of equation (6) refers to a discretized energy. Later, we show how from continuous energies using linear basis functions, we retrieve this generic formulation.

The key idea of projective dynamics is to exploit these specialized potential energies and reformulate the variational form of implicit Euler as

[TABLE]

where $g$ is the objective function whose minimum corresponds to $\bm{x}_{n+1}$ and $w_{i}=A_{i}k_{i}$ is a weighting term acting like a “stiffness” that converts the distance to an energy term ($A_{i}$ is the area of the element and $k_{i}$ is the stiffness factor). As a reminder, $||A||^{2}_{F}=tr(A^{T}A)$. The projections $\bm{p}$ are essentially functions of $\bm{x}$ and thus developing even further equation (5) into the variational form (4), we get

[TABLE]

where $\mathbf{L}=\sum w_{i}\mathbf{G}_{i}^{T}\mathbf{G}_{i}$, $\mathbf{J}=\sum w_{i}\mathbf{G}_{i}^{T}\mathbf{S}_{i}$ and $\mathbf{S}=\sum w_{i}\mathbf{S}_{i}\mathbf{S}_{i}^{T}$.

Seeking for greater material expressivity, and thus more general potential energies, we can straightforwardly extend equation (7) as

[TABLE]

where PD stands for energies that follow the quadratic formulation of projective dynamics, see equation (6), and non-PD for arbitrary energies. A remark is the absence of a weighting term in the non-PD energies, which is because the stiffness is already integrated into them. The minimization of $g$ is performed using a quasi-Newton technique with a well-suited Hessian approximation, which is deeply rooted in the projective dynamics part. The first step is to compute the gradient of the objective function $g(\bm{x})$ as defined by equation (9). An in-depth derivation of the gradient (regarding the PD part) and an explanation on why the derivative of $\frac{1}{2}tr\left(\bm{p}^{T}(\bm{x})\mathbf{S}\bm{p}(\bm{x})\right)$ vanishes can be found in [25, 33]. The derivative is given by

[TABLE]

The projective dynamics energies demand first the computation of $\bm{p}$ for a given $\bm{x}$. This step can be massively parallelized given the small computational cost of the projection per element. The objective function could be equivalently minimized by Newton’s method. Newton’s method would proceed by computing the Hessian matrix $\nabla^{2}g(\bm{x})$ and then using it to compute a descent direction as $d=-(\nabla^{2}g(\bm{x}))^{-1}\nabla g(\bm{x})$. However, Liu et al. 2017 [25] suggested replacing the Hessian by $\widetilde{\mathbf{H}}=\widetilde{\mathbf{M}}/h^{2}+\mathbf{L}$ and thus opting for a quasi-Newton approach avoiding the expensive computation of the Hessian at every time step and the possible definiteness fixes to guarantee that $d$ is a descent direction. Using only quadratic energies, the minimization of $g$ could be realized by an alternating (local/ global) solver [26], which reveals how $\widetilde{\mathbf{H}}$ is formed. A detailed presentation can be found in A. The approximated Hessian plays a major role in the fast convergence of our solver and in the overall computational efficiency of the npFEM solver without the expensive computation of the true Hessian. Furthermore, matrix $\widetilde{\mathbf{H}}$ is constant, symmetric positive definite and thus can be prefactorized once at the beginning, e.g., using Cholesky factorization, allowing for very fast computation of the descent direction $d=-(\widetilde{\mathbf{H}})^{-1}\nabla g(\bm{x})$.

The non projective dynamic energies could contribute as well to the formulation of the Hessian approximation by considering their contribution to $\sum w_{i}\mathbf{G}_{i}^{T}\mathbf{G}_{i}$, where $\mathbf{G}_{i}$ is the deformation gradient operator for the $i$-th triangular element and $w_{i}=A_{i}k_{i}$ ($A_{i}$ is the area of the element and $k_{i}$ a stiffness factor which can be computed as described in [25]). If on the other hand there is no contribution from the non-PD part, then the iterations for convergence are slightly increased.

The overall convergence is affected by the selected potential energies. The more the PD energies, the faster the convergence, which is the case for our study since out of the four different potentials (area, volume, bending, material) only the material energy does not fall into the PD category, see sections 2.1.2 & 2.1.3 for more details. The reason for this, is that the Hessian approximation is mainly built upon the PD part.

It is worth noting that the size of our matrices is $n\times n$ instead of $3n\times 3n$, typically found in other FEM solvers [5], leading to an at least 3$\times$ speedup compared to them.

Consequently, moving towards the descent direction, this iterative procedure guarantees to converge to the minimum of $g$, i.e., $\bm{x}_{n+1}$. A recapitulation of the algorithm is presented in section 2.1.5.

2.1.2 Quadratic Energies to simulate RBC membrane

In this section, we develop potential energies suitable for projective dynamics and critical for the simulation of the red blood cell membrane, namely local area conservation, bending and global volume preservation. The analysis reveals that the above-mentioned energies can be written in the quadratic form described by equation (5). The goal is to develop discrete energies and compute the auxiliary variables, i.e., the projections, in order to quantify the energy term itself. Since the current state $\bm{x}$ is known, the computation of $\bm{p}$ results in the quantification of the energy $E$ from equation (5).

Let the undeformed surface be a differentiable 2-manifold surface $S$ embedded in $\mathbb{R}^{3}$. Let us define piecewise linear coordinate functions of the undeformed and deformed surface, $\bm{g}:S\rightarrow\mathbb{R}^{3}$ and $\bm{f}:S\rightarrow\mathbb{R}^{3}$, respectively. By introducing a set $M$ of desired point-wise transformations $\mathbf{T}$, we can formulate a potential energy [26] as

[TABLE]

where $\nabla_{S}$ is the gradient operator on $S$. The set $M$ determines all allowed rest configurations $\mathbf{T}\nabla_{S}\bm{g}$. The energy $E\left(\bm{f},\mathbf{T}\right)$ is discretized over triangular elements using piecewise linear hat functions [34], and the integral is transformed to a sum over triangle potentials:

[TABLE]

where $A_{i}$ is the triangle area, $\mathbf{D_{s}}=[\bm{x}_{j}-\bm{x}_{i},\bm{x}_{k}-\bm{x}_{i}]\in\mathbb{R}^{2\times 2}$ contains the triangle edges of the deformed state embedded in 2D and $\mathbf{D_{m}}$ corresponds to the undeformed state. At this point, we have introduced a generic quadratic energy similar to what is described by equation (5). Based on the FEM reminder of B, the generic quadratic energy defined above can be rewritten as

[TABLE]

where $\mathbf{F}\in\mathbb{R}^{2\times 2}$ is the deformation gradient of the $i$-th triangular element. In this case, the deformation gradient acts like a deformation descriptor, and $\mathbf{T}$ corresponds to the projection on the manifold defined by the set $M$.

Area preservation energy. Let us define the set $M$ as the set of matrices with determinant 1, i.e., $det(\mathbf{T})=1$. In mathematical terms, this requirement is equivalent to $M\equiv SL(2)$, where the special linear group of degree n, $SL(n)$, is the set of $n\times n$ matrices with determinant 1. Let us consider the singular value decomposition (SVD) of the deformation gradient of a triangular element

[TABLE]

where $\mathbf{\Sigma}$ is a diagonal matrix storing the singular values or principal stretches $\lambda_{1},\lambda_{2}$ of $\mathbf{F}$. Area conservation of the element is equivalent to

[TABLE]

Consequently, the generic quadratic energy is rewritten as

[TABLE]

which is the distance of the current state $\mathbf{\Sigma}$ from the manifold defined by $SL(2)$. A rigorous solution to this projection $\mathbf{\Sigma}^{*}$ can be found in the Appendix of [26]. Essentially, this energy penalizes how far the element is from the state that conserves the area $\mathbf{\Sigma}^{*}$.

Bending Energy. The RBC membrane simulation demands a bending energy, commonly formulated by employing the dihedral angles across edges. The bending energy used in this study [26] measures the squared difference of mean curvatures $H$ as

[TABLE]

where $\mathbf{R}\in SO(3)$ are rotation matrices that best match the initial to the deformed state, $\bm{\Delta}_{S}$ is the Laplace-Beltrami operator on $S$. As a reminder, this operator applied on the coordinate function $\{\bm{g},\bm{f}\}$ of the surface, gives $\bm{\Delta}_{S}\{\bm{g},\bm{f}\}=-2H_{\{\bm{g},\bm{f}\}}\bm{n}$, where $H$ is the mean curvature and $\bm{n}$ is the surface normal. The discretized version [26] of the bending energy per vertex $i$, using piecewise linear hat functions, is

[TABLE]

where $A_{i}$ is the Voronoi area of the vertex $i$, $\bm{X}_{f}$ and $\bm{X}_{g}$ contain the one-ring edges of the vertex in the current and initial configuration, respectively, and $c$ stores the common cotangent weights divided by the Voronoi area. The projection $\mathbf{R}$ has a closed form solution and is presented in the Appendix of [26]. For a thorough analysis on the different bending energy algorithms, consider the work of Guckenberger et al. 2016 [35].

Global volume conservation. The volume $V$ of a body embedded in 3D can be computed employing the divergence theorem as

[TABLE]

where $\mathcal{\partial V}$ is the boundary of the 3D body $\mathcal{V}$ (RBC membrane in our case) and $\bm{n}$ is the surface normal. Given the tessellation of the membrane in triangular elements, we can further simplify [36] the volume as

[TABLE]

where $A_{i}$ is the area and $i_{1},i_{2},i_{3}$ are the vertex indices of the $i$-th element. Having defined the above toolset, we can now formulate the global volume conservation manifold as $C=V(\bm{x})-V_{0}=0$, where $V_{0}$ is the volume at rest. Therefore, the global volume conservation energy can be written as

[TABLE]

The deformation descriptor is the absolute positions of the deformed body, and the projection corresponds to a state vector which respects the volume conservation.

To find the projection $\bm{p}$, we reformulate the above energy into a minimization problem as

[TABLE]

the second constraint that restricts $\Delta\bm{x}$ to be in the direction of $\nabla C$ is a requirement for linear and angular momentum conservation [36]. The solution to the above optimization problem using Lagrange multipliers is [36, 26]

[TABLE]

while the gradient can be approximated [36] by

[TABLE]

where $\bar{\bm{n}}_{k}=\sum A_{j}\bm{n}_{j}$ is the sum of the area weighted normals of all elements which contain vertex $k$.

Global area conservation could be imposed by a similar approach. Nevertheless, after extensive testing we found out that only the local area preservation term is enough to describe the RBC behavior, without increasing more the computational cost. Essentially, since the area preservation is a property of the material, it is sufficient to encode it locally, through the elements that compose the body.

The energies developed in this section follow the quadratic formulation described by equation (5). However, to build the stacked system of projections as described by equation (8), the user has to place carefully each deformation descriptor and projection to the sparse matrices $\mathbf{G}_{i}\bm{x}$ and $\bm{p}_{i}$. Additionally, the descriptors and projections isometrically embedded in $\mathbb{R}^{2}$, have to be re-embedded back in $\mathbb{R}^{3}$, since the equations of motion refer to three dimensions.

2.1.3 Non-PD energies & internal forces

There are many energies following the above-mentioned quadratic distance formulation but the simulated materials are quite limited. Even classical models from continuum mechanics, such as Neo-Hookean, St. Venant-Kirchhoff, Mooney-Rivlin, are not supported. More importantly, the RBC falls into the latter category and for this reason, we develop the non projective dynamic energies and forces.

Consequence of elastic deformation is the accumulation of potential energy in the body. Our focus is on hyperelastic materials for which the potential energy $E[\bm{\phi}]$ is determined by the deformation function $\bm{\phi}$ (see B) and is independent of the deformation history. To describe a deformed body, we introduce the energy density function $\Psi[\bm{\phi};\bm{X}]$ which measures the potential energy per unit undeformed volume/ area on an infinitesimal domain $dV$ around the material point $\bm{X}$. The total potential energy is

[TABLE]

Every material follows a specific constitutive model which links the stresses with the strains. There is a variety of stress descriptors; for our discussion we use the $1^{st}$ Piola-Kirchhoff stress tensor $\mathbf{P}$ which for hyperelastic materials is

[TABLE]

The assumption of isotropy and rotation invariance allows us to express the energy density function $\Psi$ as a function of the principal stretches of the deformation gradient ($\lambda_{1},\lambda_{2}$ for triangular elements). To simplify the material space, we use the Valanis-Landel hypothesis [37] which assumes that

[TABLE]

where $f,g$ are scalar non-linear elastic strain energy density functions. The original formulation of the hypothesis is for solid bodies (typically rubber-like materials) and not for surfaces, thus involving the three principal stretches ($\lambda_{1},\lambda_{2},\lambda_{3}$) of the deformation gradient and three scalar functions ($f,g,h$). It covers a vast amount of materials, such as St. Venant-Kirchhoff, Neo-Hookean and Mooney-Rivlin.

The well-known energy introduced by Skalak et al. 1973 [8] for the description of the membrane of a red blood cell falls into the category of the Valanis-Landel materials. Given $I_{1}=\lambda_{1}^{2}+\lambda_{2}^{2}-2$ and $I_{2}=\lambda_{1}^{2}\lambda_{2}^{2}-1$ the energy density function is

[TABLE]

where $B$ and $C$ are membrane material properties. The above energy density function can be alternatively written in the Valanis-Landel form as

[TABLE]

Given that we enforce area conservation through the projective dynamics part, i.e., $\lambda_{1}\lambda_{2}=1$, $g$ can be neglected and only $f$ should be considered.

Since we are dealing with linear finite elements, the deformation gradient is constant over each element ($\mathcal{T}_{i}$), see B, and this observation reduces the elemental potential energy to

[TABLE]

where $W_{i}$ is the undeformed area of the element. The contribution of element $\mathcal{T}_{i}$ to the elastic forces is $\bm{F}_{int}^{non-PD}=-\nabla E_{i}$ which can be further simplified for a triangular element to [32]

[TABLE]

where $\vec{f_{1}},\vec{f_{2}},\vec{f_{3}}$ are the internal forces on the vertices of $\mathcal{T}_{i}$. Moreover, the $1^{st}$ Piola-Kirchhoff stress tensor can be further simplified to [38]

[TABLE]

where $\mathbf{P}(\mathbf{\Sigma})=\text{diag}\left(\frac{\partial\Psi}{\partial\lambda_{1}},\frac{\partial\Psi}{\partial\lambda_{2}}\right)$.

Stability. Special care is needed when designing new materials to satisfy Drucker’s first stability criterion (also known as Hill’s stability criterion), which requires a monotonic increase of strain energy density with increase in strain [39]. In practice, this stability criterion is satisfied by demanding $f^{\prime\prime},g^{\prime\prime}>0$ [38].

Modification on Skalak’s Energy. This energy behaves quite well for small deformations as we present in section 4, but for larger deformations it fails to robustly describe the mechanics of the RBC membrane. For this reason, we enhance the strain hardening behavior by adding one more term in $f$. The modified version is as follows:

[TABLE]

where $D$ is a material property like $B$. The modified Skalak’s energy satisfies Drucker’s stability criterion and has a superior behavior than the original one, see section 4. The correction term has zero contribution at rest, where $\lambda_{i}=1$, and its form was defined through fitting on the available experimental data.

2.1.4 Membrane viscoelasticity

Membrane viscoelasticity is tackled in a two-way approach. Initially, through the addition of the Rayleigh dissipation term in Newton’s equation (2) and secondly through a post-convergence step.

Typically, the Rayleigh damping matrix is a blend of the mass and stiffness matrices

[TABLE]

where $\mathbf{M}\in\mathbb{R}^{n\times n}$ is the mass matrix, $\mathbf{K}\in\mathbb{R}^{n\times n}$ is the stiffness matrix and $\alpha_{D},\beta_{D}$ are coefficients. In our analysis we do not construct explicitly the stiffness matrix (Hessian of the elastic potential energy), and thus we are bypassing this issue by using matrix $\mathbf{L}=\sum w_{i}\mathbf{G}_{i}^{T}\mathbf{G}_{i}$ instead [40]. Essentially, we make use of a part of the Hessian approximation instead of the Hessian itself. A critical point to consider is that Rayleigh damping can affect the rigid body motions (translation, rotation) if the damping ratio is high. The damping ratio ($\zeta$) is defined as $\zeta(\omega_{n})=0.5(\alpha_{D}\omega_{n}^{-1}+\beta_{D}\omega_{n})$, where $\omega_{n}$ is the $n$-th modal frequency [5]. In our study, we do an approximate modal analysis by solving $\mathbf{L}\bm{x}=\bm{\omega}^{2}\mathbf{M}\bm{x}$. For our analysis, we set $\alpha_{D}=0$ and make sure that $\zeta(\omega_{n})\ll 1$ for every mode. In this way, we guarantee that the trajectory of the simulated bodies follows the proper dynamics and that the momenta are conserved.

Due to the above restriction of Rayleigh damping, we cannot fully emulate the viscoelastic behavior of the RBC membrane. This is why we introduce the post-convergence step that complements Rayleigh damping, does not affect the rigid body dynamics, and also conserves linear & angular momenta. The key idea is to calculate the difference between the vertex velocity $\bm{v}_{i}$ and its velocity in the best-fit rigid body motion. Then, we dissipate partially the non-rigid-body part. The algorithm [41] to do so is summarized in the following lines:

$\bm{x}_{cm}=(\sum_{i}\bm{x}_{i}m_{i})/(\sum_{i}m_{i})$ $\triangleright$ center of mass (cm)

$\bm{v}_{cm}=(\sum_{i}\bm{v}_{i}m_{i})/(\sum_{i}m_{i})$ $\triangleright$ velocity of cm

$\mathbf{L}=\sum_{i}\bm{r}_{i}\times(m_{i}\bm{v}_{i})$ $\triangleright$ angular momentum

$\mathbf{I}=\sum_{i}\widetilde{\bm{r}}_{i}\widetilde{\bm{r}}_{i}^{T}m_{i}$ $\triangleright$ inertia matrix

$\bm{\omega}=\mathbf{I}^{-1}\mathbf{L}$ $\triangleright$ angular velocity

for all vertices  i do

     $\Delta\bm{v}_{i}=\bm{v}_{cm}+\bm{\omega}\times\bm{r}_{i}-\bm{v}_{i}$

     $\bm{v}_{i}\leftarrow\bm{v}_{i}+\kappa_{\text{damping}}\Delta\bm{v}_{i}$

end for

where $\bm{r}_{i}=\bm{x}_{i}-\bm{x}_{cm}$, $\widetilde{\bm{r}}_{i}\in\mathbb{R}^{3\times 3}$ is a matrix with the property $\widetilde{\bm{r}}_{i}\bm{v}=\bm{r}_{i}\times\bm{v}$, and $\kappa_{\text{damping}}\in[0,1]$ is the damping coefficient. This dissipation term is deployed after the convergence to $\bm{x}_{n+1}$ (see algorithm 1).

An important remark to avoid confusion is that in our study we are not simulating a viscoelastic material, i.e., the potential energy to be dependent on the strain rate, but we are emulating the viscoelastic behavior through viscous dissipation (proxy for viscoelasticity). This approach is proven to have equivalent rheological behavior with the generalized Maxwell model as proposed by Semblat 1997 [42]. Similar approach is followed by Fedosov et al. 2010 [19] & Hemocell [21].

2.1.5 Deformable bodies solver: npFEM algorithm

In this section, we summarize the algorithm for the deformable bodies, focused on the simulation of red blood cells (algorithm 1 describes one time step of the npFEM solver). Its extension to any other type of cell or body is based on the same principles. The advancement of the body from time $t$ to $t+1$ is done through the minimization of equation (9). The main constituents of this equation are two categories of energies. The former follows the projective dynamics setup with quadratic energies. The latter are energies that do not belong in any particular category, but they offer a more advanced description of the material space.

Equation (9) is a typical non-linear unconstrained minimization problem for which we employ a quasi-Newton technique. To further speed-up the convergence, we deploy the LBFGS algorithm [43] that takes into consideration previous iterates and gives a better Hessian approximation based on the initial one. Moreover, given the positive definiteness of $\widetilde{\mathbf{H}}$, $\bm{d}(\bm{x})$ is guaranteed to be a descent direction. However, to guarantee a sufficient and optimal decrease, we implement a popular inexact line search technique based on Armijo and curvature conditions, collectively known as Wolfe conditions. For a detailed overview on quasi-Newton techniques consider the seminal work of Nocedal & Wright [43].

Regarding the interior fluid of the red blood cell (cytoplasm), we simulate it implicitly by considering the global volume conservation (incompressibility of cytoplasm) and through the viscoelasticity of the full body as described in section 2.1.4. The decision not to simulate the interior fluid of RBCs is not because of some limitation of the computational framework but rather a matter of computational efficiency. We could have simulated the interior fluid through Palabos and subsequently to tune the viscoelastic terms to correspond to the membrane only. In this case, the global volume conservation energy has to be removed, since it would be implicitly conserved through the fluid solver, while the rest of the energies (area, bending, material) do not have to change. Moreover, from the literature there is no indication that the simulation of the interior fluid plays a crucial role on accurately reproducing the collective transport of blood cells, and for this reason we have decided to follow the more lightweight implementation path.

2.2 Lattice Boltzmann Method (LBM)

The lattice Boltzmann method (LBM) is employed to solve indirectly the Navier-Stokes equations. Densities of virtual particles, also known as populations, move on a regular grid/ lattice and collide at the lattice nodes. The LBM originally grew out of the lattice gas models [44, 45] and while lattice gases track particles, LBM tracks their statistical distribution. LBM is a second-order accurate solver for the weakly compressible Navier-Stokes equation, where the weak compressibility refers to errors that amplify as Mach number tends to one [23].

The basic quantity of the LBM is the discrete-velocity distribution function $f_{i}(\bm{x},t)$, also known as particle populations. It represents the density of particles with velocity $\left\{\bm{c}_{i}\right\}_{i=0}^{q-1}$ at position $\bm{x}$ and time $t$. The number of populations $q$ per lattice node defines the resolution of the velocity space discretization. The points $\bm{x}$ are equally spaced in a regular 3D grid with a lattice spacing at $\Delta x$. By discretising the Boltzmann equation (which describes the evolution of the continuous distribution function $f$) in velocity space, physical space and time, we get the lattice Boltzmann equation (LBE)

[TABLE]

which describes the collision and propagation of the mesoscopic particle packets. The velocities $\left\{\bm{c}_{i}\right\}_{i=0}^{q-1}$ are specially chosen such that $\bm{c}_{i}\Delta t$ point to a neighboring lattice site, with the consequence that populations are never trapped between nodes. For the collision operator $\Omega$, we are using the so-called Bhatnagar-Gross-Krook (BGK) approach [46]

[TABLE]

which relaxes the populations towards an equilibrium $f^{eq}$ at a rate determined by the relaxation time $\tau$. The dimensionless relaxation parameter $\tau$ is linked to the speed of sound $c_{s}$ and the physical kinematic viscosity $\nu$ by $\nu=c_{s}^{2}(\tau-1/2)\Delta t$, where $c_{s}=\sqrt{1/3}\Delta x/\Delta t$ ($\Delta x,\Delta t$ space and time discretization, respectively). The equilibrium populations are given by

[TABLE]

where $w_{i}$ are lattice weights and $\bm{u}$ is the macroscopic velocity at $\bm{x}$ and time $t$. In the present study, we use a 3D model with $q=19$ velocities, namely D3Q19. The lattice structure, the velocities $c_{i}$ and the lattice weights $w_{i}$ are thoroughly described in Krüger and colleagues 2017 [23]. A sketch of the D3Q19 lattice is shown in Figure 2.

The macroscopic properties of the fluid such as the density $\rho$, the velocity $\bm{u}$ and the hydrodynamic stress tensor $\bm{\sigma}$ can directly be recovered at each lattice node per time step as follows

[TABLE]

where $f_{i}^{neq}=f_{i}-f_{i}^{eq}$.

To further account for the presence of immersed objects, a forcing term $\bm{f}_{imm}$ should be added to the system. The method used in the present study is proposed in Shan & Chen 1993 [47] and applies a correction to the momentum term of the equilibrium populations. The method is summarized in algorithm 2.

2.3 Immersed Boundary Method (IBM)

The coupling between the red blood cells membrane and the blood plasma is realized through the immersed boundary method (IBM). Essentially, the IBM imposes a no-slip boundary condition, so that each point of the surface and the ambient fluid move with the same velocity. The original version of this method was developed by Peskin [24, 48] to study heart valves. The key idea is that the surface of the deformable object is viewed as a set of Lagrangian points, which do not have to conform with the fluid mesh/ lattice. The fluid feels the presence of the body only via a force field which is added to the momentum conservation equations (see algorithm 2). The interaction between the off-lattice marker points and the fluid is done via interpolation stencils, i.e., regularized Dirac delta functions. There are many variations of the IBM [49] and in this study we are using a modified version of the multi-direct forcing method introduced by Wang and colleagues 2008 [50], in a form close to the one proposed by Ota et al. 2012 [31].

Our goal is to compute a forcing term $f_{imm}$ along the immersed boundary and apply it as body force to the fluid. Let us denote with $\bm{x}$ and $\bm{X}_{k}$ the Eulerian and Lagrangian points respectively, with $\bm{U}_{k}$ the velocity of the Lagrangian point $k$ as computed by the solid solver (see algorithm 1) and with $\bm{u}^{*}(\bm{x},t)$ the fluid velocity as computed by equation (41). Using the interpolation stencil, the velocity on the Lagrangian points $\bm{X}_{k}$ is given by

[TABLE]

where the Dirac function

[TABLE]

is the product of the 1D weight function

[TABLE]

The body force $\bm{f}_{imm}$ is determined by the following iterative procedure [31]:

1. Step 0.

Initial estimation of the body force on the Lagrangian points by

[TABLE] 2. Step 1.

Compute the body force on the Eulerian grid at the $l$th iteration by considering the nearby Lagrangian points

[TABLE]

where $A_{k}$ is the corresponding surface area of the Lagrangian point $k$. 3. Step 2.

Correct the lattice velocity

[TABLE] 4. Step 3.

interpolate the corrected velocity at the Lagrangian points

[TABLE] 5. Step 4.

update the body force of Lagrangian point $\bm{X}_{k}$

[TABLE] 6. Step 5.

Repeat from Step 1.

The above iterative procedure can be applied for as many cycles as the application demands. For the case of a single red blood cell, even one cycle gives satisfactory results. The body force computed by the last cycle on the Eulerian grid is the $\bm{f}_{imm}$ of algorithm 2.

The immersed body force term $\bm{f}_{imm}$ could be used as the $\bm{F}_{ext}$ for the solid solver. However, we make use of the hydrodynamic stress tensor, equation (42), for higher accuracy. To compute the force at the Lagrangian point $k$, we project the stress tensor $\bm{\sigma}$ onto the surface normal $\bm{n}_{k}$ and we get

[TABLE]

At this point we should highlight that in the immersed boundary method the same fluid covers the whole computational domain. This contradicts with the existence of interior fluid in the RBCs, of different viscosity and density, than the exterior blood plasma. However, for our study we disregard the contribution of the interior fluid on $\bm{F}_{k}^{ext}$ by summing only the points $\bm{x}$ that are on the outside region of the bodies. For increased stability, we found out that instead of summing the contribution of the exterior Eulerian points $\bm{x}$, it is sufficient to compute the force by considering the stress tensor of the most distant point from the Lagrangian point. Regarding the interior fluid, we follow the modeling approach described in section 2.1.5, taking into account both the incompressibility and viscosity of cytoplasm, and thereby leaving the interior flow generated by the IBM as an artifact flow with no significance for our simulations.

3 Computational Framework

After the thorough presentation of each constituent of our framework, namely solid, fluid and FSI solvers, it would be appropriate to put all these pieces together and see the exact order of execution for one time step:

Compute the macroscopic fluid properties: density ($\rho$), momentum ($\rho\bm{u}$) and stress tensor ($\bm{\sigma}$) from equations (40),(41),(42) 2. 2.

Using the stress tensor ($\bm{\sigma}$) of the previous step, compute the external forces on the solids from the fluid, i.e., $\bm{F}^{ext}(t)$ from equation (51) 3. 3.

Apply the immersed boundary method (IBM) as presented in section 2.3 4. 4.

Impose the immersed boundary force ($\bm{f}_{imm}$) and any other forcing term to the fluid through the Shan-Chen forcing scheme as presented in algorithm 2 5. 5.

Collide & Stream, lattice Boltzmann (LB) step as presented in section 2.2, thus advance the fluid from $t$ to $t+1$ 6. 6.

Using $\bm{F}^{ext}(t)$ solve the dynamics of the immersed bodies (blood cells) as shown in algorithm 1, thus advance the solids from $t$ to $t+1$

The steps above advance the whole system from $t$ to $t+1$. These steps are repeated for $m$ times such that $m\cdot\Delta t=$ physical simulated time needed. As stated by MacMeccan and colleagues [5], the FE solver time step can be generally longer than the lattice time step. As a result, step 6 can be called every $n_{npFEM}$ steps. However in the current work, we set $n_{npFEM}=1$ to avoid complicating the validation of our framework. Furthermore, in the case of multiple bodies we have to consider collisions to avoid interpenetrations. The method is described in section 4.7 and this step is realized together with step 2.

4 Results

In this section, we present a series of in silico experiments and compare them with their in vitro counterparts. A complete validation is thus performed, examining thoroughly all the constituents of the proposed framework. More in details, we validate the proposed RBC model disregarding any external fluid, by testing a RBC under stretching and recovery in sections 4.1 and 4.2. Following the validated deformable body model, we test the fluid solver and the fluid-structure interaction (FSI), placing a single RBC in shear and Poiseuille flows in sections 4.4 and 4.5, while measuring suitable deformation indices. We then compare our framework against other state-of-the-art solvers in section 4.6, and finally in section 4.7 we demonstrate the capability of our framework to cope with multiple deformable bodies, which means, its ability to examine the collective transport of cells in the blood plasma for our future research endeavours.

The tuning of the deformable RBC model parameters presented in this section is minimal, since we are using well-known and tested materials to describe its mechanical and dynamical properties. Most importantly, the same set of parameters which is independent of the mesh resolution or regularity is used for all the in silico experiments.

Our framework is written in C++, it utilizes Palabos [30] for the IB-LBM parts and our own npFEM library for the solid part. The npFEM library is based on ShapeOp [51] which is an implementation of projective dynamics, freely available on the web. Currently, the complete framework is written for CPUs with extensive use of Message Passing Interface (MPI), which deals with the distribution of the load to all the available cores. In forthcoming work, we will present a hybrid framework where the IB-LB parts are executed on CPUs while the npFEM solver is fully deployed on the available general purpose Graphics Processing Units (gpGPUs). Extensive performance tests are to be expected with this work.

For an easier overview, we provide table 1 that consolidates the most critical parameters used in our numerical experiments.

4.1 Stretching experiment

The in vitro stretching of RBCs is performed using the optical tweezers technique (Dao et al. 2003/ 2006 [14, 15], Mills et al. 2004 [52] and Suresh et al. 2005 [53]), where small silica beads are attached on opposite sides of the membrane and are moved by a laser beam that traps and manipulates them. The resulting deformation leads to an extension along the longitudinal axis and compression along the transverse axis. To simulate the above stress and strain condition, we apply equal and opposite forces on selected vertices of the RBC membrane. The force is increased gradually up to $200~{}pN$ (pN: pico Newton). Force increments are applied only when the red blood cell reaches an equilibrium state. The viscous properties of the membrane cannot be tested in this experiment, since we are interested in the converged equilibrium state only. The measurements taken to assess the deformation in this kind of experiment are the axial and transverse diameters.

The simulation of the RBC includes all the energy terms presented in section 2.1, namely the ones representing area and volume conservation, as well as the bending and the material energies. There are not many free parameters to be tuned, since we adopt the values from the original model proposed by Skalak [8], which are $B=0.5~{}10^{-2}~{}dyn/cm=5~{}pN/\mu m$ and $C=5~{}dyn/cm=5000~{}pN/\mu m$. The $C$ parameter corresponds to the stiffness $k_{i}$ of the weight $w_{i}(=A_{i}k_{i})$ of the area preservation energy presented in 2.1.2. For the bending rigidity of the membrane we use an augmented value from the one presented in [19, 21], i.e., $k_{bending}=2\widetilde{\kappa}=4\kappa/\sqrt{3}\approx 1~{}pN\mu m$, where $\kappa=100\kappa_{B}T$. The weight for the volume preservation can take any arbitrary, sufficiently high value. In our simulations, we chose $k_{volume}=250$. The only free parameter to be tuned is $D$, which is the stiffness of the correction term in Skalak’s energy (see section 2.1.3). We found out that $D=35~{}pN/\mu m$ results in a close fitting to the experimental results of [53]. As shown in Figure 3, our model (npFEM - continuous line) presents a fit of better quality than other solvers, especially concerning the transverse diameter at low deformations, where most of the solvers fail. The modification of Skalak’s model adds the energy term $(\lambda_{i}-1)^{4}D/4$ in the energy density function $\Psi$. This term is mainly activated whenever the stretching ratios $\lambda_{i}$ diverge from the rest state ($\lambda_{i}=1$). Therefore, it improves the strain-hardening behavior, whereas the original model fails to fit the experimental data. The set of parameters defined here is used for any other case study presented in the results section and does not depend on the underlying mesh structure.

We represent our experimental results by a shaded region (Figure 3), which encompasses the output from all the tested cases with varying RBC mesh resolution, i.e., from 66 surface vertices up to 1026. The continuous line within the shaded region represents the case with $719$ surface vertices. It is obvious that our solver is characterized by strong mesh independence which is a direct consequence of the proper discretization of the continuous energies ( a general advantage of FEM-based solvers). Currently, all the state-of-the-art solvers [19, 21] are based on coarse-grained spectrin-link models which are essentially mass-viscoelastic spring systems. Their calibration is based on a method first introduced by Dao et al. 2006 [15] and later extended by Fedosov et al. 2010 [19], which relies on a regular two-dimensional sheet of springs, i.e., equilateral triangular elements. In this case, the physical parameters are linked to the mesh resolution. A disadvantage is that coarser meshes cannot maintain mesh regularity and therefore these methods fail to maintain mesh independence, demanding tedious tuning of the parameters set.

The time discretization is irrelevant in this kind of experiments since we are interested on the converged state of the stretched RBC. Using the backward Euler integration scheme, we guarantee unconditional stability of the simulation under any unfavorable situation, e.g., arbitrary large time step, large deformations, extreme mesh underresolution. Indeed, Figure 3 validates this statement and shows not only how stable is the npFEM solver (under-resolved case of 66 surface vertices), but also its accuracy since all the different simulations fall very close to the reference solution (experimental data).

As stated by Siguenza and colleagues [12], the axial and transverse diameters measured in optical tweezers experiments are found to be insufficient for validation purposes. They propose other quantities such as the height or the profile of the cell, which seem to be more sensitive to different membrane models. Nevertheless, they conclude and agree with the findings of Dimitrakopoulos et al. 2012 [11] that the only constitutive law able to properly match the wide variety of experimental data available in the literature is Skalak’s constitutive law, a variation of which we are using for our simulations. Additionally, we decided not to include other measured quantities as it is very difficult to find these deformation indices in the literature and compare with.

4.2 Recovery experiment

Stretching experiments help us define some of the most important mechanical properties of the RBC, but they do not reveal information about the viscoelastic nature of the membrane. There are two experiments that contribute towards this direction, namely the recovery and the tank-treading. The initial stage of the recovery experiment is the axial deformation of the RBC by optical tweezers [54], just like in the stretching experiment. When the deformation reaches the desired value, the forces are released and the time until the RBC reaches the initial state (rest shape) is tracked. To study the recovery experiment, Hochmuth et al. 1979 [55] introduced the time dependent elongation index

[TABLE]

where $\lambda=D_{A}/D_{T}$ ($D_{A}$ and $D_{T}$ are the axial and transverse diameters), $\lambda_{0}$ and $\lambda_{\infty}$ correspond to the ratios at time $t=0$ and $t=\infty$, $t_{c}$ is the characteristic recovery time, and $\delta$ is the exponent as suggested by [19]. As shown by Hochmuth et al. 1979 [55] the time rate of recovery is related mainly to the membrane viscosity and much less to the RBC cytoplasm or the outside aqueous environment. Given this observation, the numerical experiments are performed disregarding the existence of surrounding fluid.

As observed from Figure 4, the emulated viscoelastic model described in section 2.1.4 fits the curve of equation (52) very accurately, for $\delta=0$ and $t_{c}=150~{}msec$. As the stress levels increase, non-linear effects grow and a deviation from the theory is unavoidable, since it was developed for smaller deformation levels. Hochmuth suggested $t_{c}\approx 100-130~{}msec$, while Hénon et al. 1999 [54] suggested a characteristic time in the range of $100-300~{}msec$ where the lower bound corresponds to younger RBCs and the upper to older. Our tuning describes a RBC somewhere in the middle with $t_{c}=150~{}msec$, $\beta_{D}=0.01$ and $\kappa_{\text{damping}}=0.6$.

In the mass-spring solvers [19, 21], the viscosity of the RBC membrane is tackled through the addition of a dissipative force $\bm{F}_{ij}^{D}=-\eta\bm{v}_{ij}$, where $\eta$ is the membrane viscosity and $\bm{v}_{ij}$ is the relative velocity of the spring ends. However, this forcing term has the same disadvantage as the Rayleigh damping, i.e., it damps rigid body motions instead of only the relative ones. For example, a pure rotation produces $\bm{v}_{ij}\neq 0$, and thus internal forces are produced, even if they should not. This inherent drawback of $\bm{F}_{ij}^{D}$ has a severe impact on the dynamics of RBCs. A possible remedy would be to project the velocity difference $\bm{v}_{ij}$ onto the vector separating the masses.

4.3 Wheeler experiment

In this experiment, as performed by Yao et al. 2001 [56], RBCs are positioned in a shear flow with their axis of symmetry perpendicular to the flow velocity. The experiment, also known as wheeler experiment, takes its name from the characteristic wheel shape of the deformed RBCs. The red blood cells are submerged in a phosphate buffered solution (PBS) of low dynamic viscosity ($0.707cP$) and are submitted to a shear flow with shear rates in the range of $15-200~{}s^{-1}$. The deformation of the RBCs is quantified via the deformation index (DI) given by

[TABLE]

where $D_{0}$ is the initial RBC diameter ($\sim 7.82~{}\mu m$) and $D_{max}$ is the maximal diameter after deformation at constant shear rate. This is the first from a series of experiments that employs the fluid-structure interaction framework, which is composed of the lattice Boltzmann (LB) fluid solver, the nodal projective Finite Element (npFE) solid solver and the immersed boundary (IB) method for the coupling of the two solvers.

The RBC is placed at the center of the tank ($20\times 20\times 20~{}\mu m^{3}$) with periodic boundary conditions except for the walls. The wall velocity tunes the shear rate to the desired one. The spatial discretization is $\Delta x=0.5~{}\mu m$, the relaxation time $\tau\in[1,2]$, and the time step is defined through the diffusive scaling formula as

[TABLE]

where $\nu$ is the kinematic viscosity of the ambient fluid. A remark at this point is that in lattice Boltzmann simulations, if one fixes the spatial discretization and the relaxation time for a given fluid (known viscosity), then there is no freedom over the selection of the time discretization, as it is dictated by the diffusive scaling formula. If $\Delta t$ is not defined through equation (54), then physics are violated. However, in cases where the particle Reynolds number and the capillary number are very small, then a violation of $\Delta t$ can lead in an acceptable deviation of the simulated dimensionless numbers from the real ones. For a detailed overview on selecting model parameters for LBM simulations, we suggest chapter 8 from Krüger 2012 [57].

As seen from Figure 5, our framework (continuous line) fits very accurately the experimental data [56] (reference solution to investigate accuracy of the solver). The continuous line corresponds to a simulation of $514$ surface vertices, $\tau=2$ and $l=5$ iterations for the IBM algorithm. The shaded region around our results encompasses all the different case studies that we simulate with varying RBC mesh resolution (from 66 to 1026 surface vertices), different values of $\tau\in[1,2]$ and different iteration steps for the IBM iterative algorithm $l\in[1,5]$ as presented is section 2.3. This region is wider than the one of Figure 3 where we do not have fluid-solid interaction, nonetheless the property of mesh independence is validated once more, since all the different experiments fall very close to the reference solution. The accuracy of the simulations seems more sensitive when the fully coupled framework is deployed. Essentially, the immersed boundary method introduces stricter requirements on the resolution of the different solvers and how the resolutions match with each other. Nevertheless, for this single-cell case study even the under-resolved case with 66 surface vertices seems to be very close to the reference solution by Yao. At this point, we have to mention the work by Krüger 2012 [57] where the author suggests that the ratio $\bar{l}/\Delta x$, where $\bar{l}$ is the average edge length of the membrane mesh, should be somewhere between $0.5$ and $1.5$. The RBC membrane with $258$ vertices satisfies the upper bound of the above requirement and clearly the cases with $66$ or $130$ vertices violate this criterion. In spite of this violation, the framework accurately reproduces the reference solution. On the other hand, for cases of higher hematocrit ($>30\%$) we expect to face accuracy issues and probable instabilities. More in details, IBM demands $\bar{l}/\Delta x\sim 1$ to avoid fluid penetrating the solid surfaces and also mesh under-resolution leads to poorly resolved collisions between bodies. This topic will be thoroughly investigated in upcoming work.

The proposed solver presents also great robustness, since it withstands shear rates as high as $2000~{}s^{-1}$ without any stability issues (see Figure 5). Once more, Skalak’s original formulation presents a correct behavior for small deformations but fails for higher shear rates. The modification term, which is significant for larger deformations, repairs the RBC behavior within the whole deformation range. It is important to point out that a modification like the one presented by our framework cannot be easily realized by the coarse-grained spectrin link solvers, because the spring parameters are not intuitively linked to the material properties. On the contrary, the npFEM solver allows great material expressivity with minimum effort from the side of the user.

4.4 Tank-Treading experiment

Experimental observations of RBCs in shear flows [58, 59, 60, 61] reveal a tumbling (T) motion at low shear rates and tank-treading (TT) at higher shear rates. In this experiment, the RBC’s axis of symmetry is oriented along the flow velocity and its behavior is examined for different shear rates in the tank-treading regime. The ambient fluid dynamic viscosity ($5~{}cP$) is much larger than in the wheeler experiment. We are following the same experimental setup as Fedosov et al. 2010 [19] and Tran-Son-Tay et al. 1984 [60]. The RBC (258 vertices mesh resolution) is positioned at the center of a domain $20\times 20\times 20~{}\mu m^{3}$, the spatial resolution is $0.5~{}\mu m$ and the relaxation time $\tau=1$. The IBM iterations can vary from one to five without any difference on the resulting deformation.

Tran-Son-Tay et al. 1984 [60] studied this setup and defined two curves of the TT frequency (TTF) for young and old RBCs. In the present study, we are averaging the two curves and the “error” bounds correspond to the maximum and minimum values. The TTF is calculated by tracking a bead attached on the RBC surface and counting the cycles per second (see Figure 7). Figure 6 shows an excellent agreement of our framework compared to the experimental data of [60]. Figure 7 shows that the tank treading motion is very precisely captured, while the continuous curve is qualitatively the same as compared to Fischer et al. 1978 [58]. The TT motion at higher shear rates results in a “breathing” of the RBC, i.e., in an alternating compression-extension in the longitudinal axis as shown in Figure 7.

4.5 Poiseuille flow experiment

The ability of RBCs to travel through microcapillaries much smaller than their size is paramount for the transfer of oxygen to the organs. Thereby, the ability of a solver to simulate confined flows is a great advantage. Several researchers did experimental studies on the deformation of RBCs in small tubes and channels, e.g., Tomaiuolo et al. 2011 [62], Tsuakada et al. 2001 [63], and Suzuki et al. 1996 [64]. The main feature of this flow is the transition of the RBC biconcave shape to a parachute-like shape, which intensifies as the velocity increases. A common way to quantify the deformation of the RBC in the parachute-like regime is via a deformation index defined as $\Gamma=L/D$, where $L$ is the length of the parachute (length of the dome) and $D$ is the deformed diameter, with an approximate initial value of $7.82~{}\mu m$. Tsukada et al. 2001 [63] provide a graph where they correlate $\Gamma$ with the RBC velocity. Generally for tube diameters close to the RBC diameter, $\Gamma$ adopts values larger than 1. In our numerical setup, the tube diameter is $9.2~{}\mu m$, the RBC velocity varies from $4~{}\text{to}~{}9~{}mm/s$, and for all cases we achieved $\Gamma>1$. In Figure 8, we present the transition from the biconcave to the parachute-like shape. The resulted shape is very close to the experimentally observed one by Tomaiuolo et al. 2011 [62]. Simulations in confined spaces demand finer resolutions and smaller time steps to maintain stability, and our spatial discretization consquently is chosen to be $\Delta x=1/3~{}\mu m$ with a relaxation time $\tau=1$. These stability issues are not to be confused with the unconditional stability of the npFEM solver, since they refer to the FSI coupling. The explicit coupling of the fluid solver with the solid npFEM solver introduces some limitations at extreme flow cases like the investigated one. With explicit coupling, we mean that the external forces on the RBC come from the positions at time $t$ and not at $t+1$, see equation (2) for more details.

The parachute-like shape is a challenging deformation to capture. Most of the solvers [21, 19] present a parachute with $\Gamma\ll 1$, where the dome is very small. Our solver on the contrary captures the expected shape even if the employed bending rigidity is higher than the suggested one.

4.6 Computational efficiency

Main focus of our work is the solid solver, and thus the computational efficiency concerns the npFEM solver with respect to other well-known implementations. Currently, the state-of-the-art open-source solver for the simulation of the collective transport of cells in blood plasma is the software Hemocell [21, 22]. The solver uses Palabos [30] for the blood plasma, the RBCs are represented as discrete element membranes (mass-viscoelastic spring systems based on the model introduced by [19]), and the FSI is realized by a modified IBM introduced in Mountrakis 2015 [1]. Their framework has a proven capability of simulating multiple deformable bodies, thanks in part to the computationally efficient solid structure solver. Moreover, other state-of-the-art high-performing solvers, e.g., Blumers et al. 2017 [65], Rossinelli et al. 2015 [66], use almost the same methods for the solid solver, i.e., coarse grained spectrin-link models, and thus we focus on Hemocell which is available to the public and has many similarities on the fluid & FSI parts with our computational framework. Additionally, we avoid comparisons with GPU-accelerated solid solvers (see Blumers et al. 2017 [65]) since it would not be fair to compare a CPU with a GPU implementation. In an upcoming project, we will present a GPU version of the npFEM solver that achieves execution time per iteration at the order of nanoseconds in contrast to the current efficiency which is at the order of milliseconds.

The reference simulation is chosen to be the wheeler experiment, presented in section 4.3 and both solvers are executed on the same workstation. As expected, comparing an explicit mass-viscoelastic spring system with an implicit FEM-based solver, the efficiency comparison favors the former, as seen in table 2. Nevertheless, and in spite of the higher complexity of the approach, the npFEM framework appears to exhibit a computational efficiency quite close to Hemocell and consequently to any other mass-spring system solver.

It is also worth noting the benchmark of the linear explicit FE solver of MacMeccan et al. 2009 [5] for which the average computational time for an update is $2.3~{}msec/step$ for a RBC with $254$ surface vertices. A follow-up work from the same group, Reasor et al. 2011 [20], stated that the aforementioned solver was abandoned because of accuracy issues, since it could not capture tank-treading and the parachute-like shape in confined flows.

4.7 Proof of Capability

Recapitulating the performance tests of our solver, one can observe that the extension to flows of multiple RBCs should not be a problem. Indeed, with a computational time per update very close to the one of mass-spring systems, and with proven capability of the latter to scale well with increasing number of deformable bodies [67], we present Poiseuille and shear flows with O($10^{2}$) RBCs. The deliberately small number of RBCs helps us present the rich dynamics of individual RBCs produced by our framework. Studies with larger numbers of RBCs and higher hematocrit will be part of upcoming work. Furthermore, dealing with vessels of few micrometers, the hematocrit is by definition low. As stated by Zhao & Shaqfeh 2011 [68], although the normal range of human body hematocrit is 40%-50%, the local hematocrit drops significantly in small vessels and is only about 40% of the body average at $30~{}\mu m$ vessel diameter. Hence the mean hematocrit in this study is chosen to be around 20%.

To run a case with multiple RBCs in a robust manner, we need on top of our framework a collision detection/ handling system. Regarding the collision detection, there are two stages. In the first, Palabos [30] uses an efficient data structure for identifying, for each vertex of the examined RBC, neighboring vertices on other RBCs. The algorithm is based on the idea that all RBC vertices are attached to a closely located fluid cell, allowing to efficiently locate them on ground of their spatial position. At the second stage, npFEM solver builds a local kd-tree [69] for an even faster collision detection system, reducing the colliding data received from Palabos. The collision handling fits naturally in our solver in the sense that whenever a collision is detected, a collision energy/ force is added to the non-PD part of our solver (see equations (9) and (10)) acting upon the vertex under collision. Every colliding point has a desired position ($\bm{p}_{i}$) which is defined from the colliding neighbors and their normals. Thus the collision energy and force are defined as

[TABLE]

where $w_{col}$ is the weight of the collision term in the optimization problem. An important precaution is that at every iteration of the solver, the algorithm checks if $\bm{x}_{i}$ still penetrates the neighboring body. If it does not, the corresponding collision term is cancelled. This procedure avoids sticking of the bodies, which would be unrealistic. The collision energies and forces have no contribution to the Hessian approximation matrix ($\widetilde{\mathbf{M}}/h^{2}+\mathbf{L}$) and thus they may increase the convergence steps. However, after extensive testing, we found out that because of the small time steps (limitation introduced by the explicit nature of LBM), the collision handling process has no measurable impact on the total computational time. While the suggested collision energies favor simplicity over accuracy, any other type of energy can replace our simple spring-like equation. The efficient handling of collisions is one of the striking advantages of our method, as other frameworks such as Hemocell [21] skip collision detection, and rely on a sufficient fluid resolution to avoid collisions thanks to lubrication effects.

Both Poiseuille and shear flows are executed on a regular lattice of size $32\times 32\times 80~{}\mu m^{3}$ with periodic boundary conditions except from the walls. The spatial discretization for the blood plasma is $0.5~{}\mu m$, the LBM relaxation time $\tau=1$, the dynamic viscosity of the blood plasma $\sim 1.2~{}cP$, and the time discretization is computed from equation (54).

The Poiseuille flow as seen in Figure 9 consists of 123 RBCs (258 surface vertices each) inside a tube of $15.75~{}\mu m$ radius at a hematocrit of $\sim 20\%$. The shear flow as seen in Figure 10 consists of 123 RBCs (258 surface vertices each) at a hematocrit of $\sim 14\%$. The shear rate of the simulation is $500~{}s^{-1}$.

The npFEM solver has the capability to be executed independently of the fluid solver, using a larger time step, and thus speeding up substantially the overall simulation time. This useful capability is backed by the increased stability of our deformable bodies solver given its implicit nature. Additionally, the implicit Euler introduces a small numerical dissipation which is beneficial for the coupling. The upper limit on the solid solver time step is indirectly imposed by the collision resolution. Large time steps of the npFEM solver may lead to a potential under-resolution of the collisions detection, and thus interpenetrations may be observed, although the overall stability of the collective transport is not affected.

5 Conclusion

We presented a novel, highly accurate and robust model for the description of the mechanical, rheological and dynamical properties of red blood cells. The model is characterized by generality, robustness, accuracy and performance, fulfilling all the criteria we consider to be critical aspects for a universal model. To our knowledge, the npFEM solver is a novel approach in the field of computational biophysics to simulate viscoelastic bodies. The solver is independent of mesh resolution or mesh regularity and just one set of parameters can accurately describe the behavior of the body for our full range of case studies. Its extension from a thin-shell model to fully 3D bodies with internal structure and to any type of cells (platelets, white blood cells) is trivial based on the principles discussed at section 2.1. The performance of the new continuum-based solver is very close to other state-of-the-art mesoscopic solvers, and thus we managed to bridge the computational gap which was the main obstacle for using the former models.

The framework that combines the npFEM solver, the fluid solver and the FSI is versatile, highly modular and reliable, allowing the user to extend it based on other preferred numerical methods, as well.

As the number of the coupled models is increased, it is not uncommon that the quality of the output decreases. In our case however, extensive validation for many single RBC test cases showed a close match with the experimental data for multiple scenarios.

In the future, based on the proof of capability (see section 4.7), we plan to validate the framework for multiple RBCs test cases and further investigate its scaling with an increasing number of various deformable bodies. Of special interest would be to investigate reduced basis finite elements techniques, that would lead to an even faster deformable bodies solver. Finally, our npFEM solver with its modular interface on evaluating energies and forces could allow an efficient porting of the code to a general-purpose GPU platform.

Funding

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 675451 (CompBioMed project).

Acknowledgements

We acknowledge support from the PASC project 2017-2020: Virtual Physiological Blood: an HPC framework for blood flow simulations in vasculature and in medical devices.

Appendix A Alternating solver

The minimization of the objective function (7) with strictly quadratic potential energies is performed using an alternating (local/ global) solver. In the local step, $\bm{x}$ is fixed while the projections $\bm{p}$ per element are computed. In the global step, $\bm{p}$ is fixed and the minimum of $g$ is found by equating its gradient to zero:

[TABLE]

so

[TABLE]

The matrix $\widetilde{\mathbf{M}}/h^{2}+\mathbf{L}$ is symmetric positive definite and thus $\bm{x}_{n+1}^{*}$ is a global minimum for fixed projections. The above scheme is guaranteed to converge monotonically to a minimum. The local step decreases each constraint keeping fixed the vector $\bm{x}$. The global step, keeping $\bm{p}$ fixed, finds a compromise of the new state towards the projected positions and thus minimizes the objective function. Consequently, this sequence is non-increasing and bounded from below (minimization of Frobenius norms). The local/ global steps are repeated until convergence is satisfied (typically below 10 iterations).

This alternating solver is a special case of the quasi-Newton method presented in section 2.1. Using the same Hessian approximation $\widetilde{\mathbf{H}}=\widetilde{\mathbf{M}}/h^{2}+\mathbf{L}$, it is revealed that the descent direction $\bm{d}$:

[TABLE]

Thus, projective dynamics can be viewed as a quasi-Newton method that computes the next iterate as $\bm{x}+\bm{d}$. For normal quasi-Newton methods, we have to use line search techniques to find $\alpha$ such that $\bm{x}+\alpha\bm{d}$ reduces the objective function as much as possible [43]. The particularity of PD is that $\alpha=1$ [25] and thus there is no need for line search techniques. At this point, we must highlight that the choice of $\widetilde{\mathbf{H}}$ is inspired by the local/ global solver, as seen from equation (58), and thus it is deeply rooted in projective dynamics.

Appendix B FEM reminder

When an object deforms, every material point $\bm{X}$ (undeformed state) is displaced to a new deformed location $\bm{x}$. The deformation function $\bm{\phi}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ maps every material point to its respective deformed location, such that $\bm{x}=\bm{\phi}(\bm{X})$. The Jacobian matrix of the deformation function $\partial\bm{\phi}/\partial\bm{X}$ is an important physical quantity for describing elastic bodies and is called deformation gradient $\mathbf{F}$. Let us focus on a triangular element isometrically embedded in $\mathbb{R}^{2}$. The deformation function can be defined to be a piecewise linear function $\hat{\bm{\phi}}$ (approximation of $\bm{\phi}$) over the element. In more details, for the triangular elements we have [32]

[TABLE]

where $\mathbf{A}$ & $\bm{b}$ are element specifics. Differentiating equation (60) with respect to $\bm{X}$, we conclude that the deformation gradient $\mathbf{F}=\partial\hat{\phi}/\partial\bm{X}=\mathbf{A}$ is constant for the element. Let us denote with $\vec{X}_{1},\vec{X}_{2},\vec{X}_{3}$ the undeformed vertex positions of an arbitrary triangular element, and let $\vec{x}_{1},\vec{x}_{2},\vec{x}_{3}$ be the deformed counterparts. Each vertex must satisfy the deformation function

[TABLE]

a more compact form reveals that

[TABLE]

Consequently, the deformation gradient $\mathbf{F}\in\mathbb{R}^{2\times 2}$ of an arbitrary triangular element embedded in $\mathbb{R}^{2}$ is

[TABLE]

Of special interest is the derivation of various strain tensors from the deformation gradient. The Green strain tensor is given by

[TABLE]

while for small deformations, the infinitesimal strain tensor is given by

[TABLE]