UEF//eRepository
DSpace https://erepo.uef.fi
Rinnakkaistallenteet Luonnontieteiden ja metsätieteiden tiedekunta
2020
Discrete element and finite element
methods provide similar estimations for hip joint contact mechanics during
walking gait
Li, Mao
Elsevier BV
Tieteelliset aikakauslehtiartikkelit
© 2020 Elsevier Ltd.
CC BY-NC-ND https://creativecommons.org/licenses/by-nc-nd/4.0/
http://dx.doi.org/10.1016/j.jbiomech.2020.110163
https://erepo.uef.fi/handle/123456789/24435
Downloaded from University of Eastern Finland's eRepository
Journal Pre-proofs
Discrete Element and Finite Element Methods Provide Similar Estimations for Hip Joint Contact Mechanics during Walking Gait
Mao Li, Mikko S. Venäläinen, Shekhar S. Chandra, Rushabh Patel, Jurgen Fripp, Craig Engstrom, Rami K. Korhonen, Juha Töyräs, Stuart Crozier
PII: S0021-9290(20)30587-X
DOI: https://doi.org/10.1016/j.jbiomech.2020.110163
Reference: BM 110163
To appear in: Journal of Biomechanics Received Date: 7 August 2020
Revised Date: 7 November 2020 Accepted Date: 25 November 2020
Please cite this article as: M. Li, M.S. Venäläinen, S.S. Chandra, R. Patel, J. Fripp, C. Engstrom, R.K. Korhonen, J. Töyräs, S. Crozier, Discrete Element and Finite Element Methods Provide Similar Estimations for Hip Joint Contact Mechanics during Walking Gait, Journal of Biomechanics (2020), doi: https://doi.org/10.1016/
j.jbiomech.2020.110163
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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.
© 2020 Published by Elsevier Ltd.
Discrete Element and Finite Element Methods Provide Similar Estimations for Hip Joint Contact Mechanics during Walking
Gait
Mao Lia*, Mikko S. Venäläinena, b*, Shekhar S. Chandraa, Rushabh Patelc, Jurgen Frippa, d, Craig Engstrome, Rami K. Korhonenf, Juha Töyräsa, Stuart Croziera
aSchool of Information Technology and Electrical Engineering, University of Queensland, Brisbane, Australia
bTurku Bioscience Centre, University of Turku and Åbo Akademi University, Turku, Finland
cSchool of Mechanical and Mining Engineering, University of Queensland, Brisbane, Australia
dThe Australian e-Health Research Centre, CSIRO Health and Biosecurity, Brisbane, Australia
eSchool of Human Movement Studies, University of Queensland, Brisbane, Australia
fDepartment of Applied Physics, University of Eastern Finland, Kuopio, Finland
*Contributed equally
Correspondence: Mikko S. Venäläinen
School of Information Technology and Electrical Engineering, University of Queensland, Brisbane, Queensland 4072, Australia. Email: m.venalainen@uq.edu.au
Original Article
Number of words: 3, 897
Abstract
Finite element analysis (FEA) provides a powerful approach for estimating the in-vivo loading characteristics of the hip joint during various locomotory and functional activities. However, time-consuming procedures, such as the generation of high-quality FE meshes and setup of FE simulation, typically make the method impractical for rapid applications which could be used in clinical routine. Alternatively, discrete element analysis (DEA) has been developed to quantify mechanical conditions of the hip joint in a fraction of time compared to FEA.
Although DEA has proven effective in the estimation of contact stresses and areas in various complex applications, it has not yet been well characterised by its ability to evaluate contact mechanics for the hip joint during gait cycle loading using data from several individuals. The objective of this work was to compare DEA modelling against well-established FEA for analysing contact mechanics of the hip joint during walking gait. Subject-specific models were generated from magnetic resonance images of the hip joints in five asymptomatic subjects. The DEA and FEA models were then simulated for 13 loading time-points extracted from a full gait cycle. Computationally, DEA was substantially more efficient compared to FEA (simulation times of seconds vs. hours). The DEA and FEA methods had similar predictions for contact pressure distribution for the hip joint during normal walking. In all 13 simulated loading time-points across five subjects, the maximum difference in average contact pressures between DEA and FEA was within ±0.06 MPa. Furthermore, the difference in contact area ratio computed using DEA and FEA was less than ±6%.
Keywords
Hip joint, Contact mechanics, Discrete element, Finite element, Walking gait
1. Introduction
The hip joint is commonly affected by osteoarthritis (OA), a condition associated with variable levels of chronic pain, reduced mobility and physical impairment (Lespasio et al., 2018; March et al., 2014). Clinical and epidemiological studies have indicated that abnormal mechanical conditions within the hip joint can predispose to the development and progression of OA (Felson, 2013; Guilak, 2011). Research involving quantitative analyses of contact mechanics for the hip joint, such as contact pressures and areas, can provide important insights into the understanding and prediction of joint degeneration and OA (Abraham et al., 2017; Ateshian et al., 2015).
Extensive experimental studies have enhanced understanding of the mechanical function of the hip joint using either cadavers and pressure-sensitive film or piezo-resistive sensors to measure static contact forces (Anderson et al., 2008; Dienst et al., 2002; Townsend et al., 2018;
von Eisenhart et al., 1999), or in-vivo measurements with instrumented prostheses to directly measure contact pressures (Bergmann et al., 2016; Bergmann et al., 2018). Although cadaver experiments provide valuable information on contact pressures and areas these are only available in limited areas of the hip joint (e.g., large coverage leads to excessive overlap and crinkle artefacts) (Anderson et al., 2008). In contrast, in-vivo approaches allow detailed measurement of subject-specific mechanical parameters over the entire surface of the instrumented prostheses under various locomotory and functional activities (Bergmann et al., 2016). However, these in-vivo measurements are highly invasive and only realistic in patients, under extremely controlled conditions, scheduled for hip joint surgery (e.g., arthroplasty).
Finite element analysis (FEA) has been widely used for non-invasive, in-vivo analyses of contact mechanics of the hip joint (Abraham et al., 2017; Anderson et al., 2008; Ateshian et al., 2015; Bachtar et al., 2006; Henak et al., 2013; Vafaeian et al., 2017). The implementation of
subject-specific FEA, however, involves specialised procedures related to segmentation of radiological images and mesh generation to reach a stable numerical solution (Li et al., 2015).
The rapid evolution of computing capabilities and the development of novel approaches, i.e.
total Lagrangian explicit integration dynamic relaxation and meshless, have improved the efficiency of computational solutions (Horton et al., 2010; Joldes et al., 2010; Li et al., 2014a;
Li et al., 2016; Li et al., 2014b; Miller et al., 2007). However, analyses of hip biomechanics using FE procedures are typically still labour-intensive and computationally expensive (Vafaeian et al., 2017).
Discrete element analysis (DEA) (Abraham et al., 2013; Anderson et al., 2008; Chao et al., 2010; Genda et al., 2001; Townsend et al., 2018) has been used successfully in analyses of joint contact mechanics and, like FEA, provides an in-vivo computational technique for dedicated investigations of human joints (Abraham et al., 2013; Thomas-Aitken et al., 2018; Townsend et al., 2018; Vafaeian et al., 2017). Whilst being a compromise in accounting for the joint sliding contact mechanisms between articular surfaces, the DEA method has advantages such as the computational efficiency and easy generation of computational grid making this approach an attractive tool for clinical use (Townsend et al., 2018).
The DEA method and its variants have been shown to provide accurate estimations for trends of contact pressures and areas, and contact pressure distributions during daily activities, such as level walking, ascending and descending stairs (Abraham et al., 2013; Anderson et al., 2010; Benemerito et al., 2020; Chao et al., 2010). Despite the computationally efficient estimation of contact parameters and the generally good agreement between FEA (Anderson et al., 2010; Chao et al., 2010), DEA has been also reported to result in slight overestimation of absolute magnitudes of contact stresses in the hip joint compared to FEA (Abraham et al., 2013). Furthermore, the ability of DEA to evaluate contact parameters in the hip joint during gait cycle loading has not yet been well characterised using data from several individuals. The
aim of the this study was to analyse contact mechanics for five asymptomatic hip joints under loading conditions of walking gait using both DEA and FEA methods. The subject-specific models were based on automated 3D segmentation from magnetic resonance (MR) images.
The modelling results at 13 different loading configurations during the gait cycle were compared in terms of maximum and average contact pressures, contact pressure distributions, and contact area ratio.
2. Methods
2.1. Hip joint images and cartilage segmentation
To demonstrate the potential of our in-house DEA implementation, MR images of the hip joint from five asymptomatic subjects were acquired (3T Magnetom whole-body system, Siemens Healthcare, Germany) with a true fast imaging with steady-state free precession sequence (Xia et al., 2014). All images were obtained with a 150 mm field of view and voxel sizes ranged from 0.63 mm×0.63 mm×0.63 mm (subject 1), 0.23 mm×0.23 mm×0.49 mm (subject 2), and 0.78 mm×0.78 mm×1.00 mm (subjects 3-5). Informed consent was obtained from all participants. The research was approved by the University of Queensland Human Research Committee (approval number: 2018000405).
The hip joint cartilage was segmented from the MR images using the automated methods developed previously in our group based on 3D statistical shape modelling approaches (Chandra et al., 2016; Chandra et al., 2014; Xia et al., 2013). The windowed sinc algorithm was adopted for further smoothing the surfaces of automated segmentation results (Taubin et al., 1996).
The FEA method requires the geometry (i.e., the 3D cartilage surface model created from segmentation of the MR images) to be spatially discretised using volume elements. In this work, FE models were constructed from automated segmentation results using a custom MATLAB (Mathworks, Natick, MA) script, in which the segmented cartilage geometries were meshed with quadratic tetrahedral elements using the built-in generateMesh function from the Partial Differential Equation Toolbox. This element type has been reported to perform similarly as hexahedral elements typically preferred in contact analysis (Tadepalli et al., 2011). To quantitatively compare capabilities of DEA and FEA methods for analysing hip contact mechanics, the nodes and elements (i.e., triangular cells on articular cartilage surfaces) used by DEA were extracted from corresponding FEA models. As a result, comparisons could be conducted based on the same articular cartilage geometries, loading data and boundary conditions. The DEA method, however, does not inherently require the generation of FE meshes as the 3D surface model with triangular or rectangular elements is sufficient. An example 3D surface model (smoothed) of acetabular cartilage constructed using automated segmentation results from MR images of the (right) hip joint is shown in Fig. 1, along with sub-divisions into anterior, middle and posterior regions.
Suggested position of Fig. 1 2.3. Subject-specific DEA model
In DEA modelling, a spring array over the subchondral surfaces of the joint is created to represent articular cartilage. For each spring in the hip joint, one terminal is attached at the centre of each triangular cell of the femoral surface, and the other is then determined by projecting a point along the normal of the attached cell onto the opposing acetabular surface.
The spring length is calculated by summing the total cartilage (i.e., both acetabular and femoral) thickness and the gap between these opposed articular surfaces. For a specific loading condition, the developed search and optimisation algorithms are utilised to 1) determine the attachment
point of each spring on the triangular cells of acetabular surface, 2) measure the length for each spring and 3) calculate compression forces for those spring elements being compressed.
The spring stiffness 𝐾 is a function of the material properties and geometry of articular cartilage of the hip joint, which can be determined as follows:
𝐾𝑖 = 𝐸(1 − 𝜐)𝑆𝑖
(1 − 2𝜐)(1 + 𝜐)ℎ𝑖𝑐𝑖 (1)
where 𝑖 represents the 𝑖𝑡ℎ spring, 𝐸 is Young’s modulus, 𝜐 is Poisson’s ratio, 𝑆𝑖 is the area of the triangular cell where the spring attached, ℎ𝑖 is the sum of thickness of both acetabular and femoral cartilage, and 𝑐𝑖 is a non-linear factor (Genda et al., 2001; Yoshida et al., 2006). In previous studies, the non-linear behaviour of articular cartilage was approximated using a set of parameters such as the axial Cauchy stress, axial stretch and material properties (Ateshian et al., 1997; Huang et al., 2005; Volokh et al., 2007). In these studies, however, a (pre-set) uniform thickness was assumed for the hip cartilage layers rather than subject-specific geometries. In the current study, we modified the non-linear factor by taking subject-specific cartilage geometry information into account as follows:
𝑐𝑖 = 𝛼𝑑𝑖/ℎ𝑖 (2)
where 𝑑𝑖 is the compression distance of the 𝑖𝑡ℎ spring, 𝛼 is an empirical constant (𝛼 = 1.25) in this work. To address the convergence issue, the non-linear factor 𝑐𝑖 defined in Eq. (2) will impose a penalty function to penetration by exponentially increasing the stiffness of individual spring with the increment of compression distance 𝑑𝑖.
For a specific loading condition, the generated DEA model was solved iteratively using Newton's method to search global optima by updating the position of the femoral head towards the articular surface of acetabular cartilage. Eventually, the acetabulum-femur contact system
was in mechanical equilibrium and residual forces (i.e., differences between the applied forces and the sum of spring forces) were minimised.
The spring forces were calculated at projection points on the acetabular surface. Once the acetabulum-femur contact system was in equilibrium, the nearest neighbour interpolation algorithm was adopted to calculate nodal forces using the computed spring forces on the neighbouring triangular facets (Ni and Nguyen, 2009). The DEA simulation was an in-house implementation in MATLAB.
2.4. FE simulation
In FE simulation, cartilage to cartilage contact was assumed to be frictionless and modelled using a hard pressure-overclosure relationship and surface-to-surface discretisation. For each time point, an Abaqus input file with appropriate contact definitions, material properties, boundary and loading conditions was generated, and the created FE models were solved using Abaqus/Standard 3DEXPERIENCE R2019x (Dassault Systèmes SIMULIA Corp., Johnston, RI, USA). A diagram illustrating workflows for both DEA and FEA in subject-specific biomechanical analyses is shown in Fig. 2. The time-consuming and computation-intensive procedures involved in FEA, such as FE mesh generation and FE simulation, are avoided in DEA.
Suggested position of Fig. 2
2.5. Loading and boundary conditions
The hip joint loading data were extracted from published in-vivo experimental results (Bergmann et al., 2016). In a previous study (Abraham et al., 2013), the acetabulum was assumed fixed; while the femoral head was allowed to translate in all axes (i.e., x, y and z) with
restrained rotation, which could describe joint postures in static loading scenarios. However, this work aims to analyse hip contact biomechanics under dynamic loading (walking gait), whereby the femoral head is free to rotate (i.e., flexion or extension) in the sagittal plane with respect to a transverse axis (Lewis et al., 2010). The full gait cycle extracted from continuous walking was divided into 13 time-points in the current work, covering gait events such as heel- strike (time-point 2), mid-stance (time-point 4), toe-off (time-point 7) and terminal-swing (time-point 11) (Krebs et al., 1998).
2.6. Material properties
Since the main objective of this paper is to evaluate the capability of the subject-specific DEA method for quantifying mechanical responses of the hip joint cartilage during walking in comparison with well-established FEA, a simple elastic material model was utilised. The Young’s modulus and Poisson’s ratio for the hip articular cartilage used in this examination were 10 MPa and 0.43 (Abraham et al., 2013; Richard et al., 2013; Shepherd and Seedhom, 1999), respectively, to represent the dynamic mechanical properties of cartilage under short- term loading conditions.
3. Results
3.1. Contact pressure
Overall, the contact pressures estimated using DEA and FEA were highly similar throughout the gait for all subjects. Illustrative examples of obtained contact pressure distributions and their element-wise differences for subjects 1 and 3 are provided in Figs. 3 and 4, respectively, whereas the results for remaining subjects (2, 4 and 5) can be found in Supplementary Figs. 1-
3. Despite the good agreement in contact pressures in the more central contact regions, some local differences between DEA and FEA were observed particularly at the periphery of contact.
Suggested position of Fig. 3 Suggested position of Fig. 4
The estimates of maximum contact pressures obtained using DEA and FEA were also well in line with each other (Fig. 5). The greatest differences between the two methods were observed for subjects 1 - 3 during the second peak load (from Mid-stance to Toe-off) where the maximum contact pressures obtained using DEA were approximately 1.6 MPa higher compared to FEA.
Furthermore, for subject 3, the maximum contact pressures estimated by DEA were slightly higher than FEA predictions from Heel-strike to Toe-off of the gait cycle, and slightly lower than the FEA results in the Terminal-swing phase. In other gait phases, the differences in maximum contact pressures between the two methods were only modest for all subjects and on an average of 0.28 MPa, see Fig. 5.
Suggested position of Fig. 5
For all subjects, DEA and FEA produced similar spatial averages of contact pressures throughout the gait cycle (Fig. 6). The spatial averages were calculated as the average of contact pressures on all the surface faces in contact (i.e. pressure >0) at each individual time point.
Variations in the average contact pressures between DEA and FEA among the five subjects were further analysed at several phases of the walking gait cycle, i.e., heel-strike, mid-stance, toe-off and terminal-swing (Figs. 7). Differences in average contact pressures for the DEA and FEA data for the 13 loading time-points were quantitatively evaluated in terms of absolute value and percentage (Fig. 8). The maximum (absolute and percentage) difference in average contact pressure among the five subjects between DEA and FEA ranged from -0.05 MPa to 0.06 MPa, and -6% to 10%, respectively.
Suggested position of Fig. 6 Suggested position of Fig. 7 Suggested position of Fig. 8
3.2. Contact area ratio
Contact area ratio was calculated as a percentage of the total area of the acetabulum cartilage surface. For all 13 loading time-points, all triangular cells on the articular surface being contacted (i.e., contact pressure > 0) were extracted for calculating contact area ratio. Fig. 9 compares the contact area ratio among the 13 time-points of a gait cycle for the five subjects, and the results show good agreement between DEA and FEA. Differences in the contact area ratio were further assessed and shown as box plots (Fig. 10), which indicate that the maximum difference in contact area ratio computed for the 13 time-points of gait using DEA and FEA is
±6%.
Suggested position of Fig. 9 Suggested position of Fig. 10 3.3. Computation speed
To compare the computation speed, the computational times for the two methods were measured and scaled to per CPU core. Overall, the DEA method achieved ~300 times higher computational efficiency over the conventional FEA approach, i.e., the average computation time for each DEA model was ~30 seconds (per CPU core), while FEA took about 150 minutes (per CPU core).
4. Discussion
A comparison study between DEA and FEA methods for analysing contact mechanics of the hip joint during walking gait was conducted. Subject-specific biomechanical models were generated from MR images of hip joint cartilage in asymptomatic subjects. The results showed a consistent agreement between these two methods in the context of local contact pressure (i.e., peak magnitude and distribution pattern), average contact pressure, and contact area ratio.
The local contact pressure on the articular surface of the hip joint plays an important role on assessment of risk in the development of OA (Daniel et al., 2006). Therefore, the capability to predict the peak magnitude and distribution pattern of contact pressures in joints is critical for DEA to be routinely used. The results (Figs. 3 - 5), indicate that DEA provides estimates of hip contact mechanics during gait cycle loading similar to FEA. A recent cadaver experiment on the hip joint reported that the highest contact pressure during gait was found in the posterior- middle and anterior areas (of the acetabulum) at heel-strike gait phase (Townsend et al., 2018).
This is in good agreement with present simulation results (Figs. 3 - 5).
The maximum differences in contact pressures estimated using DEA and FEA were observed at the periphery of contact (i.e., up to 5.5 MPa in Fig. 4). This might be explained by differences in contact detection between the two methods as considerable disagreement in estimated contact pressures may be generated if one method detects contact whereas the other does not. However, this affects only a small number of peripheral elements and the differences in absolute terms are still below the maximum contact pressures in the loading scenario, and therefore does not lead substantial disagreement in the estimation of overall contact mechanics.
Furthermore, similarly, it has been also reported before that the largest discrepancies between DEA and FEA can be found at the periphery of contact (Kern and Anderson, 2015).
Based on literature the contact pressure varies more between individuals in a specific loading condition than between different loading conditions for a single subject (Harris et al., 2012). In the present study, however, contact pressure distributions were found to be relatively
similar for all subjects in a specific loading condition, and different loading for an individual produced greater variations in contact pressure distribution than that between different subjects.
This discrepancy might arise from how biomechanical models were generated. Harris et al used computed tomography (CT) images (Harris et al., 2012), while this paper utilised MR imaging enabling much better soft tissue contrast. Therefore, more accurate geometrical information, in particular for cartilage, could be extracted from MR images to generate subject-specific models.
However, the quantitative influences of geometrical deviations (i.e., caused by image segmentation and mesh generation errors) on analyses of hip joint contact mechanics have not been sufficiently investigated, and will be further quantified in our future work.
The average contact pressures obtained using DEA and FEA were similar (Fig. 6), e.g., heel-strike (DEA: 1.17 ± 0.23 MPa; FEA: 1.20 ± 0.21 MPa) and mid-stance (DEA: 0.97 ± 0.18 MPa; FEA: 0.98 ± 0.16 MPa). These results correspond well to a previous study in which 16 subject-specific FEA models were generated from a variety of healthy (hip) joints, and the average contact pressures were reported to be 1.08 ± 0.32 MPa (heel-strike) and 0.94 ± 0.23 MPa (mid-stance), respectively (Harris et al., 2012).
Absolute average contact pressure estimated with DEA was slightly lower than the estimation with FEA at heel-strike in subjects 1, 2 and 5, and slightly higher than that estimated with FEA from mid-stance to toe-off in subjects 2 and 4. In the following loading scenarios (i.e., from toe-off to terminal-swing), however, differences in absolute average contact pressures were very small (< 0.02 MPa) between these two methods (Fig. 3). This observation agrees with previous research reporting that the difference in contact pressures between DEA and FEA increased with the magnitudes of contact pressures (i.e., up to 40%) (Abraham et al., 2013). However, this difference was markedly lower in the current study and remained within 10% (Fig. 5), which might be attributed to the introduction of the nonlinear factor (Eq. (2)) that
contact pressures between DEA and FEA might be partially explained as different model representations were used in these two methods (Abraham et al., 2013), where DEA represented deformable cartilage tissues using one-dimensional springs rather than a continuum modelled explicitly in FEA.
During the full gait cycle, the average contact area for the five subjects estimated by DEA and FEA were similar, i.e., 2493 ± 147 𝑚𝑚2 and 2502 ± 128 𝑚𝑚2, respectively. These average values reported in the present study are similar to those reported in a 10-year follow-up (computational) study on 12 subjects after peri-acetabular osteotomy with an average contact area of 2381 ± 376 𝑚𝑚2 (Armiger et al., 2009). These average contact areas, however, are higher than those reported by Abraham et al., where contact areas computed using DEA and FEA were 1200 𝑚𝑚2 and 1300 𝑚𝑚2, respectively (Abraham et al., 2013). This, however, might be explained by the fact that contact mechanics varies between individuals based on subject-specific joint geometry and only one subject was evaluated in the previous study.
The variation of contact area ratio on the hip joint articular surface during the gait cycle ranged from 52% (mid-swing) to 98% (mid-stance) across the five subjects in the current study.
This is in line with earlier research reporting that contact area ratio of the hip joint during walking gait varied between 72% and 100% (Yoshida et al., 2006). However, it has also been reported that contact area ratio at heel-strike of the gait cycle was 36% (Harris et al., 2012), which was much lower compared to those reported by Yoshida et al (~75% at heel-strike) (Yoshida et al., 2006), and based on our computational analyses, i.e., 86 ± 2.5% (DEA) and 87
± 4.6% (FEA) (Fig. 8). Hence, further experimental and computational studies are necessary.
Although a consistent agreement between presented DEA and well-established FEA for analysing hip contact mechanics in a full walking gait cycle was achieved, this work still has limitations. 1) Whilst DEA and FEA models were generated based on individual geometries, in particular subject-specific cartilage segmented from MR images, the lack of exact body mass
for each subject limited the computation accuracy, especially for absolute values of contact pressures and contact areas. 2) During gait, the hip flexion-extension angle varies continuously.
However, subject-specific hip angles were not available for the present examination. Instead, the same flexion-extension angle related to a specific loading condition was applied to all subjects (Lewis et al., 2010). Such generic loading data (i.e., body mass and joint angles mentioned above) cannot accurately represent the subject-specific mechanical factors. Harris et al have pointed out that subject-specific biomechanical models loaded by generic data will lead to more uniform contact pressures and area distribution patterns (Harris et al., 2012). 3) The five hip joints from asymptomatic subjects in 13 loading scenarios extracted from an entire gait cycle enabled us to conclude the key trends in hip contact mechanics computed using DEA and FEA methods. The comparison results in the present study are important for implementing the DEA method for analysing contact mechanics of the pathological joints, where changes in volume and/or thickness of cartilage caused by joint diseases (i.e., OA) might introduce larger difference between DEA and FEA.
Despite the fact that in general FEA is substantially more time-consuming approach than DEA, it has been demonstrated that with atlas-based FEA, the time needed to generate and run the models can be reduced down to a few minutes (Mononen et al., 2019). Contrary to DEA, atlas-based FEA also enables the estimation of local stresses and strains acting on different tissue constituents and can therefore be used for investigating the mechanisms leading to the development of OA. By using the results of DEA as an input for FEA, however, it might be possible to perform these analyses with minimum computational cost even without the need for an atlas of tissue geometries, therefore bringing the computational evaluation of joint function even closer to clinical applications.
DEA has proven successful for assessing hip contact mechanics and providing comparable
presented a subject-specific DEA method with the capability for predicting not only trends but also absolute values for contact pressures and contact areas. Importantly, DEA modelling is computationally efficient compared to the conventional FEA approach. This together with the present promising results indicate potential for implementing subject-specific DEA biomechanical analyses to routine diagnosis and treatment assessment.
Acknowledgements
This research was supported under Australian National Health and Medical Research Council Development Grant: APP1139868. Dr. Venäläinen received funding from the Academy of Finland (grant number 322123) during the course of the study. Centre for Science, Finland, is acknowledged for providing computational resources and software for FEA simulations.
Conflict of interest statement
No.
Reference:
Abraham, C.L., Knight, S.J., Peters, C.L., Weiss, J.A., Anderson, A.E., 2017. Patient-specific chondrolabral contact mechanics in patients with acetabular dysplasia following treatment with peri-acetabular osteotomy. Osteoarthr Cartilage 25, 676-684.
Abraham, C.L., Maas, S.A., Weiss, J.A., Ellis, B.J., Peters, C.L., Anderson, A.E., 2013. A new discrete element analysis method for predicting hip joint contact stresses. J Biomech 46, 1121-1127.
Anderson, A.E., Ellis, B.J., Maas, S.A., Peters, C.L., Weiss, J.A., 2008. Validation of finite element predictions of cartilage contact pressure in the human hip joint. J Biomech Eng-T Asme 130.
Anderson, D.D., Iyer, K.S., Segal, N.A., Lynch, J.A., Brown, T.D., 2010. Implementation of Discrete Element Analysis for Subject-Specific, Population-Wide Investigations of Habitual Contact Stress Exposure. J Appl Biomech 26, 215-223.
Armiger, R.S., Armand, M., Tallroth, K., Lepisto, J., Mears, S.C., 2009. Three-dimensional mechanical evaluation of joint contact pressure in 12 periacetabular osteotomy patients with 10-year follow- up. Acta Orthop 80, 155-161.
Ateshian, G.A., Henak, C.R., Weiss, J.A., 2015. Toward patient-specific articular contact mechanics. J Biomech 48, 779-786.
Ateshian, G.A., Warden, W.H., Kim, J.J., Grelsamer, R.P., Mow, V.C., 1997. Finite deformation biphasic material properties of bovine articular cartilage from confined compression experiments.
J Biomech 30, 1157-1164.
Bachtar, F., Chen, X., Hisada, T., 2006. Finite element contact analysis of the hip joint. Medical &
Biological Engineering & Computing 44, 643-651.
Benemerito, I., Modenese, L., Montefiori, E., Mazza, C., Viceconti, M., Lacroix, D., Guo, L.Z., 2020.
An extended discrete element method for the estimation of contact pressure at the ankle joint during stance phase. P I Mech Eng H 234, 507-516.
Bergmann, G., Bender, A., Dymke, J., Duda, G., Damm, P., 2016. Standardized Loads Acting in Hip Implants. PloS one 11.
Bergmann, G., Kutzner, I., Bender, A., Dymke, J., Trepczynski, A., Duda, G.N., Felsenberg, D., Damm, P., 2018. Loading of the hip and knee joints during whole body vibration training. PloS one 13.
Chandra, S.S., Surowiec, R., Ho, C., Xia, Y., Engstrom, C., Crozier, S., Fripp, J., 2016. Automated Analysis of Hip Joint Cartilage Combining MR T2 and Three-Dimensional Fast-Spin-Echo Images.
Magnet Reson Med 75, 403-413.
Chandra, S.S., Xia, Y., Engstrom, C., Crozier, S., Schwarz, R., Fripp, J., 2014. Focused shape models for hip joint segmentation in 3D magnetic resonance images. Med Image Anal 18, 567-578.
Chao, E.Y., Volokh, K.Y., Yoshida, H., Shiba, N., Ide, T., 2010. Discrete element analysis in musculoskeletal biomechanics. Mol Cell Biomech 7, 175-192.
Daniel, M., Herman, S., Dolinar, D., Iglic, A., Sochor, M., Kralj-Iglic, V., 2006. Contact stress in hips with osteonecrosis of the femoral head. Clinical orthopaedics and related research, 92-99.
Dienst, M., Seil, R., Godde, S., Brang, M., Becker, K., Georg, T., Kohn, D., 2002. Effects of traction, distension, and joint position on distraction of the hip joint: An experimental study in cadavers.
Arthroscopy-the Journal of Arthroscopic and Related Surgery 18, 865-871.
Felson, D.T., 2013. Osteoarthritis as a disease of mechanics. Osteoarthr Cartilage 21, 10-15.
Genda, E., Iwasaki, N., Li, G.A., MacWilliams, B.A., Barrance, P.J., Chao, E.Y.S., 2001. Normal hip joint contact pressure distribution in single-leg standing - effect of gender and anatomic parameters.
J Biomech 34, 895-905.
Guilak, F., 2011. Biomechanical factors in osteoarthritis. Best Pract Res Cl Rh 25, 815-823.
Harris, M.D., Anderson, A.E., Henak, C.R., Ellis, B.J., Peters, C.L., Weiss, J.A., 2012. Finite element prediction of cartilage contact stresses in normal human hips. J Orthop Res 30, 1133-1139.
Henak, C.R., Carruth, E.D., Anderson, A.E., Harris, M.D., Ellis, B.J., Peters, C.L., Weiss, J.A., 2013.
Finite element predictions of cartilage contact mechanics in hips with retroverted acetabula.
Osteoarthr Cartilage 21, 1522-1529.
Horton, A., Wittek, A., Joldes, G.R., Miller, K., 2010. A meshless Total Lagrangian explicit dynamics algorithm for surgical simulation. Int J Numer Meth Bio 26, 977-998.
Huang, C.Y., Stankiewicz, A., Ateshian, G.A., Mow, V.C., 2005. Anisotropy, inhomogeneity, and tension-compression nonlinearity of human glenohumeral cartilage in finite deformation. J Biomech 38, 799-809.
Joldes, G.R., Wittek, A., Miller, K., 2010. Real-time nonlinear finite element computations on GPU - Application to neurosurgical simulation. Comput Method Appl M 199, 3305-3314.
Kern, A.M., Anderson, D.D., 2015. Expedited patient-specific assessment of contact stress exposure in the ankle joint following definitive articular fracture reduction. J Biomech 48, 3427-3432.
Krebs, D.E., Robbins, C.E., Lavine, L., Mann, R.W., 1998. Hip biomechanics during gait. J Orthop Sports Phys Ther 28, 51-59.
Lespasio, M.J., Sultan, A.A., Piuzzi, N.S., Khlopas, A., Husni, M.E., Muschler, G.F., Mont, M.A., 2018.
Hip Osteoarthritis: A Primer. Perm J 22, 17-084.
Lewis, C.L., Sahrmann, S.A., Moran, D.W., 2010. Effect of hip angle on anterior hip joint force during gait. Gait Posture 32, 603-607.
Li, M., Miller, K., Joldes, G., Kikinis, R., Wittek, A., 2014a. Patient-Specific Meshless Model for Whole-Body Image Registration. Lect Notes Comput Sc 8789, 50-57.
Li, M., Miller, K., Joldes, G.R., Doyle, B., Garlapati, R.R., Kikinis, R., Wittek, A., 2015. Patient- specific biomechanical model as whole-body CT image registration tool. Med Image Anal 22, 22- 34.
Li, M., Miller, K., Joldes, G.R., Kikinis, R., Wittek, A., 2016. Biomechanical model for computing deformations for whole-body image registration: A meshless approach. Int J Numer Meth Bio 32.
Li, M., Wittek, A., Joldes, G., Zhang, G., Dong, F., Kikinis, R., Miller, K., Year Whole-Body Image Registration Using Patient-Specific Nonlinear Finite Element Model. In Computational Biomechanics for Medicine. New York, NY.
March, L., Smith, E.U.R., Hoy, D.G., Cross, M.J., Sanchez-Riera, L., Blyth, F., Buchbinder, R., Vos, T., Woolf, A.D., 2014. Burden of disability due to musculoskeletal (MSK) disorders. Best Pract Res Cl Rh 28, 353-366.
Miller, K., Joldes, G., Lance, D., Wittek, A., 2007. Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation. Commun Numer Meth En 23, 121-134.
Mononen, M.E., Liukkonen, M.K., Korhonen, R.K., 2019. Utilizing Atlas-Based Modeling to Predict Knee Joint Cartilage Degeneration: Data from the Osteoarthritis Initiative. Ann Biomed Eng 47, 813-825.
Ni, K.S., Nguyen, T.Q., 2009. An Adaptable k-Nearest Neighbors Algorithm for MMSE Image Interpolation. Ieee T Image Process 18, 1976-1987.
Richard, F., Villars, M., Thibaud, S., 2013. Viscoelastic modeling and quantitative experimental characterization of normal and osteoarthritic human articular cartilage using indentation. J Mech Behav Biomed 24, 41-52.
Shepherd, D.E.T., Seedhom, B.B., 1999. The 'instantaneous' compressive modulus of human articular cartilage in joints of the lower limb. Rheumatology 38, 124-132.
Tadepalli, S.C., Erdemir, A., Cavanagh, P.R., 2011. Comparison of hexahedral and tetrahedral elements in finite element analysis of the foot and footwear. J Biomech 44, 2337-2343.
Taubin, G., Zhang, T., Golub, G., Year Optimal surface smoothing as filter design. Berlin, Heidelberg.
Thomas-Aitken, H.D., Willey, M.C., Goetz, J.E., 2018. Joint contact stresses calculated for acetabular dysplasia patients using discrete element analysis are significantly influenced by the applied gait pattern. J Biomech 79, 45-53.
Townsend, K.C., Thomas-Aitken, H.D., Rudert, M.J., Kern, A.M., Willey, M.C., Anderson, D.D., Goetz, J.E., 2018. Discrete element analysis is a valid method for computing joint contact stress in the hip before and after acetabular fracture. J Biomech 67, 9-17.
Vafaeian, B., Zonoobi, D., Mabee, M., Hareendranathan, A.R., El-Rich, M., Adeeb, S., Jaremko, J.L., 2017. Finite element analysis of mechanical behavior of human dysplastic hip joints: a systematic review. Osteoarthr Cartilage 25, 438-447.
Volokh, K.Y., Chao, E.Y., Armand, M., 2007. On foundations of discrete element analysis of contact in diarthrodial joints. Mol Cell Biomech 4, 67-73.
von Eisenhart, R., Adam, C., Steinlechner, M., Muller-Gerbl, M., Eckstein, F., 1999. Quantitative determination of joint incongruity and pressure distribution during simulated gait and cartilage thickness in the human hip joint. J Orthop Res 17, 532-539.
Xia, Y., Chandra, S.S., Engstrom, C., Strudwick, M.W., Crozier, S., Fripp, J., 2014. Automatic hip cartilage segmentation from 3D MR images using arc-weighted graph searching. Physics in medicine and biology 59, 7245-7266.
Xia, Y., Fripp, J., Chandra, S.S., Schwarz, R., Engstrom, C., Crozier, S., 2013. Automated bone segmentation from large field of view 3D MR images of the hip joint. Physics in medicine and biology 58, 7375-7390.
Yoshida, H., Faust, A., Wilckens, J., Kitagawa, M., Fetto, J., Chao, E.Y.S., 2006. Three-dimensional dynamic hip contact area and pressure distribution during activities of daily living. J Biomech 39, 1996-2004.
Figure legends
Figure 1. An example 3D surface model of the acetabulum cartilage generated using automated segmentation of MR images of the (right) hip joint with regional definition of cartilage partitions.
Figure 2. Workflow of FEA and DEA approaches for subject-specific analysis of hip joint biomechanical function
Figure 3. Contact pressure distributions on the articular surface of the hip joint (subject 1) computed using DEA (left column) and FEA (middle column) methods. Patterns of contact pressure distribution are compared at heel-strike, mid-stance, toe-off and terminal-swing of a walking gait cycle. Differences in absolute values of contact pressures are evaluated on an element-wise basis in the right column, demonstrating an excellent agreement between the methods.
Figure 4. Contact pressure distributions on the articular surface of the hip joint (subject 3) computed using DEA (left column) and FEA (middle column) methods. Patterns of contact pressure distribution are compared at heel-strike, mid-stance, toe-off and terminal-swing of a walking gait cycle. Differences in absolute values of contact pressures are evaluated on an element-wise basis in the right column, demonstrating an excellent agreement between the methods.
Figure 5. A comparison of absolute values of peak contact pressures computed using DEA and FEA methods for each of the 13 time-points of a walking gait cycle among the five subjects (a – e).
Figure 6. A comparison of absolute values of average contact pressures computed using DEA and FEA methods for each of the 13 time-points of a walking gait cycle among the five subjects (a – e).
Figure 7. Variations of average contact pressure magnitudes determined using DEA (blue) and FEA (red) at heel-strike, mid-stance, toe-off and terminal-swing phases across the five subjects. The DEA shows slightly greater variability in contact pressure values than FEA only between heel-strike and mid-stance.
Figure 8. Differences in average contact pressures computed using DEA and FEA during walking gait for the five subjects. Both median (i.e., red solid line in each box) and average (i.e., diamond in each box) differences in computed average contact pressure values are within ±0.02 MPa (a) and ±4% (b).
Figure 9. A comparison of contact area ratio computed using DEA and FEA methods in the 13 time-points of a walking gait cycle for five subjects (a - e). Whilst DEA computation has some differences in contact area ratio magnitudes (subjects 2, 3 and 4) compared to FEA, trends in contact area variation are successfully predicted for all subjects.
Figure 10. Box plots of differences in contact area ratio between DEA and FEA across all 13 time-points among the five subjects. Both median (i.e. red solid line in box) and average (i.e., diamond in box) difference varies within ±2%, and the maximum difference (subject 4) ranges from -6.0% to 4.0%.
Figure 1 Click here to access/download;Figure;Fig. 1.tiff
Figure 2 Click here to access/download;Figure;Fig. 2.tiff
Figure 3 Click here to access/download;Figure;Fig. 3.tiff
Figure 4 Click here to access/download;Figure;Fig. 4.tiff
Figure 5 Click here to access/download;Figure;Fig. 5.tiff
Figure 6 Click here to access/download;Figure;Fig. 6.tiff
Figure 7 Click here to access/download;Figure;Fig. 7.tiff
Figure 8 Click here to access/download;Figure;Fig. 8.tiff
Figure 9 Click here to access/download;Figure;Fig. 9.tiff
Figure 10 Click here to access/download;Figure;Fig. 10.tiff
