Fluid-structure modeling of flow-induced alveolar epithelial cell deformation

Computers and Structures 85 (2007) 1066–1071 www.elsevier.com/locate/compstruc

Fluid-structure modeling of flow-induced alveolar epithelial cell deformation H.L. Dailey a, H.C. Yalcin a, S.N. Ghadiali a

a,b,*

Mechanical Engineering and Mechanics, Lehigh University, Rm. 265 Packard Lab, 19 Memorial Drive West, Bethlehem, PA 18015, USA b Bio Engineering Program, Lehigh University, Bethlehem, PA 18015, USA Received 31 May 2006; accepted 20 November 2006 Available online 11 January 2007

Abstract Fluid buildup in the small pulmonary airways can exert injurious stresses on the epithelial cells which line airway walls. Under these conditions, the amount of deformation-induced cell injury may depend on the magnitude of the hydrodynamic stresses and the cells’ mechanical properties. In this study, we present 2D and 3D fluid-structure interaction models of flow-induced cell deformation. We report wall shear stress on the cells and cell membrane strain for a range of cell and membrane stiffness values. Results indicate that more compliant cells experience lower wall shear stresses but higher membrane strains. We correlate our results to experimental studies of cellular injury under airway reopening conditions and validate our computational method by comparison to an analytical solution of flow over a rigid protrusion in a channel. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Lung mechanics; Epithelial cell; Shear stress; Membrane strain; Deformation

1. Introduction The deformation of biological cells exposed to flow is a fluid-structure interaction problem with important clinical implications. For example, small pulmonary airways can fill with fluid during pathological conditions such as pneumonia. The buildup of fluid can lead to a cascade of deteriorating lung function known as acute respiratory distress syndrome (ARDS). Approximately 40% of patients with ARDS die due to lung injury, blood infections, and multi-organ failure [1]. In many cases, ARDS patients require mechanical ventilation. Although ventilation therapy can be life saving, it can also damage lung epithelial cells and exacerbate lung trauma, resulting in a condition *

Corresponding author. Address: Department of Mechanical Engineering and Mechanics, Lehigh University, Rm. 265 Packard Lab, 19 Memorial Drive West, Bethlehem, PA 18015, USA. Tel.: +1 610 758 5168; fax: +1 610 758 6224. E-mail address: [email protected] (S.N. Ghadiali). 0045-7949/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.11.012

called ventilator-induced lung injury [2]. During the ventilation of fluid-filled lungs, the flow of air-liquid interfaces and microbubbles in small pulmonary airways induces shear and normal stresses that may produce significant cellular deformation. Although the bubble-induced stress problem has been investigated for rigid channels [3,4], cell deformation may be critical for understanding the mechanisms of lung injury. The purpose of this study is to examine the relationship between cell deformability and the flow-induced stresses on cells. 2. Model development 2.1. Cell mapping In this study, we present 2D and 3D fluid-structure interaction (FSI) models of epithelial cells (EC) on a channel wall. We define the cell morphology from laser-scanning confocal microscopy (LSCM) images of EC cultured in-vitro. LSCM is a fluorescent microscopy technique that

H.L. Dailey et al. / Computers and Structures 85 (2007) 1066–1071

eliminates out-of-plane light to obtain high resolution images within a given cross-sectional plane. Sequential cross-sectional images obtained from LSCM can be used to visualize the 3D morphology of the epithelial cells. The images obtained from the confocal microscope were imported into the Rhinoceros (Seattle, WA) CAD package and used to generate 3D surface maps of the EC. To create a 2D model from this 3D mapping, we selected a typical cell cross section and scaled the slice to generate a representative cell geometry with length, L = 50 lm, and height, h = 5 lm. 2.2. Fluid-structure interaction (FSI) model setup For the 2D FSI model, we created an array of four representative cells as shown in Fig. 1. The cells are attached to a deformable substrate (Esub = 1E6 dyn/cm2) in a channel with height, H = 250 lm.The steady-state velocity profile at the inlet is specified based on the desired wall shear stress, sfw, for Poiseuille flow in an equivalent flat-walled channel. The flat wall shear values, sfw = 5 and 30 dyn/ cm2 were selected to match the range of conditions observed in experimental studies of bubble-induced cell damage [5]. In the solid domain (i.e. within the cells), the finite element mesh consists of 7-node triangular planestrain elements. Additional 3-node plane-strain isobeam elements along the cell surfaces represent the cell membrane where the membrane thickness is d = 5 nm. A range of the cell elastic modulus, 500 < Ecell < 4500 dyn/cm2, was investigated based on experimental microrheology measurements in cultured human alveolar epithelial cells [6– 8]. The maximum membrane stiffness, Emem = 1E15 dyn/ cm2, was chosen such that changes in Ecell had no effect on the solution. The minimum stiffness, Emem = 1E9 dyn/ cm2, was the value below which changing Emem had no effect on the solution.

