Prosenjit Bagchi

Department of Mechanical & Aerospace

Engineering,

Rutgers,The State University of New Jersey,

98 Brett Road,

Piscataway,NJ 08854

Paul C.Johnson

Department of Bioengineering,

University of California,San Diego,

La Jolla,CA 92093

Aleksander S.Popel

Department of Biomedical Engineering,School of

Medicine,

Johns Hopkins University,

720 Rutland Avenue,611 Traylor Bldg.,

Baltimore,MD 21205

Computational Fluid Dynamic

Simulation of Aggregation of

Deformable Cells in a Shear Flow

We present computational ﬂuid dynamic (CFD) simulation of aggregation of two deform-

able cells in a shear ﬂow.This work is motivated by an attempt to develop computational

models of aggregation of red blood cells (RBCs).Aggregation of RBCs is a major deter-

minant of blood viscosity in microcirculation under physiological and pathological con-

ditions.Deformability of the RBCs plays a major role in determining their aggregability.

Deformability depends on the viscosity of the cytoplasmic ﬂuid and on the rigidity of the

cell membrane,in a macroscopic sense.This paper presents a computational study of

RBC aggregation that takes into account the rheology of the cells as well as cell-cell

adhesion kinetics.The simulation technique considered here is two dimensional and

based on the front tracking/immersed boundary method for multiple ﬂuids.Results pre-

sented here are on the dynamic events of aggregate formation between two cells,and its

subsequent motion,rolling,deformation,and breakage.We show that the rheological

properties of the cells have signiﬁcant effects on the dynamics of the aggregate.A stable

aggregate is formed at higher cytoplasmic viscosity and membrane rigidity.We also show

that the bonds formed between the cells change in a cyclic manner as the aggregate rolls

in a shear ﬂow.The cyclic behavior is related to the rolling orientation of the aggregate.

The frequency and amplitude of oscillation in the number of bonds also depend on the

rheological properties.DOI:10.1115/1.2112907

1 Introduction

Aggregation of cells plays a key role in many important bio-

logical processes.In the blood of humans and many other mam-

mals,red blood cells RBCs aggregate at low shear rates forming

a linear or fractal-like structure called “rouleaux.” During me-

tastasis,cancer cells often ﬂow through the blood stream as ag-

gregates.Platelets form aggregates to prevent blood clotting.The

present work is primarily motivated by an attempt to develop

computational models and simulations of aggregation of RBCs in

a shear ﬂow.RBCs under normal conditions are highly deform-

able.Under pathological conditions,such as sickle cell anemia

and bacterial infections,deformability of RBCs is greatly reduced.

Clinical studies have shown that under such conditions aggrega-

bility of RBCs is also increased 1,2.It is however not estab-

lished whether there is a connection between reduced deformabil-

ity and increased aggregability of RBCs.We hypothesize that

deformability of cells,and hence their rheological property,plays

a key role in determining their aggregability.To that end,we

develop a computational ﬂuid dynamic simulation to study aggre-

gation of deformable cells.The simulations considered here are

two dimensional2D,though the methodology can be extended to

3D.This is a simulation of aggregation of cells that takes into

account cell rheology.The methodology is general and can be

applied to any deformable cell,though the speciﬁc rheological

values used here correspond to RBCs.

The mechanism of RBC aggregation in blood ﬂow is as fol-

lows.RBCs are randomly drawn close to each other by the ﬂow of

plasma.If the shearing forces by ﬂuid motion are small,the cells

adhere to each other and form aggregates.Currently there are two

theories that describe the mechanism of aggregation:bridging be-

tween cells by cross-linking molecules 3,and osmotic force gen-

erated by the depletion of molecules in the intercellular space 4.

RBC aggregation leads to increased blood viscosity,and hence

elevated resistance to blood ﬂow.Aggregation is common in pa-

tients with peripheral vascular disease.Elevated aggregation is

often associated with a higher risk of cardiovascular disease.RBC

aggregation is elevated after myocardial infarction,ischemic brain

infarcts,in diabetes,and during sepsis.Thus,understanding the

mechanics of RBC aggregation may lead to a better understanding

of cardiovascular diseases and the pathology of blood.Studies on

the effect of RBC aggregation in microcirculation are summarized

in recent reviews 1,5.

Theoretical approaches to the study of aggregation can be

broadly classiﬁed into two categories.The ﬁrst one is the quasi-

steady approach which predicts the equilibrium shape of the

RBCs forming the aggregate through a balance between the adhe-

sive surface energy and the elastic energy stored in the RBC mem-

brane 6,7.The motion and breakage of an aggregate in a ﬂuid

ﬂow cannot be modeled in this approach.The second approach is

based on the theory of kinetic modeling of colloidal suspension

8–10.In the kinetic approach deformation of a cell is not con-

sidered.Further,the details of the ﬂow pattern around each de-

formable cell and aggregate are also neglected which can have a

signiﬁcant effect on the stability and motion of the aggregate.It

has been shown earlier by in vitro 11,and recently by in vivo

experiments that aggregation has a signiﬁcant effect on the veloc-

ity proﬁles in microvessels 12–15.

Clearly,a theoretical/computational model of cell aggregation

that would take into account deformability of the cells and their

dynamic motion in a ﬂow is lacking.Aggregation is a multiscale

process.At the nanoscale the molecular bridging between adjacent

cells initiates aggregation.At the microscale,deformability of in-

dividual cells dictates the process.At the macroscale,the ﬂow of

plasma imparts a shear force causing the rolling and breakage of

the aggregates.Here we present a computational ﬂuid dynamic

model that combines the three scales of the problem.Speciﬁcally

we study the effect of cell deformation,strength of the adhesion

molecule,and the shear rate of the bulk ﬂow on the motion and

Contributed by the Bioengineering Division of ASME for publication in the J

OUR-

NAL OF

B

IOMECHANICAL

E

NGINEERING

.Manuscript received November 12,2004;revi-

sion received July 27,2005.Review conducted by:C.Ross Ethier.Technical Editor:

Frank Yin.Approved August 15,2005.

1070/Vol.127,DECEMBER 2005 Copyright © 2005 by ASME Transactions of the ASME

breakage of an aggregate formed between two cells placed in a

linear shear ﬂow.Simulations presented here are two dimensional,

but the method can be extended to 3D.

2 Simulation Technique

2.1 Front-tracking Method for Deformable Cells.The

simulation technique considered here is the immersed boundary

method developed by Peskin 16,and later extended by Univerdi

and Tryggvason 17 and Tryggvason et al.18 as the front-

tracking method for multiple ﬂuids with different properties.The

