Numerical simulations of elastic capsules with nucleus in shear flow

TL;DR

This paper presents a numerical method to simulate the deformation of elastic capsules with a nucleus in shear flow, revealing how nuclear properties affect cell deformation and flexibility.

Contribution

The study introduces a spectral method for simulating elastic capsules with a nucleus, analyzing the impact of various parameters on cell deformation in shear flow.

Findings

Eukaryotic cells deform less than prokaryotic cells under shear.

Deformation decreases with higher membrane bending stiffness.

Cells with a nucleus are less flexible passing through narrow capillaries.

Abstract

The shear-induced deformation of a capsule with a stiff nucleus, a model of eukaryotic cells, is studied numerically. The membrane of the cell and of its nucleus are modelled as a thin and impermeable elastic material obeying a Neo-Hookean constitutive law. The membranes are discretised by a Lagrangian mesh and their governing equations are solved in spectral space using spherical harmonics, while the fluid equations are solved on a staggered grid using a second-order finite differences scheme. The fluid-structure coupling is obtained using an immersed boundary method. The numerical approach is presented and validated for the case of a single capsule in a shear flow. The variations induced by the presence of the nucleus on the cell deformation are investigated when varying the viscosity ratio between the inner and outer fluids, the membrane elasticity and its bending stiffness. The deformation of the eukaryotic cell is smaller than that of the prokaryotic one. The reduction in deformation increases for larger values of the capillary number. The eukaryotic cell remains thicker in its middle part compared to the prokaryotic one, thus making it less flexible to pass through narrow capillaries. For a viscosity ratio of 5, the deformation of the cell is smaller than in the case of uniform viscosity. In addition, for non-zero bending stiffness of the membrane, the deformation decreases and the shape is closer to an ellipsoid. Finally, we compare the results obtained modeling the nucleus as an inner stiffer membrane with those obtained using a rigid particle.

Equations

$$\nabla \cdot \text{u}=0,$$

$$\frac{\partial \text{u}}{\partial t}+\text{u}\cdot \nabla \text{u}=-\nabla P+\frac{1}{Re}\nabla \cdot \left[{\mu }^{\ast }(\nabla \text{u}+\nabla {\text{u}}^T)\right]+\text{f},$$

$$Re=\frac{\rho \dot{\gamma }{R}^2}{{\mu }_o}.$$

$${\text{a}}_1=\frac{\partial \text{x}}{\partial \theta },{\text{a}}_2=\frac{\partial \text{x}}{\partial \varphi },{\text{a}}_3=\frac{{\text{a}}_1\times {\text{a}}_2}{\left|{\text{a}}_1\times {\text{a}}_2\right|},$$

$${\text{a}}_{\alpha \beta }={\text{a}}_{\alpha }\cdot {\text{a}}_{\beta },{\text{a}}^{\alpha \beta }={\text{a}}^{\alpha }\cdot {\text{a}}^{\beta },$$

$$I_1={\text{A}}^{\alpha \beta }{\text{a}}_{\alpha \beta }-2,I_2=\left|{\text{A}}^{\alpha \beta }\right|\left|{\text{a}}_{\alpha \beta }\right|-1.$$

$$I_1={\lambda }_1^2+{\lambda }_2^2-2,$$

$$I_2={\lambda }_1^2{\lambda }_2^2-1=J_s^2-1.$$

$$\text{T}=\frac{\text{1}}{{\text{J}}_s}\text{F}\cdot \frac{\partial {\text{W}}_s}{\partial \text{e}}\cdot {\text{F}}^T$$

$$T^{\alpha \beta }=\frac{2}{J_s}\frac{\partial W_s}{\partial I_1}{\text{A}}^{\alpha \beta }+2J_s\frac{\partial W_s}{\partial I_2}{\text{a}}^{\alpha \beta }.$$

$$W_s^{NH}=\frac{1}{2We}\left(I_1-1+\frac{1}{I_2+1}\right),$$

$${\nabla }_s\cdot \text{T}+\text{q}=0,$$

$$\frac{\partial T^{\alpha \beta }}{\partial {\xi }^{\alpha }}+{\Gamma }_{\alpha \lambda }^{\alpha }{T}^{\lambda \beta }+{\Gamma }_{\alpha \lambda }^{\beta }{T}^{\alpha \lambda }+q^{\beta }=0,\beta =1,2$$

$$T^{\alpha \beta }{b}_{\alpha \beta }+q^n=0$$

$${\text{M}}_{\beta }^{\alpha }=-B\left(b_{\beta }^{\alpha }-B_{\beta }^{\alpha }\right),$$

$$M_{|\alpha }^{\alpha \beta }-Q^{\beta }=0,$$

$${\text{boldmath}\epsilon }_{\alpha \beta }\left(T^{\alpha \beta }-b_{\gamma }^{\alpha }{M}^{\gamma \beta }\right)=0,$$

$$\frac{\partial T^{\alpha \beta }}{\partial {\xi }^{\alpha }}+{\Gamma }_{\alpha \lambda }^{\alpha }{T}^{\lambda \beta }+{\Gamma }_{\alpha \lambda }^{\beta }{T}^{\alpha \lambda }-b_{\alpha }^{\beta }{Q}^{\alpha }+q^{\beta }=0,\beta =1,2$$

