Optimization of nonlinear hyperelastic coefficients for foot tissues using a magnetic resonance imaging deformation experiment

Abstract

Accurate prediction of plantar shear stress and internal stress in the soft tissue layers of the foot using finite element models would provide valuable insight into the mechanical etiology of neuropathic foot ulcers. Accurate prediction of the internal stress distribution using finite element models requires that realistic descriptions of the material properties of the soft tissues are incorporated into the model. Our investigation focused on the creation of a novel three-dimensional (3D) finite element model of the forefoot with multiple soft tissue layers (skin, fat pad, and muscle) and the development of an inverse finite element procedure that would allow for the optimization of the nonlinear elastic coefficients used to define the material properties of the skin muscle and fat pad tissue layers of the forefoot based on a Ogden hyperelastic constitutive model. Optimization was achieved by comparing deformations predicted by finite element models to those measured during an experiment in which magnetic resonance imaging (MRI) images were acquired while the plantar surface forefoot was compressed. The optimization procedure was performed for both a model incorporating all three soft tissue layers and one in which all soft tissue layers were modeled as a single layer. The results indicated that the inclusion of multiple tissue layers affected the deformation and stresses predicted by the model. Sensitivity analysis performed on the optimized coefficients indicated that small changes in the coefficient values (±10%) can have rather large impacts on the predicted nominal strain (differences up to 14%) in a given tissue layer.

Fig. 1

The foot-loading device used during collection of MRI deformation data. For a more thorough description, see Petre et al. [21].

Fig. 2

Top: MRI image sets collected in the unloaded state and with 12.5% and 50% of the subject′s body weight applied to the forefoot. Middle: The 3D nature of the data allows for visualization and analysis of the tissue deformation along any plane. Bottom: Frontal plane sections through the approximate metatarsal head region. Note that most of the deformation has already occurred when only one quarter of the total load (12.5% BW) was applied.

Fig. 3

(a) The outer surface of the 3D finite element forefoot mesh. (b) The FE mesh with the toes, skin, and muscle removed to display the bones and plantar fat. Note that the phalanges of each toe were fused for this study. (c) Joint definitions used to describe model bone kinematics. Open circles represent the locations of joint centers of rotation. The midfoot node was fixed in space throughout simulation of the loading experiment (see Fig. 1). (d) FE model of the MRI loading experiment incorporating the 3D, layered tissue FE mesh. The midfoot and cast were fixed during the simulation and forces were applied to the rigid plate, as in the MRI experiment.

Fig. 4

An exemplar slice from the undeformed MRI image set (a). The same slice showing the segmentation of structures; red denoting skin, blue, fat pad, and green, bone, where all remaining tissue was classified as muscle (b). The recolored ternary image based on the segmentation of structures (c). It should be noted that the tissue on the dorsal side of the foot was considered to be either muscle or skin only to ease in the construction of the FE mesh as this simplification did not impact the results of the FEM model prediction of plantar pressure.

Fig. 5

A schematic description of the mesh overlay process in the region of interest shown at the top. (a) Each node of the undeformed FE mesh is assigned a pixel intensity value corresponding to its location in the undeformed MRI image set. This pixel value is fixed to the node as it undergoes deformation during simulation of the experiment. The deformed mesh can be used to generate a (relatively coarse) mode-predicted deformed image that can be compared to the experimental MRI images on a node by node basis. (b) If the current material coefficients are; for example, stiffer than the real material, then the mesh will not deform as much as the tissue and the node intensity values will fail to correlate with the deformed MRI image set. (c) If the material properties governing node motion correctly match the tissue, then the deformed mesh will correlate well with the deformed image set. Although this process is shown here in 2D with 8 nodes, the technique was applied in 3D allowing for out-of-plane deformations of more than 17,000 nodes induced by nearly incompressible materials.

Fig. 6

Schematic diagram of the inverse finite element technique used to determine material coefficients. Inputs to the technique (at the left) are deformed and undeformed MRI image sets, undeformed FE mesh, and an initial guess for material coefficients. The mesh overlay process, used to determine the pixel value at each node, is further described in Fig. 5. The error calculation and its inputs are described with Eq. (3). Briefly, each vector I_n,f contains the corresponding image intensities for each of the n nodes at loading level f. Error is calculated by comparing predicted intensities at each of the loading levels to the intensities in the unloaded state I_n,0.

Fig. 7

Qualitative comparison of the deformations predicted in the region of the metatarsal heads when a compressive load equal to 50% body weight was applied to the model versus deformations observed in MRI images when the equivalent load was applied to the imaged forefoot

Fig. 8

Overall deformation (nominal strain, bottom row) of the forefoot at full 50% BW loading using optimized coefficients for layered tissue representation (best fit, see Table 1 for coefficients), along with hydrostatic pressure (middle row), and von Mises stress distribution (top row) within tissue layers. The loading and boundary conditions were representative of the MRI experimentation. The forefoot model was sectioned at an approximate level of metatarsal heads to emphasize a potential region of clinical interest. Distributions of all strain components illustrate the complex nature of the combined deformation modes in the foot. See text for further discussion.

Fig. 9

Clinically relevant mechanical variables were not highly sensitive to increased mesh density. While a coarser mesh was used for the identification of material coefficients (left: 26,335 elements, 29,697 nodes), a higher mesh density may be desirable for convergence (right: 39,675 elements, 43,893 nodes). Plantar pressures and von Mises stress distribution are plotted. Both the magnitude and location of peaks were similar between the predictions of coarse and fine meshes as obtained by using the best material coefficient fit (see Table 1 for coefficients). Loading and boundary conditions of these simulations approximated those of peak forefoot loading instant of the push-off phase.

Fig. 10

Comparison of clinically relevant mechanical variables predicted when tissue layer material properties were incorporated in the model (left) and when lumped material properties were used (right). Top row illustrates similar predicted plantar pressure distributions in both models. However, the lumped material model predicted higher magnitudes of pressure. Middle row illustrates the comparison of von Mises stresses in the region under the metatarsal heads. While stresses are concentrated under the third metatarsal head in both models, stress distribution is more uniform in the lumped model (right) than in the layered model (left) where reduced levels of stress are observed in the highly compliant fat pad tissue layer. Bottom row illustrates the comparison of maximum principal strains in the region under the metatarsal heads. Peak strain values were concentrated under the first metatarsal head in the layered model (Left) and under the third metatarsal head in the lumped model (right). Strain distribution was more uniform in the lumped model (right) than in the layered model where high strain values were observed throughout the highly compliant fat pad layer.