cells are modeled as viscous liquid drops surrounded by elastic

membranes.The viscosity of the liquid inside the cell cytoplasm

can differ from that of the liquid outside plasma.The main idea

of the front-tracking method is to use a single set of equations for

all ﬂuids,inside and outside.The governing equations for the ﬂuid

ﬂow are solved on a ﬁxed Eulerian grid,and the cell-plasma in-

terface is tracked in a Lagrangian manner by a set of moving grid

used to discretize the interface Fig.1.This way the interface

tracking is decoupled from the ﬂuid solver.Thus the method is a

mixed Eulerian-Lagrangian method.A body force term Fx,t is

introduced in the governing equations that accounts for the inter-

face.It is zero everywhere in the ﬂow except at the interface

Fx,t =

S

fx

,tx − x

dx

1

where x is a point in the ﬂow domain,x

is a point on the inter-

face,S is the entire interface,and is the Delta function.Here

f=f

s

+f

m

is an interfacial force that arises due to the elastic force f

s

generated in the membrane,and the intercellular bridging force f

m

due to cell-cell aggregation.For incompressible and immiscible

ﬂuids,the governing equations are the continuity and the Navier-

Stokes equations:

u = 0,and

u

t

+ u u

= − p + + F 2

Here ux,t is the ﬂuid velocity anywhere in the ﬂow, is the

density,p pressure,=u+u

T

is the viscous stress,and

x,t is the viscosity in the entire ﬂuid.For any point within a

blood cell =

c

and for any point outside =

p

,where

c

is the

viscosity of the liquid interior of the cell,and

p

is the viscosity

of the extracellular ﬂuid.As the cells move and deform,x,t

needs to be updated in the ﬂow ﬁeld.The details about updating

x,t can be found in 17 and are not repeated here.The gov-

erning equations are discretized spatially using a ﬁnite difference

scheme,and temporally using a two-step time-split or projection

scheme.A typical grid used in the simulation is 128128 for the

ﬂow domain and 512 for the cell membrane.Details of the time-

split scheme can be found in 19.

Once the ﬂuid velocity is obtained by the method described

above,the cells are advected in the Lagrangian manner.The ve-

locity of the cell membrane is obtained by interpolating the veloc-

ity of the ﬂuid as

ux

=

S

uxx − x

dx 3

where S denotes the entire ﬂow domain.The membrane is then

advected by

dx

dt

= ux

4

The multidimensional function used above is constructed by

multiplying one-dimensional delta functions,such as x−x

=x−x

y−y

,in two dimension.The key to the numerical

implementation of the above formulation of the membrane force

is a smooth representation of the function as 16,17

Dx =

1

2x

1 + cos

x

2x

for x 2x 5

Dx = 0 for x 2x 6

where x is the grid size.The above representation approaches

the actual delta function as the grid size approaches zero.In the

discrete form,the sharp interface is replaced by a rapid variation

over four grid points.The integrals in Eqs.1 and 3 can be

written in the discrete form as

Fx

j

=

i

Dx

j

− x

i

fx

i

7

ux

i

=

j

Dx

j

− x

i

ux

j

8

where i and j represent grid points along the interface and in the

interior of the domain,respectively.

2.2 Cell Model.The cells are modeled as viscous liquid drops

surrounded by elastic membranes.The viscosity of the cytoplas-

mic ﬂuid can be different from that of the extracellular ﬂuid.This

difference is taken into account in the front tracking method de-

scribed above.The elastic force f

s

in the membrane can be ob-

tained from a strain energy function.In the present model we

assume that the membrane behaves like a neo-Hookean material.

We note that an RBC membrane is strongly resistant to area dila-

tation.The neo-Hookean model employed here does allow area

dilatation.Membrane models that restrict area dilatation and

hence are more accurate for the RBC have been developed

20,21.The present methodology does allow incorporation of

these models.The neo-Hookean model is chosen because of its

simplicity.An accurate representation of cell membrane is not

central to this study.For the present purpose,a cell model that

takes deformability into account is sufﬁcient.For a 2D neo-

Hookean membrane,the strain energy function is given by see,

e.g.,22

W= E

s

h

1

2

+

2

2

+

1

−2

2

−2

9

where E

s

is the shear modulus of elasticity of the membrane,h is

the thickness,and

1

and

2

are the principal stretch ratios.The

tensions T

1

and T

2

in the principal directions are then given by

22

T

1

=

E

s

h

1

2

1

2

−

1

1

2

2

and T

2

=

E

s

h

1

2

2

2

−

1

1

2

2

10

For the 2D simulations considered here,the membrane is a closed

curve.The 2D cell is then equivalent to an actual 3D cell subject

to a stretching in one direction only.That is,T

1

0,T

2

=0,where

“1” indicates the “in plane” direction along the membrane,and

“2” indicates the “out of plane” direction normal to Fig.1.The

deformation

2

in the out of plane direction is not zero.But since

T

2

=0,we can express

2

in terms of

1

.Then,for a 2D cell,we

have

Fig.1 Computational domain and grids

Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1071

T =

E

s

h

3/2

3

− 1 11

where T=T

1

and =

1

.For a discretized cell,T is the tension

acting along a line segment connecting two adjacent Lagrangian

grid points on the membrane,and is the stretch ratio unde-

formed length by deformed length of the line segment.At any

Lagrangian grid point on the membrane two line segments meet.

The membrane elastic force f

s

at the Lagrangian point is then the

resultant vector of the tensions in the two adjacent segments:f

s

=T

i

e

i

−T

j

e

j

where i and j denotes two adjacent line segments,and

e

i

and e

j

are the unit tangent vectors along them.

2.3 Aggregation Model.Even though the molecular mecha-

nism of aggregation has not been fully resolved bridging versus

depletion,here we use a formalism of bond formation that has

been applied to a wide range of cell-cell adhesion phenomena

23,24.This is a convenient way to simulate cell-cell adhesion

and it does not imply a speciﬁc molecular mechanism.The inter-

membrane force f

m

due to molecular crosslinking is governed by

reaction equations with the reaction rates which are functions of

the local distance between the membranes.The reaction term for a

speciﬁc segment of a cell membrane is computed by visiting all

other cells and obtaining the minimum distance between the

neighboring cells.If the minimum distance is less than a threshold

length,the segment pair is given an opportunity to form bonds.

Assuming that both cells contribute an equal number of molecules

to bond formation,the reaction equation for bond density n

b

is

given by

n

b

t

= 2

k

+

n −

n

b

2

2

− k

−

n

b

2

12

where n is the density of the cross-linking molecules on each cell,

and k

+

and k

−

