J. theor. Biol. (1992) 158, 173-193

Leukocyte Deformability: Finite Element Modeling of Large Viscoelastic Deformation CHENG

DONGt AND

RICHARD SKALAK

Department of Applied Mechanics and Engineering Sciences, Bioengineering, R-012, University of California, San Diego, La Jolla, CA 92093, U.S.A. (Receioed on 7 May 1991, Accepted on 6 February 1992) An axisymmetric deformation of a viscoelastic sphere bounded by a prestressed elastic thin shell in response to external pressure is studied by a finite element method. The research is motivated by the need for understanding the passive behavior of human leukocytes (white blood cells) and interpreting extensive experimental data in terms of the mechanical properties. The cell at rest is modeled as a sphere consisting of a cortical prestressed shell with incompressible Maxwell fluid interior. A largestrain deformation theory is developed based on the proposed model. General nonlinear, large strain constitutive relations for the cortical shell are derived by neglecting the bending stiffness. A representation of the constitutive equations in the form of an integral of strain history for the incompressible Maxwell interior is used in the formulation of numerical scheme. A finite element program is developed, in which a sliding boundary condition is imposed on all contact surfaces. The mathematical model developed is applied to evaluate experimental data of pipette tests and observations of blood flow. Introduction It is the object of this study to develop quantitative theoretical and computational analyses of rheological behavior of human leukocytes in a passive (non-activated) state. A method of computing finite deformation of an axisymmetric shell-like body containing incompressible viscoelastic material under external stress is developed as a model of the leukocyte. This study will attempt to understand how a leukocyte deforms in response to the pressure applied during entry into small capillary vessels. Although the present theories and computations with finite element formulation are motivated by our interest in the problem of modeling leukocyte deformability, the general method of analysis should be applicable to other problems in cell biology. Leukocytes are usually classified as granulocytes, lymphocytes and monocytes, in which granulocytes include neutrophils, eosinophils and basophils. All of these cells are morphologically distinct and circulate in the blood. Leukocytes constitute a very small volumetric fraction of the total formed elements of the normal blood. Leukocytes occupy 1/600th of the total cellular volume, which is very small comparing to

t Current address: BioengineeringProgram, Collegeof Engineering,The PennsylvaniaState University, 229 HallowellBuilding, UniversityPark, PA 16802-6804, U.S.A. 173

0022-5193/92/180173 + 23 $08.00/0

© 1992 AcademicPress Limited

174

C.

DONG

AND

R.

SKALAK

erythrocytes which occupy almost one half of the blood volume according to Caro et al. (1978) and Fung ( 1981 ). Leukocytes have generally been neglected in consideration of blood rheology or hemorheology in bulk because of their small volumetric percentage. However, it has been found that leukocytes are larger and much stiffer than erythrocytes. They are about twice as big (in volume per cell) and 2000 times more viscous than red blood cells (Schmid-Sch6nbein et al., 1980). Hence, the effect of their presence can be very important for blood flow in the narrow capillaries where the vessels and cells are comparable in diameter. There is a need, therefore, to study leukocyte deformability, especially large deformations of the cell during entry into small capillary vessels. In general, the leukocyte ultrastructure consists of the nucleus, granules, cytoplasm, cytoskeleton and membrane (Fig. 1). More than half of the cell volume is occupied by cytoplasm, the main cell continent between the nucleus and the cell membrane. Cytoplasm itself contains many different types of granules or organelles and forms a cortical layer (cytoskeleton) beneath the membrane (Alberts et al., 1983). The whole cell is bounded by a bilayer lipid membrane which separates the cell from its environment and controls transport of material into and out of the cell. There are many fine folds or wrinkles on the cell membrane surface which provide about 80-120% excess area over that needed to envelope the cell volume by a smooth spherical surface (Fig. 2). Micropipette techniques have been used to study a single cell's rheological behavior (Chien et al., 1984; Evans, 1984). We will focus on the deformability of non-activated neutrophils under pipette testing in the present paper. Experiments show that when aspiration pressure is applied, a cell near the pipette is drawn toward the pipette tip by fluid flow. At the moment the cell seals the pipette, an initial deformation is observed and then a continuous deformation of the cell body inside the pipette follows with time in a non-linear manner (Chien et al., 1984). The passive deformation of leukocytes under applied pressure in the micropipette experiments has been studied with several analytical models. These models have given a description of either the time course of entire cell deformation under the assumption of small strain (Schmid-Sch6nbein et al., 1981 ; Dong et al., 1988) and a simplified structural approximation of cell behavior (Evans, 1984; Yeung & Evans, 1989). It is useful to develop a more detailed biomechanical model with appropriate physical assumptions about the cell structure and its intracellular components, especially when cells are subjected to a finite deformation. A numerical formulation using the finite element method is needed with the aid of a super-computer to carry out useful computations. A finite element program is developed, in which the leukocyte model consists of a Maxwell fluid (Fig. 3) bounded by a prestressed cortical layer. This layer undergoes elastic stretching and shearing (Fig. 4) in large deformations. The interior Maxwell material is assumed to be incompressible and represents the heterogenous interior of the cell in an average sense. The prestressed elastic cortical shell is considered to represent an actin-rich layer near the cell surface. The effects of cell surface folds or wrinkles on cortical elasticity is taken into account in a general non-linear constitutive relation. The membrane theory is adopted for this cortical layer, neglecting the