$$T^{\alpha \beta }{b}_{\alpha \beta }-Q_{|\alpha }^{\alpha }+q^n=0.$$

$$\frac{{\mathbf{u}}^{\ast }-{\mathbf{u}}^n}{\Delta t}=\mathbf{RHS}(\mathbf{u},p),$$

$${\nabla }^{2}p=-\frac{1}{\Delta t}\nabla \cdot {\mathbf{u}}^{\ast }.$$

$${\mathbf{u}}^{n+1}={\mathbf{u}}^{\ast }+\Delta t\nabla p,$$

$$(\mathbf{I}-\Delta t\mathbf{D}){\mathbf{u}}^{n+1}={\mathbf{RHS}}^n$$

$${\mu }^{\ast }(\mathbf{x})$$

$$\mathbf{Du}=\underset{{\mathbf{D}}_1\mathbf{u}}{\underbrace{\frac{1}{Re}{\nabla }^2\mathbf{u}}}+\stackrel{{\mathbf{D}}_2\mathbf{u}}{\overbrace{\frac{1}{Re}\nabla \cdot \left({\mu }^{\prime }(\mathbf{x})\left[\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right]\right)}}.$$

$$(\mathbf{I}-\Delta t{\mathbf{D}}_1){\mathbf{u}}^{n+1}={\mathbf{RHS}}^n+\Delta t{\mathbf{D}}_2{\mathbf{u}}^n.$$

$$x(\theta ,\varphi )=\sum _{n=0}^{N_{SH}-1}\sum _{m=0}^n{\bar{P}}_n^m(\cos \theta )(a_{nm}\cos m\varphi +b_{nm}\sin m\varphi ),$$

$${\bar{P}}_n^m(x)=\frac{1}{2^{n}n!}\sqrt{\frac{(2n+1)(n-m)!}{2(n+m)!}}(1-x^2{)}^{m/2}\frac{d^{n+m}}{d{x}^{n+m}}(x^2-1{)}^n.$$

$${\mu }^{\ast }(\mathbf{x})=(1-I(\mathbf{x}))+I(\mathbf{x})\lambda .$$

$${\nabla }^{2}I=\nabla \cdot \text{G}$$

$${\text{U}}_{ib}(X,t)={\int }_{\Omega }\text{u}(x,t)\delta (\text{X}-\text{x})dx,$$

$$\text{f}(\text{x},t)=\int \text{q}(\text{X},t)\delta (\text{x-X})ds.$$

$${\text{X}}^{n+1}={\text{X}}^n+{\int }_0^t{\text{U}}_{ib}dt.$$

Full text

1 Introduction

The deformation of a cell in shear flows is one of the fundamental mechanical problems in cell biology. A living cell is subjected to mechanical forces of various magnitude, direction and distribution throughout its life. The cell response to those forces reflects its biological function [35]. As an example, red blood cells (RBC) have a diameter of about $7\mu m$ in the underformed state. Their ability to deform quite significantly allow them to pass through narrow capillaries having a diameter of $3\mu m$ [55]. This high deformability enables them to reach various parts of the human body and to distribute oxygen and nutrient to cells. RBC can however be affected by protozoan Plamodium falciparum, the parasites that cause malaria. These parasites change the red blood cell chemical and structural composition [3, 9], thus inducing a stiffening of the cell membrane [56, 22]. Such variations in the mechanical properties of the cells affects the blood rheology which may help to diagnose diseases.

Many early experimental studies have addressed the interaction between tiny deformable particles and an external flow. Several interesting types of motion have been discovered such as tumbling and tank-treading in shear flow [19, 11], the zipper flow pattern [15] or parachute cell shapes [48]. More recent studies focused on cells that exhibit very large deformations at high shear rates, which can cause breaking [7], just to mention few examples. Most of these studies are of experimental nature. Such investigations can however be quite expensive since they require dedicated facilities not easy to fabricate. In addition, experimentally measuring the exact deformation and stresses can be rather complicated. Developing robust and reliable numerical platforms is thus of increasing importance in order to perform high-fidelity simulations beside laboratory experiments.

Many cells, including red blood cells, can be modeled as capsules. Capsules consist of a droplet enclosed by a thin membrane: the membrane area can vary while the enclosed volume is constant. Nowadays, several numerical studies on the deformation of a capsule in shear flow have been reported in literature. At certain shear rates, the capsule reaches a steady shape while its membrane exhibits a rotation known as tank-treading motion [23]. This tank-treading motion disappears when the viscosity or shear rate of the external fluid becomes low enough and instead a flipping or tumbling motion similar to that of a rigid body appears [46, 12]. Membranes can also undergo buckling or folding for high elastic moduli or at low and high shear rates in absence of bending rigidity [53, 23]. A solution to this problem is proposed by introducing a stress on undeformed membrane, the so-called pre-stressed capsule [30]. As regards the motion of non-spherical capsules in shear flow, different types of motion occur when changing the fluid viscosity, the membrane elasticity, the geometry of the problem or the applied shear rate. In [49], a phase diagram is presented for biconcave shaped capsule in which the transition from tank-treading to tumbling motion is identified when decreasing the shear rate.