are the forward and reverse reaction rate coefﬁ-

cients,respectively.Note that the above equation is slightly dif-

ferent from that for a receptor-ligand interaction,but consistent

with the models of Zhu 24.The reaction rates are computed as

k

+

= k

+

0

exp

−

k

ts

l − l

o

2

2K

B

T

,if x = l l

t

13

k

−

= k

−

0

exp

−

k

b

− k

ts

l − l

o

2

2K

B

T

14

where x is the distance between two facing segments of the cells,

l

t

is the threshold distance below which the bond formation is

initiated,l =x,and l

o

are the stretched and unstretched bond

length,k

b

is the spring-constant force per stretched length,k

ts

is

the transition spring constant,K

B

is the Boltzmann constant,T is

the absolute temperature,and k

+

0

,and k

−

0

are the ratecoefﬁcients in

equilibrium.The bonds behave like stretched springs,and the

force per bond is given by

f

b

= k

b

l − l

0

.15

The aggregation force f

m

per unit length of the cell membrane is

then obtained as

f

m

= f

b

n

b

x

x

.16

Note that in this formulation the bond density,rate coefﬁcients,

bond length,and the intermolecular force are not constant and

vary along the contact area of the cells.The values of the param-

eters,such as the spring constants and the reaction rates used in

the present study are given in Table 1.

In the present computations,the threshold distance l

t

is taken

to be about 1.2 times the unstretched length l

0

;thus the bonds do

not stretch signiﬁcantly.We further emphasize that the detailed

molecular mechanism of aggregation is not fully understood,and

the “stretched spring” model assumed here may not be the best

model for RBC aggregation.The theory of polymer bridging is

poorly understood 25.There is,however,evidence in support of

the “stretched spring” model used here.In experiments,RBC ag-

gregation is mediated by Dextran molecules 3,12–15.Many ex-

perimental studies,such as those by Chien and Jan 3 and Brooks

26,concluded that the Dextran polymers behave like “ﬂexible

linear chains with random coiling.” The more recent works of

Swenson et al.25 and Barshtein,Tamir,and Yedgar 27 also

concluded that the polymer bridges are “partly stretched” chains

between the cells.The theoretical work of Skalak and Zhu 28

also supports the “stretched spring” model for the RBC aggrega-

tion,similar to the one described here.

In the present study we consider two-dimensional simulations.

The ﬂow domain considered here is a Couette ﬂow in a channel

driven by two walls at the top and bottom that move in directions

opposite to each other.In the absence of any cell,the undisturbed

ﬂow is thus a linear shear ﬂow with shear rate denoted by

˙

.The

distance between the top and bottom wall is about 25 m,typical

of a blood vessel dimension in microcirculation.The length of the

channel is also taken as 25 m,and the ﬂow at the left and right

boundaries of the channel is assumed to be periodic.The govern-

ing equations are made dimensionless using a characteristic length

a such as cell radius,the applied shear rate

˙

,and the plasma

viscosity

p

.The dimensionless time is then t˙

where t is the

physical time.The characteristic scale for velocity is a˙

.This

results in the following major dimensionless parameters:

Reynolds number,Re =

˙

a

2

p

Membrane rigidity,E

*

=

p

˙

a

E

s

h

Table 1 Physical and dimensionless parameters used in the simulations.Here a

0

is a characteristic length scale which is often

the radius of an undeformed cell.For an elliptic cell studied here,the semimajor axis of the undeformed cell is taken as the cell

radius.The number of molecules per cell is in the range 50–300.The number of bonds formed in equilibrium is in the range

30–100.The number density n is obtained assuming the cell surface area 100 m

2

.

Physical parameters Dimensionless parameters

Cell radius a 4 m a/a

0

1–0.667

Vessel diameter D 25 m D/a 2

˙

25-225 s

−1

Re 0.0003–0.003

E

s

1.2-3010

−3

dyn/cm

E

*

0.02–0.5

k

b

10

−6

–10

−8

N/m

K

b

0.125–3.0

p

=1.2 cP

c

=1.2–6 cP

1–5

k

ts

10

−10

–10

−12

N/m

C

1

*

0.3–0.003

k

+

0

10

−7

–10

−9

s

−1

cm

−2

k

+

0

*

0.1–0.001

k

−

0

10

−11

–10

−13

s

−1

cm

−2

k

−

0

*

10

−3

–10

−5

1072/Vol.127,DECEMBER 2005 Transactions of the ASME

Viscosity ratio, =

c

p

Spring constant,K

b

=

k

b

p

˙

a

Values of the dimensionless parameters used in this study are

given in Table 1.The “membrane rigidity” E

*

is the inverse of the

shear modulus E

s

.For a normal RBC,E

s

=610

−3

dyn/cm,and

=5.However,under pathological conditions,such as sickle cell

disease and malaria,signiﬁcant deviation from these values may

occur.Further,the parameters governing the kinetics of the RBC

aggregation under pathological conditions are not clearly known.

Clinical studies do show a tendency of high aggregation under

such conditions.Thus considering a range of the parameters as

used in our simulations Table 1 is of biological interest.

3 Results and Discussion

3.1 Simulation of Single Cell.We ﬁrst discuss some results

on the simulation of deformation of a single cell subject to a linear

shear ﬂow.The geometry of the ﬂow domain is the same as that

described in Fig.1.A single,isolated cell is ﬁrst placed at the

center of the channel.The initial cell shape is assumed to be

elliptic with an aspect ratio of 0.375.As the ﬂow starts,the cell

deforms,and the cell surface and the cytoplasmic ﬂuid rotate in a

tank-treading fashion.The trajectory of a ﬁxed point on the cell

surface and streamlines obtained from the velocity ﬁeld are shown

in Figs.2a and 2b,respectively.In Fig.2c,we plot the steady

shape of the cell at various shear rates

˙

=25,75,and 225 s

−1

.

The cell elongates and aligns itself with the direction of the ﬂow

as the shear rate is increased.The effect of varying shear modulus

is described in Fig.2d which also shows that the cell becomes

more elongated with decreasing rigidity.The elongation of the cell

is accompanied by a slight increase in the surface area which is

found to be at most 2%.In Fig.2e,the result for =5 is shown.

A rotation of the cell like a rigid body at higher cytoplasmic

viscosity is observed.The direction of rotation is consistent with

the direction of the applied shear.In the ﬁgure only one instanta-

neous orientation is shown after the cell starts rotating.There is

almost no visible change in the cell shape.Though the cell mem-

brane undergoes tank-treading motion,its velocity is much less

than the earlier cases of =1 shown in Figs.2b–2d.These