FIG. 1. Transmission electron micrograph of human neutrophil in the resting state (courtesy of Dr K.-L. P. Sung).

C. DONG and R. SKALAK

(Fa~:ing p. 174)

J. theor. Biol.

FIG. 2. Scanning electron micrograph of human neutrophils showing pronounced ruffles and folds of the membrane while the cell is approximately spherical in shape (courtesy of Dr G. W. Schmid-SchSnbein).

LEUKOCYTE

DEFORMABILITY

175

FIG. 3. Schematic 1-D Maxwell fluid for the cell interior; k and p are elastic and viscous c o m p o n e n t s respectively. Constitutive equations are 3-D (see text).

layer thickness and bending stiffness. A contact sliding boundary condition has been imposed on the cortical layer surface where the cell contacts a rigid wall. In our present model, this rigid wall is used to model a micropipette in which the aspiration pressure is applied to the cell surface. The formulations are based on the principle of virtual work and satisfy the governing equations of the interior viscoelastic fluid and the exterior elastic thin shell, as well as all of the corresponding stress boundary conditions. Finite Element Formulation MAXWELL FLUID INTERIOR

Consider the leukocyte interior as an incompressible Maxwell fluid. Assume at some time t, the body has been deformed from the initial configuration x" into the current configuration zi, where x ~ (a = !, 2, 3) designate a set of initial curvilinear co-ordinates and zi (i= 1, 2, 3) define another set of curvilinear co-ordinates for the deformed configuration. Let ~'J be the contravariant components of Cauchy stress (or spatial stress) referred to the current state in the z~ co-ordinate system at time t. To obtain the stressdeformation relation applicable to the case of incompressible, isotropic, viscoelastic material, Pipkin (1964) introduced a general constitutive representation in the form of a series of multiple integrals of strain history, in which the first order approximation for a Maxwell fluid can be described by the following equation at time t:

Oz' Oz' .... b,,f 2k2 f ' expf_kt,~dt,~ rr~'(t)=-P(t)G'~+dx" Ox hg g 1 2 k ~ ' " " ( t ) - ~ - Jo Z,,,.(t-t') . Ia /.

(I)

176

C. D O N G

A N D R. S K A L A K

NI

I I

' N2

I !

! I

I I

I I

! I

,y \

\

FIG. 4. An axisymmetric thin shell with principal tensions N~ and N2 in the meridional and circumferential directions respectively.

where ~,,,, is the Green strain tensor whose covariant components in the x ° configuration are defined as (m, n = 1, 2, 3) :

r~.=] \dx"

dx"