For eukaryotic cells, the overall mechanical properties of a cell are not only determined by its membrane but also by other cell organelle such as the cell nucleus [43]. Typically, the nucleus is stiffer than the surrounding cytoplasm which results in lower deformation when subject to the external stimuli [6, 21]. To model and predict the cell behavior, the mechanical properties of the nucleus need to be quantified. To this end, both experimental tests and numerical simulations have been carried out in [6]. The elastic modulus of the nucleus in round and spread cells was found to be around $5000N/m^{2}$, roughly ten times larger than for the cytoplasm. As further example, the nucleus of bovine cells is nine times stiffer than the cytoplasm [36], yet small deformations of the nucleus may occur when a cell is subjected to flow [16]. Though it can exhibit large deformation on a substrate when highly compressed, stretched or flattened [20, 5, 24], the nucleus may be assumed as a rigid particle for an intermediate range of the applied forces (external shear).

The objective of present research is to quantify the deformation of a cell with nucleus under simple shear flow at low but finite Reynolds numbers (fluid inertia) using numerical simulations. For the present model, the cell is made of two thin (two-dimensional) elastic membranes as originally proposed by [25]: the outer membrane separating the cell from the ambient fluid and an inner membrane acting as the boundary of the nucleus. The neo-Hookean hyperelastic model is chosen for the strain energy of the membranes, while the inner and outer fluids are assumed to be Newtonian and can have different viscosities. The viscosity inside nucleus is assumed to be the same as the surrounding cytoplasm. As the deformation of the nucleus is assumed to be negligible, this is numerically treated either with a second inner stiffer membrane or, alternatively, as a rigid particle using a different numerical approach [31, 2]. The presence of a nucleus and the effect of viscosity ratio and bending stiffness on the cell deformation are investigated to document the potential of the numerical method here developed. The results are expected to provide an interesting comparison between mechanical properties of cells with and without nucleus.

This manuscript is organised as follows. In §2, the geometry of the problem and the governing equations for the membranes and flow dynamics are presented. Section §3 provides a brief introduction to the numerical methods used in order to simulate the problem. Validations of the present implementations are presented in §4 while section §5 reports the main results of the present study. Finally, conclusions and perspectives are given in section §6.

2 Problem statement

Figure 1 depicts the flow configuration considered and the coordinate system adopted. An initially spherical cell located at the geometrical center of a rectangular computation box is considered. The upper and lower walls of the domain move at opposite velocities in the streamwise direction while periodicity is assumed in the other two directions. As boundary conditions we therefore impose no-slip at the walls and periodicity in the steamwise and spanwise directions .

2.1 Navier-Stokes equations

The dynamics of the incompressible flow of a Newtonian fluid are governed by the Navier-Stokes equations,

[TABLE]

where $\textbf{u}=(u,v,w)^{T}$ is the velocity vector, $P$ the hydrodynamic pressure, and ${\bf f}$ indicates the fluid-solid interaction force. The Reynolds number $Re$ is defined as

[TABLE]

In the expression above, $\dot{\gamma}$ is shear rate, R the reference cell radius, $\rho$ the reference density (assumed here to be the same for the fluid inside and outside the cell) and $\mu_{o}$ the reference viscosity. In the present work, the reference viscosity is set to be the viscosity of the fluid outside the cell and the ratio between inner and outer viscosity is defined as $\lambda=\frac{\mu_{i}}{\mu_{o}}$, with $\mu_{i}$ the inner viscosity. In the expression above $\mu^{\ast}=\frac{\mu(\bf{x})}{\mu_{o}}$ indicates the ratio of viscosity at each point to the reference viscosity.

2.2 Membrane dynamics

Cells are surrounded by a deformable membrane known as the plasma or cytoplasmic membrane. Along with a number of biological functions, its main purpose is to separate the interior of each cell from the external environment. It can moreover deform quite significantly, as in the case of red blood cells traveling through capillary vessels. In the present numerical study, such membrane is modelled using a hyper-elastic model.

A point on the surface of the cell is expressed by using the curvilinear coordinates $(\xi^{1},\xi^{2})$. To define the cell, two different coordinate bases are used, see figure 2. The first is a fixed cartesian base, $\left({\bf e_{1},e_{2},e_{3}}\right)$ corresponding to position $\textbf{x}\left(\xi^{1},\xi^{2}\right)$. The second coordinate is a local covariant base $({\bf a}_{1},{\bf a}_{2},{\bf a}_{3})$ which follows the local deformation of the membrane. The unit vectors of the local base are

[TABLE]

where $\theta$ and $\phi$ are latitudinal and longitudinal angles on the cell surface. The co-variant and contra-variant metric tensors are defined as

[TABLE]

where $\alpha,\beta=1,2$. The basis vectors and metric tensors in the undeformed (reference) state are hereafter denoted by capital letters $(\textbf{A}^{\alpha},\ \textbf{A}^{\alpha\beta})$. The invariants of the transformation $I_{1}$ and $I_{2}$ are defined as

[TABLE]

Equivalently, they can also be determined from the principal stretching ratios $\lambda_{1}$ and $\lambda_{2}$ as

[TABLE]