results agree qualitatively with those of previous researchers,e.g.,

Pozrikidis 29 and Ramanujan and Pozrikidis 30.These studies

considered 3D cell deformation at zero inertia,whereas the

present study is 2D and inertia is not zero,though small.

It should be noted that the smoothing of the cell membrane

using the discrete delta function in the immersed boundary

method introduces a numerical viscosity.While such a numerical

viscosity is extremely bad for ﬂow computations at high Reynolds

numbers,it is not detrimental for the low Reynolds number cases

considered here Re0.001.The projection method that is used

here to solve the Navier-Stokes equations satisﬁes the incompress-

ibility condition up to the precision of the computer 10

−14

.Thus

mass conservation is exactly satisﬁed on a local basis.We also

keep track of the cell volume during the simulations.The change

in the cell volume is less than ±0.1% for a grid resolution of

128128.The cell does not show any net increase/decrease in

volume over a long time.Thus the small error is probably due to

the numerical error in computing the cell volume,rather than an

error in the front-tracking method.Since the simulations presented

here are 2D,there is no direct way of validating the results on cell

deformation.However,we have been able to extend the method to

3D,and thus been able to compare the deformation of a 3D

spherical capsule obtained from our immersed boundary simula-

tion with the boundary integral calculation of Pozrikidis 29.The

comparison is shown in Fig.2f for a low Reynolds number of

0.001.Note that Eggleton and Popel 31 had also used a similar

immersed boundary method for large deformation of capsules.A

detailed comparison of the immersed boundary results,boundary

integral results,and asymptotic results in the limit of small defor-

mation,is given in their paper.

3.2 Simulation of Aggregation.Simulations are performed

for aggregation of two elliptic cells in Couette ﬂow at shear rates

typical of physiological conditions.Simulations with biconcave

shape of an RBC are also performed and described later.Table 2

lists the range of shear rate,spring constant,membrane rigidity,

and viscosity ratio used in the simulations.The simulation strat-

egy is as follows.First we let the cells aggregate under a no-ﬂow

condition as shown in Fig.3a.The static,no-ﬂow condition is

chosen as the starting condition since aggregation is believed to be

the highest under such a condition.To initiate aggregation,we

ﬁrst place the cells in close proximity such that the minimum

distance between the adjacent cells is below the threshold distance

l

t

.The number of bonds and the contact area between two cells

increase until they reach steady values.The cells reach ﬁnal

steady shapes,and an aggregate is formed.Once the equilibrium

is reached under no-ﬂow condition,the shear ﬂow is imposed.

Under the action of the shear ﬂow,the aggregates roll,deform,

and depending on the spring constant,shear rate and other param-

eters,they either remain intact or break apart.In the following

section we describe such dynamic events as revealed by our

simulations.

Figures 3b–3e show the instantaneous orientations of the

rolling aggregate at various ,and K

b

.Figure 3b corresponds to

˙

=75 s

−1

,K

b

=2,=1.The aggregate rolls almost like a rigid

body.Consequently,the tank-treading motion of an individual cell

in the aggregate is not present any more.The trace of a point on

the membrane surface is shown which clearly indicates no tank

treading.The individual cells undergo some degree of deforma-

tion by the ﬂuid.The contact area and the total number of bonds

remain nearly steady.Consequently,the aggregate does not break

and remains stable.Interestingly the contact area does not remain

ﬂat;it changes shape during each turn of the aggregate.In Fig.4,

the instantaneous ﬂow ﬁeld around the rolling aggregate is pre-

sented.We also note that the rolling of the aggregate induces

ﬂuctuations in the surrounding ﬂow as evident from the color

contours in the ﬁgure.Such disturbed ﬂow ﬁeld can have an effect

on the cells,such as leukocytes and platelets adhering to the wall,

and can either dislodge them or stabilize the cell-to-wall adhesion.

The folded shape of the cell interface in these ﬁgures is some-

what peculiar.There is little experimental work reporting the

shape of a rotating aggregate.Most previous works have reported

aggregate shapes under static conditions 3,28,32.However,

Skalak et al.6 reported the shape Fig.8 in their paper of an

aggregate obtained via electron microscopy which shows the sig-

moidal shape of the interface and thus is very similar to the ones

found here.It is thus possible to see such shapes in 3D.The

sequence of four snapshots given in Fig.3b suggests that the

interface is rotating like a ﬂexible but nonextensible element or,

thread.This is expected since the membrane can easily deform

but not extend.Simulations of deformable and nonextensible ele-

ments in shear ﬂow have shown similar nature of folding see,

e.g.,32.Our simulations also show strong curvatures near the

edge of the interface,particularly at high K

b

.Skalak and Chien

33 showed that as the attractive force increases,the curvature

near the edge also increases.Their work considered a two-cell

rouleaux similar to the present work,and the semicircular shape

of individual cell obtained in their work is similar to our result

corresponding to the no-ﬂow condition Fig.3a.It is shown

later Sec.3.3 that inclusion of bending resistance in the model

would prevent the folding of the interface and the membrane.

Results at a lower spring constant K

b

=0.875 are shown in Fig.

3c.The contact area is now less than before as many bonds are

broken by the ﬂuid force.However,the bond breaking rate is

slow,and the aggregate does not break instantly.We have per-

formed this simulation up to t˙

=1101.5 s at which point

Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1073

about 35% of bonds still exist.At K

b

=0.175,shown in Fig.3d,

the aggregate deforms signiﬁcantly from its static shape.The as-

sembly very quickly breaks,and the cells are completely detached

from each other.The effect of increasing is shown in Fig.3e.

Increased cytoplasmic viscosity reduces the deformability of the

cells as seen before in Fig.2e for a single cell.The aggregate

shape is now nearly circular,and individual cells almost retain

their semicircular shape that they obtain under no-ﬂow conditions.

Increased viscosity makes the bonds very stable,and hence no

disaggregation is observed.

Fig.2 Simulation of single cell deformation.„a… An elliptic cell in a shear ﬂow:The trajectory of a point on the cell

surface is indicated to show tank treading.„b… Streamlines constructed fromthe velocity ﬁeld in and around the cell.

The cytoplasmic ﬂuid rotates with the cell surface.Also shown are the effects of varying shear rate „c…,membrane

rigidity „d…,and cell-to-plasma viscosity ratio „e….Initial cell shape is shown by the solid line.In „c… Œ ˙

=25 s

−1

,¯˙

=75 s

−1

,---

˙

=225 s

−1

.In „d… Œ E

*

=0.02,¯E

*

=0.1,--- E

*

=0.5.In „e… =5.In „f…,the Taylor deformation parameter for