1067

2.3. Computational FSI formulation We employ the ADINA 8.2 (Watertown, MA) finite element (FE) package to solve our FSI problem [9,10]. The ADINA solver employs a mixed discretization with an arbitrary Lagrangian–Eularian (ALE) formulation. The deformable walls have a Lagrangian formulation and the fluid domain has an Eularian formulation. The fluid and solid solutions are coupled through continuity boundary conditions for displacement, traction, and velocity at the fluid-structure boundary. Fluid flow in the ADINA model is governed by the incompressible form of the Navier–Stokes equations. ovi ¼0 oxi ovi ovi op o2 v i þ qvj q ¼ þl oxi ot oxj oxj oxj

q

ð1Þ ð2Þ

where q is the fluid density, l is the fluid viscosity, vi is the velocity vector, xi is the position vector, t is time, and p is the fluid pressure. Note that all equations are written in standard indicial notation. Cellular deformation in the solid domain is governed by a standard stress balance relationship and a linear elastic solid constitutive relationship. orsij ¼0 oxj E mE eij þ ekk dij rsij ¼ 1þm ð1 þ mÞð1  2mÞ   1 od i od j where eij ¼ þ 2 oxj oxi

ð3Þ

ð4Þ

Here rsij is the stress tensor, E is the Young’s modulus, m = 0.49 is the Poisson ratio, eij is the strain tensor, ekk = e11 + e22 + e33 is the strain invariant, and di is the solid displacement in a given direction. In the FSI models,

Fig. 1. FSI model setup (2D).

1068

H.L. Dailey et al. / Computers and Structures 85 (2007) 1066–1071

the fluid and solid domains are coupled by satisfying the following three boundary conditions at the interface between the two domains. d fi ¼ d si nj rfij ¼ nj rsij ni v i ¼ ni

where rfij ¼ pdij þ l

od si ot



ovi ovj þ oxj oxi



ð5Þ ð6Þ ð7Þ

where d fi and d si are the fluid and solid displacements at the fluid-structure interface, nj is the interface normal vector component in the j-direction, dij is the Kronecker delta function, and rfij and rsij are components of the fluid and solid stress tensors. The kinematic condition in Eq. (5) enforces displacement continuity between the fluid and solid domains and the dynamic condition in Eq. (6) specifies stress continuity between the fluid and solid domains. Finally, the no-slip fluid velocity continuity condition in Eq. (7) requires equal fluid and solid velocities at the fluid–solid boundary. We use an iterative solution scheme in ADINA for computing the two-way coupling between the fluid and solid domain. First, Eqs. (1) and (2) are solved given the no-slip condition in Eq. (7) for the flow velocities in the fluid domain and fluid tractions, n Æ rf, at the fluid-structure boundary. The value of n Æ rf calculated at the fluid-structure interface is then used as a boundary condition to solve Eqs. (3) and (4) for deformations in the solid domain. These deformations result in a new shape for the fluid domain and the iteration technique checks whether the kinematic and dynamic conditions (Eqs. (5) and (6)) are satisfied. If these conditions are not satisfied, the solver repeats the fluid and solid computations with the updated traction and displacement information. To maintain stability and to speed the approach to convergence, the solver employs displacement and stress relaxation factors when passing boundary information between the fluid and solid domains.

[11]. Fig. 2 shows shear stress amplification due to the presence of cells on the channel wall when the flat-wall shear stress is 30 dyn/cm2. When the membrane and cell stiffness are both low (e.g. Emem = 1E9 dyn/cm2, Ecell = 500 dyn/ cm2), the cell experiences larger deformations and this deformation mitigates the shear stress amplification. As the cell becomes more rigid (increasing Ecell and Emem), cell deformation decreases and the shear stress amplification approaches a constant value of sw, max/sfw = 2.13. Cell and membrane mechanical properties may also affect deformation-induced membrane failure. As a measurement of the potential for cellular injury, we calculate membrane strains along the beam elements used to model the cell membrane. Fig. 3 shows the strain profiles for all four cells when Ecell = 500 dyn/cm2 and Emem = 1E9 dyn/ cm2. Note that the positive (tension) and negative (compression) strain values correspond to the upstream and downstream edges of the cell. Fig. 4 compares the maximum positive and negative membrane strains, emem, max, for a range of membrane stiffness values when sfw = 30 dyn/cm2 and Ecell = 500 and 4500 dyn/cm2. As expected, cell deformation and membrane strain decrease as the cell becomes more rigid (increasing Emem or Ecell). For the stiffest cells, the results converge to a load-independent strain magnitude of approximately 0%. Interestingly, the positive and negative membrane strains are asymmetric, with the membrane experiencing larger tensile (positive) strains than compressive (negative) strains. In the preceding models, we used poiseuille flow conditions to generate physiologically-appropriate wall shear stress levels. The pressure gradient range for these simulations was 0.02 < dP/dz < 0.12 dyn/cm2/lm. However, when a bubble tip passes over cells, it can exert steep pressure gradients in the range 1 < dP/dz < 7 dyn/cm2/lm [5,13]. These higher pressure gradients have been implicated as a primary mechanism of cell injury [13]. Therefore, we conducted a separate set of simulations to investigate cellular deformation due to pressure gradients. In these simulations, a pre-defined pressure gradient load was applied to the apical surface of a single-cell. Fig. 5 shows