The ratio of the deformed to the undeformed surface area is defined by the Jacobian $J_{s}=\lambda_{1}\lambda_{2}$. The two dimensional Cauchy stress tensor, ${\bf T}$, is computed from the strain energy function per unit area $W_{s}(I_{1},I_{2})$ of the undeformed membrane as

[TABLE]

where $\textbf{F}=\textbf{a}_{\alpha}\bigotimes\textbf{A}^{\alpha}$ and $\textbf{e}=(\textbf{F}^{T}\cdot\ \textbf{F}-\textbf{I})/2$ is the Green-Lagrange strain tensor. Equation (7) can be further expressed component-wise as

[TABLE]

In the rest of this work, the strain energy function $W_{s}$ is modeled using the neo-Hookean (NH) law [40, 41, 58]. Using this model, the strain energy function is expressed as

[TABLE]

where $We=\frac{\rho R^{3}\dot{\gamma}^{2}}{G_{s}}$ is the Weber number (or non-dimensional surface shear modulus). The local equilibrium relates the tensor T to the external elastic load q according to

[TABLE]

with $\nabla_{s}\cdot$ the surface divergence operator. In curvilinear coordinates, the load vector can be written as ${\bf q}={q}^{\beta}\textbf{a}_{\beta}+{q}^{n}\textbf{n}$. The force balance in equation (10) is further decomposed into tangential and normal components,

[TABLE]

where $\Gamma_{\alpha\lambda}^{\alpha}$ and $\Gamma^{\beta}_{\alpha\lambda}$ are the Christoffel symbols.

In some cases, due to noninfinitesimal membrane thickness or a preferred configuration of an interfacial molecular network, bending moments accompanied by transverse shear tensions play an important role on cell deformation [39]. Bending stiffness is incorporated into the model using a linear isotropic model for the bending moment [40, 39]:

[TABLE]

where $B=\frac{G_{b}}{\rho R^{5}\dot{\gamma}^{2}}$ is the non dimensional bending modulus, and $b^{\alpha}_{\beta}$ is the Gaussian curvature ($B^{\alpha}_{\beta}$ corresponds to that of the reference configuration). According to the local torque balance, including bending moments on the membrane, we obtain the transverse shear vector Q and in-plane stress tensor T,

[TABLE]

where ${}^{,}\ _{|\alpha}\ {}^{,}$ represents the covariant derivative and $\varepsilon$ is the two-dimensional Levi-Civita tensor. The left hand side of equation (13) identifies the antisymmetric part of the in-plane stress tensor, which is always zero as proved in [57]. Including the transverse shear stress Q, the local stress equilibrium, including bending finally gives

[TABLE]

The non-dimensional numbers in equations (3-14) are obtained using the radius of the cell $R$ as length scale, the shear rate $1/\dot{\gamma}$ as time scale, thus coupling membrane deformation and flow dynamics, $\rho_{o}\dot{\gamma}^{2}R^{3}$ as reference surface shear modulus and $\rho_{o}\dot{\gamma}^{2}R^{5}$ for the bending stiffness. The Weber number can be re-written as $We=Re\cdot Ca$ where $Ca=\mu_{o}R\dot{\gamma}/G_{s}$ is the Capillary number. In the present study, we shall consider different stiffnesses of the cell membrane by varying the Capillary number and keeping the Reynolds number constant.

3 Numerical methods

3.1 The Navier-Stokes solver

The Navier-Stokes equations are discretized using a staggered uniform grid to prevent checkerboard numerical instability, while the time integration relies on the classical projection method [8]. This method is a three-step procedure: first, a non-solenoidal velocity field ${\bf u}^{*}$ is computed as

[TABLE]

where ${\bf RHS}({\bf u},p)$ is the right-hand side of the discretised Navier-Stokes equations and contains the fluid-structure interaction (FSI) forces. In a second step, the pressure field is obtained as the solution to the following Poisson equation

[TABLE]

Finally, the corrected velocity field ${\bf u}^{n+1}$ is obtained as

[TABLE]

where $\nabla p$ is the pressure gradient required for the velocity field ${\bf u}^{n+1}$ to be divergence-free. Second-order central differences are used for the spatial discretization of the convective terms, while their temporal integration relies on the Adams-Bashforth explicit method.

As we allow for viscosity contrast between the fluid inside and outside the cells, the classical Fast Fourier spectral method cannot be readily used to evaluate the diffusive term ${\bf Du}=\nabla\cdot\left(\mu\left[\nabla{\bf u}+\nabla{\bf u}^{T}\right]\right)$. Indeed, the viscosity field being a function of space, this operator cannot be reduced to a constant coefficient Laplace operator. However, Dodd & Ferrante [10] have recently introduced a splitting operator technique able to overcome this drawback. Though it has initially been derived for the pressure Poisson equation, this splitting approach can easily be extended to the Helmholtz equation resulting from an implicit (or semi-implicit) integration of the diffusive terms,

[TABLE]

where ${\bf I}$ is the identify matrix, and ${\bf RHS}^{n}$ the discretized right-hand side including the non-linear advection terms. Given the viscosity field

[TABLE]

where $1$ is the constant part and $\mu^{\prime}({\bf x})$ the space-varying component, the diffusive term ${\bf Du}$ can be re-written as