an initially spherical capsule deforming in a shear ﬂow is shown.The solid line corresponds to the present front-

tracking method extended to 3D,and the symbols represent the boundary integral results †29‡.

1074/Vol.127,DECEMBER 2005 Transactions of the ASME

In the absence of bending resistance,as in the above ﬁgures,the

cell shape is not unique at zero transmural pressure,i.e.,the pres-

sure difference across the RBC membrane,inside pressure - out-

side pressure under hydrostatic conditions.A detailed discussion

on the effect of transmural pressure on the resting shape of a cell

was given by Pozrikidis 34,35.Thus in the simulations presented

here,an initially elliptic cell that undergoes deformation in a shear

ﬂow may not return to the same initial shape if the ﬂow is

stopped.Note that in our simulations,the transmural pressure is

not a parameter that we specify.The pressure inside and outside

the cell evolves naturally.Furthermore,the shear ﬂow considered

in our case is not a pressure driven ﬂow.Thus in absence of a cell,

the pressure is constant everywhere in the ﬂow ﬁeld and can be

taken to be zero.Under the static no-ﬂow condition,the transmu-

ral pressure is zero.As the ﬂow starts,the cell deforms,and an

internal pressure higher than the external pressure is established.

Thus the transmural pressure is positive for a deformed cell in a

ﬂow.The cell shape and the transmural pressure under the inﬂu-

ence of the external shear is then unique for a given shear rate and

membrane properties.When two cells form an aggregate,the

transmural pressure is again found to be positive in both static and

rolling conditions.If the aggregate breaks and the cells detach

completely from each other with the shear ﬂow still acting the

cells will return to the deformed shape that they assumed when

placed in a shear ﬂow without aggregation.If the aggregate does

not break at all,but the shear ﬂow is switched off,the shape of the

aggregate returns to the initial steady shape that was obtained

under static conditions.A negative transmural pressure is not ob-

served in our simulations.The conclusions remain the same even

when bending resistance is included.

In Fig.5,we study the effect of varying K

b

,,

˙

and E

*

on the

evolution of total number of bonds between the cells.The increase

in the number of bonds under no-ﬂow condition is shown in Fig.

5a.The equilibrium is reached at around t˙

95.At this point

the shear ﬂow is started,and the number of bonds starts changing

as the aggregate begins to roll.The number of bonds decreases

very rapidly for K

b

=0.175.In this case 95% of the bonds are

broken in 0.25 s t˙

19.At K

b

=0.875 and 1,the decrease is

quite slow.Finally,at K

b

=2 and 3 the aggregate is stable,and the

number does not decrease at all.The effect of viscosity ratio is

shown in Fig.5b,where results for =1 and 5 are plotted.

Higher cytoplasmic viscosity increases the stability of the aggre-

Fig.3 Simulation of aggregate formation between two cells.

„a… Evolution of cell shape under the no-ﬂow condition.Three

time instants are t˙

=0,50,and 97.„b… ˙

=75 s

−1

,K

b

=2,=1.

The four instants are approximately at equal time interval dur-

ing a rotation of the aggregate.„c… Effect of reduced bond

strength at K

b

=0.875,=1,

˙

=75 s

−1

.„d… K

b

=0.175,=1,

˙

=75 s

−1

.„e… Effect of cell viscosity:K

b

=0.875,=5,E

*

=0.096,˙

=75 s

−1

.

Fig.4 Fluid velocity vectors and color contours showing the

ﬂuid velocity component in the x direction during aggregate

rolling

Table 2 Simulation runs for two cell aggregate

Simulation runs ˙

s

−1

K

b

E

*

1 75 0.175 1 0.096

2 75 0.875 1 0.096

3 75 1.0 1 0.096

4 75 2.0 1 0.096

5 75 3.0 1 0.096

6 75 0.175 5 0.096

7 75 0.875 5 0.096

8 75 0.875 1 0.02

9 75 0.875 1 0.5

10 25 2.625 1 0.096

11 225 0.292 1 0.288

Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1075

Fig.5 Evolution of the number of bonds with time.„A… Effect of spring constant

K

b

.Curves „a… to „e… represent K

b

=0.175,0.875,1.0,2.0,and 3.0.„B… Effect of

viscosity ratio .The four curves are for „a… K

b

=0.175,=1,„b… K

b

=0.875,=1,

„c… K

b

=0.175,=5,„d… K

b

=0.875,=5.„C… Effect of varying shear rate „a…,„b…,and

„c… in the plot… and varying membrane rigidity „b…,„d…,and „e….„a…

˙

=225 s

−1

,„b…

˙

=75 s

−1

,and „c… ˙

=25 s

−1

.„b… E

*

=0.096,„d… E

*

=0.5,and „e… E

*

=0.02.

1076/Vol.127,DECEMBER 2005 Transactions of the ASME

gate.For K

b

=0.175 at which a rapid aggregate breakage was ob-

served earlier at =1,now results in a much slower bond break-

ing at =5.

The effects of shear rate and membrane rigidity are considered

in Fig.5c.Curves a,b,and c are the results for

˙

=225,75,

and 25 s

−1

,respectively.At the highest shear rate the aggregate is

quickly broken.At the intermediate shear rate,bond breaking

takes place over a longer time.At the lowest rate,the aggregate is

completely stable.The effect of membrane rigidity is elucidated

by the curves b,d,and e in the same ﬁgure.The curves

represent E

*

=0.096,0.5,and 0.02 with 0.096 corresponding to

the normal RBCs.The bond breaking decreases for both cases

when E

*

increases or decreases from this value.Thus it appears

that E

*

=0.096 is near an optimum value at which the aggregate

breaks at the highest rate.We also note that the aggregate breaking

does not happen instantly.For different combinations of param-

eters,the breakage can take place over a few seconds.

Figures 5a–5c also suggest that the number of bonds oscil-

lates as the aggregate rolls in a shear ﬂow.The shape and orien-

tation of the aggregate at ﬁve different time instants are shown in

Fig.6 together with the number of bonds.The oscillations are

related to the rolling orientation of the aggregate.A rapid bond

breakage is initiated when the contact area is along the ﬂow di-

rection.On the contrary,the number of bonds increases when the

contact area is normal to the ﬂow direction.In the ﬁrst case,the

bonds break primarily by the shearing force of the ﬂuid acting

parallel to the contact area.In the latter orientation,the bond

breakage is by pulling the cells in opposite directions by the ﬂuid

forces normal to the contact area.The results are in agreement

with the experimental observation by Chien and Jan 3 and Chien

et al.36 that the shearing force is more effective to break an

aggregate than the normal force.Figures 5a–5c also suggest