3.1. Influence of cell deformability Cellular deformability can play a significant role in determining the normal stress, shear stress and stress gradients experienced by EC. Other investigators have examined the flow-induced stresses on rigid cell-like protrusions from channel walls [4,11,12]. In this study, we extend this concept to include the effects of cell deformability for realistic cell geometries. To investigate the effect of cell deformability on the shear stress experienced by the cells, we calculated a shear stress amplification factor, sw, max/sfw, where sw, max is the maximum wall shear stress value and sfw is the flat wall shear stress for the same inlet conditions. In these models, the maximum shear stress occurs at the apex of the cell, which agrees with the results of Gaver and Kute

Shear Stress Amplification, τw,max/τ fw

3. Results (2D) 2.15

2.10

Emem [dyn/cm2] 2.05

1E15 1E12 1E9

2.00 500

1500

2500

Cell Stiffness, Ecell

3500

4500

[dyn /cm2]

Fig. 2. Cell deformability (low Ecell and Emem) mitigates shear stress amplification.

H.L. Dailey et al. / Computers and Structures 85 (2007) 1066–1071

Membrane Strain, εmem [ %]

Cell Height, h [μm]

6

1069

Compression

4 2 Tension

0 0

20

40

60

80

100

120

140

160

20% Tension

10% 0% -10%

Compression

-20% 0

20

40

60

80

100

120

140

160

Axial coordinate, z [μ m] Fig. 3. Cells deform in the direction of flow. Each cell experiences positive (compressive) and negative (tensile) membrane strains. Positive membrane strains correspond to the upstream side of the cells. Data presented for Ecell = 500 dyn/cm2, Emem = 1E9 dyn/cm2.

Ecell = 500

Ecell = 4500

the cellular deformations and membrane strains resulting from two different pressure gradient loads when Ecell = 500 and 4500 dyn/cm2. Note that the compression/tension asymmetry is more significant than in the shear flow case (Fig. 4), especially for higher pressure gradients (dP/dz = 7 dyn/cm2/lm).

Max Membrane Strain, ε mem,max [%]

20%

Tension

15% 10% 5%

3.2. Comparison to experiments

0% -5% -10%

Compression -15% 1E+09 1E+10 1E+11 1E+12 1E+13 1E+14 1E+15

Membrane Stiffness, log(Emem) [dyn /cm2] Fig. 4. Both positive and negative membrane strain magnitudes decrease with increasing Ecell and Emem. Membrane strain is asymmetric, with larger positive (tensile) strains than negative (compressive) strains. Results presented for sfw = 30 dyn/cm2.

Recent experiments by Ghadiali and colleagues [14,15] indicate that cell morphology may influence the risk of cellular injury or death during airway reopening. In these studies, alveolar epithelial cells were cultured to achieve subconfluent and confluent monolayers. The subconfluent cultures were characterized by taller isolated cells while confluent cells were shorter and more tightly packed. Microbubbles were passed over the cells to mimic airway reopening conditions and the cells were stained with fluorescent markers to distinguish between living and dead cells. For constant hydrodynamic conditions (i.e. bubble velocity) these studies indicate that cells in subconfluent

Height, h [μm]

6 4 2 0 0

10

20

30

40

Length, z [μm]

dP/dz = 7 dyn/cm2/µm dP/dz = 1 dyn/cm2/µm

Max Membrane Strain, ε mem,max [%]

40%

Deformed Shape

30% 20%

dP/dz = 7 dyn /cm2/µm dP/dz = 1 dyn /cm2/µm Tension

10% 0% -10%

Compression -20% 1E+09 1E+10 1E+11 1E+12 1E+13 1E+14 1E+15