[TABLE]

The constant coefficients operator ${\bf D}_{1}$ can then be treated implicitly while the variable coefficients operator ${\bf D}_{2}$ is treated explicitly. The resulting Helmholtz equation then reads

[TABLE]

Since ${\bf D}_{1}$ is now a constant coefficient Laplace operator, equation (21) can be solved using a classical Helmholtz solver based on Fast Fourier transforms. Note that a similar Fast Fourier-based solver is used to solve the Poisson equation (16) for the pressure.

3.2 Membrane representation: Spherical harmonics

The membrane shape has been modeled as linear piece-wise functions on triangular meshes by Pozrikidis [38], Ramanujan & Pozrikidis [42] and Li & Sarkar [34] among others. The finite element method has also been employed by Walter et al. [54] for its generality and versatility. Another interesting method is the global spectral method. Fourier spectral interpolation and spherical harmonics have been used for two-dimensional [14] and three-dimensional simulations [26, 57]. Here, we follow the approach of Zhao et al. [57] , previously implemented in [58, 60, 59, 45]. This is briefly outlined below.

The capsule surface is mapped onto the surface of the unit reference sphere $S^{2}$, using the angles in spherical coordinates $(\theta,\phi)$ for the parametrization. The parameter space $\{(\theta,\phi)\mid 0\leq\theta\leq\pi,0\leq\phi\leq 2\pi\}$ is discretized by a quadrilateral grid using Gauss-Legendre (GL) quadrature intervals in $\theta$ and uniform spacing in the $\phi$ direction. All other surface quantities are stored on the same mesh, i.e. the grid is co-located. The surface coordinates $x(\theta,\phi)$ are expressed by a truncated series of spherical harmonic functions,

[TABLE]

yielding $N^{2}_{SH}$ spherical harmonic modes. The corresponding normalized Legendre polynomials are given by

[TABLE]

The SPHEREPACK library [1, 50] is employed for the forward and backward transformations. To deal with aliasing errors arising due to the nonlinearities in the membrane equations (products, roots and inverse operations needed to calculate the geometric quantities), we use an approximate de-aliasing by performing the nonlinear operations on $M_{SH}>N_{SH}$ points and filtering the result back to $N_{SH}$ points. A detailed discussion on this issue is provided in Freund & Zhao [13]. For most of the results presented in the present study, 576 spherical harmonics with 24 modes are used to define the cell shape.

Considereing different viscosity inside and outside the cell, a space and time dependent viscosity field is defined by an indicator function $I({\bf x},t)$ related to the membrane location,

[TABLE]

Here we follow Unverdi and Tryggvason [52] for the definition of the indicator function as the solution to the following Poisson equation

[TABLE]

where the Green’s function $\textbf{G}=\int\delta(X-x)\textbf{n}\ \mathrm{d}s$, and n is the unit normal vector to the cell surface. Using the smooth Dirac delta function introduced below in the computation of G makes the indicator function smoother near the boundary [28]. Such indicator function is similar to the regularised Heaviside function used in the levelset framework.

3.3 Immersed boundary method

3.3.1 Immersed boundary method for deformable partices

The immersed boundary method [37] is commonly used to solve fluid-solid interaction problems. In this method, two distinct sets of grid points are used: (i) an Eulerian mesh to solve fluid flow, see section 3.1, and (ii) a Lagrangian mesh for solving the particle motion, section 3.2. In the present low Reynolds number framework, we start from the original approach of Peskin [37]. At each time step, the fluid velocity defined on the Eulerian mesh is first interpolated onto the Lagrangian mesh,

[TABLE]

where $x$ and $X$ are the Eulerian and Lagrangian coordinates and $\delta$ is a smooth Dirac delta function, here that proposed in Ref. [44]. The elastic force per area q and surface normal vectors ${\bf n}$ are then computed from the membrane equations described above. As next step, the normal vectors are used to compute the indicator function $I({\bf x},t)$ on the Eulerian mesh. The force is then spread to Eulerian mesh and added to the momentum equations as

[TABLE]

Thereafter, the positions of the Lagrangian points are updated according to

[TABLE]

Note that equation (28) assumes an over-damped regime, i.e. the Lagrangian points go to their equilibrium position immediately after the FSI force is applied. Finally, the fluid flow is solved in the Eulerian framework as explained in section 3.1. A flowchart for computational procedure at each iteration is depicted in figure 3 .

The method described above is not particularly efficient when the Reynolds number increases since it requires very small time steps, which increases the computational time. At moderate and high Reynolds numbers, a modification of the method by Kim et al. [28], is employed to be consistent with the assumption of inertialess membrane. In our approach, in addition to the Lagrangian coordinates $\bf{X}$, we introduce the additional immersed boundary points $\bf{X}_{ib}$ whose motion is governed by equation (28). Since the total force exerted on each element on the membrane surface is equal to the difference between its acceleration and the acceleration of fluid element at the same location, the motion of the real Lagrangian points is governed by

[TABLE]

where $\rho_{os}$ is the surface density of the base fluid. The two sets of Lagrangian points, X and $\textbf{X}_{ib}$, are connected to each other by a spring and damper, i.e. a fluid-solid interaction force $\textbf{F}_{FSI}$ computed using the following feedback law