that the frequency and amplitude of the oscillations in bond num-

bers are not the same for all cases.For easily deformable cells the

amplitude of oscillation is higher and the frequency is lower.With

decreasing deformability,e.g.,at higher and K

b

the amplitude

of oscillation is lower but the frequency is higher.

It should be noted that the total number of bonds depends both

on the change in bond density and change in contact area.If the

aggregate is stable,the contact area is found to remain nearly the

same,but the bond density varies with time giving rise to the

observed oscillatory behavior of the total number of bonds.If,on

the other hand,the aggregate is breaking,the contact area and the

bond density both are decreasing.At any given time,the bond

density is found to be nearly constant over the contact area,except

near the two edges where it smoothly goes to zero.Thus it seems

that instead of using Eq.12,one can assume a constant bond

density and contact area to get the similar result with an appropri-

ate choice of the spring constant.Such an approach would be

correct for quasi-steady simulations.However,the bond density

and contact area evolve in time,and they are not known a priori.

Thus Eq.12 is a convenient way to model the bond kinetics.

Figures 5a–5c show that with different choices of the spring

constant one can either get a stable or breaking aggregate for a

constant shear rate.Alternatively one can ﬁx the spring constant to

a speciﬁc value and get either a breaking or a stable aggregate

under varying shear rates.

The instantaneous orientation of an aggregate deﬁned by an

angle that a line joining the center of mass of each cell forms

with the direction of the shear ﬂow x axis is shown in Fig.7a.

The location of the center of mass of each cell can be obtained

from the simulation data.As the aggregate rolls,varies periodi-

cally from 0 to 2.When =0 and ,the contact area is normal

to the ﬂow direction,and when =/2 and 3/2,it is parallel.

The results show that as the spring constant or cell viscosity in-

creases the period of oscillation decreases.This can be understood

from the theory of rotation of an ellipsoidal particle in a shear

ﬂow 37.In the limit of a 2D circular particle,the angular orien-

tation changes linearly with time,and the period of rotation is

minimum.A departure from the circular shape means nonmono-

tonic behavior in angular orientation and angular velocity,and an

increase in period of rotation.At =5 at which the aggregate

shape is nearly circular,the result is close to Jeffery’s solution.

However,at reduced or K

b

,is nonmonotonic,and the period

of rotation is higher.The angular velocity

˙

shown in Fig.7b is

also periodic with the period half that of ,as expected from the

Jeffery’ solution.The angular velocity is low when the contact

Fig.6 Explanation of oscillations in the total number of bonds as the aggregate rolls in

shear ﬂow.The shape and orientation of the aggregate are shown at a few time instants.

Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1077

area is normal to the ﬂow direction,and high when it is parallel.

The simulations presented here are 2D.Quantitative differences

may arise if 3D simulations are performed.However,the major

conclusion of the paper,that aggregability increases with reduced

deformability,is expected to hold even in 3D.

3.3 Effect of Viscoelasticity and Bending Resistance.Vis-

coelastic nature of the membrane of a red blood cell has been

reported in the literature e.g.,38 and measured to be

about 0.001 dyn s/cm.The results described in the previous sec-

tion do not include the effect of membrane viscoelasticity.In order

to include this effect,Eq.11 is modiﬁed by adding a term

m

/d/dt,where

m

is the membrane viscoelasticity.Several

dimensionless parameters then can be formed,such as

m

/

p

a

and

m

˙

/E

s

,where a is the characteristic length of the cell.The

parameter

m

/

p

a is about 100 for a normal RBC implying a

strong effect of viscoelasticity.The second parameter,also called

Deborah number 39,40,is in the range 0.01 and 1 for shear rates

between 1 s

−1

and 100 s

−1

.The Deborah number compares the

membrane viscous and elastic forces,and could be more relevant

for the present problem.Since RBC aggregation is dominant at a

low shear rate,a corresponding low value of the Deborah number

indicates that the membrane viscous force may be small compared

to the elastic force justifying the exclusion of the viscoelastic

effect in the results presented so far.Nevertheless,we have in-

cluded the viscoelasticity in our model and repeated several simu-

lations for elliptic cells at

m

/

p

a=20.Sample results are shown

in Figs.8 and 11 for one representative case.The results suggest

that the rotation rate of the rolling aggregate and the bond break-

ing rate are reduced when viscoelasticity is included.However,

the shape of the cells,the oscillating nature of the bond breaking,

and the overall decrease in the number of bonds remain nearly the

same compared to the case without viscoelasticity.The decrease

in rotation rate and the bond breaking rate is probably due to the

memory effect of the membrane viscosity.Adetailed investigation

of this effect is beyond the scope of this paper.It should be noted

that there has been very little theoretical work done on the defor-

mation of viscoelastic membranes in shear ﬂow.

The RBC membrane also has a bending modulus E

B

of about

1.810

−12

dyn cm.The bending resistance prevents formation of

high curvature during deformation.The RBC also has a biconcave

shape under resting condition.The shape of the cell is expected to

have an effect on aggregation.The bending resistance is also re-

sponsible for the biconcave shape of the cell.To include the bend-

ing resistance in our simulation,we follow the approach by

Pozrikidis 41.Equation 11 is modiﬁed as T=t+qn,where

is the membrane tension obtained by differentiating the strain-

energy function with respect to the extension ratio,q=dm/dl is

the transverse shear tension,m=E

B

−

r

is the bending moment,

is the local curvature,

r

is the reference curvature in resting

conﬁguration,l is the arc length along the membrane,and t and n

are the unit tangent and normal vectors.The relevant dimension-

less parameter is E

B

/a

2

E

s

;for a normal RBC and a4 m,this

is about 10

−3

.First we have validated our simulations by repeating

some results given in Pozrikidis 41 which are shown in Fig.9.

This ﬁgure shows the stages of evolution of an initially elliptic

cell and a biconcave cell toward the reference state of a circle.The

2D immersed boundary simulations show qualitatively similar be-

havior of the 3D boundary integral simulations by Pozrikidis 41.

We have repeated several simulations of cell aggregation pre-

sented in the last sections with the bending moment included.Two

representative cases are shown in Figs.10 and 11.Both the static

Fig.7 „a… Angular orientation and „b… velocity of an aggregate

at ˙

=75 s

−1

.The dash line is Jeffery’s †37‡ solution.-- K

b

=0.875,=5;-- K

b

=0.875,=1;-- K

b

=3,=1.

Fig.8 Effect of membrane viscoelasticity on the instanta-

neous shape of the rolling aggregate.„a… Viscoelasticity in-

cluded.Four different time instants are shown:t˙