g.,.

(2)

and p is arbitrary hydrostatic pressure in the Maxwell material; gob and G u are the contravariant components of the metric tensor referred to x a and z i co-ordinates, respectively (a, b, i,j= 1, 2, 3); g,,b are the covariant components o f the metric tensor referred to x a co-ordinates (a, b = 1, 2, 3). In eqn (1), the smoothness properties o f the Maxwell relaxation function has been applied for the time interval 0 ÷ < t < oo. Also, 7/°b(t) is assumed to be zero for t < 0 +. We define the reference (undeformed) state as a sphere since it is observed that the cell is roughly spherical before it is drawn into the micropipette in the aspiration experiment.

LEUKOCYTE

DEFORMABILITY

177

Let s ab be the contravariant components of Piola-Kirchhoff stress (or material stress) tensor referred to the initial x ° configuration, which are defined by: Sab

OXa OXb ij

=J-:=oz"oZj

(3)

where J is the Jacobian determinant. The value of J will be unity for incompressible material. The equilibrium equation in terms of the material stress for a Maxwell fluid at time t can be written as ( l = 1, 2, 3) :

(: Oxb/l oz31=0,

(4)

in which the vertical bar is used to indicate covariant differentiation with respect to x ~ co-ordinates. All gravitational and inertial forces are assumed to be negligible. A material form of the principle of virtual work which is equivalent to the equilibrium eqn (4) is given by (l= 1, 2, 3) :

fvo sabSY"b d Vo= f Aott°SutdAo,

(5)

where 6 Y.b, 6ut are variations of Green strain and displacement from an equilibrium state; Vo, Ao are the volume and the area of the body in the x ° configuration; t~0is the surface traction measured per unit initial area. It can be shown that eqn (5) also leads to satisfaction of the associated stress boundary conditions (1= 1,2, 3) : /o= s °b 0z t -noa, Oxb

(6)

in which n0o are the components of initial unit normal vector to area A0 in the x ° coordinate system. The Maxwell fluid representing the leukocyte interior can be modeled by an assembly of finite elements which are defined by nodal points. Introducing the nodal displacement components u,~ associated with the nodal point a, the displacement field u,,, in every finite element for the Maxwell fluid is approximated by the following interpolation at time t (m = 1, 2, 3): Ne

u,.(t)= E

(7)

tt=!

where N e is the total number of nodes in the finite element and ~a is the shape function associated with the node a (a = 1, 2 . . . . . N e) defined in the initial x ° coordinate system.

178

c. DONG AND R. SKALAK

In terms of nodal displacement u,~,,from eqn (7), the components of the Green strain tensor in eqn (2) can be written as (a, b, c, l = 1, 2, 3 ; a, fl = 1, 2 . . . . . N ~) :

=

,~,j

f

~

o

+do •,,,u~

i

a

["Its .m T M . n ~

IdllA c - - W

di,~q." I o n ~

UlU b

P ~ ,'lh~

,~harc a.,p "-'~ .,,'*" ~ h,,~~ b t .,t,,,~

--

a

[3

b

al

a

/3

dO do Fh,,g U,,Ut,] j.k

(8)

F~.. is the Christoffel symbol of the second kind associated with x" configuration and the comma in do~. denotes the partial derivative of shape function with respect to the co-ordinate x". From eqns (I-3), it is seen that the state of stress depends on the deformation history and the current isotropic pressure. In the case of incompressible viscoelastic material, the pressure is an unknown scalar function of location and time. For the finite element formulation in the present analysis, this pressure term is assumed to be equal to the element pressure pe, i.e. the pressure is considered to be uniform within each element e ( e = i, 2 . . . . . NP): p=pe,

(9)

where N p is the total number of finite elements used in representing the Maxwell fluid. Introducing all of these approximations, the Piola-Kirchhoff material stresses s "b can be expressed in terms of the nodal displacements u,~ and the element pressure pC. We apply the principle of virtual work [eqn (5)] to a typical finite element e, utilizing the eqns (7-8). Then the virtual work equations for the Maxwell fluid become (/, c, m, n = 1,2, 3) :

(

,agmb--

doa F /.~g,,,/

[O ..O ,h6,.-do ..do F,.b+do~doP,bF,"..--...~ . . . . .¢hju.~t,u

dV~

O

= f,4 to,,, ~, O " 6 u " " dA~, 6 a

(1o)

where V~, d~ are the initial volume and surface area of the typical finite element e in x" configuration: 6u .... is the variation of the contravariant components of displacement referred to the node a for the typical finite element e; 6',',, is the Kronecker symbol; a, fl = I, 2 . . . . . N".

LEUKOCYTE

DEFORMABILITY

179

Since the material is incompressible, the possible deformations of the Maxwell fluid will be restricted to the volume preserving deformation. The continuity equation for the typical finite element e in the Maxwell fluid can be expressed by (Oden, 1974): f

v~

(x~3-- l) d V ; = 0 .

(ll)

The invariant/3 can be defined in terms of Green strain components (/, m, n = 1, 2, 3): I3 = det 16,",+ 2y,,,/g~nl.

(12)

It will be seen that eqns (l 1-12) can be expressed in terms of finite element nodal displacement components if eqn (8) is employed. Then the eqns (10-12) constitute a system of governing equations for the Maxwell fluid, which are in the form of finite element discretization. The system of equations is written for each element and finally assembled into a single system. By performing the integrations over the initial volume V~, and initial surface area At of the finite element in these equations, an algebraic system of equations involving all nodal displacements and unknown pressures over the whole Maxwell fluid interior can be derived. The unknown displacements and pressures will depend on the surface tractions which are applied on the Maxwell fluid by the exterior cortical layer. Hence solving the coupled system of equations is necessary to complete the solution.

CORTICAL SIIELL

The leukocyte exterior is modeled as a prestressed cortical layer (thin shell) under elastic stretching and shearing during large deformation. The initial (undeformed) configuration is described in a system of curvilinear co-ordinates, x" (a = 1, 2, 3), and a second curvilinear co-ordinate system, z; (i= 1, 2, 3), is used to describe the current configuration at time t. An elastic potential (or membrane strain energy) IF,,, measured per unit initial area A0 of the cortical shell, is defined as the function of the principal stretch ratios ~.~ (n = 1,2): W,,, = W,,(A,1, A.2).

(13)

For the case of axisymmetric deformation, it is convenient to define the principal stretch ratios, At and ~-2, for the thin elastic shell in the meridional and circumferential directions, respectively (Fig. 4). A.,, ( n = l, 2) are defined by: 21

ds ds0

and

dr A.2=-7- , dr0

(14)

where So, s are lengths along the meridian in the initial x ~ and the current z i configurations, respectively; r0, r are radii of the revolution with respect to the axis of symmetry for the shell in the initial x ~ and the current z i configurations, respectively.

180

c.

DONG

AND

R. S K A L A K

In terms of the strain energy function W,,,, the material form of the principle of virtual work in eqn (5) can be expressed for the cortical shell as (1= 1, 2, 3) :

f

d.o=f /OSu,dgo,

(15)

in which 6 W,, is the variations of strain energy function; Ao is the initial surface area of the cortical shell in the x a co-ordinate system. The cortical layer tensions (or stress resultant), Ni (i= 1, 2) measured in dyn cm-' along the principal directions (Fig. 4) can be expressed in terms of the strain energy function (Dong, 1988): Nj-

1 OW,,, ;h 0~.1

and

N2-

1 014",,, ~10Z2

(16)

Now let displacement components ul in the cortical shell be interpolated by the finite element approximation at time t: ~e

ut(t) = E ~tVu~(t),

(17)

v=l

then a discrete representation of the principle of virtual work for the cortical shell in terms of the finite element nodal displacement uy (/= 1, 2, 3) is:

fj ~o{~.2Nl~-q-~.lN20~2vl~u~dAeo=~ Ito.~.~t~u~dA~,

(18)

in which ~0 is the shape function associated with the nodal point o (v = 1, 2 , . . . , ~e). N" is defined as the total number of nodes used in the typical finite element ~ (~= 1, 2 . . . . . /VP), where .~P is the total number of finite elements used in representing the cortical shell. For the present model of a leukocyte in a large aspiration deformation, the cortical shell is assumed to be prestressed, i.e. before any deformation there is a n isotropic tension To which persists and is additive to any stresses due to deformation. The strain energy function for the cortical shell will be based on the assumptions of an isotropic material. A stress-strain constitutive relation for the thin shell is written in terms of the principal cortical layer tensions and the principal stretch ratios. In particular, it is proposed that these tensions Ni (i= 1, 2) are comprised of both the isotropic part (due to prestressed tension and elastic stretching) and the anisotropic part (elastic shearing). Thus

N,= To + E,,A(A.,, A.2)+ EsB,(~.,, A,2),

(19)

where Ea and E, are, respectively, the elastic moduli for area dilatation and membrane shear; A and B; are functions of the principal stretch ratios which give the isotropic and anisotropic elastic tensions accordingly. Due to the effect of leukocyte surface unfolding on the cortical layer elasticity, we expect a non-linear behavior of the

LEUKOCYTE

181

DEFORMABILITY

prestressed elastic layer to be qualitatively as shown in Fig. 5 for a case in which the elastic layer is under isotropic tension N [assuming the stretch ratio Z~= Z~= Z in eqn (19)]. The critical point at which the cell membrane is unfolded is assumed to be at e=Z,~.2- 1 ~ 1.0 (i.e. a maximum 100% excess area can be provided). At this point membrane tensions suddenly increase in Fig. 5. The smooth portion of the membrane may have a large resistance to any further extension of area due to the tension increase. To express this effect based on the existence of strain energy function, a set of functions A and Bi in eqn (19) is suggested for the cortical shell in the present leukocyte model as follows (i = I, 2) : A(;t,, ~2) = (;L,;~2 - I) ~, (20) 1

2

1

in which p is a positive finite number, p > 1. In Fig. 5, p = 29 is chosen in the present cortical shell model to represent a steep change in tension after the cell membrane is unfolded. Equations (19-20) represent the unfolding characteristics of the neutrophil membrane in the experiments described by Schmid-Sch6nbein et al. (1980). /

,.0] I

/ / "o

/

/

I

I

I

/ jd J

.- ~

0 . .~

• ( --------)c

~o . ='~[

••

!

f

// wl

_j.

( .......

)p=9

( .......

)p:15

(

)p:29

// i

1.00

I 1.25

I

I 1.50

I

I 1-75

Stretch rotio ,~" ,~ : A z FiG. 5. Stress-strain and area strain curves in isotropic tension for the proposed leukocyte membrane model with parameters 7"o=0 • 12 dyn cm -t, Eo=0.01 dyn cm -~ and E,=O"14 dyn c m - L

The eqns (18-20) constitute a system of governing equations for the cortical shell, which can be written for each element and then assembled in a global matrix. During the finite element discretization, the principal tensions and stretch ratios are also expressed in terms of the finite element nodal displacement in the specific curvilinear co-ordinate system chosen. The equations finally derived are in the form of algebraic

182

C. DONG AND R. SKALAK

system involving unknown nodal displacements. The equations for the cortical shell must be coupled to those derived from the Maxwell fluid developed in the previous section via the common surface traction and continuity of displacement. Then assembling all of the individual contributions of each element yields a global system for the composite structure. The surface traction components {0 have been assumed to be specified along the initial surface of the cortical shell. However, the current values of surface tractions applied by the external pressure loading and the internal Maxwell stress depend on the complete solutions from the system of coupled equations governing the cortical shell and Maxwell fluid during large deformation. Therefore, the applied stress vector will be modified by using the updated stresses including the effects of the changes in the surface configuration.

Numerical Procedures

The finite element formulation presented in the preceding section leads to an assembly of finite elements representing both the Maxwell fluid and the cortical shell during deformation. The finite element matrices involved will be derived in this section. Since the deformation of the present leukocyte model is axisymmetric, global circular cylindrical co-ordinates (r, z, 0) are introduced for the finite element matrix equations. A typical finite element grid used for the numerical computation is sketched in Fig. 6 by using MAZE, an input generator for finite element mesh developed by Lawrence Livermore Laboratory (Hallquist, 1983). The elements used are rings in three-dimensional space, but the computation deals only with the typical meridional cross-section in which the elements are plane. Rectangular linear elements with four nodes are used for the Maxwell fluid interior while two node one-dimensional elements are used for the cortical shell.

FIG. 6. Finiteelementgrid used in numericalcomputationwith 150 elementsfor the interior Maxwell fluid and 20 elementsfor the exteriorcortical layer.

LEUKOCYTE MAXWELL

DEFORMABILITY ELEMENT

183

MATRICES

Let the element displacement vector be U e for the typical finite element e (e= 1, 2 . . . . . NP), such that

re=

{u;, u'_ . . . . .

u~ ~, u~°} T,

(21)

where uT, u~ are displacements of the nodal point a (a = 1, 2 . . . . , N ~) in the r and z directions, respectively. Starting with the virtual work form of equilibrium equation for the element e in the Maxwell body, eqn (I0) can be expressed in terms of physical components in the cylindrical co-ordinate system. For the axisymmetric case, a simple expression can be achieved by taking m, n = 1,2 in eqn (10), which gives ( a = 1 , 2 , . . . ,Ne):

y~ { rc~6u ~,+ 7r~ au~ } = 0,

(22)

Ct

where 7r7, Jr-~ are functions involving volume integrals and area integrals, which will be in terms of finite element shape functions ~ , nodal displacements u~, uP and element pressure pC. Let the vector H ~ be defined in the element e (e = 1, 2 . . . . . NP), such that /7 e= {TrJ, rr~ . . . . . trY', Jr_Y'}r

(23)

and it should be noted that the vector/7 e will be the function of the nodal displacement vector U e and the pressure pC. Moreover, the v e c t o r / 7 e also depends on the complete deformation history of the Maxwell material up to the current time t through the vector U~(t'), 0 < t' < t. Therefore, a system of simultaneous non-linear equations, involving 2N e unknown nodal displacement components in Ue(t) and one unknown uniform pressure p(t) for each element e at current time t, follows from eqn (22). Assuming that the variations of element nodal displacement components gut and cSu~ can be arbitrary gives: I V = [U~(t),pe(t), U~(t')] =0.

(24)

Equation (24) is a vector form that is equivalent to 2N ~ individual equations on its components, i.e. ~r~, 7 r ~ , . . . , trY', rr_-u'=0. In order to complete the system, the continuity eqn (11) needs to be added for each element e: A~[U~(0] =0,

(25)

where eqns (8), (12), and (21) have been employed. Considering Taylor expansions of eqns (24-25) with respect to the element displacement vector and pressure at a reference time to (0 < to < t) gives (no sum on e) : ll~(t) = H~(to) -~ aH~(t°~) A U~(t) + Oil~(t°) Ap~(t) = O, OU ~ Op~

A'(t) = A~(t0) + OA~(t°) AUg(t) = O, 8U ~

(26)

184

C. D O N G

AND

R. S K A L A K

where the terms of higher order in the residue in the expansions are neglected and A U ' ( t ) = U~(t) - U~(to); Ap~(t)=pe(t)-p~(to). Considering eqn (26), the element matrix for the Maxwell fluid in the present leukocyte model can be written at time go a s "

OU" OA~(to) OU ~

OH

Leukocyte deformability: finite element modeling of large viscoelastic deformation.

An axisymmetric deformation of a viscoelastic sphere bounded by a prestressed elastic thin shell in response to external pressure is studied by a fini...
3MB Sizes 0 Downloads 0 Views