[TABLE]

The final modified procedure is therefore as follows. At each time step, we first compute $\textbf{U}_{ib}$ and the fluid-solid interaction force $\textbf{F}_{FSI}$ from equation (30). The indicator function $I({\bf x},t)$ is then computed to identify the interior of the cell and impose viscosity contrasts, and the momentum equation solved to obtain the flow velocity ${\bf u}$. Finally, the positions of the Lagrangian points $\bf{X}$ are updated using equation (29). This additional equation is made non-dimensional as above, in particular with $\rho_{o}R$ for the surface density of the base fluid and $\rho_{o}R^{2}\dot{\gamma}^{2}$ for the elastic and fluid-solid interaction forces per unit area. For completeness, we report the non-dimensional form of equation (29),

[TABLE]

where $d^{\ast}=d/R$ is ratio between the membrane thickness and initial radius of the cell, assumed in the present study to be $d^{\ast}=0.01$. In the above, $\textbf{F}_{A}$ is the penalty force used to enforce volume conservation, calculated as in [28]:

[TABLE]

Here $\Delta p$ represents the pressure generated by the volume change, $\eta(\theta,\phi)$ is the surface area of each element and $\textbf{e}_{n}$ the local unit normal vector. This force is also added to the elastic force q before spreading it to the Eulerian mesh according to equation (27).

3.3.2 Immersed boundary method for rigid particles

In the present study, two different immersed boundary methods are considered to simulate the stiff nucleus. First, the method described above is used with high surface shear modulus for the inner capsule. In the second approach we follow the implementation by Breugem [4], which has been widely used in the framework of rigid particles, see e.g. [31]. In this method, a moving Lagrangian mesh is adopted to impose no-slip and no-penetration on the surface of a rigid object. The numerical procedure is as follows: first, the prediction velocity $\textbf{u}^{\ast}$ is computed from the Navier-Stokes equations neglecting the fluid-solid interaction force. This fluid velocity is then interpolated onto the Lagrangian mesh $(\textbf{U}^{\ast})$ and the fluid-solid interaction force computed, based on the difference between the fluid and the solid body velocity at each Lagrangian point,

[TABLE]

This force is spread to the Eulerian grid and the second prediction velocity $\textbf{u}^{\ast\ast}$ obtained by solving the Navier-Stokes equations with the fluid-solid interaction force. The divergence-free constraint is then imposed on the velocity field by solving the pressure Poisson equation and correcting the velocity field appropriately. Finally, the total force and torque on each particle is computed, and the translational and rotational velocities of the particle obtained by integrating the Newton-Euler equations. Readers are also referred to [51] for further details.

3.4 Notes on the parallelisation and implementation

The Eulerian mesh is decomposed using a 2D-pencil domain decomposition in the streamwise and spanwise direction. For that purpose, the library 2DECOMP & FFT [33] is used. The same library handles all of the transpose operations required for the Helmholtz and Poisson solvers based on the fast Fourier transforms. Regarding the parallelization of each particle/capsule, each processor can either be master or slave. A processor is labelled master if it contains most of the Lagrangian points describing the given particle, while those containing only some of these points are labelled as slaves. The rest of the processors do not have any rule for the considered particle. Only the master processor has the information of the particle (e.g. Lagrangian points and their velocities) in its memory, though the slaves can access it for interpolation and spreading operations, which might require information from the neighbours. For the particle equations, the master is responsible for all the numerical procedures, e.g. transformation using spherical harmonics. Such parallelization saves memory usage but requires communication between cores at each time step.

In order to obtain accurate results, the density of Eulerian and Lagrangian grid points should be similar. In some cases, it is nonetheless necessary to have a very fine Eulerian mesh, thus requiring very fine Lagrangian mesh. However, since the spherical harmonic calculations are costly, the overall computational time increases significantly. To make the code more efficient, two different sets of Lagrangian points are therefore considered: forcing points that are used for interpolation-spreading and the spherical-harmonic points that are used to define the shape of the cell and the elastic stresses. While the density of the forcing points has to be similar to that of the Eulerian points, fewer points are required for the spherical harmonics representation of the cell, especially in the case of stiffer membranes deforming less. At each time step, before computing the elastic forces, the spherical-harmonic points are obtained using spectral interpolation. These points are then used to compute the elastic forces and surface normal vectors. Once computed, the elastic forces and surface normal vectors are interpolated onto the finer mesh so that the ealstic forces are spread on the Eulerian mesh. All spectral interpolations are done with the SPHEREPACK library [1, 50].

4 Validations

In order to validate the current implementation, two different cases are considered for a simple capsule: deformation of a capsule in shear flow at low Reynolds number and equilibrium position of a capsule in Poiseuille flow at finite Reynolds number. The domain is a box with the size $10\times 10\times 10$ in units of cell radius, and the cell is located at its center. The Eulerian grid is $128^{3}$ whereas for the Lagrangian mesh $24\times 48$ points are chosen in the latitudinal and longitudinal directions respectively. The non-dimensional bending stiffness $\frac{G_{B}}{\rho_{o}R^{5}\dot{\gamma}^{2}}$ is zero unless otherwise mentioned. The dealiasing ratio is kept at $M_{SH}/N_{SH}=2$.