=42,48,53,

and 58.„b… Viscoelasticity not included.Shapes at t˙

=42 and

48 are shown.

Fig.9 Relaxation of an initially elliptic „left… cell and a biconcave „right… cell toward their resting circular

shape under the action of bending resistance at E

B

/„a

2

E

s

…=0.003,=1

1078/Vol.127,DECEMBER 2005 Transactions of the ASME

and rolling cases for elliptic and discoid shapes are considered.

For the elliptic cells,the results under the no-ﬂow condition show

that the inclusion of bending moment removes any high curvature

near the edge of the contact area.However the contact area re-

mains nearly the same in the two cases.The shape of the rolling

aggregate is similar to the one observed without bending moment,

except that the folding of the contact line is reduced.The variation

in the number of bonds Fig.11 also shows similar trends in

terms of bond breaking and oscillations for the two cases.Similar

results were found for all cases that were repeated.Thus we con-

clude that for an aggregate rolling in a shear ﬂow,bending resis-

tance may not have a signiﬁcant effect.Of course,under static

condition,as seen in Fig.10,bending resistance has an effect on

cell shape.We have also performed one simulation using bicon-

cave shape Fig.10.The overall trends of aggregation,aggregate

rolling,and bond breakage are similar to those for elliptic cells.

As expected,the bond breaking rate for the biconcave cells is

slower due to higher contact area.

4 Conclusion

In this paper we present 2D numerical simulations on the ag-

gregation of deformable cells in a Couette ﬂow.The simulation

technique couples the bulk hydrodynamics with the rheology of

the cells and the kinetics of intercellular bond formation.The

focus of the paper is to describe 1 the effect of cell rheological

Fig.10 Effect of bending resistance and biconcave shape on the shape evolution of an aggre-

gate.Both elliptic and biconcave cells are considered.Figures at left and right show the formation

of the aggregate under static and rolling conditions,respectively.In „a…,the solid line is for bend-

ing moment included,and the dash line is without bending moment.

Fig.11 Effect of viscoelasticity and bending resistance on the

number of bonds.Bold line is for elliptic cells without bending

resistance and viscoelasticity.…,viscoelastic membrane;—,

elliptic cell with bending resistance;–––––,biconcave discoid

cell with bending resistance.

Journal of Biomechanical Engineering DECEMBER 2005,Vol.127/1079

properties on the stability of an aggregate and 2 the dynamics of

the rolling motion of the aggregate in the shear ﬂow.

The major conclusion from this work is that the deformability

of the cells as determined by cell-to-plasma viscosity ratio and the

shear modulus of the cell membrane,has a signiﬁcant effect on

the stability and motion of the aggregate.Aggregates made of

deformable cells are easily breakable by a shear ﬂow,while those

made by less deformable cells are not.These observations are

consistent with many clinical studies which show that increased

aggregation of red blood cells is often accompanied by reduced

cell deformability.Such is the case in sickle cell disease,malaria,

and hypertension.The present study,however,is the ﬁrst theoret-

ical study to establish a direct link between deformability and

aggregability.

We also observe that in a shear ﬂow,the number of bonds

oscillates with time.The oscillations are related to the rolling

orientation of the aggregate.The bond breakage is maximum

when the contact area is parallel to the ﬂow direction,and mini-

mum when the contact area is normal.This result implies that the

shearing force is more efﬁcient to break the bonds than the normal

pulling force.The frequency and amplitude of the oscillations are

dictated by the rheology of the cells and the bond strength.A

deformable aggregate has higher amplitude and longer period of

rotation than a stable aggregate.This is evident when angular

orientation and velocity are taken into consideration which are

found to match qualitatively with the analytical prediction of Jef-

fery 37.

Acknowledgment

This work was supported by a grant from the National Heart,

Lung,and Blood Institute HL 52684.We also thank the referees

for their comments which have signiﬁcantly improved the quality

of the paper.The computations were performed on IBM p690 at

NCSA,UIUC and HP V2500 at CACR,Caltech.

References

1 Rampling,M.W.,Meiselman,H.J.,Neu,B.,Baskurt,O.K.,2004,“Inﬂuence

of Cell Speciﬁc Factors on Red Blood Cell Aggregation,” Biorheology 41,pp

91–112.

2 Ami,B.,Barshtein,G.,Zeltser,D.,Goldberg,Y.,Shapira,I.,Roth,A.,Keren,

G.,Miller,H.,Yedgar,S.,2001,“Parameters of Red Blood Cell Aggregation

as Correlates of the Inﬂammatory state,” Am.J.Physiol.,280,pp.H1982–

H1988.

3 Chien,S.and Jan,K.-M.,1973,“Ultrastructural Basis of the Mechanism of

Rouleaux Formation,” Microvasc.Res.5,pp.155–166.

4 Neu B.,and Meiselman H.J.2002,“Depletion-Mediated Red Blood Cell

Aggregation in Polymer Solutions,” Biophys.J.,83,pp.2482–2490.

5 Popel,A.S.and Johnson,P.C.,2005,“Microcirculation and Hemorheology,”

Annu.Rev.Fluid Mech.37,pp.43–69.

6 Skalak,R.,Zarda,P.R.,Jan,K.-M.,and Chien,S.,1981,“Mechanics of

Rouleau Formation,” Biophys.J.,35,pp 771–781.

7 Evans,E.A.,1985,“Detailed Mechanics of Membrane-Membrane Adhesion

and Separation.” Biophys.J.48,pp 175–183.

8 Murata,T.,and Secomb,T.W.,1988,“Effects of Shear Rate on Rouleau

Formation in Simple Shear Flow,” Biorheology,25,pp.113–122.

9 Chen,J.,and Huang,Z.,1996,“Analytical Model for Effects of Shear Rate on

Rouleau Size and Blood Viscosity,” Biophys.Chem.96,pp.273–279.

10 Lim,B.,Bascom,P.A.J.,and Cobbold,R.S.C.,1997,“Simulation of Red

Blood Cell Aggregation in Shear Flow,” Biorheology,34,pp.423–441.

11 Goldsmith,H.L.,and Karino,T.,1980,“Physical and Mathematical Models of

Blood Flow:Experimental Studies, Erythrocyte Mechanics and Blood Flow,

Cokelet,G.R.,Meiselman,H.J.,and Brooks,D.E.,eds.,John Willey & Sons.

12 Bishop,J.J.,Nance,P.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2001a,

“Effect of Erythrocyte Aggregation on Velocity Proﬁles in Venules,” Am.J.

Physiol.280,pp.H222–H236.

13 Bishop,J.J.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2001b,“Effects

