|Home | About | Journals | Submit | Contact Us | Français|
Shape plays an important role in determining the biomechanical response of a structure. Specimen-specific finite element (FE) models have been developed to capture the details of the shape of biological structures and predict their biomechanics. Shape, however, can vary considerably across individuals or change due to aging or disease, and analysis of the sensitivity of specimen-specific models to these variations has proven challenging. An alternative to specimen-specific representation has been to develop generic models with simplified geometries whose shape is relatively easy to parameterize, and can therefore be readily used in sensitivity studies. Despite many successful applications, generic models are limited in that they cannot make predictions for individual specimens.
We propose that it is possible to harness the detail available in specimen-specific models while leveraging the power of the parameterization techniques common in generic models. In this work we show that this can be accomplished by using morphing techniques to parameterize the geometry of specimen-specific FE models such that the model shape can be varied in a controlled and systematic way suitable for sensitivity analysis. We demonstrate three morphing techniques by using them on a model of the load-bearing tissues of the posterior pole of the eye. We show that using relatively straightforward procedures these morphing techniques can be combined, which allows the study of factor interactions. Finally, we illustrate that the techniques can be used in other systems by applying them to morph a femur. Morphing techniques provide an exciting new possibility for the analysis of the biomechanical role of shape, independently or interaction with loading and material properties.
Finite element (FE) modeling has been employed successfully to study the mechanics of biological structures. One essential ingredient of FE modeling is the geometry or shape of the structure to model (the other essential ingredients being the mechanical properties, loading and boundary conditions). Traditionally there have been two approaches to defining the model geometry for FE modeling in biomechanics: generic or specimen-specific. Generic models are developed from general dimensions and mechanical properties of a population. A particularly powerful feature of generic models is that they can be parameterized, such that aspects of the model geometry, mechanical properties and loading are defined in a way that they can be varied by specifying a parameter value. Parametric studies using generic models have shown that interactions between parameterized factors can be substantial and often non-intuitive (Rekow, Harsono et al., 2006; Sigal, 2009). In many cases, however, generic models cannot represent the complex details of the 3D architecture of biologic structures, and thus cannot make predictions about the biomechanics of a particular specimen.
Specimen-specific models incorporate more of the details that make a specimen unique, and therefore may predict individual biomechanics more accurately. However, specimen-specific models also have limitations: by nature, specimen-specific models account for a limited set of combinations of geometry, mechanical properties and boundary conditions, those combinations particular to the specimens included in the study. Hence, studies using specimen-specific models cannot control parametric variations in the models to the extent possible with generic models, and rely instead on statistical analysis. Unfortunately, despite recent advances, preparation and analysis of specimen-specific models is often time consuming and difficult, which limits the number of models that can be used in a study, thereby reducing the statistical power restricting the results and conclusions (Viceconti and Taddei, 2003; Taddei, Pancanti et al., 2004).
When generic and specimen-specific models have been combined, most often it has been accomplished using specimen-specific geometries and parameterizing the mechanical properties and boundary conditions (Brolin and Halldin, 2004; Anderson, Peters et al., 2005; Laz, Stowe et al., 2007; Rissland, Alemu et al., 2009; Sigal, Flanagan et al., 2009). The results are useful, but without parameterizing the geometry these studies cannot determine the effects of variations in shape, or of interactions between geometry, mechanical properties and boundary conditions.
In this study we show that it is possible to apply morphing techniques to parameterize the geometry of specimen-specific FE models such that the model shape can be varied in a controlled and systematic way suitable for sensitivity analysis. We demonstrate three morphing techniques by using them on a model of the load-bearing tissues of the optic nerve head (Figure 1). We show that these morphing techniques can be combined for the study of factor interactions, in a way which allows pre-processing and simulation to be accomplished with relatively straightforward procedures. Finally, we also illustrate that the techniques can be used in other systems by applying them to morph a femur into alternate configurations.
Morphing is accomplished by extracting the nodes of the baseline model surface, displacing these nodes (the actual morphing), and then remeshing the interior of the morphed geometry. The total displacement applied to the nodes is the linear combination of partial displacements, each one the result of parameterizing a factor of interest (Figure 2). We demonstrate three techniques for defining the partial displacements: a scaling-based (SB) technique to vary the thickness of the scleral shell; an analytic-based (AB)technique to change the diameter of the scleral canal; and a landmark-based (LB) technique to change the anterior-posterior position of the central lamina cribrosa (LC).
Each partial displacement is the product of three components: a deformation field, a smoothing function, and a scaling constant. The deformation field defines the overall nature of the shape change. The smoothing function eliminates discontinuities and steep gradients in the deformation field that distort the surface mesh during morphing and hinder remeshing. The scaling constant sets the reference for the magnitude of the changes in shape and allows definition of the morphing in simple terms, such as double the thickness, or half the diameter. It is convenient to define partial displacements that are independent of each other. A simple way to achieve this is to base all the displacements on the baseline geometry. If the displacements are not independent, the scaling constants may become scaling functions, which are difficult to calibrate, complicating and reducing the generality of the analysis. In the end, there are two numbers setting the scale and sense of the morphing, the scaling portion of the definition of the partial displacement, and the linear coefficient of the parameterization. The techniques described below were implemented using C++ (Visual Studio v6.0, Microsoft, Redmond, WA), and Amira (vDev4.1.1, Visage Imaging, Carlsbad, CA). See the appendix for the equations and constants.
The techniques are demonstrated using a baseline specimen-specific FE model of the eye of a normal monkey reconstructed following a procedure we have described previously (Burgoyne, Downs et al., 2004; Roberts, Liang et al., In Press (Accepted July 2009)) (Figure 1). The hex-based models of the lamina cribrosa and scleral shell used in our previous work were imported into Amira, converted to 4-node tetrahedral elements and refined twice. The triangulated surfaces were then extracted and smoothed to remove sharp angles. The baseline triangulated surface consisted of 28,398 nodes forming 57,254 triangles.
The model was translated and rotated so that the anterior lamina cribrosa surface centroid coincided with the origin (0,0,0), and the first, second and third moments of inertia coincided with the X, Z and -Y axes. This meant that the lamina cribrosa was roughly aligned with the XZ plane, the Y axis increased in the anterior direction. The origin was just anterior to the lamina cribrosa due to the natural curvature of its anterior surface.
In SB morphing the deformation field, USB, is defined from a vector field implied by the geometry(Figure 3). In the example we parameterized the thickness of the scleral shell using a vector field, t, representing the scleral shell thickness at the nodes on the exterior scleral shell surface. The smoothing was a step function, f(r), where r was the distance from a node to the posterior lamina-sclera interface (the posterior insertion of the lamina cribrosa (Yang, Downs et al., 2009)). This procedure was used to assign a deformation vector to the nodes on the exterior surface of the sclera. All other nodes were assigned a null deformation. Hence, deformations preserved the distribution of scleral thickness variations as well as the shapes of the interior scleral surface, the equator, the scleral canal and the lamina cribrosa. The scleral thickness vectors had a natural scale, the thickness of the baseline model, and therefore no further scaling of the deformation was necessary. When parameterized, the magnitude and sense of the change in scleral shell thickness was controlled by the parameter β, such that β=0.5 reduced the scleral shell thickness to half and β=−1 doubled the thickness.
In AB morphing, the deformation field, UAB, is defined by analytical functions (Figure 4). By combining relatively simple functions, it is possible to define complex deformations while maintaining control of the direction, sense and magnitude of the prescribed nodal displacements. In the example, to parameterize the size of the scleral canal we first computed a function g(a) for all the nodes. The vectors defined by g(a) were radial from the origin on the XZ plane, and had a zero component in the Y-direction. Smoothing was defined as a function of the elliptical distance to the origin, r, normalized using the vertical and horizontal radii of the scleral canal Rx and Rz. This function was chosen to prescribe an expansion which preserves canal wall shape and mesh quality. Unlike SB morphing methods, the scleral canal expansion deformation does not have an implied scale, so it was necessary to define one; in this case k1=200 μm. When parameterized, the magnitude and sense of the changes in scleral canal expansion were controlled by the parameter λ such that λ=1 expanded the canal by 200 μm and λ=−1 shrunk the canal by 200 μm.
In LB morphing, the deformation field, ULB, is defined as the product of a warping field h(r), a smoothing function f(r), and a scaling constant k2 (Figure 5). The warping field is defined through corresponding sets of “source” and “target” landmarks located in space. The classical Bookstein thin-plate spline method is then used to compute a 3D warping of space matching the source to the target landmarks (Bookstein, 1987; Zelditch, Swiderski et al., 2004). In the example, we used nine landmarks on each of the anterior and posterior surfaces of the lamina cribrosa to parameterize the anterior-posterior position of the lamina cribrosa. The target landmarks were placed using a combination of manual displacements and scaling (Sigal, Hardisty et al., 2008). The method produced a smooth deformation field h(r), which was restricted, also smoothly, by f(r), to deform only the portions of the structure intended. We used separate smoothing functions for the anterior and posterior surfaces of the lamina to allow more precise control of the smoothing at the edges and limit the deformations to the lamina cribrosa. The scaling constant k2=0.66 was chosen so that the negative of the displacement produced by the deformation field output of the Bookstein method produced a flat lamina cribrosa (i.e. an anterior lamina cribrosa with zero average depth relative to its periphery, Figure 5, left panel, middle row). When parameterized, the magnitude of the changes in lamina cribrosa position were controlled by the parameter η such that η<0 produced a shallower lamina (a flat lamina when η=−1) and η>0 produced a deeper lamina.
To demonstrate the techniques described above we produced a set of models with various combinations of the parameters β, λ and η, and used these models in FE simulations to predict their biomechanical response to an increase in intraocular pressure (IOP). A total of 125 models were produced, representing all combinations of the three aspects of shape parameterized (laminar position, scleral thickness, and scleral canal size), with each factor varied among five levels. After the triangulated surfaces were deformed the model volume was remeshed with 4-node tetrahedra using Amira. The mesh interior was then iteratively smoothed and relaxed using a Laplacian algorithm until the largest change in nodal locations was smaller than 0.1 μm. For simulation the meshes were converted to 10-node tetrahedra by adding mid-side nodes to the element edges. Simulations were carried out using Abaqus (v6.8, Dassault Systems, Providence, RI, USA). Pre and post-processing were carried out on a Windows XP-based desktop workstation with 2 GB of RAM and two Intel XEON X5365 3GHz CPUs. Models were solved on a Linux server with 32 GB of RAM and four Intel Itanium2 CPUs. From definition through analysis, each model required less than 7 minutes of wall clock time on average. Morphing itself required less than a second for each model. Remeshing (including smoothing-relaxing) and FE solving required, each, about 45% of the wall clock time. Another 5% of the time was spent moving models to and from the server.
Tissues were assumed to be homogeneous, practically incompressible (ν= 0.45), and linearly elastic, with Young’s moduli of 7.56 MPa for the sclera, and 2.04 MPa for the lamina cribrosa (Roberts, Liang et al., In Press (Accepted July 2009)). The effects of IOP were modeled as a homogeneous force of 5 mmHg acting on the element faces exposed to the interior of the eye. The nodes on the equator were restricted to deform on the plane of the equator. Defining boundary conditions and mechanical properties was especially simple because all the models were related to the baseline. Before morphing we identified the surfaces on which boundary conditions had to be defined, and these definitions were preserved through the morphing and remeshing.
We evaluated the quality of the mesh of the morphed models by determining the number distorted elements in each volume mesh using Abaqus’ internal routines (Couteau, Payan et al., 2000; Abaqus, 2003). Adequate element size could also potentially depend on the model geometry, and therefore vary between models. We also selected some cases with particularly large strains or displacements, refined them and checked that their meshes were sufficiently refined based on the maximum principal strain. In all cases they were.
All three morphing techniques were used successfully, independently and in combination, for parameterizing a model of the lamina cribrosa and sclera (Figure 6). The techniques produced smooth transformations in space that preserved the quality of the surface mesh, and were simple to remesh as a result. Model volume meshes were good, with less than 1% of the elements distorted (Table 1). This proportion is on the same order as what has been observed with model morphing using other algorithms (Krause and Sander, 2006; Bade, Haase et al., 2007; Sigal, Hardisty et al., 2008). Interestingly, the number of distorted elements was not lowest in the baseline model. In no case was mesh quality an impediment for FE analysis.
The biomechanical response of the lamina cribrosa and peripapillary sclera was strongly affected by all three parameterized geometric factors (Figure 6). This demonstrates that variations in geometry produce meaningful and strong effects worthy of consideration. Moreover, the effects of the three parameterized geometric factors were not independent, i.e. the factors interacted.
In this study, we introduced three morphing techniques and demonstrated their application by parameterizing a specimen-specific model of the load-bearing tissues of the posterior pole of the eye. The parameterizations can be applied independently, or in combination, allowing rapid production of families of models with carefully controlled variations in geometry that are useful for sensitivity analysis, including analysis of factor interactions. This is useful because it harnesses the detail available in specimen-specific models, while leveraging the power of the parameterization techniques common in generic models.
The morphing techniques also enable a more thorough exploration of the biomechanical consequences of changes in the shape of a specimen, such as may occur during aging or disease, for example the age-related changes in the human rib cage and forearm bones (Bouxsein, Myburgh et al., 1994; Gayzik, Yu et al., 2008) and the remodeling of the lamina cribrosa in early glaucoma (Roberts, Grau et al., 2009). We have focused on parameterizing model geometry in this study because methods to parameterize mechanical properties are well developed for the relatively simple materials we used (Brolin and Halldin, 2004; Anderson, Peters et al., 2005; Laz, Stowe et al., 2007; Rissland, Alemu et al., 2009; Sigal, 2009). For more complex mechanical properties the techniques for parameterizing geometry and mechanical properties must be generalized, and may have to be considered simultaneously. An example would be the anisotropy of the sclera and lamina cribrosa. As the geometry is morphed, the directions of anisotropy should change as well. Deepening of the lamina cribrosa will require adapting the orientation of the mechanical anisotropy of elements near the edge of the lamina to properly represent the tethering of the peripheral laminar trabeculae into the scleral canal wall (Roberts, Grau et al., 2009; Roberts, Liang et al., In Press (Accepted July 2009)). Further, for some systems, accurate predictions of mechanical behavior requires specimen-specific mechanical properties (Humphrey and Na, 2002; Viceconti, Davinelli et al., 2004; Vorp and Vande Geest, 2005; Vande Geest, Sacks et al., 2006; Girard, Suh et al., 2009). Generalizing and parameterizing these properties remains a challenge.
A useful property of the morphing techniques in sensitivity analysis is that the models produced are relatively easy to compare with one another. All the models are related to the baseline, sharing nodes and connectivity on the surface, and it is straightforward to transfer the boundary conditions from one model to another. This may be particularly useful for studies with more complex boundary conditions. In addition, models produced through morphing are related to each other by a known transformation. This transformation may be used to identify corresponding regions for analysis, which in turn allows a more direct quantitative comparison of the mechanical response. These can be considerable time savers in pre- and post-processing.
Successful morphing required a smoothing function to reduce discontinuities in the deformation field. The sensitivity of a particular model to these discontinuities depends on the details of the geometry, the quality of the mesh, and the magnitude and direction of the deformation field. For example, when elements pare highly elongated (high aspect ratio), the mesh is more sensitive to deformations that increase the aspect ratio such as compressive deformations on the short axis. Hence, it is not possible to prescribe a priori a degree of smoothness that will be satisfactory in all cases. Therefore, all deformation fields need to be tested over the ranges of parameters and models on which they will be used. We acknowledge that the smoothing functions used in the example are somewhat arbitrary and that their applicability may be limited to the particular cases shown. We explored other functions such as normal distributions, but the long tails were problematic. The trigonometric functions we chose were relatively simple because their magnitudes, and those of their derivatives, were clear and easy to scale. From an implementation standpoint, considerable effort was initially required to identify simple and useful morphing techniques and develop the scripts and modules which support them. Once the scripts were completed, however, the process was almost completely automated, requiring only minimal user intervention. The difficulties and time requirements of applying the techniques to other systems will depend on their complexity, but it is generally much faster than producing new specimen-specific models. As an example, we have applied the morphing techniques to a femur (Figure 7), which required approximately 3 hours of setup time, after which the time required to generate morphed surfaces was negligible.
We have demonstrated morphing applied to model surfaces, which required the internal volume of each geometry to be remeshed. When the morphing deformation vector fields are smooth, and the baseline mesh is of high quality it is possible to morph the baseline volume mesh directly (Sigal, Hardisty et al., 2008). Direct morphing of the volume mesh simplifies pre- and post-processing even further.
The methods described herein have some limitations that deserve consideration. First, there is some degree of arbitrariness in the deformation vector fields, even though our choices for variation were informed by an understanding of the anatomy and biomechanics of the structure of interest. The resulting morphed models agree with the variability of the anatomy and biomechanics of the posterior pole of the eye (Sigal, 2009; Sigal, Flanagan et al., 2009; Yang, Downs et al., 2009). Recent reports have shown that the three geometric properties varied in this study can vary considerably between individuals and also change with aging and disease (Roberts, Grau et al., 2009; Yang, Downs et al., 2009; Yang, Downs et al., 2009). Ideally, the shape variations introduced through morphing would be informed by a population study establishing the nature and magnitude of the physiologic variations (Van Essen, 2005; Laliberte, Meunier et al., 2007). However, we believe that the relationship between morphometry and analysis of sensitivity to variations in shape works both ways. Morphometry informs sensitivity analysis to keep the results relevant, and in turn, sensitivity analysis helps focus morphometry by identifying the shape variations of biomechanical consequence. A second potential limitation is that the morphing was applied to the baseline model. Hence, any inherent problems with the original geometry would propagate to the morphed models. In addition, the deformations obtained with morphing were not always ideal. For example, reducing the canal size produced a small distortion of the posterior lamina cribrosa surface near the insertion into the sclera (Figure 4, middle left panel). Similarly, zero displacement in the Y direction when morphing the size of the canal changes slightly the shape of the peripapillary sclera. Distortions such as these can be avoided using more complex functions. We evaluated several such functions and found their contributions to the biomechanics to be minimal (results not shown). However, these complicate the description of the method, and are therefore not included in this proof-of-concept work. Morphing of independent factors lends itself naturally to multivariate analysis. We are currently developing the scripts to couple the morphing techniques with variations in tissue material properties and loading for use within factorial and response surface experimental designs, which will be the subject of future reports.
Some authors refer to morphing as warping (Zöckler, Stalling et al., 2000). Variations of these techniques have been used for rapid reconstruction of specimen-specific models (Fernandez, Ho et al., 2005; Brock, Dawson et al., 2006; Sigal, Hardisty et al., 2008), nonlinear strain computation (Veress, Weber et al., 2002; Phatak, Sun et al., 2007), medical image registration (Todd-Pokropek, 2002) and segmentation (Bowden, Rabbitt et al., 1998), and to make up for the sparsity of data in low quality datasets (Blanz, Mehl et al., 2004; Shim, Pitto et al., 2007). Morphing is also popular in the animation and computer graphics community (Sederberg and Parry, 1986).
We have introduced morphing methods to parameterize specimen-specific models, and demonstrated their application on models of the posterior pole of the eye and the femur. The originality of this work lies in the application of morphing techniques for parametric analysis suitable for FE modeling and sensitivity analysis. While none of the concepts of morphing, specimen-specific modeling, or sensitivity analysis are novel by themselves, the integration of the three techniques shows promise for the study of biomechanics.
Supported in part by NIH-BRIN/INBRE grant P20 RR16456 and NIH grant R01 EY018926
We thank Dr. Claude Burgoyne for access to the data from which the baseline geometry was reconstructed. We also thank Jonathan Grimm, Juan Reynaud and Dr Cari Whyne.
Disclosures: - All authors were fully involved in the study and preparation of the manuscript.
- The material has not been, and will not be, published or submitted for publication elsewhere.
- Proprietary Interest: None
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.