A grid independence study has been carried out for $Ca=0.6$ by increasing the number of the Eulerian grid points and of spherical-harmonic grids by a factor of 1.5. In this case, the change in the deformation parameter is found to be less than $2\%$. In addition, we also increased the box size by a factor 1.5 and measured a difference in the deformation parameter below $0.3\%$. To check time step independency of the results, we decreased the time step by $50\%$ and obtain a change in the deformation parameter of only $0.002\%$. Finally, the effect of the ratio between the forcing points and the number of spherical harmonics describing the membrane deformation is investigated using 1152 force points. When the ratio between the two is $4$, the difference measured by the deformation parameter is less than $0.5\%$ with respect to the case with ratio $1$. Increasing the ratio from $1$ to $2.25$, the change is of about $0.1\%$ . Here we used the more demanding ratio of 1 to be sure to capture all features of the cell shape for large deformations.

4.1 Single cell in shear flow

In this section, the deformation of a cell in a simple shear flow obtained with the present implementation is compared to the results by Pranay et al. [41] and Zhu et al. [58]. These authors have used boundary integral approaches to solve for the Stokes flow, so the Reynolds number is chosen here small enough, $Re=0.1$, to ensure low inertia of the flow. The cell deformation is quantified by the parameter

[TABLE]

where $l_{1}$ and $l_{2}$ denote the major and minor semi-axes of the equivalent ellipsoid in the middle plane, respectively. The inertia tensor of the cell is used to obtain the equivalent ellipsoid as in [42]. As shown in figure 4, a good agreement is obtained with the different implementations.

4.2 Capsule in channel flow at finite inertia

The second validation case is provided by the equilibrium position of a cell in Poiseuille flow at finite Reynolds numbers. Indeed, after some transients, a single capsule reaches a steady state characterised by a constant wall-normal position, constant velocity and deformation. This equilibrium position is located between the wall and the channel centerline and depends on both the Reynolds and Capillary numbers [27]. In the present validation, the Capillary number is set to $Ca=0.174$ and the Reynolds number is varied. Note that since the present bending model is different from that in the work of Kilimnik et al [27], we chose the value of the bending stiffness that best fit their data, B=.02. The dependence of the wall-normal position and deformation of the capsule on the Reynolds number are reported on figures 5(a) and 5(b), respectively. As shown by these figures, a good agreement with the results in [27] is obtained. In figure 5(b), the difference appearing at $Re=100$ may be attributed to the model considered in [27]. These authors used a Hookean law with a lattice-spring model (LSM) where the membrane has a finite thickness whereas in present work we use the hyperelastic neo-Hookian model with a infinitely thin membrane. At small Reynolds numbers, when the deformation is relatively small, the last term on the right hand side of equation (9) vanishes and our model practically reduces to a Hookean law: the two methods give thus very similar results. In the case of large deformations, the term introducing non-linearity in the constitutive model is not negligible any longer and the results start to deviate from each other.

5 Results

5.1 Stiff nucleus

As a starting point, the two-membrane model, based on the original IBM by Peskin [37], is employed. In order to have a stiff nucleus, the capillary number of the nucleus is chosen to be 300 times smaller than that of the outer membrane. The cell is subject to homogeneous shear as in the validation cases presented above, i.e. same configuration and numerical parameters.

Figure 6 depicts the deformation parameter as a function of the Capillary number for viscosity ratios $\lambda=1$ and 5 in the absence of bending resistance, both for capsules with and without a stiff nucleus. Note that the deformation parameter is computed on the outer membrane, the inner one not being noticeably deformed. As shown in this figure, the presence of a nucleus reduces the deformation, and this reduction is larger the higher the Capillary number. The stiff nucleus reduces the outer membrane deformation since the minimum radius cannot be smaller than the radius of the nucleus. At higher Capillary numbers, the membrane would tend to deform more thus making the effect of the nucleus becomes more evident. It can be inferred from figure 6 that the deformation is smaller for the cases with a more viscous fluid between the outer membrane and the nucleus, $\lambda=5$. In this case, viscous forces appear to work together with elastic forces to reduce the cell deformation.

The transient evolution of the deformation parameter to reach the final state is demonstrated in figure 7 for three different capillary numbers. The figure shows that larger capillary numbers require longer time to reach the final steady state. As for the steady state, the deformation is larger for higher capillary numbers.

The first two rows of figure 8(a) depict the steady shape of the cell for three different capillary numbers of the outer membrane and zero beding stiffness. Note that, in the first row, the cell considered has no nucleus. In the absence of nucleus, the cell assumes an ellipsoidal shape, while it has a thicker middle part in the presence of the nucleus. Cells with nucleus thus have a lower flexibility and may encounter more difficulties to pass through narrow vessels.

The effect of bending stiffness on the deformation parameter is presented in figure 6(b). It can be observed that cells with bending stiffness deform less and the reduction measured by the deformation parameter increases with the Capillary number. This effect is documented by the shape of capsules with bending stiffness reported in the lowest panels of figure 8. Here, one can see that the deformation is reduced mainly on the edges of the capsule. Indeed, figure 8 shows that the cell shape is closer to an ellipsoid when adding bending rigidity. The difference between the different cases shown in the figures demonstrates that the effect of the bending rigidity is not negligible in such conditions and should be accounted for to obtain more accurate predictions.

