Numerical Simulation of Biological Structures Session Organizer: Wilkins AQUINO (Cornell University) Keynote Lecture Computational modeling of glucose distribution in hollow fiber membrane bioreactors G.U. UNIKRISHNAN, V.U. UNIKRISHNAN, J.N. REDDY* (Texas A&M University) Characterization of viscoelastic properties of cylindrical vessels using the velocity response produced by an impulsive force Daniel E. ROSARIO, Wilkins AQUINO* (Cornell University) Solid versus membrane finite elements in analysis of the mitral valve: A case study Victorien PROT*, Bjorn SKALLERUD (Norwegian University of Science and Technology) An inverse problem approach for elasticity imaging through vibroacoustics Miguel AGUILO, Wilkins AQUINO* (Cornell University) Proper orthogonal decomposition model reduction for inverse problems in acoustic-structure interaction John C. BRIGHAM, Wilkins AQUINO* (Cornell University) Shell analysis of elliptical guard cells in higher plants: A review J. Robert COOKE*, Richard RAND (Cornell University), Herbert MANG (TU Vienna), Josse DeBAERDEMAEKER (Katholieke Universiteit Leuven), Jae Young LEE (Chonbuk National University) In vivo ultrasound bone property determination through inverse finite element modeling Mija HUBLER*, Wilkins AQUINO, Christopher EARLS (Cornell University) Finite element analyses of palm leaf petiole-sheath junctions in simple bending and twisting and in dynamic (oscillatory) flexure Karl NIKLAS*, J. Robert COOKE (Cornell University), Jae Young LEE (Chonbuk National University) Modeling pipette aspiration of biological membranes Philip BUSKOHL (Cornell University) A new model for nucleation in two-phase lipid bilayer membrane vesicles Sanjay DHARMAVARAM*, Timothy HEALEY (Cornell University) For multiple-author papers: Contact author designated by * Presenting author designated by underscore Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) Computational modeling of glucose distribution in hollow fiber membrane bioreactors G.U. UNNIKRISHNAN, V.U. UNNIKRISHNAN and J.N. REDDY* *Advanced Computational Mechanics Laboratory Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, USA jnreddy@tamu.edu Abstract Bioreactors play an important role in tissue engineering, as they simulate physiological environment required for the development of tissue substitutes. The nourishments supplied by the capillaries in a bioreactor reach the cells through diffusion and convection from the scaffold boundaries. This limits the size of cultured tissue, thus making it clinically nonusable. Hollow fiber membrane bioreactors (HFMB) however are capable of overcoming the size restrictions induced due to mass transport limitations. In this work, a fluid-biphasic interface model is developed to analyze the glucose transport in a hollow fiber membrane bioreactor. The scaffold-fiber-fluid domain is treated as a single biphasic continuum and the interface is coupled automatically into the finite element model. A biphasic finite element representation of the bioreactor provides a fluid-scaffold coupling, thus making it possible to study fluid flow in lumen and fiber and scaffold domains. The novelty of this biphasic-fluid interface model is that no additional boundary conditions at the interface are necessary. The influence of bioreactor parameters like porosity, stiffness, and permeability of the scaffold, and also inlet flow velocity are studied in this work, as these optimum parameters help in preventing a nonhomogenous distribution of nutrients. The developed finite element model would thus be a better tool to aid in the design and analysis of HFMB and would help in the optimization of the bioreactor working conditions. Acknowledgement. The research is carried out with the support of the Oscar S. Wyatt Endowed Chair account at Texas A&M University. The authors are grateful for the support. Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) Characterization of viscoelastic properties of cylindrical vessels using the velocity response produced by an impulsive force Daniel E. ROSARIO and Wilkins AQUINO* *Cornell University, School of Civil and Environmental Engineering 313 Hollister Hall, Ithaca NY 14853, United States wa27@cornell.edu Abstract An inverse methodology for non destructive characterization of viscoelastic material properties will be presented in this work. This study shows the characterization of viscoelastic material properties of an isotropic cylindrical vessel using the surface velocity response produced by an impulsive force. Sensitivity analyses and an inverse problem are shown to quantify the potential for material characterization through the use of the velocity response obtained from an experimental setup. The experiment consisted of a silicone cylinder filled with water embedded in gelatin. An impulsive force was generated on the surface of the cylinder using the radiation force of ultrasound, and the resulting velocity response was measured. The viscoelastic material parameters considered in the analyses were the equivalent modulus and the mechanical loss. It is shown that this methodology have great promises for the estimation of viscoelastic material properties using the cylinder surface velocity response. 1. Introduction Characterization of mechanical properties of biological structures has been shown to pose particular difficulty due to the rate dependent material behavior (see Fung[4]). In order to characterize these properties, it is necessary to perform tests which contain sufficient information relating to this rate dependent behavior without disturbing the current state of the biological structure. Fortunately, several methods which were initially developed for medical imaging have been shown to provide the information necessary to characterize biological material properties. Some examples are the work by Catheline et al.[2] and Zhang et al.[5]. In this study, an inverse methodology is applied to the characterization of viscoelastic properties of an isotropic cylindrical vessel using the surface velocity response produced by an impulsive force. The experiment was conducted by the Ultrasound Research Laboratory at the Mayo Clinic. In this experiment, an impulsive ultrasound radiation force was applied to a silicone rubber cylinder and the resulting propagating velocity wave was measured using a laser vibrometer. An example is shown in where the velocity response obtained from the experiment is used to quantify the cylindrical vessel viscoelastic material parameters (i.e. mechanical loss and equivalent modulus) through an inverse problem approach. 2. Experiment A silicone cylindrical vessel was filled with water and embedded in a gelatin block to simulate in-vivo characterization of tissue. An impulsive radiation force of 500 Hz was generated by the excitation transducer at a point on the cylinder surface. The excitation transducer is moved along the longitudinal axis of the cylinder to produce an impact at varying distances from the measurement point, which is always fixed. These distances vary from 2 mm to 20 mm in intervals of 2 mm. A laser vibrometer was used to measure the velocity response at the measurement point on the cylinder surface. A schematic representation of the experimental setup is shown in Figure 1. 1 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca The cylindrical vessel had an outer diameter of 5 mm, an inner diameter of 3 mm, and a length of 25 cm. The gelatin block had a width of 6.2 cm, a height of 6 cm, and a length of 25 cm. The density of the cylinder and the gelatin was 1130 kg/m3 and 1000 kg/m3, respectively. Cylinder Gelatin Block Measurement Point Excitation Transducer Laser Vibrometer Figure 1: Schematic of the experimental setup. 3. Foward problem A finite element fluid-structure interaction formulation presented in Cook et al.[3] was used to model the dynamic response of the experimental setup presented in Section 2 for a given set of material parameters. For this work, the cylinder and gelatin material behavior was taken to be linear isotropic viscoelastic and nearly incompressible (Poisson’s ratio of 0.48). The strains and displacements in both solids, the cylinder and gelatin, were assumed to be small due to the force excitation. The stress-strain constitutive relationship for the linear viscoelastic solids was considered to obey a Prony series relaxation law in the time domain and is given by ∑ ∫σ ( G x,t ) = E∞ε ( G x, t ) + N ⎡t ⎢ m=1 ⎢⎣ 0 −(t−ξ ) Eme τi ∂ε ( xG,ξ ∂ξ ) ⎤ dξ ⎥ ⎥⎦ , (1) where G x is the spatial position vector, t is the total time, σ ( G x, t ) is the stress tensor, ε ( G x,t ) is the strain tensor, E∞ is the long- term elastic moduli, Em are the viscous moduli, τ m are the relaxation times, and N is the number of terms for the Prony series. The fluid was assumed to be compressible with no net flow, small pressure amplitudes, and negligible viscosity (i.e. acoustic medium). In addition, the only excitation in the fluid was assumed to be due to the motion of the cylinder. 4. Inverse problem The inverse problem considered was to determine the viscoelastic material parameters of a cylindrical vessel using the measured velocity response obtained from the experiment previously describe. To obtain a solution, this inverse problem was cast as a constrained optimization problem that minimize an error functional which quantifies the distance between the experimental response and the response computed using the Finite Element method for a given set of parameter estimates. For this work, the error functional was defined as J (αG ) = vexp (t ) − vsim (αG,t ) , L∞ (2) 2 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca where αG is the vector of unknown parameters, vexp (t ) and vsim (αG,t ) are the experimental and computed velocity responses, respectively, and L∞ is the L∞ -norm. The optimization problem was solved using the Surrogate-Model Accelerated Random Search (SMARS) algorithm introduced by Brigham and Aquino[1]. This optimization algorithm combine an iterative application of a random search algorithm and a surrogate-model optimization approach to minimize the error functional defined in eq. (2). 5. Example One example was considered to show the potential of the proposed approach to inversely characterize the viscoelastic material behavior of a cylindrical vessel using the velocity response obtained from the experiment showed in Section 2. For this example, only the velocity response produced by the impulsive force applied at 20mm from the measurement point was used. Due to the linearity in the representative numerical models, the velocity response could be normalized to remove the effect of any uncertainty in the load magnitude, and as such, the amplitude of the impact force is irrelevant. In addition, the responses for the experiment and the optimization simulations were filtered through a low-pass filter with a cutoff frequency of 800 Hz. The parameters to be identified by the inverse problem were the long- term elastic moduli ( E∞ ), the viscous moduli ( Em ), and the relaxation times (τ m ). Characterization of viscoelastic material behavior using these material parameters can lead to problems of nonuniqueness in the solution of the optimization algorithm since different combinations of the parameters or different representations of rheological models can outcome in the same material behavior. For this reason, in the example, the material behavior is represented using the mechanical loss ( tan[δ (w)] ) and the equivalent modulus ( Eeq (w) ) of the material over the frequency domain. These quantities can represent in a unique way the viscoelastic material behavior and are derived from the complex Young’s modulus which is given by ∑E* (w) = E∞ + N Em m=1 ⎛ ⎜⎜⎝ w2τ 2 m + iwτ m 1 + w2τ 2 m ⎞ ⎟⎟⎠ , (3) where w is the circular frequency. The equivalent modulus ( Eeq (w) ) is the magnitude of the complex modulus and the mechanical loss ( tan[δ (w)] ), which is associated with viscous damping of the material, is defined as the ratio of the imaginary to the real part of the complex modulus. Several forward simulations were performed before solving the inverse problem to obtain a coarse sensitivity estimate of the velocity response to changes in the viscoelastic properties of the gelatin and cylinder. If the velocity response were found to have a weak sensitivity to the viscoelastic properties of a particular material (i.e. gelatin or cylinder), then the solution of the inverse problem would likely be non-unique. In other words, a viscoelastic material with weak sensitivity will not be identifiable in the inverse problem and the viscoelastic parameters can be assumed to be known since they do not affect the velocity response. 5.1 Results The sensitivity results for the viscoelastic properties of the gelatin block show a weak sensitivity to the velocity response. In contrast, the sensitivity results for the viscoelastic properties of the cylindrical vessel demonstrate a great sensitivity to the velocity response. Therefore, only the viscoelastic properties of the cylinder were considered for the inverse problem. The results shown were using four terms in the Prony series relaxation law ( N = 4 ) which means that the optimization algorithm had to search for nine viscoelastic parameters. In Figure 2, the velocity response for different excitation points is shown for the experiment and the optimization algorithm solution, it is clear that the optimization algorithm found a solution that matched the experimental response with high accuracy. 3 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca Experiment Optimization algorithm Figure 2: Normalized velocity response for varying distance of excitation point from measurement point. Acknowledgement This work was partially supported by The National Institute of Biomedical Imaging and Bioengineering and Cornell University. References [1] Brigham JC and Aquino W. Surrugate-Model Accelerated Random Search algorithm for global optimization with applications to inverse material identification. Computer Methods in Applied Mechanics and Engineering 2007; 196: 4561-4576. [2] Catheline S, Gennisson JL, Delon G, Fink M, Sinkus R, Abouelkaram S, and Culioli J. Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: An inverse problem approach. Journal of the Acoustical Society of America 2004; 116: 3734-3741. [3] Cook RD, Malkus DS, Plesha ME, and Witt RJ, Concepts and Applications of Finite Element Analysis,(Fourth edn). John Wiley & Sons, 2002. [4] Fung YC, Biomechanics: mechanical properties of living tissues,(2nd edn). Springer-Verlag, 1993. [5] Zhang XM, Kinnick RR, Fatemi M, and Greenleaf JF. Noninvasive method for estimation of complex elastic modulus of arterial vessels. Ieee Transactions on Ultrasonics Ferroelectrics and Frequency Control 2005; 52: 642-652. 4 Proceedings of the 6th International Conferense on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) Solid versus membrane finite elements in analysis of the mitral valve: a case study Victorien PROT∗, Bjørn SKALLERUD∗ ∗Department of Structural Engineering, Norwegian University of Science and Technology Richard Birkelands vei 1a, 7041 Trondheim, Norway victorien.prot@ntnu.no Abstract The present work addresses an incompressible transversely isotropic hyperelastic material used in three–dimensional solid finite element analysis of a porcine mitral valve. The aims of the study are to investigate the influence of the layered structure of the mitral valve on stress distributions and global responses. First, the implementation of the material model is described. Then, three three-dimensional finite element analyses of the mitral valve, using three different collagen layer arrangements, are conducted. Our simulations show that when the leaflets are arranged in a two layer structure with the collagen fibres on the ventricular side, the stress in the fibre direction through the thickness in the central part of the anterior leaflet is homogenized and the peak stress is reduced. Finally, the global response of the mitral valve finite element model is compared to corresponding echocardiographic measurements and a simulation using membrane elements. It appears that our model bulges too much in the left atrium. This may be due to the fact that there is evidence of active muscle fibres in some parts of the anterior leaflet, whereas our constitutive modelling is based on passive material. 1 Introduction The mitral valve is a thin walled complex connective tissue structure located between the left atrium and left ventricle of the heart, preventing the blood from flowing back into the atrium when the ventricle contracts. The mitral apparatus consists of two leaflets (anterior and posterior) attached to the annulus and the chordae tendinae. The chordae are further attached to the papillary muscles. According to previous studies (Kunzelman et al. [1], Grande-Allen et al. [2], Jensen et al. [3]) the mitral leaflets can be considered as three–layered laminated structures: the atrialis on the atrium side, the ventricularis on the ventricular side and the inner fibrosa layer, with each layer having different mechanical properties. The present study employs a three dimensional finite element model of the mitral apparatus using solid elements in order to take into account different material properties for each layer. With the membrane approach these finer structural details are not accounted for. Earlier studies have presented general anisotropic hyperelasticity formulations including incompressibility (Weiss et al. [4], Ru¨ter and Stein [5], Holzapfel et al. [6]). In this work we develop a transversely isotropic hyperelastic material model for solid elements using an uncoupled strain energy function based on the model described by Prot et al. [7] for the mitral leaflets. In Section 2 we present the constitutive equations of a transversely isotropic hyperelastic material model for solid elements using a penalty method to treat incompressibility. Finite element analyses of the mitral valve are conducted in order to investigate the influence of the collagen structure on the mitral valve response and particularly on the resistance of leaflets to bending. The global response from simulations is compared to corresponding ultrasound recordings. The comparison illustrates some shortcomings of the simulations. 1 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca 2 Hyperelastic transversely isotropic material model for the mitral leaflets 2.1 Kinematics Let Ω0 and Ω be the reference and current configurations, respectively. The deformation map ϕ(X ) : Ω0 → R3 transforms a material point X ∈ Ω0 into the related current position x = ϕ(X ) ∈ Ω. Hence, the deformation gradient F is defined as F = ∂ϕ(X )/∂X = ∂x/∂X , with the volume ratio J = detF > 0 (J = 1 for an incompressible material). We consider the multiplicative decomposition of the deformation gradient F , i.e. F = (J1/31)F¯ and C = (J2/31)C¯ . F¯ and C¯ = F¯T F¯ are the modified deformation gradient and modified right Cauchy-Green tensor. Note that det(F¯ ) = 1 and det(C¯ ) = det(F¯ )2 = 1. The modified left Cauchy-Green tensor then reads B¯ = F¯ F¯T . We assume that the only anisotropic property arises from a fibre family embedded in the continuum and that the direction of the fibre at point X in the reference configuration Ω0 is defined by a unit vector a0(X). During deformation this fibre moves with the material points of the continuum body and arrives at the deformed configuration Ω. We define the vector a¯ = F¯ a0. 2.2 Strain energy function In order to describe the anisotropic hyperelastic response of the mitral valve leaflets, we employ the following form for the strain energy function Ψ: Ψ(I¯1, I¯4, J) = c0 ec1(I¯1−3)2+c2(I¯4−1)2 − 1 + κ(J − 1)2 , Ψ¯ (I¯1,I¯4): isochoric part U(J): volumetric part (1) fbwuyhneucrtseiioncng0a,ancld1a,rI¯g1ce2=vaatrrelCu¯me=aotfterrκB¯i.aalΨn¯pdairsI¯a4am=neatC¯eltre:sr,an0κa⊗tiisvaea0s.ptTroashiinetiviedeneepareginsyatflotuyanpcpptairrooanxmitmeoteathrteeatnohdnee(mJpar−toep1roi)as2leaidssbksylniogMwhatnlyy-aNcsoetmwhepmrpeaesnnsiaablntldye Yin [8]. 2.3 Stress and elasticity tensors The material model was implemented in the finite element program ABAQUS using the user subroutine UMAT. The implementation of the six Cauchy stress components and the tangent stiffness is required. Now we provide the expressions of the Cauchy stress and the spatial elasticity tensors. The Cauchy stress tensor σ is expressed as follows: σ = 2κ(J − 1)1 + 1 J devσ¯ , where σ¯ = 2ψ1B¯ + 2ψ4a¯ ⊗ a¯ and dev[·] = (I − 1 3 1 ⊗ 1) : (·), ψi = ∂Ψ¯ ∂I¯i , i = 1, 4. (2) (3) In tensor notation the spatial elasticity tensor C is expressed as: C = Cvol + Ciso, Cvol = −4κ(J − 1)I + 4κ(J − 1 2 )1 ⊗ 1, JCiso = 4ψ11devB¯ ⊗ devB¯ + 4ψ14 devB¯ ⊗ dev(a¯ ⊗ a¯ ) + dev(a¯ ⊗ a¯ ) ⊗ devB¯ +4ψ44dev(a¯ ⊗ a¯ ) ⊗ dev(a¯ ⊗ a¯ ) − 2 3 (1 ⊗ σ¯ + σ¯ ⊗ 1) + 2 9 (trσ¯ )1 ⊗ 1 + 2 3 trσ¯ I. ψi j = ∂2Ψ¯ ∂I¯i∂I¯j , i, j = 1, 4 (4) (5) (6) (7) 2 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca 3 Mitral valve finite element simulations In this section a three dimensional finite element model of a porcine mitral valve is presented. The closure of the valve between early systole (p = 0) and peak systole (p = 120 mmHg) is considered. The geometry of the model shown in Figure 1 and the modelling of the chordae tendinae is described by Prot et al. [9]. The chordae were modeled as an incompressible hyperelastic material. The leaflets were meshed with eight Region 1 Annulus Anterior leaflet 23 mm Chordae tendinae Papillary muscles Posterior leaflet 11 mm Figure 1: Initial geometry of the mitral valve at beginning of systole. noded brick hybrid elements and the chordae with truss elements. The mesh in the initial configuration is shown in Figure 1. The displacements of two nodes representing the papillary muscles were constrained. The annulus was assumed to be fixed, with translations constrained but not the rotations. In order to capture coaptation (i.e., apposition of the posterior and anterior leaflets) a contact condition was prescribed on the atrium surface of the leaflets. The blood pressure measured in the left ventricle of the pig between beginning of systole and the maximum pressure in the ventricle was applied on the ventricular surface of the leaflets as load history. Four analyses were performed: the first one using membrane elements [9] with the same material properties through the thickness of the leaflet, in the second one, the leaflets were modeled with two layers of solid elements (0.5 mm each) using the same material properties, the third one with two layers with solid elements (0.5 mm each), the ventricular layer containing the collagen and the atrial layer without collagen and the fourth one with three lFaiygeurrseo2f(sao) liildluesltermateens ttsh(e13stmremssedaicshtr)ibthuetiionnnethr rloauygehr containing the collagen. All simulations were quasi–static. the thickness of the anterior leaflet at about one third of the distance between the annulus and the coaptation line (region 1 on Figure 1). As can be seen in Figure 2(a), the stresses in the atrium surface and in the ventricular surface are reduced and increased, respectively, with the leaflet composed of two different layers compared to the case with constant material properties through the thickness. In region 1, the model using two different layers homogenizes the stress distribution through the thickness in the fibre direction. The analyses show that the anterior leaflet experiences greater stresses on the atrium side compared to a membrane approach, but no compressive stress on the ventricular side. In Figure2(b) the solid line representing the deformed shape of the leaflet obtained from the finite element analysis with the same material properties in the solid elements layers is compared with the ultrasound measurements carried out on the pig for several pressure levels. As can be seen, the anterior leaflet bulges too much into the left atrium compared to the echocardiographic measurement. 4 Conclusion The present study suggests that the use of membrane-shell finite elements to model the leaflets may be sufficient to capture the global response of the mitral valve. However, for finer details such as the stress distribution through the thickness of the leaflets solid finite elements are needed. Our results show that the stress in the leaflets are affected by using different material properties in the layers. We show that using a passive transversely isotropic hyperelastic material model for the leaflets the valve bulges too much in the left atrium. We conclude that the flat shape of the 3 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca atrium side one layer two layers three layers membrane elements 0 A20 40 LVOT 60 80 LA P LV thickness ventricle side 0 Cauchy stress in the fiber direction 100 0 20 40 60 80 (a) Cauchy stress in the fibre direction, region 1 (see, Figure 1) (b) peak systole (120 mmHg) Figure 2: (a) Stress distribution along the fiber direction through the thickness of the anterior leaflet in region 1, (b) Comparison with ultrasound measurements at 120 mmHg: the solid dashed lines represent the solid and membrane element simulations, respectively. (LVOT = Left Ventricular Outflow Tract, LA = left atrium, LV = left ventricle, A = anterior, P = posterior) leaflets shown by the ultrasound recording may be due to the presence of active muscle fibres. Hence, for improved material modelling, the account of muscle fibres need further studies and layer–specific material parameters need also to be established. REFERENCES [1] Kunzelman KS, Cochran RP, Murphree SS, Ring WS, Verrier ED, Eberhart RC. Differential collagen distribution in the mitral valve and its influence on biomechanical behavior. J. Heart Valve Dis. 1993; 2:236-244. [2] Grande-Allen KJ, Calabro A, Gupta V, Wight TN, Hascall VC, Vesely I. Glycosaminoglycans and proteoglycans in normal mitral valve leaflets and chordae: association with regions of tensile and compressive loading. Glycobiology 2004 14(7):621–633. [3] Jensen A, Baandrup U, Hasenkam J, Kundu T, Jørgensen C. Distribution of the microelastic properties within the human anterior mitral leaflet. Ultrasound in Medicine & Biology 2006; 2(12):1943–1948, 2006. [4] Weiss JA, Maker BN, and Govindjee S. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput. Meth. Appl. Mech. Eng. 1996; 135:107–128. [5] Ru¨ter M and Stein E. Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity. Comput. Meth. Appl. Mech. Eng. 2000; 190(5): 519–541. [6] Holzapfel GA, Gasser TC, Ogden RW. A new constitutive framework for arterial wall mechanics and a comparative study of material models. Journal of Elasticity 2000; 61:1-48. [7] Prot V, Skallerud B, Holzapfel GA. Transversely isotropic membrane shells with application to mitral valve mechanics. Constitutive modeling and finite element implementation. Int. J. Numer. Meth. Eng. 2007; 71(8), 987–1008. [8] May-Newman K and Yin FCP. A constitutive law for mitral valve tissue. J. Biomech. Eng. 1998; 120:38–47. [9] Prot V, Haaverstad R, Skallerud B. Finite element analysis of the mitral apparatus: annulus shape effect and chordal force distribution. J. Biomechanics and Modeling in Mechanobiology. Accepted, 2007. 4 Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) An inverse problem approach for elasticity imaging through vibroacoustics Miguel AGUILO, Wilkins AQUINO* *Cornell University, School of Civil and Environmental Engineering 313 Hollister Hall, Ithaca NY 14853-3501, United States wa27@cornell.edu Abstract A new methodology for estimating the spatial distribution of elastic moduli using the steady-state dynamic response of solids immersed in fluids is presented. The direct problem consists of a coupled acoustic-structure interaction boundary value problem solved in the frequency domain using the finite element method. The ensuing acoustic pressure and velocity fields are used to inversely estimate the spatial distribution of Young’s modulus. This work proposes the use of Gaussian Radial Basis Functions (GRBF) to represent the spatial variation of elastic moduli in soft tissue with localized stiff regions. In this work, it is shown that GRBF posses the advantage of representing smooth functions with quasi-compact support and can efficiently represent elastic moduli distributions such as those that occur in soft biological tissue in the presence of tumors. The inverse problem is then cast as an optimization problem in which the objective function is defined as a measure of discrepancy between an experimental response and a finite element model. Non-gradient based algorithms in combination with a divide and conquer strategy are used to solve the resulting optimization problem. The feasibility of the proposed approach is demonstrated through a series of numerical and a physical experiment. 1. Introduction Noninvasive techniques to characterize material properties are of particular interest in many science and engineering fields. These techniques are of great significance in the medical field since material properties provide valuable insight about the onset and progression of certain diseases. In recent years, researchers have tried to develop techniques that could quantitatively reconstruct the spatial distribution of mechanical properties in soft tissue; the fact that the elastic moduli of human tissue can be used as a metric to differentiate between malignant tissue and healthy tissue has served as their primary motivation. Different methodologies have been employed to image the distribution of elastic moduli in soft tissue. The main premise in determining the mechanical properties of soft tissues is that current medical imaging techniques such as ultrasound and magnetic resonance imaging can provide information about complex deformation processes in human tissue. If information is available about the tractions acting on a given region, it is possible to inversely estimate the distribution of the mechanical properties (Barbone et al. [1]). Furthermore, medical imaging techniques not only provide a means to measure deformation, but some techniques can also provide information about material properties. One such technique is vibroacoustography (VA) (Fatemi et al. [2]). Vibroacoustography VA is based on tissue vibration induced by the radiation force of ultrasound (US). In this method, the radiation force (a nonlinear, second-order phenomenon in wave theory) of ultrasound is used to vibrate tissue at low (kHz) frequencies at the focal point of the US beam. The resulting vibrations produce an acoustic field that is detected by a sensitive microphone (or hydrophone). Stiffer tissues normally produce a significantly different acoustic field than soft tissues. Conceptually, VA resembles palpation; i.e., it detects tissue response to a force. The force used in VA is highly localized and taps on a point of the object at a time. Thus, this method can provide detail information on tissue mechanics at high resolution that is not available from conventional palpation or US. 1 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca It has been shown that the radiation force of ultrasound and the ensuing acoustic emissions can be used to characterize the viscoelastic properties of solids (Brigham et al. [3]). However, in their work, it was assumed that the boundary of different material regions was known a priori. The present work extends the work of Brigham et al. by using the vibration response of bodies excited with the radiation force of ultrasound in combination with computational inverse problem techniques to estimate spatially-varying elastic properties. The inverse estimation of spatially-varying mechanical properties in this work is cast as an optimization problem in which an error functional that measures the misfit between experimental and computed responses is minimized by searching over a space of admissible functions that describe the distribution of elastic moduli. Different techniques have been proposed in the literature for solving inverse problems in which the unknown is a function. For instance, the adjoint method was used by Oberai and coworkers (Oberai et al. [4]) for solving elasticity imaging problems. This method has the advantage that only two solutions of a forward model are required in each iteration of the minimization process. The disadvantage of the method is that only local minima are guaranteed to be found. An alternate approach to the adjoint method is to approximate the unknown function using a finite dimensional basis and search for the unknown coefficients using global optimization methods. The advantage of this approach is that global minima can then be found irrespective of initial guesses, but with the disadvantage that, depending on the basis that is selected, a large number of parameters may be needed, increasing the computational cost significantly. Despite this, global optimization methods are favored in this work for their robustness and lack of sensitivity to initial guesses. Since a finite dimensional representation will be used to approximate spatially varying elastic moduli, it is very important that an adequate basis be selected. The quality of a basis for the problem at hand can be judged based on the dimensionality of the subspace needed to “satisfactorily” approximate the desired elastic modulus distribution. Intuitively speaking, a good basis should intrinsically carry the nature of the solution. One of the main applications of the present work is the detection and characterization of tumors in breast tissue. Tumors are expected to be stiff, localized regions embedded in a relative homogeneous soft matrix (e.g. soft tissue). It is proposed in this work that this type of property distribution can be suitably represented with Gaussian radial basis functions (GRBF). It is important to point out that the notion of representing tissue stiffness distribution by using quasi-compactly supported functions has also been suggested (Skovoroda et al. [5]). However, in their work GRBF were not used as a finite dimensional basis for inverse elasticity imaging as proposed herein. In addition, a reduction of the number of parameters needed is accomplish by using a divide and conquer strategy (DCS). The DCS consist in dividing the spatial domain in halves at the end of each optimization run. The optimal solution of the previous run is used as reference, allowing the generation of unknown coefficients around a feasible solution. 2. Forward problem The acoustic-structure response of the systems studied in this work can be described by a coupled system of partial differential equations derived from the conservation of linear momentum. The governing equation for the solid behavior is given by ∇ ⋅ σ ( G x, t ) = ρuG( G x, t ) on Ω. (1) For a compressible fluid with no net flow, negligible viscosity, and small pressure amplitudes, the governing PDE from conservation of linear momentum is given by In the above equations G x is the ∇p( K x, t ) + ρ f position vector, uK f t K (x, t ) = 0 represents on time, Φ σ . ( G x, t ) is the stress tensor, ρ (2) is the mass density (assumed constant in this work), and uG( G x, t ) is the acceleration vector, Ω is the solid domain, K p(x,t) is unf the ( xK, scalar acoustic fluid pressure in excess of t) is the acceleration of the fluid boundary in hydrostatic pressure, ρ f is the mass density the direction of the normal, and Φ is the fluid of the fluid, domain. 2 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca 3. Inverse problem The inverse problem considered was the estimation of spatially-varying elastic properties from the system response due to a radiation force of US. An error functional that measures the misfit between experimental and computed responses was defined as ∫ ∫ ( )G J (E(x)) = ω2 G r Exp ( G x, ω ) − G r FEA ( G x, ω , E ( G x), ...) 2 dxdΩ , Ω ω1 (3) where J is the error functional, E ( G x ) is the spatial distribution of elastic moduli, and ω represents frequency. The quantities G r Exp and G r FEA represent the experimental and computed responses respectively. These quantities can denote pressure or velocity fields depending on the experiment performed. The spatial distribution of elastic moduli was represented as a weighted sum of N radially symmetric radial basis functions augmented by a constant term ( )∑E ( G x ) = Eo + N λiθi GG x − xi 2 , i =1 (4) where θi is the defined as the GRBF θ .= e( )G x − G xi 2 2 ci i (5) In the above formulation Eo is the matrix Young’s modulus, λ is the Young’s modulus at a point in space, and c provides the function quasi-compact support. After parametrization of the spatial distribution of elastic moduli, the error functional is express as ( )m n J ({c},{λ}, Π) = ∑ ∑ ⎡⎣rijExp G xi , ω j − r FEA ({c} , {λ} , Π, G x, ω , ....)⎤⎦ 2 . i =1 j =1 (6) The inverse problem was cast as an optimization problem in which the above error functional is minimized by searching over a space of admissible functions that describe the spatial distribution of elastic moduli minimize J ({c},{λ},Π) {λ},{c}⊂ RN Π ⊂ R3 (7) 4. Results Numerical and physical experiments were used to demonstrate the potential of the proposed technique to be applied in elasticity imaging. A steady-state dynamic FE analysis was performed and the acoustic response at several measurement points was used to solve the inverse problem. In this section, relevant results pertaining to the one inclusion numerical experiment is shown. Case 1 Inclusion Target Eo (Pa) 1.00x109 Mean Eo (Pa) 1.03x109 Standard Deviation 1.01x107 Table 1: Results from 5 optimization runs for the matrix elastic modulus using the pressure field. 3 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca a) b) Figure 2: (a) Schematic of the distribution of elastic moduli obtained from the pressure field for the numerical experiment. (b) Schematic of the distribution of elastic moduli for the target solution. Acknowledgement I want to acknowledge Cornell University and the Coleman Foundation for supporting this work. References [1] Barbone, P.E. and J.C. Bamber, Quantitative elasticity imaging: what can and cannot be inferred from strain images. Physics in Medicine and Biology, 2002; 47: p. 2147-2164. [2] Fatemi, M. and J.F. Greenleaf, Ultrasound-stimulated vibro-acoustic spectrography. Science, 1998; 280(5360): p. 82-85. [3] Brigham, J.C., et al., Inverse estimation of viscoelastic material properties for solids immersed in fluids using vibroacoustic techniques. Journal of Applied Physics, 2007; 101(2): p. 023509-1 - 023509-14. [4] Oberai, A.A., N.H. Gokhale, and G.R. Feijoo, Solution of inverse problems in elasticity imaging using the adjoint method. Inverse Problems, 2003; 19(2): p. 297-313. [5] Skovoroda, A., S. Emelianov, and M. O'Donnell, Tissue elasticity reconstruction based on ultrasonic displacement and strain images. Ieee Transactions of Ultrasonic Ferroelectrics Frequency Control, 1995; 42: p. 747-765. 4 Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) Proper orthogonal decomposition model reduction for inverse problems in acoustic-structure interaction John C. BRIGHAM and Wilkins AQUINO* *Cornell University, School of Civil and Environmental Engineering 313 Hollister Hall, Ithaca NY 14853, United States. wa27@cornell.edu Abstract This presentation will analyze the capabilities of the proper orthogonal decomposition (POD) technique for model reduction in model-updating solution strategies for acoustic-structure interaction inverse problems. POD is a technique which can be used with either experimental or simulated field data to derive a reduced order set of basis functions capable of being used in a numerical representation of the system (e.g. finite element analysis). These basis functions are the optimal set, in an average sense, for representing the field data for the given order. Therefore, the basis is expected to produce an accurate and efficient numerical representation of the system, provided sufficient information exists in the field data. Strategies will be presented for generating the field data which provides the most suitable information for the reduced order model for the solution of inverse acoustic-structure interaction problems. Then the proposed methodology will be illustrated through a numerical and a laboratory experiment, in which a vibroacoustic testing approach is applied to evaluate the viscoelastic behavior of structures. Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) Shell analysis of elliptical guard cells in higher plants: a review J. Robert COOKE*, Richard H. RAND, Herbert A. MANG, Josse G. DeBAERDEMAEKER, Jae Young LEE *Professor Emeritus, Department of Biological and Environmental Engineering 214 Riley-Robb Hall, Cornell University, Ithaca, NY 14853 USA jrc7@cornell.edu and j_robert_cooke@mac.com Abstract This presentation reviews finite element shell analyses (linear and nonlinear, isotropic and anisotropic) of microscopic stomata modeled as doubly-elliptic toroidal shells. Stomata are the regulating valve in the seminal issue of water and carbon dioxide transport in plant biology. The finite element method allows more realistic modeling of complex geometries and material properties that are common in biology. 1. Introduction The physical sciences and engineering are increasingly broadening their traditional focus to include topics formerly within the exclusive domain of the biological sciences. Biological organisms, of course, are governed by and must function within the constraints of physics, so this expansion is direct. As the tools of engineering evolve and mature, especially the digital computer and software implementations such as the finite element method, the opportunities for cross-fertilization between the physical and biological domains increase. We hope that this presentation will encourage other attendees of this conference to turn their attention to the study of biology and to the transfer of the insights gained and methods developed in the physical sciences. This presentation will review an application in plant physiology – an area that has attracted far less attention to date than has been the case for biomedical studies and animal physiology. The finite element method makes it feasible to deal effectively with the typical biological circumstances of increased geometrical complexity and highly anisotropic and nonlinear material properties and composites. Our focus will be on the surprisingly deep insights derived from studies of 30 and 20 years ago (Cooke et al. [1, 2] and Lee [4]). Our efforts utilized software, FESIA, being developed at the time by Richard Gallagher and colleagues (Thomas and Gallagher [5]) here at Cornell in structural engineering. Their interest had been motivated by some rather dramatic collapses in multi-story tall concrete structures coming into usage at that time as cooling towers. They had developed a shell element that would accommodate a doubly curved surface such as occurs in a hyperbolic paraboloid. As an illustration of the theme of this conference, they were interested in a megastructure and we were interested in a structure of microscopic dimensions, i.e., invisible to the unaided eye. One other important attribute deserves mention at the outset: The system we are studying already exists and is functioning. Our interest primarily is one of reverse engineering, i.e. analysis rather than design. 2. Background 2.1 Plant biology considerations: Our project sought to understand better the processes of carbon dioxide and water exchange between plants and the environment, especially for agronomically important crops, that through evolution developed an effective feedback control system for regulating this gas exchange process. Higher plants provide an internal environment in which the energy received from the sun is converted into energy usable by the plant. That process, called photosynthesis, requires a source of carbon, which is obtained from the carbon dioxide in the atmosphere. Terrestrial plants are able to avoid desiccation in a relatively dry environment as a consequence of a relatively impervious outer surface. In many respects this is analogous to the thermal 1 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca environment provided by a building’s outer surfaces. Passage of carbon dioxide into the plant’s interior is facilitated by specialized passages analogous to windows that are opened and closed by the plant. Specialized epidermal cells on plants, called guard cells, form pores called stoma that open and close (Figure 1) acting as an actuator valve for the passage of gases between the plant’s interior and its environment. Usually these pores are open during the daylight when photosynthesis can proceed, and closed at night, blocking gas exchange. Because the plant’s interior is water vapor saturated and therefore at higher concentration than the ambient environment, any time stomata are open for carbon dioxide uptake, water vapor moves in the opposite direction through the same pathway. Figure 1. Fully opened stomata on lower surface of a cucumber leaf. (from Troughton and Donaldson [6]) One of the two most common geometrical shapes for stomata consists of a pair of kidney-shaped cells which we modeled as an elliptical torus, i.e., as a torus that is elliptical in plan view as seen from the outer surface of the leaf, and also elliptical in elevation view. The opening and closing of the pore was conjectured as early as 1856 to result from opposing forces caused by hydrostatic pressures inside the guard cell and from the adjacent cells. Before the emergence of digital computers and the finite element method, only less realistic modeling of the behavior of the mechanics of the opening and closing of the pore was possible. Various beam theories were applied, but offered little understanding of the process by which pores opened. Two of the present authors (Rand and Cooke) applied Flügge’s (then) recent thin shell theory of a circular torus with disappointing results. Namely, a circular torus of biologically realistic dimensions, when inflated, causes the pore size to decrease, rather than increase. Such behavior would be catastrophic for a plant – the pore must open when water is in ample supply (when some water loss can be tolerated) and close when water is in short supply (flaccid guard cells). (Note: An unphysiological garden hose shape when joined end-to-end and inflated would exhibit an increasing ‘pore’.) This suggested that the observed non-circular shape of stomata played an essential role in their functioning as valves. However, since no analytical solution was available, it was natural to utilize the finite element method to study their behavior. 2.2 Shell considerations Figure 2 shows two views of the doubly elliptical shell studied, which includes a plate at each end of the guard cell pair. This plate does limit the outward expansion of the shell at each end of the pair, but is not required for the pore to open when the internal pressure is increased. Geometrical nonlinear aspects reveal some additional properties. One of the authors (Mang) generalized the FESIA software to handle geometric nonlinearities. Unsurprisingly, the nonlinear studies revealed that there are limits to the size of the pore opening as a result of increases in guard cell pressure. That is, the pore cannot be made arbitrarily large. Indeed, as the circular model showed, a contrary behavior would be expected as a more circular geometry is approached. Mang also generalized the FESIA software to include orthogonally anisotropic material properties so that we could model the role of the micellae (radially oriented microfibrils). A new finite element for shell analysis was introduced by Lee in 1986 and used in the subsequent calculations. Figure 2. Toroidal shell – elliptical in plan and elevation views. Bottom view at right shows the plate at the end. From Lee [4]. 2 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca 3. Results The elliptical geometry is the fundamental attribute governing the mechanism of pore opening! Nowhere in the plant kingdom does a circular stoma exist. That shape simply would not operate correctly. Only half of a guard cell is shown in Figure 3. Elevation views: The top row has circular and the bottom elliptical cross sections. Plan views: The left column has circular and the right one elliptical top views. Figure 3. Four geometrical shapes: a) circular torus (plan and elevations) with plate in end, b) circular torus (plan view) elliptical in elevation with plate in end, c) elliptical torus (plan view) circular in elevation with plate at end, d) elliptical in plan view and elliptical in elevation (without end plate) From Lee [4]. a c A circular torus with a circular cross section (a) would disturb the adjacent epidermal surface and the pore would not open. When the elevation cross section is made elliptical (b) the disturbance to the adjacent leaf surface is less, but the pore does not open sufficiently when inflated. If the top (or plan) cross section is made elliptical (c), but the elevation cross section remains circular, the adjacent cells would be disturbed and the pore would be little influenced by pressurization of the guard cell. Finally, when the torus (d) is doubly elliptical – plan and elevation – the guard cell does not bulge into the adjacent cells of the leaf surface and the pore width increases (moving in a direction opposite to the internally applied pressure). b d Figure 4 shows the deformed shape of a doubly-elliptical toroidal shell. The shape of the deformed elliptical torus reveals another attribute of importance that has previously been ignored. Namely, when inflated, the shell expands mostly vertically (out of the plane of the leaf surface), minimizing the stretching of the leaf surface. Figure 4 shows the deformed shape for an elliptical torus having an elliptical elevation cross section and end plate. Top-right shows the deformed shape superimposed on the mesh. Bottom-right shows the major principal stresses at the middle surface. From Lee [4]. The shell model did not impose a constraint on the length of the pore but that movement is inherently negligible, as reported in the experimental literature and confirmed by both the linear and nonlinear analyses. 3 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca Prior to our work, the two prevailing theories of the mechanics of pore opening attributed the response to 1) a thicker ventral wall of the guard cell (facing the pore) than occurs for the opposite guard cell wall (facing the surrounding cells) and 2) the presence of radial cellulose microfibrils in the guard cell, resulting in anisotropic material properties. Our finite element analyses showed that a) the dominant consideration is geometry, b) a perfectly uniform wall thickness for a doubly-elliptical shell will work correctly, and c) a perfectly isotropic elliptical torus will open under static pressure when inflated BUT the anisoptropic property becomes a decided advantage when the dynamical response of the feedback control system is considered. Specifically, the internal volume (lumen) of the guard cell is much smaller than the enclosed volume of the surrounding cells so when a unit of water diffuses from the larger volume into the smaller volume, the pressure drop in the surrounding volume is less than the pressure increase in the smaller volume. Hence, the radial stiffening allows a smaller external pressure change to offset a larger pressure change inside the guard cell. That is, in general, an increase in internal pressure opens the pore while a smaller increase in the surrounding pressure is able to counteract and close the pore. Indeed, the pore width can be represented as a multilinear function of the internal pressure and the pressure in the adjacent cell. This result is consistent with two very different experimental approaches: 1) pressurized inflation of the guard cell with the surrounding cells ruptured and pressuration of the surrounding cells with the guard cell ruptured, and 2) a plasmolytic approach in which osmotic conditions were used to alter the pressures. (A subsequent study of the dynamics of pore opening (Delwiche and Cooke [3]) revealed that this relationship must be clipped to exclude negative pore widths and that this is the key issue in causing a stable oscillatory behavior (a limit cycle) for stomata under drought stress.) (Double-click to see the animation at right.) 1.50 1.25 B 1.00 C D 0.75 0.50 The pore width rises in a nonlinear manner and plateaus A with increasing pressure, i.e., the width cannot be increased arbitrarily. Decreasing the wall thickness increases the pore size for a given pressure. The greater the stiffness provided by the micellae, the greater the maximum pore width achievable. However, most of the diffusive control is achieved with smaller pore widths, making larger widths unnecessary. 0.25 Figure 5 Pore width vs normalized internal pressure for the nonlinear model (Cooke et al. [2].) 0.00 References P/E [1] Cooke JR, DeBaerdemaeker JG, Rand RH and Mang HA. A Finite Element Shell Analysis of Guard Cell Deformations. Transactions of the ASAE, 1976; 19(6):1107-1121. [2] Cooke JR, Rand RH, Mang HA, DeBaerdemaeker JG. A Nonlinear Finite Element Analysis of Stomatal Guard Cells. Unpublished ASAE Paper No. 77-5511. 1977 Winter Meeting, Chicago, IL. 19 pages. ASAE, St. Joseph, MI 49085 microfiche. [3] Delwiche MJ, Cooke JR. An analytical model of the hydraulic aspects of stomatal dynamics. J Theoretical Biology. 1977. 69:113-141 [4] Lee JY. A Finite Element for Shell Analysis and its Application to Biological Objects. 1986. 249 pages. Ph. D. Thesis, A. R. Mann Library. Cornell University, Ithaca, NY 14853. [5] Thomas GR & Gallagher RH., A triangular element based on generalized potential energy concepts in finite elements for thin shells and curved members. Ashwell and Gallagher (eds) John Wiley & Sons. London 1976. pp155-169 [6] Troughton J and Donaldson L. Probing plant structure. Chapman and Hall, London. 1972. 116 pages. 4 Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) In vivo ultrasound bone property determination through inverse finite element modeling Mija HUBLER*, Wilkins AQUINO, Christopher EARLS *Civil and Environmental Eng. Cornell University 220 Hollister Hall mhh33@cornell.edu Abstract Numerical modeling of the femur has been used as a means of better understanding loading conditions of bones and in the design of prosthetics. Recent advances have been made in using numerical models of bone and inverse solution techniques which exploit the properties of surface waves [1]. In an effort to advance dynamic structural modeling and disease detection capabilities below the bone surface, finite element models of vibroacoustic tests are used in this paper in conjunction with inverse solution techniques to determine bone properties. Preliminary results show the finite element model is sensitive to changes in the inner layers of the bone. 1. Introduction Bone fracture and bone diseases pose a national health threat, affecting over an estimated 44 million Americans according to the National Osteoporosis Foundation in 2007. Previous work using numerical modeling techniques to better understand bone structure and composition has focused on in vitro conditions [2], dynamic loading conditions of the bone [3], or the development of prosthetics. Yet scientific studies correlating in vitro bone strength and density have not successfully been able to extrapolate out of their experimental population, or to in vivo conditions [4]. Recently, it has been demonstrated that ultrasonic guided waves through the surface of the outer cortical layer of bone provide a feasible method for assessing structural and material properties of cortical bone [1]. This work develops a method to work with data from in vivo samples to determine properties in the depth of the bone. Using finite element modeling allows for the consideration of complicated wave dynamics which closed form based inverse solvers cannot consider. This is an advantage since the disease develops gradually in the lower bone layer of trabecular structure, and later causes the thinning of the outer cortical bone, where it can be detected currently. Dispersion and Biot theories have attempted to describe ultrasound wave propagation in porous structures such as bone, yet consistently underestimate the experimentally observed results [5]. 2. Method A forward finite element model will first be constructed of the bone and surrounding fluids. This model is built using continuum and fluid potential elements discretized enough in space and time to allow the required wave modes to manifest themselves. Careful convergence studies have been carried out to ensure correct fluidstructure interactions. This finite element model will serve as an environment to reproduce experimental results done by using vibroacoustograhy on patients at a clinic. Vibroacoustograhy was developed at the Mayo Clinic as an imaging modality based on ultrasound [6]. The vibroacoustography method uses two transducers to produce ultrasound beams at slightly shifted frequencies to intersect in the fluid-encased body as shown in Figure 1. The ultrasound radiation force created at the intersection causes the body to dynamically respond. This response is then detected using a hydrophone and processed to obtain a clean signal. Using this data and the forward finite element model, inverse solvers, such as Monte Carlo or genetic algorithms, may be applied to solve for material properties of the given human femur. 1 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca Figure 1: The vibroacoustography method for bone uses two ultrasound transducers of slightly different frequencies which intersect on the specimen and cause an oscillating stress field to form. The resulting dynamic response of the structure is measured using a hydrophone. 3. Sensitivity Studies Sensitivity studies are a crucial element to the inverse solution process. They provide the link between the system response and material properties. Since no particular wave is excited in our experimental setup, the response parameters are not obvious. Our goal is to detect material variations within the lower layers of the bone from the dynamic response, and to determine the limits of our forward model. All of the following sensitivity studies were done on a two dimensional continuum finite element model that neglects fluid effects, which have been shown to be minor. Sensitivity plots are shown in Figure 2. The intersecting ultrasound beams create a force perpendicular and parallel to the beam at the focal point. Thus, it is possible to angle the excitation in order to obtain a response that may be more sensitive to changes in the properties we would like to detect. In our first study the forward finite element model was excited with a pulse load at the center of the bone section with various angles with respect to the surface. A maximum angle of 30 degrees was selected since this approaches the critical angle due to Snell’s law. Sensitivity was characterized by comparing the norm of the vertical surface displacements through time as given by equation 1. (1) where n refers to the current configuration, x the surface location, t the time, and y the measured vertical surface displacement. The units on the norm were removed by division of the largest norm. Next, the sensitivity of the bone response to different material properties of the trabecular bone layer was investigated. According to literature the Young’s modulus has been estimated as an average of 0.8Mpa for bone of trabecular structure, with a vartiation from 0.05GPa to 1.05Gpa [7]. Knowing the sensitivity of responses to this variation would show that an inverse solution can detect this material property. The forward model was again excited with a pulse load at the center of the bone section while the modulus of the trabecular bone was varied. Sensitivity to this variation was characterized using the same norm approach. Once a forward model has been developed to be used in the inverse solution algorithm it is a challenge to vary boundary conditions. As a result, it is of importance to investigate the effect of varying boudary conditions on the forward model as prescribed by the natural disease progression. In the last study the trabecular bone layer was kept at it’s average material property values while studying the response sensitivity to a thickness variation from 3.41mm to 1.45mm. 2 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca (a) (b) (c) Figure 2: Sensitivity studies of variations in loading angle (a), trabecular bone Young’s modulus (b), and thickness (c) were done on a two dimensional continuum finite element model. 4. Discussion In the sensitivity study of loading angle, the response becomes more pronounced when the load is less perpendicular to the bone surface. This enhanced response, combined with other data analysis techniques will allow for a more refined conclusion to the inverse problem. As a result, the other two sensitivity studies are shown at a 30 degree loading angle. The forward model displays sensitivity to both variations in Young’s modulus and thickness of the trabecular bone layer. In the case of sensitivity to the variation in Young’s modulus, this is a positive result, suggesting it will be possible to determine such changes in patients. Yet, sensitivity to thickness describes a further challenge in constructing a finite element model for a system with undetermined boundaries. Correlations between bone thickness and strength may provide a means to work around this challenge. 5. Conclusion The design of an inverse solution method, in conjunction with a forward finite element model, to work with vibroacoustic data of in vivo bone would allow for osteoporosis diagnosis at early stages of the disease. Sensitivity studies presented here on the forward finite element model suggest the inverse solution is feasible if a method of considering sensitivity to bone thickness may be devised. Further work will focus on the development of a more refined three dimensional forward model in conjunction with an inverse solution scheme. 3 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca References [1] Moilanen P. et. al. Assessment of the cortical bone thickness using ultrasonic guided waves: Modeling and in vitro study. Ultrasound in Medicine and Biology 2007; 33:254-262. [2] Barkmann R. et. al. A method for the estimation of femoral bone mineral density from variables of ultrasound transmission through the human femur. Bone 2007; 40:37-44. [3] Taddei F. et. al. Subject-specific finite element models of long bones: An in vitro evaluation of the overall accuracy. Journal of Biomechanics 2006; 39:2457-2467. [4] Lotz J, Cheal E, and Hayes W. Fracture Prediction for the proximal femur using finite element models: Part 1 – Linear Analysis. Journal of Biomechanical Engineering 1991; 113:353-360. [5] Haire T and Langton C. Biot Theory: A review of its application to ultrasound propagation through cancellous bone. Bone 1999; 24:291-295. [6] Fatemi M and Greenleaf J. Vibro-acoustography: An imaging modality based on ultrasound-simulated acoustic emission. Proceedings of the National Academy of Science USA 1999; 96:6603-6608. [7] Wirtz D. et. al. Critical evaluation of known bone material properties to realize anisotropic FE-simulation of the proximal femur. Journal of Biomechanics 2000; 33:1325-1330. 4 Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) Finite element analyses of palm leaf petiole-sheath junctions in simple bending and twisting and in dynamic (oscillatory) flexure Karl J. NIKLAS*, J. Robert COOKE, Jae Young LEE *Department of Plant Biology 208 Plant Science Building Cornell University kjn2@cornell.edu Abstract Finite element analyses were used to estimate qualitatively stress and strain distributions in a stereotypical palm leaf at the petiole-sheath junction when subjected to simple bending and twisting and to simple linear dynamic wind-inducted flexure. This qualitative approach was used because (1) palm leaf morphology and mechanical behavior are unquestionably the most complex in the plant kingdom (e.g., leaf length varies between 0.25 m and 25.1 m), (2) the geometry and material properties of palm leaves vary ontogenetically as well as among different species, and (3) plant morphology and anatomy in general manifests extreme heterogeneity in cross-sectional shape and size and extreme anisotropic material composition. Our analyses reveal heterogeneous stress and strain distributions within the petiole-sheath junction when leaves experience simple bending or twisting and morphological distortions; predicted morphological distortions corresponding to those observed when real palm leaves were bent and twisted. Likewise, simulations of dynamically loaded petiole-sheath junctions reveal strains that are consistent with morphological changes when real leaves are subjected to strong gusts of wind. Finite element analyses indicate that maximum stresses occur in anatomical regions occupied by plant tissues that are strong in tension and relatively extensible compared to neighboring tissues. From these analyses, we conclude that the stereotypical palm leaf base is highly adapted to mechanically support large leaves experiencing large static and dynamic loadings. 1. Introduction In a seminal paper dealing with the leaf base of palms, Tomlinson wrote “. . . the construction and behavior of the leaf base of palms is [sic] a complicated subject which is of considerable morphological interest. It is a problem which can only be understood by considering mechanical aspects of leaf development in relation to the growth of the stem” (Tomlinson, 1962). This assessment remains unchallenged even today. The mechanical behavior and complexity of the palm leaf have no parallels anywhere in the plant kingdom, because in addition to their often great size (which can reach a record length of 25.11 m in the case of Raphia regalis), palm leaves possess a tubular leaf sheath (which completely encircles the stem) that supports a cantilevered petiole bearing the leaf lamina (Fig. 1). This tubular sheath, which is anatomically complex owing to a vascular infrastructure of intersecting fibrous strands (Esau, 1967; Gifford and Foster, 1996), persists in some palms as the leaf matures, grows in size, and as the cantilevered petiole bearing the leaf lamina develops (e.g., Cocos, Phoenix, and Trachycarpus). In terms of its biomechanics, the stereotypical palm leaf base normally experiences four kinds of mechanical stresses: stresses resulting from the growth of younger leaves it Figure 1: Representative palm petiole-sheath junctions. The fleshy tissues of the tubular leaf sheaths of older leaves have broken down to reveal the infrastructures of fibrous vascular strands that provide mechanical support. Younger leaves are enveloped by the sheaths of older leaves. 1 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca envelopes, stresses due to the expansion of the stem, and stresses resulting from the weight of the leaf (self- loading), and those resulting from dynamic wind-loading. In addition, the leaf base may experience stresses resulting from the growth of reproductive structures bearing flowers (inflorescences). Because of this morphological and anatomical complexity, traditional biomechanical analyses of surgically isolated tissues or portions of the palm leaf have been largely uninformative, particularly in helping to evaluate stress-strain distributions in the petiole-sheath junction. To remedy this situation, we used finite element analyses to assess the distribution of the stress magnitudes in this portion of a stereotypical leaf base when it is subjected to simple bending or twisting analyses. We also used this approach to explore the flexural distortions in the petiole-sheath when a leaf is subjected to dynamic wind-induced displacements. Because the size, shape, and tissue material properties of palm leaves change ontogenetically over the course of leaf development and differ among the mature leaves of different species, our analyses necessarily were qualitative rather than quantitative, i.e., we focused on general patterns in stress distributions and gross displacements (stresses). Nevertheless, our analyses point to a number of broad patterns that are consistent with empirical observations and highlight the remarkable functional (adaptive) nature of palm leaf morphology and anisotropic tissue composition. 2. Materials and Methods 2.1 Simple bending and twisting simulations The nodal displacements were computed by solving the system equations and the strains and stresses were recovered at each integration points using the formulas circumferential direction length direction thickness direction rs t ε = BΔ σ = Eε (1) (2) thickness direction length direction length direction circumferential direction thickness direction circumferential direction € € Figure 2. Mesh and anisotropic properties used NOTE: See the Extended Abstract CD for enlarged figures. Figure 2 lists the numerical values used to evaluate the stresses and strains. However, it is important to note that our analyses provide only qualitative assessments of stress and strain distributions because the foregoing numerical values are simply indicative, and do not apply to all palm species. 2.2 Dynamic simulations Dynamic equilibrium equation MΔ˙˙+ CΔ˙ + KΔ = F where M=mass matrix, C=damping matrix, K =stiffness matrix, Δ˙ =velocity, and Δ =displacement. (3) F =force vector, Δ˙˙ =acceleration, € € € € 2€ € 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca The vibration modes are extracted by solving the eigenvalue problem MΔ˙˙+ KΔ = 0 (4) for free vibrations without damping and external forces. The square roots of eigenvalues represent the frequency of vibration (radian/sec), and corresponding eigenvectors represent the vibration mode shapes. The nodal displacements at a given time step is computed by linear combination of the vibration mode shapes. Ienigeonuvralauneaslyasreis,ext€htreactmedodferomsuptehreposmsitailolnestmveatlhuoed, was applied with Newmark integration scheme. The i.e., lowest frequency. The modes with lower frequency usually have more effect to the dynamic motion than those with higher frequency. Only the lowest 7 modes were used in the analysis, because the rest of the modes were assessed to have insignificant effects. 3. Results 3.1 Simple bending and twisting simulations -W -W W Figure 3. Simple bending (left) and twisting (right) loading Figure 4: A finite element predicted stress distributions for a stereotypical palm leaf petiole (at right)- sheath (left) junction subjected to simple bending. Top: Length direction; middle: thickness direction; bottom: radial direction Figure 5: A finite element predicted stress distributions for a stereotypical palm leaf petiole (at right)- sheath (left) junction subjected to twisting. Top: Length direction; middle: thickness direction; bottom: radial direction 3 6th International Conference on Computation of Shell and Spatial Structures 3.2 Dynamic wind-induced flexure simulations IASS-IACM 2008, Ithaca Figure 6: Seven lowest modal shapes for a stereotypical palm leaf petiole (at right)-sheath (left) junction. Figure 7: Transient: Left column (wind horizontal: middle max; bottom min); right column (wind upward: middle at max; bottom at minimum). Double-click the top row of images for animations. Figure 6 displays the mode shapes for the seven lowest natural frequencies because the rest of the modes turned out to have an insignificant effect. The modes with lower frequency usually have more effect to the dynamic motion than those with higher frequency. Figure 7 shows snapshots of the response to a single half-sine-wave when directed horizontally (left column) and upward (right column). NOTE: Flash animations (.SWF and .MOV files) of wind-induced transient oscillations are contained in the CD version of this abstract. The embedded .mov files require that a QuickTime player be installed. References [1] Hughes TJR. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis (2nd edn). Dover, 2000. [2] Simo JC and Armero F. Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes. International Journal for Numerical Methods in Engineering 1992; 33:1413-1449. [3] Wilson EL, Taylor RL, Doherty WP, Ghaboussi J. Incompatible displacement models. In Numerical and Computer Models in Structural Mechanics, Robinson AR, Fenves SJ, Perrone N, Schnobrich WC (eds). Academic Press: New York, 1973; 43-57. 4 Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) Micropipette aspiration of lipid vesicles: A 2-D approach Philip BUSKOHL*, Sovan L. DAS, James T. JENKINS *Theoretical & Applied Mechanics, Cornell University Kimball Hall 315, Ithaca, NY 14853-1503 prb28@cornell.edu Abstract Biological systems extensively utilize lipid membranes to carry out the various functions for life. The boundaries of the cell are defined by lipid bilayers that protect the cytoplasm from intruders and also regulate the traffic of nutrients and by-products to and from the cell. In regions of shallow curvature or substantial lateral tension, the effect of bending is typically neglected due to the thinness of the membrane. However, when the curvature of the membrane becomes large, the bending resistance plays an integral role. The formation of tubules, budding of vacuoles, and mitosis are all vital processes that sustain cellular life and involve significant changes in membrane curvature. Micropipette aspiration of synthetic lipid membranes is often employed to study their bending resistance. Here, a boundary layer analysis on a 2-D outline of a vesicle under aspiration is carried out. For comparison, the Finite Element Method (FEM) will be applied to the governing non-linear differential equation and used to numerically investigate the influence of the bending stiffness parameter. The aim of this work is to develop a method to quantitatively calculate the bending stiffness of the membrane from quantities measured during the aspiration experiment. 1. Introduction The function of lipid membranes in sustaining life motivates the study of their mechanical properties. Elasticity models of lipid vesicles have been proposed by Helfrich[8], Jenkins[10,11], and others as a special case of thin shell theory. Helfrich [8] proposed a free energy function dependent on the Gaussian curvature and the square of the mean curvature based on symmetry arguments. The general equilibrium equations were then derived by Jenkins [11] for a vesicle with a constant pressure loading and an area constraint. In the axisymmetric case these partial differential equations reduce to a system of first order non-linear ordinary differential equations (ODEs). Taking advantage of this reduced form, Das and Jenkins [3] performed a boundary-layer analysis for a vesicle with two different phases of lipid composition and determined the line tension at the interface between the phases. A natural experimental extension of this work is the aspiration of two-phase vesicles. Pipette aspiration of two-phase lipid bilayers has recently been studied by Tian et al. [12] and Das et al. [4]. Figure 1 shows images of two-phase and single-phase vesicles undergoing aspiration. Figure 1: Two-phase lipid vesicle and single phase vesicle undergoing aspiration. (Tian et al. [12]) The initial analysis considers only the two dimensional projection, or outline, of the vesicle as shown in Figure 1. In the following, boundary-layer and FEM techniques are formulated to model the “truncated-ring” shape of this projection. 1 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca 2. Aspiration Analysis The vesicle is divided into three regions that differ either by loading condition or material composition. Boundary and continuity conditions are imposed at the interface of each region. Figure 2a outlines these three regions, and Figure 2c defines the Cartesian and local frames used in this analysis. The solution curve is parameterized by arc length and, due to symmetry, only half of the vesicle is analyzed. Figure 2: (a) Vesicle description (b) FBD of ring element (c) Cartesian and local frames The equilibrium equations are derived from linear and angular momentum balance, neglecting higher order terms, using the free body diagram (FBD) of the ring element shown in Figure 2b. Evans[6] derived the same system of ODEs for the analysis of membrane adhesion. The resulting system of ODEs, combined with a linear constitutive relationship between moment and curvature, are made dimensionless using pipette pressure, total arc length, and bending stiffness. After some manipulation, the system of first-order equations are consolidated into one second-order differential equation for the curvature, h, as a function of dimensionless arc length, s .  2h  1  2h3  h 1  0 2 (1) In equation (1), the lateral tension  is constant and has been normalized by p1  p. The bending stiffness,  , acts as the small parameter   ( p1  p)l3 (2) because the pipette pressure p1  p and total arc length, l are assumed to be order one quantities. As seen in equation (3a), continuity of curvature is imposed at the mouth of the pipette, s  l1 , along with a jump condition in transverse shear force due to the pipette contact force, f1 at the mouth. Pipette Mouth: s  l1 h  h   2 h  h  f1  0 (3a) Separation of Phases: s  l2 h  h   2 h  h  f2  0 (3b) At the interface between the two phases, continuity is not enforced, but instead a jump in both shear and curvature occurs(equation 3b). The shear force f2 directly corresponds to the line tension between lipid phases. The curvature is assumed known at this interface, and at the boundary of region (-) that is not in contact with the pipette mouth(refer to Figure 2a). 2 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca 2.1 Perturbation Method Equation (1) represents a nonlinear “beam-string” problem. Far away from the mouth of the pipette, the curvature approaches the solution of a string which, when under a uniform load, is equal to a constant. The string solution is denoted as the “outer” layer of the region (-). In a narrow regime on either side of the pipette mouth, bending significantly influences the solution because of large curvature changes. By appropriate scaling, the arc length of this narrow regime is magnified, so that the higher order terms of equation (1) act as order one quantities. The solution of this scaled equation forms the “inner” layer. Finally the layers are matched in the limit as the two solutions approach each other. For a more detailed discussion of perturbation methods and lipid vesicle boundary-layer analysis see Hinch [9], Das and Jenkins [3], Das [2]. The preliminary results of this analysis, up to first order in epsilon, are outlined below:        h(s, )   2 0 sech  0 ~s  B  C1 sec h 0 ~s  B tanh 0 ~s  B 2 tanh2  0 ~s  B 0 1 0 (4)        h(s, )   2 0 sech  0 ~s  B  C1 sec h 0 ~s  B tanh 0 ~s  B 2 tanh2  0 ~s  B 0  0 (5)    C1  C3  (1  ) 1 2 tanh2B 0 sech(B) tanh(B) C1  C3   4 tanh(B)sech(B) 1 o sec h2 (B)  tanh2 (B) sinh2 (B)  40 sinh(B) 1  0 f1 (6) where ~s  s  is the stretched arc length ,   p0  p , 0 is the first term in the expansion of lateral tension, p1  p and C1, C3, B are constants determined by equations (6). To analyze the second phase region (++) a new small parameter  must be introduced because the bending stiffness is different in this region. An expression for the ratio of the two small parameters will potentially yield the desired method to quantitatively measure bending stiffness. The boundary-layer calculation of an aspirated single phase vesicle is compared with the numerical solution from Matlab’s boundary value problem solver bvp4c in Figure 3. There is good agreement between the order one analytical solution and the numerical solution. Figure 3: Comparison of boundary-layer and numerical solutions: Single Phase Aspiration [Dark Blue: Aspirated Region(-); Light Cyan: Non-Aspirated Region(+)] 2.2 Finite Element Method The finite element method (Belytschko [1]) has been formulated to solve the non-linear ODE governing vesicle curvature. The model uses linear bar elements to approximate the curvature on a non-uniform grid. The number of nodal points increases with proximity to the interface between region (-) & (+). The trial functions belong to the space H 1 0 0,1 defined as the space of all functions that vanish on boundaries with defined curvature, and 3 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008, Ithaca also have square integrable first derivatives. Using this trial function, the general and discretized weak forms are:  w 2h b  a  2   2 wh3  2wh  wh  ds  w ds  0 (8) n h  Ni hie i 1  h  n dNi i 1 ds hie n  Bi hie i 1 n w  Ni wie i 1  w  n dNi i 1 ds wie n  Bi wie i 1  {w}T    2[N ]T [B]{he} b a      2 2 [ N ]T [ N ]{h}2 [ N ]   2 [ B]T [ B]   [ N ]T [ N ]ds  {h  e }   [N ]T ds   0  (9) The FEM will be implemented using a Newton-Raphson iteration scheme for each region of the vesicle. The boundary condition at the mouth, equation (3), will be enforced by modifying the residual of both regions simultaneously during the iteration. Through numerical integration and appropriate boundary conditions, the tangent angle and position of the vesicle will then be determined from the curvature solution. 3. Future Work The 2-D analysis will serve as a stepping-stone toward understanding the 3-D aspiration boundary-layer calculation. The relative transparency of the 2-D calculation will offer insight into the nature of the full solution. Numerically, the goal is to model the 3-D vesicle aspiration with finite elements. Dong and Skalak [5] pursued this problem by implementing sliding boundary condition on the aspirated portion of the mesh. More recently Feng and Klug[7] developed an FEM model of lipid bilayer membranes using specific C1 -conforming triangular loop subdivision elements. An understanding of bending effects in pipette aspiration will not only give researchers a means to quantitatively measure membrane bending stiffness, but also offer insight into the mechanics of many complex biological processes. References [1] Belytschko T, Liu WK, Moran B. Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons Ltd, 2000. [2] Das SL. “Studies of axisymetric lipid bilayer vesicles: parameter estimation, micropipette aspiration, and phase transition”. Ph.D dissertation Cornell University. 2007; 21-40, 53-85. [3] Das SL and Jenkins JT. A higher order boundary-layer analysis for lipid vesicles with two fluid domains. Journal of Fluid Mechanics. 2007; 597: 429-448 [4] Das SL, Tian A, Baumgart T. Mechanical stability of micropipette aspirated giant vesicles with fluid phase coexistence. Journal of Physical Chemistry B. 2008; (submitted in revision) [5] Dong C and Skalak R. Leukocyte deformability: finite element modeling of large visco-elastic deformation. Journal of Theoretical Biology. 1992; 158: 173-193. [6] Evans EA. Detailed mechanics of membrane-membrane adhesion and separation. Biophysical Journal. 1985; 48: 175-183. [7] Feng F, Klug W. Finite element modeling of lipid bilayer membranes. Journal of Computational Physics. 2006; 220: 394-408. [8] Helfrich W. Elastic Properties of Lipid Bilayers: Theory and Possible Experiments. Z. Naturforsch. 1973; 28c; 693-703. [9] Hinch EJ. Perturbation Methods. Cambridge University Press, 1991. [10] Jenkins JT. The equations of mechanical equilibrium of a model membrane. SIAM Journal of Applied Math. 1977; 32; 755-764. [11] Jenkins JT. Static equilibrium configurations of the model red blood cell. Journal of Mathematical Biology. 1977; 4: 149-169. [12] Tian A, Johnson C, Wang W., Baumgart T. Line tension at fluid membrane domain boundaries measured by micropipette aspiration. Physical Review Letters. 2007; PRL 98: 208102. 4 Proceedings of the 6th International Conference on Computation of Shell and Spatial Structures IASS-IACM 2008: “Spanning Nano to Mega” 28-31 May 2008, Cornell University, Ithaca, NY, USA John F. ABEL and J. Robert COOKE (eds.) A new model for nucleation in two-phase lipid bilayer membrane vesicles Sanjay DHARMAVARAM* and Timothy HEALEY Department of Theoretical and Applied Mechanics Cornell University Ithaca NY 14853, USA sd282@cornell.edu; tjh10@cornell.edu Abstract This work is motivated by the experiments of Baumgart, et. al. (Nature, 425, (2003) p. 821) on giant unilamellar vesicles (GUVs) in which they observe remarkable two-phase patterns on GUVs when these are subjected to varying osmotic pressure or temperature. These man-made vesicles are of interest because they are easier to work with (when compared to real cell membranes) but at the same time exhibit many features of real cell membranes. GUVs also have a great potential in drug delivery systems and other biomedical applications, and hence, understanding their mechanical behaviour is of great interest. The vesicles used exhibit two phases - ordered (Lo) and disordered (Ld) - which are both fluid-like. Their molecular structure suggests that lipid layers in the ordered phase have a greater thickness than the ones in the disordered phase. The aim of our work is to propose a new model for two-phase lipid bilayers and use the GUVs as an example to test our theory. Although Ginzburg-Landau type theories have been proposed to explain phase separation in multi-phase lipid bilayers. These theories usually introduce an ad-hoc curvature-phase coupling energy. Here, we instead propose a pressureinduced phase transition model using a van der Waals non-convex potential for in-plane effects and a bending potential energy which takes into consideration the change in thickness associated with each phase. The latter turns out to be crucial and introduces a natural coupling between curvature and phase. We derive the equations of equilibrium and consider a nominally spherical vesicle under inflationary pressure (so that we do not confuse things with standard shell buckling under compression). We then linearize the equations about a spherically symmetric state and using rigorous nonlinear bifurcation analysis and group-theoretic strategies show the existence of non-spherical bifurcated equilibria. These equilibria represent phase-nucleated states from the homogeneous vesicle and bear striking similarity to many of the non-spherical states observed in the experiments mentioned above.