Membrane Stiffness, log(Emem ) [dyn/cm2] Fig. 5. Directly-applied pressure gradients produce significant membrane strain asymmetry in a solid-only cell model. The deformed cell shapes and membrane strain profiles are given for dP/dz = 1 dyn/cm2/lm and 7 dyn/cm2/lm. Results presented for Ecell = 500 dyn/cm2.

1070

H.L. Dailey et al. / Computers and Structures 85 (2007) 1066–1071

Morphology Cases A B C

Shear Stress Amplification, τ w,max / τ fw

monolayers exhibit a higher rate of death than the cells in confluent monolayers. To compare with these experiments, we used our computational models to calculate the shear stress amplification factor for three geometry cases: tall isolated cells (‘‘subconfluent’’ case), short isolated cells, and short neighboring cells (‘‘confluent’’ case). In all three cases, we used rigid cells (Ecell = 4500 dyn/cm2, Emem = 1E15 dyn/cm2). Fig. 6 shows that the subconfluent morphology experiences the highest shear stress amplification, while the confluent morphology has the lowest shear stress. These results are consistent with the experimental observation of higher death rates in subconfluent cells. Specifically, the taller,

more isolated subconfluent cells experience larger shear stress amplifications and these larger shear stresses may be responsible for the higher death rates. 3.3. Comparison to lubrication theory Although no analytical solution exists for a deformable cell in channel flow, we compared the results of our computational methods to the lubrication theory predictions offered by Gaver and Kute [11] for Stokes flow over a rigid semi-circular cell adhering to a channel wall. The channel has height H, length L, and a semi-circular obstruction of radius R. The flow field can be described in dimensionless form by the following relations: uðx; yÞ ¼ AðxÞy 2 þ BðxÞy þ CðxÞ

2.75

0

2.50

vðx; yÞ ¼ 

2.25

where

A ðxÞ 3 B ðxÞ 2 y  y þ C 0 ðxÞy þ DðxÞ 3 2

4 dP 4 dP ; BðxÞ ¼  2 ðh1 ðxÞ þ h2 ðxÞÞ; b2 dx b dx 0 0 A ðxÞ B ðxÞ h2 ðxÞ3 þ h2 ðxÞ2 þ C 0 ðxÞh2 ðxÞ DðxÞ ¼ 3 2 AðxÞ ¼

2.00 1.75 1.50

A

B

C

Fig. 6. Shear stress amplification for three cell morphology cases: (A) tall isolated ‘‘subconfluent’’ cells, (B) short isolated cells, and (C) short neighboring ‘‘subconfluent’’ cells. Subconfluent cells (case A) experience higher shear stresses than confluent cells (case C). Results presented for rigid cells (Ecell = 4500 dyn/cm2, Emem = 1E15 dyn/cm2).

ð8Þ

0

CðxÞ ¼

ð9Þ

4 dP h1 ðxÞh2 ðxÞ b2 dx

ð10Þ

Here b = H/L is the channel aspect ratio and h1(x) and h2(x) describe the upper and lower channel surfaces. The local pressure gradient, dP/dx can be written as follows: " #1 Z x¼1 dP dx 3 ¼  ðh1 ðxÞ  h2 ðxÞÞ ð11Þ 3 dx x¼0 ðh1 ðxÞ  h2 ðxÞÞ

Fig. 7. (a) 3D cell model, (b) nodal displacements when Ecell = 500 dyn/cm2, sfw = 30 dyn/cm2, (c) and (d) wall stress maps for sfw = 5 dyn/cm2 and 30 dyn/cm2. Shear stress and displacement maxima occur at the upstream face and downstream bump.

H.L. Dailey et al. / Computers and Structures 85 (2007) 1066–1071

We calculated shear stress on cell surface from the lubrication theory formulation and compared the results to ADINA models of the fluid domain. We considered cases in which 0.6 < R/H < 0.85, where R/H is the ratio of cell height to channel height. Note that lubrication theory assumptions are only valid when the cell-channel gap is small or R/H > 0.5. We calculated the shear stress amplification factor, sw, max/sw, fw, and for all cases, the average percent difference between the ADINA models and the lubrication theory was 3.2%, which indicates good correspondence between the computational results and analytical solution.

1071

deformation for 3D models of confluent and subconfluent cells in response to the transient stress waves generated during airway reopening. These models may help elucidate the link between hydrodynamic forces and cell deformation responses in the progression of ARDS and VILI. Acknowledgements H.L.D. is a National Science Foundation Graduate Research Fellow. S.N.G. is a Parker B. Francis Fellow in Pulmonary Research. This work was partially support by a Beginning Grant-in-Aid from the American Heart Association.