For a number of microfluidic applications, it may be important to understand the effect of flow inertia on the deformation of the transported cells. We therefore report in figure 9 the effect of increasing the Reynolds number on the deformation parameter. To prevent buckling, we considered small bending stiffness in the simulations. It can be observed that when increasing the Reynolds number the steady state deformation parameter first decreases for $Re=1$ and then increases ($Re=5$). The initial deformation rate is faster when increasing inertia. Note also some oscillations in the deformation for $Re=5$, as observed in previous studies. These can be attributed to the formation of a pair of vortices inside the cell, on the two sides of the nucleus, in the transient stage (See figure 9b). Such vortices disappear at steady state but their formation and breakup results in oscillations of the cell membrane.

5.2 Rigid nucleus

As mentioned previously, two different models of nucleus are considered. For the results in the present section, the nucleus is modeled as a rigid, not deformable, particle following a different implementation of the immersed boundary method, see above. Figure 10 reports a comparison of the deformation for cells whose nucleus is modelled as a stiff membrane or as a rigid spherical particle, for viscosity ratio $\lambda=1$. The two methods produce approximately similar results, however with slightly larger deformation for cells with a rigid spherical nucleus. This fact can be related to the different nucleus rotation rates obtained with the two models. Indeed, for a rigid nucleus we assume no slip at the interface, whereas some slip is present at the surface of the nucleus if this is represented by an elastic membrane.

Finally, we report that the computational time needed for the case of a capsule with rigid nucleus is about 1.2 times that for a cell with a nucleus represented by a stiff membrane. Note that the simulations assuming a rigid nucleus have been performed by coupling together two different approaches, able to model deformable and rigid object. This implementation opens possibility of modeling new more complicated structures, which will be investigated in the future.

6 Conclusions

The deformation of a capsule containing a stiff nucleus in homogeneous shear flow is studied numerically using an immersed boundary method to account for the fluid-solid interaction. The Neo-Hookean hyperelastic constitutive model is used to describe the cell membrane deformation while the fluid inside and outside each capsule is assumed to be Newtonian. The cell nucleus is modeled in two different ways, first as a second inner capsule with a significantly stiffer membrane and then as spherical rigid particle using a different implementation of the immersed boundary method, most suited to solid objects [32].

In the immersed boundary method, a Lagrangian mesh is used to follow the deformation of the elastic membrane defining the capsule and an Eulerian mesh to solve the momentum equations. The shape of the membrane, its deformation and internal stresses are represented by means of spherical harmonics in order to have an accurate computation of the high-order derivatives of the membrane geometry. To save computational time in cases with very fine underlying Eulerian meshes, we have implemented the possibility of using a coarser Lagrangian mesh for the computation of the cell shape and stresses and a finer mesh for the communication of forces exchanged with the fluid. Spectral interpolation is employed to link the two Lagrangian representations of the geometry of the cell. Finally, to have the possibility to consider different viscosities inside and outside the cell, an indicator function is computed on the Eulerian mesh as the solution of a Poisson equation. The right hand side is obtained by spreading the normal vectors to cell surface known at the Lagrangian grid points onto the Eulerian mesh. The viscosity is then taken to be a function of this indicator function. The accuracy of the code is validated against results pertaining to the deformation of capsules without nucleus. In particular, we consider the inertialess flow of a capsule in homogeneous shear and transport in pressure-driven Poiseuille flow at moderate Reynolds numbers (finite inertia).

The behavior of cells with a stiff nucleus in homogeneous shear has then been investigated. The cell deformation parameter is reported for different Capillary numbers, two values of the viscosity ratio and with or without bending rigidity. We observe that the deformation is smaller for cells with a nucleus. Examining the shape of the cell, that with nucleus is thicker in the middle part making it less flexible to pass through narrow vessels. When also considering bending stiffness, we observe an even smaller deformation and the shape of the cell is more regular and closer to an ellipsoid. Finally, we have compared the results obtained by modeling the nucleus as a rigid particle, reporting small differences. We show also that the numerical approaches for rigid and deformable objects can coexist, which opens the possibility of modelling more complicated structures, e.g. small rigid cavities and obstacles in the flow [60].

The method presented here can be employed and extended in the future to study the behavior of cells in different and more complicated configurations, enabling us to extract qualitative and quantitative data about the maximum stress on the membrane. Possible extensions include the possibility to consider a vesicle [47], i.e. an inextensible membrane, and a density contrast. Including cell-cell and cell-wall interactions in the numerical platform would allow us to study pair interactions of cells with nucleus in shear flow [41] and ultimately investigate dense suspensions of deformable objects [17, 18, 29].

Acknowledgements

This work was supported by the European Research Council Grant No. ERC-2013-CoG-616186, TRITOS and by the Swedish Research Council (VR). The authors acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing) and the support from the COST Action MP1305: Flowing matter. Dr. Lailai Zhu is thanked for the help with the implementation of the membrane module.