of Erythrocyte Aggregation and Venous Network Geometry on Red Blood Cell

Axial Migration,” Am.J.Physiol.281,pp.H939–H950.

14 Bishop,J.J.,Nance,P.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2001c,

“Erythrocyte Margination and Sedimentation in Skeletal Muscle Venules,”

Am.J.Physiol.281,pp.H951–H958.

15 Bishop,J.J.,Popel,A.S.,Intaglietta,M.,and Johnson,P.C.,2002,“Effect of

Aggregation and Shear Rate on the Dispersion of Red Blood Cells Flowing in

Venules,” Am.J.Physiol.283,pp.H1985–H1996.

16 Peskin,C.S.1977,“Numerical Analysis of Blood Flow in the Heart,” J.

Comput.Phys.,25,pp.220–233.

17 Unverdi,S.O.and Tryggvason,G.1992,“A Front-Tracking Method for Vis-

cous Incompressible Multi-Fluid Flows,” J.Comput.Phys.100,pp.25–37.

18 Tryggvason,G.,Bunner,B.,Esmaeeli,A.,Juric,Al-Rawahi,N.,Tauber,W.,

Han,J.,Nas,S.,and Jan.,Y.-J.,2001,“A Front Tracking Method for the

Computations of Multiphase Flow,” J.Comput.Phys.169,pp 708–759.

19 Bagchi,P.,and Balachandar,S.,2002,“Steady Planar Straining Flow Past a

Rigid Sphere at Moderate Reynolds Number,” J.Fluid Mech.,466,pp.365–

407.

20 Skalak,R.,Tozeren,A.,Zarda,P.R.,and Chien,S.,1973,“Strain Energy

Function of Red Blood Cell Membranes,” Biophys.J.,13,pp.245–264.

21 Evans,E.A.,and Skalak,R.1980,Mechanics and Thermodynamics of

Biomembranes,CRC Press,Boca Raton,FL.

22 Barthes-Biesel,D.,Diaz,A.,and Dhenin,E.,2002,“Effect of Constitutive

Laws for Two-dimensional Membranes on Flow-induced Capsule Deforma-

tion,” J.Fluid Mech.,460,pp.211–222.

23 N’Dri,N.A.,Shyy,W.,and Tran-Son-Tay,R.,2003,“Computational Model-

ing of Cell Adhesion and Movement Using a Continuum-Kinetics Approach,”

Biophys.J.,85,pp.2273–2286.

24 Zhu,C.,1991,“A Thermodynamic and Biomechanical Theory of Cell Adhe-

sion.Part 1:General Formulism,” J.Theor.Biol.150,pp 27–50.

25 Swenson,J.,Smalley,M.V.,and Hatharasinghe,H.L.M.,1998,“Mechanism

and Strength of Polymer Bridging Flocculation,” Phys.Rev.Lett.,81,pp

5840–5843.

26 Brooks,D.E.,1973,“The Effect of Neutral Polymers on the Electrokinetic

Potential of Cells and other Charged Particles.III:Experimental Studies on the

Dextran/Erythrocyte System,” J.Colloid Interface Sci.,43,pp.701–713.

27 Barshtein,G.,Tamir,H.,and Yedgar,S.,1998,“Red Blood Cell Rouleaux

Formation in Dextran Solution:Dependence on Polymer Conformation,” Eur.

Biophys.J.,27,pp.177–181.

28 Skalak,R.,and Zhu,C.,1990,“Rheological Aspects of Red Blood Cell Ag-

gregation,” Biorheology,27,pp 309–325.

29 Pozrikidis,C.,1995,“Finite Deformation of Liquid Capsules Enclosed by

Elastic Membranes in Simple Shear Flow,” J.Fluid Mech.,297,pp.123–152.

30 Ramanujan,S.,and Pozrikidis,C.,1998,“Deformation of Liquid Capsules

Enclosed by Elastic Membranes in Simple Shear Flow:Large Deformations

and the Effect of Fluid Viscosities,” J.Fluid Mech.361,pp.117–143.

31 Eggleton,C.D.,and Popel,A.S.,1998,“Large Deformation of Red Blood

Cell Ghosts in a Simple Shear Flow,” Phys.Fluids 10,pp.1834–1845.

32 Derganc,J.,Bozic,B.,Svetina,S.,and Zeks,B.,2003,“Equilibrium Shapes of

Erythrocytes in Rouleau Formation.” Biophys.J.,84,pp.1486–1492.;Torn-

berg,A.-K.,and Shelley,M.,2004,“Simulating the Dynamics and Interactions

of Flexible Fibers in Stokes Flows,” J.Comput.Phys.,196,pp.8–40.

33 Skalak,R.,and Chien,S.,1983,“Theoretical Models of Rouleau Formation

and Disaggregation,” Ann.N.Y.Acad.Sci.,416,pp.138–148.

34 Pozrikidis,C.,2003a,“Membrane Theory for Capsules and Cells, in Model-

ing and Simulation of Capsules and Biological Cells,edited by C.Pozrikidis,

Chapman and Hall/CRC Press Mathematical Biology and Medicine Series.

Boca Raton,Florida.

35 Pozrikidis,C.,2003b,“Deformed Shapes of Axisymmetric Capsules Enclosed

by Elastic Membranes,” J.Eng.Math.,45,pp.169–182.

36 Chien,S.,Sung,L.A.,Sinchon,S.,Lee,M.M.L.,and Skalak,R.,1984,“En-

ergy Balance in Red Cell Interactions.” Ann.N.Y.Acad.Sci.416,pp.138–

148.

37 Jeffery,G.B.1922,“The Motion of Ellipsoidal Particles Immersed in a Vis-

cous Fluid,” Proc.R.Soc.London,Ser.A,102,161–170.

38 Evans,E.A.,and Hochmuth,R.M.,1976,“Membrane Viscoelasticity,” Bio-

phys.J.,16,pp 1–11.

39 Barthes-Biesel,D.,1991,“Role of Interfacial Properties on the Motion and

Deformation of Capsules in Shear Flow,” Physica A,172,pp.103–124.

40 Barthes-Biesel,D.,and Sgaier,H.,1985,“Role of Membrane Viscosity in the

Orientation and Deformation of a Spherical Capsule Suspended in Shear

Flow,” J.Fluid Mech.,160,119–136.

41 Pozrikidis,C.,2001,“Effect of Membrane Bending Stiffness on the Deforma-

tion of Capsules in Simple Shear Flow,” J.Fluid Mech.,440,pp.269–291.

1080/Vol.127,DECEMBER 2005 Transactions of the ASME

## Comments 0

Log in to post a comment