4. Results (3D) References We have also developed 3D FSI models of cell deformation in shear flow. Unlike the 2D model, this 3D cell model does not include any elements that model the plasma membrane. Fig. 7a shows the 3D cell model and Fig. 7b shows the displacement magnitudes generated when sfw = 30 dyn/ cm2 and Ecell = 500 dyn/cm2. The shear stress amplification is similar to the results given for the 2D case with sw, max/ sfw = 1.86 and 2.01 when sfw is 5 and 30 dyn/cm2. Shear stress maps for the whole cell are given in Figs. 7c (sw = 5 dyn/cm2) and Figs. 7d (sw = 30 dyn/cm2). The highest shear stresses occur on the upstream face and the downstream bump (circled in Fig. 7a). As expected, the regions of high shear correspond to the largest displacements in the solid model. 5. Conclusions We have developed 2D and 3D fluid-structure interaction models of deformable epithelial cells exposed to flow. The results indicate that cell deformability can reduce shear stress amplification, especially for compliant cells at high flow rates (Ecell = 500 dyn/cm2, Emem = 1E9 dyn/cm2, sfw = 30 dyn/cm2). The results also indicate that increasing either cell or membrane stiffness can decrease membrane strain. Cells exposed to flow experience asymmetrical membrane tension and compression, with the positive (tensile) strains being larger than the negative (compressive) strains. We also investigated how the large pressure gradients experienced by cells during microbubble flows influence deformation. Results indicate that high pressure gradients increase the asymmetry of the membrane strain when cells are compliant (low Ecell and Emem). These observations may be useful for developing pharmacological treatments that combat flow-induced cell injury during ARDS and VILI by altering the cells mechanical properties. Future studies will compare shear stress amplification and cell

[1] Ware LB, Matthay MA. The acute respiratory distress syndrome. New Eng J Med 2000;342(18):1334–49. [2] Ricard JD, Dreyfuss D, Saumon G. Ventilator-induced lung injury. Eur Respir J 2003;22(Suppl. 42):2s–9s. [3] Ghadiali SN, Gaver III DP. An investigation of pulmonary surfactant physicochemical behavior under airway reopening conditions. J Appl Physiol 2000;88(2):493–506. [4] Jacob AM, Gaver III DP. An investigation of the influence of cell topography on epithelial mechanical stresses during pulmonary airway reopening. Phys Fluids 2005;17(3):031502. [5] Kay SS, Bilek AM, Dee KS, Gaver III DP. Pressure gradient, not exposure duration, determines the extent of epithelial cell damage in a model of pulmonary airway reopening. J Appl Physiol 2004;97:269–76. [6] Trepat X, Grabulosa M, Puig F, Maksym GN, Navajas D, Farre R. Viscoelasticity of human alveolar epithelial cells subjected to stretch. Am J Physiol-Lung C 2004;287:L1025–34. [7] Berrios JC, Schroeder MA, Hubmayr RD. Mechanical properties of alveolar epithelial cells in culture. J Appl Physiol 2001;91:65–73. [8] Alcaraz J, Buscemi L, Grabulosa M, Trepat Z, Fabry B, Farre R, Navajas D. Microrheology of human lung epithelial cells measured by atomic force microscopy. Biophys J 2003;84:2071–9. [9] Bathe K-J, Zhang H, Ji S. Finite element analysis of fluid flows fully coupled with structural interactions. Comput Struct 1999;72:1–16. [10] Bathe K-J, Zhang H. Finite element developments for general flows with structural interactions. Int J Numer Met Eng 2004;60:213–32. [11] Gaver III DP, Kute SM. A theoretical model study of the influence of fluid stresses on a cell adhering to a microchannel wall. Biophys J 1998;75:721–33. [12] Lu H, Koo LY, Wang WM, Lauffenburger DA, Griffith LG, Jensen KF. Microfluidic shear devices for quantitative analysis of cell adhesion. Anal Chem 2004;76:5257–64. [13] Bilek AM, Dee KC, Gaver III DP. Mechanisms of surface-tensioninduced epithelial cell damage in a model of pulmonary airway reopening. J Appl Physiol 2003;94:770–83. [14] Ghadiali SN, Yalcin HC. Injury and repair of epithelial cells during the reopening of fluid-filled airways. In: Proceedings of the 5th world congress of biomechanics, Munich, Germany; 2006. [15] Yalcin HC, Ghadiali SN. Influence of cellular morphology and mechanics on injury patterns during airway reopening. In: Proceedings of the biomedical engineering society fall meeting, Baltimore, MD; 2005.