Clinical Contrast-Enhanced Computed Tomography With Semi-Automatic Segmentation Provides Feasible Input for Computational Models of the Knee Joint

Osteoarthritis (OA) is a common joint disease; the prevalence of knee OA has been estimated to be 3.8% of the world population [1], burdening the lives of over 260 million people. This painful and immobilizing joint disease may be initiated and worsened by abnormal joint loading, i.e., high impact loads [2] or repetitive joint loads generated during daily physical activities such as walking or running [3,4]. Unfortunately, the routine diagnostics of OA has mostly focused on the evaluation of symptoms or anatomic features, such as chronic pain or the decrease in joint space width, which are more characteristic of later stages of the disease. Therefore, if we wish to prevent the progression of the disease or to predict its onset, it would be advantageous to have more sophisticated imaging and diagnostic methods.

Biomechanical modeling can be used for the evaluation of joint mechanics under physiologically relevant loading conditions [5,6], providing diagnostically valuable information on the function of the joint in terms of tissue stresses and strains. Computational finite element (FE) models enable simulation of strains in the collagen fibrils and assessing progression of damage; these models have also been used to estimate cell death and changes in tissue structure around cartilage lesions as well as to predict the progression of OA [7–15]. Biomechanical analysis techniques are also relevant for the investigation of other normal and altered states in a knee joint, not only for investigation of cartilage damage and OA [16,17].

A major obstacle to the clinical application of modeling approaches and the quantitative analysis of knee joint tissues is that they demand information on subject-specific joint geometries. In most cases, the knee joint geometry for computational modeling is obtained via manual segmentation of the joint structures [18], which is a very time-consuming task to perform for each individual subject. Sophisticated techniques, such as statistical shape modeling and neural networks, have been exploited in the development of automatic segmentation methods [19–23]. These methods have proven to be superior to less complex techniques, for instance, those based on simple thresholding. Several automatic and semi-automatic segmentation methods have been introduced for segmentation of knee joint connective tissues from magnetic resonance and computed tomography (CT) images [20,22–25]. However, the feasibility of these methods has not typically been tested in computational modeling, possibly due to the laborious work required in transferring the generated segmentations into computationally sufficient and compatible forms.

Previously, a semi-automatic method was developed to segment contrast-enhanced CT images [25]. In contrast-enhanced CT imaging, contrast agent is administered into the joint space, enabling visualization of articular cartilage [26,27]. The previous method was designed to segment healthy and degenerated cartilage from contrast-enhanced CT images of knee joints. The method was based on cortical bone determination and subsequent automatic segmentation of articular cartilage. It was shown that semi-automatic segmentations of articular cartilage were reproducible and accurate as compared with manual segmentations. However, the applicability of this method to provide useful geometries for computational FE modeling of the knee joint has not yet been evaluated. This would be important as automated segmentation of joint tissues from clinical images and FE modeling of joint biomechanics would ideally enable clinical diagnostics of joint function to predict the onset and progression of OA.

This study aims to investigate the suitability of a recently developed semi-automatic segmentation method for contrast-enhanced CT images [25] in biomechanical FE modeling of knee joint function. Manually generated segmentations were used for comparison. FE model predictions were compared between the models generated by manual and semi-automatic segmentation methods during physiologically relevant loading. We hypothesize that the FE models generated based on the semi-automatic segmentation method would yield similar values for cartilage stresses and strains as the FE models generated based on manual segmentation.

In this study, both manual and semi-automatic segmentations of femoral and tibial cartilages were used. The segmentations were based on the contrast-enhanced CT data (105 mM, Hexabrix™ 320, Guerbet, Roissy, France; diluted to 0.9% saline for intra-articular injection) from seven patients (Age = 57.7±5.1 years, BMI = 27.9±4.9 kg/m²) with persistent knee pain and arthroscopically confirmed cartilage degeneration. The knees had International Cartilage Repair Society (ICRS) scored lesions (grades 1–3) evaluated by an experienced orthopedic surgeon. The subjects provided written informed consent, and the study protocol was approved by the Ethical Committee of the Northern Ostrobothnia Hospital District (Decision No. 33/2010). Moreover, the study adhered to the Declaration of Helsinki.

The proximal tibia and distal femur bones were segmented from the CT images by drawing axial contours sparsely using Stradwin (v5.2, Department of Engineering, University of Cambridge, UK). 3D surfaces were generated from these contours. Subsequently, cortical bone thickness and periosteal and endocortical surfaces, i.e., outer and inner cortical surfaces, were determined. The optimization method utilizes deconvolution of intensity profiles perpendicular to cortical surface. The algorithm captures both periosteal and endocortical surfaces simultaneously. This was done because optimization of cortical thickness is required when bone–cartilage interface is obtained accurately [28]. Regions where cartilage covers bone epiphyses were determined by registering bone surface templates, which included the bone-cartilage region information, on the generated periosteal surfaces of femur and tibia (Matlab, R2015a, MathWorks, Inc., Natick, MA). First, surfaces were registered rigidly, and then the templates were registered affinely to the periosteal bone surfaces. Subsequently, the articular cartilages were segmented. Intensity profiles were captured along surface normals for each vertex point at the bone–cartilage interface to determine the articular cartilage surface. Due to contrast enhancement, the joint space had high intensity whereas cartilage had low intensity; these local minima and maxima were used to determine the cartilage surface. A detailed explanation of the semi-automatic segmentation method was presented in a previous study [25].

Since the manual segmentation is widely regarded as the gold standard method [18,21,29], our models were compared against models generated by manual segmentation. Cartilages were segmented manually using Seg3D (v2.3, University of Utah, UT). Finally, all of the stereolithography (STL) surfaces were converted to a solid standard ACIS text format in matlab, and these geometries were meshed in Abaqus (v6.14, Dassault Systèmes, Providence, RI).

The surface mesh processing differed between semi-automatic and manual approaches and therefore the effect of the surface mesh processing was also examined. The surface meshes (STL) of semi-automatic segmentations were created directly from segmentations in matlab [25]. Additionally, the surfaces of semi-automatic segmentations were postprocessed similarly as those of the manual segmentations by using mimics with customized parameters (smoothness: 0.6, triangle reduction: 0.1). To analyze the effect of this specific postprocessing, the results from the FE models generated from these “modified semi-automatic” segmentations were compared with those generated from manual and semi-automatic segmentations.

Finite element meshes were created in abaqus (v6.14, Dassault Systèmes, Providence, RI) using first-order four-node porous continuum tetrahedral elements (type = C3D4P, one integration point in each element). The global distance between the nodes was set to 2 mm and 1 mm for the femoral and tibial cartilages, respectively, with a tolerance of 20% in size using both curvature control and minimum size control. Since the implemented loading caused high compressive strains, the stability of the meshes between the master (femoral cartilage) and slave (tibial cartilage) surfaces was ensured by doubling the density of the mesh of the slave surface. A sensitivity study was conducted using a mesh with six times the original element density (Fig. S1, which is available in the Supplemental Materials on the ASME Digital Collection).

To obtain sufficient contact interaction, the femoral cartilage was defined as a surface and the tibial cartilage was defined as a set of nodes. The discretization method between the master and slave surfaces was defined as a “surface-to-surface” contact with finite sliding; pressure-overclosure was evaluated as hard. Furthermore, tangential movements were defined as frictionless. Fluid flow was assumed to be negligible due to the application of high-rate dynamic loading [30] and hence free fluid flow was not allowed through the cartilage surfaces. The tibial bone–cartilage interface was fixed in all directions. The femoral bone–cartilage interface was coupled with a reference point, which was set in the middle-central between the medial and lateral epicondyles of the femur, similarly as in a previous study [7].

Instead of applying a simple axial loading, the simplified gait loading was applied in the models to enable a more extensive evaluation of the geometrical differences between the segmentation approaches. Similarly as in the earlier studies [7,31], joint motion and loading were run over 0.8 s, which covered one simplified stance phase. Movement determined in literature [32,33] was controlled using time-dependent boundary conditions. The force was scaled to match the weight of each subject by adjusting the maximum force to >2 times body weight. The flexion angle followed the simplified gait obtained from previous studies (Fig. 1(c)) [7,31]. Similarly as in those studies, a free varus-valgus alignment was allowed to maintain sufficient tibiofemoral contact, i.e., constant contact at the medial and lateral compartments, during the introduced loading. Further motions were fixed: anterior–posterior and medial–lateral translations due to their subject-specific variation, and internal–external rotations due to their nonsystematic behavior between subjects [34,35]. This is a reasonable assumption, since the aim of this study was to compare the models generated by two different segmentation techniques.

where I is the unit tensor. The chosen material parameters (Table 1) are based on an earlier experimental study [30]. The nonfibrillar matrix was modeled by the non-fibrillar matrix modulus (E_m) and Poisson's ratio (ν_m) using the neo-Hookean hyperelastic model, whereas the fibrillar network (17 fibrils in each integration point) was modeled by the viscoelastic primary and secondary fibrils with the initial fibril network modulus (E₀), strain-dependent fibril network modulus (E_ε) and viscoelastic damping coefficient (η). The effect of fluid support during dynamic loading was considered by Darcy's law using constant permeability (k) (see Supplemental Materials on the ASME Digital Collection for more details).

Based on previous studies [36–39], the primary collagen fibrils (four fibrils in each integration point) were oriented in specific split-line orientations, whereas the secondary fibrils (13 fibrils in each integration point) were randomly oriented. Fibril orientations for the secondary fibrils were coded directly in the user-defined material model (UMAT), model subroutine in abaqus, whereas the primary fibril orientations were first calculated in a global coordinate system separately for each integration point in each element using a custom matlab script and saved into a separate file. When solving models in abaqus, the primary fibril orientations were read by the UMAT subroutine. The definition for fibril implementation does not differ whether using hexahedral or tetrahedral meshes [15,40]. Since implementation of identical depth-dependent material properties of cartilage for each model would require identical tetrahedral element meshes, which was not the case between all compared models, depth dependency of the collagen fibril orientation was ignored. We wanted to ensure that differences between the models would be produced only by the geometry and not by the implementation of the fibril orientation. Thus, collagen fibrils were orientated parallel to the surface throughout the cartilage depth.

Model simulations were run implicitly with the abaqus/standard solver using consolidation analysis including three subsequent steps [31]. The first step included only an axial translation of femur to ensure an initial contact between the femoral and tibial cartilages (step duration = 0.1 s). The second step included an application of the initial force and flexion angle matching with the utilized data from the simplified gait loading (step duration = 0.1 s). The third step included a gait loading with time-dependent boundary conditions for the flexion angle and forces through the tibiofemoral joint (Fig. 1(c), step duration = 0.8 s).

Contact force and contact area of the tibiofemoral cartilage–cartilage contact were obtained from all the models. In addition, fibril strain, pore pressure, maximum principal strain (tensile strain), and maximum principal stress (tensile stress) distributions were averaged over the tibiofemoral cartilage–cartilage contact area (Fig. 1(c)). All these parameters were then analyzed as a function of the stance.

The three models generated by manual and semi-automatic segmentation methods were compared against each other (manual versus semi-automatic; manual versus modified semi-automatic; semi-automatic versus modified semi-automatic). For each pair of the models, statistical differences of each parameter were examined over the stance phase by using 1D statistical parametric mapping implemented in matlab [41]. This method was chosen since the traditional methods, such as the parametric t-test or nonparametric Wilcoxon signed rank test, do not account for multiple comparisons on smooth and random 1D trajectories. Due to the low number of subjects, we used a nonparametric statistical parametric mapping approach that corresponds to the two-tailed dependent-samples Wilcoxon signed rank test. For comparison, we also conducted traditional Wilcoxon signed-rank tests (IBM^® SPSS^® Statistics, v21, IBM Corp., Armonk, NY) (see Supplemental Materials on the ASME Digital Collection).

When the FE models based on manually and semi-automatically generated segmentations were compared at the first peak loading, i.e., at ∼25% of the stance phase, the distributions of maximum principal stresses were similar between the models for both the femoral and tibial cartilage surfaces (Fig. 2). Furthermore, neither of the models displayed systematically higher maximum stress values than the other; instead, the model with the highest stress values varied from subject to subject. (Fig. 2). Mean maximum principal strains on the contact area along the stance phase revealed only small and nonsystematic differences between the models (Fig. 3).

Differences in contact forces between the manual and semi-automatic FE models were minor, showing no systematic variation in the varus-valgus orientation (Fig. 4(a)). On the contrary, standard deviations (SDs) of contact area differences were relatively high. No statistically significant difference was found in the contact force or area between the models (p ≥0.05). However, the maximum principal stress differed statistically significantly both in the medial and lateral compartments (Fig. 4(a)).

When the FE models generated from manual and semi-automatic segmentations were compared, pore pressure was the parameter, which exhibited the greatest variation (SD ∼11%). Nevertheless, the differences pore pressure values were statistically insignificant in both models (p ≥0.05) (Fig. 5(a)). Differences were smallest in maximum principal strain and fibril strain. There were no significant differences in these parameters (p ≥0.05) (Fig. 5(a)).

There was a better agreement when the values of fibril strain and maximum principal stress were obtained from the models generated by the manual and semi-automatic segmentations (Figs. 4(a) and 5(a)) as compared to the comparison of the models generated by the manual and modified semi-automatic segmentations (Figs. 4(b) and 5(b)). Moreover, statistical differences were detected between the models generated from semi-automatic and modified semi-automatic segmentations (Fig. 5(c)).

No systematic differences were observed in the parameters between the manual and semi-automatic approaches in a subjectwise manner (Fig. 6). These comparisons revealed that differences in the parameter values between the models were slightly higher in the medial compartment than in the lateral compartment. When comparing the models generated by the manual and semi-automatic segmentations, the mean absolute difference in the contact area was <25 mm² and <16 N in the contact force. The mean absolute differences in the maximum principal stress were <1 MPa and in the pore pressure <1.5 MPa. For the values of maximum principal strain and fibril strain, the mean absolute differences were <0.72% and <0.40%, respectively.

The time required for manual segmentation of femoral cartilage was approximately 290 min and that of tibial cartilage was about 160 min. semi-automatic segmentation of femoral/tibial cartilage took approximately 20 min.

In this study, we compared FE models generated from manually and semi-automatically segmented joint geometries of seven osteoarthritic knees. The semi-automated segmentation method enabled more rapid FE model generation while still achieving similar results as those obtained with the more time-consuming manual method.

There were only minor differences in maximum principal strain and fibril strain between the models generated from semi-automatic and manual segmentations. This suggests that cartilage strains may represent the most useful parameters that can be analyzed from FE models generated from semi-automatic segmentation. This would be advantageous in clinical practice since excessive tissue strains may be risk factors for cell death and proteoglycan loss [9,42,43]. Moreover, collagen fibril strains may be used in the prediction of collagen failure and the progression of fibril damage [10,11,44].

Another parameter, maximum principal stress, is often analyzed from FE models of the knee joint to reflect the point of cartilage failure [6,45,46] and the initiation and progression of OA [7,8]. In this study, even though maximum principal stresses showed more variation between the models when compared to strains, the mean difference in the absolute values of this parameter was less than 1 MPa. This small difference introduces very little uncertainty in the estimation of cartilage failure since cartilage failure stress in tension can vary from approximately ∼5 MPa to ∼15 MPa and is dependent on several factors, such as age, location, strain-rate, and OA grade of cartilage [7,15,47–50].

There has been a rapid development of automatic segmentation methods in the past few years. Neural networks, statistical shape models, and atlas-based methods provide a good basis for determining the geometries for soft tissues [19–23]. This study utilized a previously presented segmentation method [25] which reveals articular cartilage geometries; the method was proven to be accurate in terms of dice similarity coefficient, specificity, and sensitivity. The utilized technique has also been shown to enable an accurate determination of cortical bone thickness from CT images [28], which is beneficial when evaluating alterations related to OA. Importantly, articular cartilage surfaces were segmented automatically with the custom matlab script, leveraging the fact that administrated contrast agent provides a high degree of image contrast. Even though the segmentation method used in this study was not fully automatic, one important aspect was that the surfaces were generated in 3D. Most of the previously introduced segmentation methods rely on slice-by-slice 2D segmentations prior to the 3D construction of surfaces. Furthermore, the surfaces required no postprocessing for the presented modeling purposes.

Few semi-automatical segmentation methods have been combined with computational modeling [24,51]. These previously applied segmentation methods require modifications of the surfaces before they can be implemented into computational models. Furthermore, in these prior studies, simulations were conducted with a simple axial loading instead of applying more physiologically relevant loading conditions. Here, the segmented surface mesh of cartilage was directly used for both tetrahedral meshing and then applied in the subsequent steps in the modeling procedure. The semi-automatic segmentation approach to generate FE models described herein shows promise for analysis of joint function. However, mesh generation, assignment of boundary and loading conditions, and model analysis can be time-consuming, and require substantial domain expertise (i.e., expertise in computational biomechanics). Thus, future research should focus on automation of all steps of the FE modeling pipeline to support clinical application. For instance, automated or template-based FE meshing could also be included in the procedure [52].

In this study, the segmentation method was utilized for contrast-enhanced CT images. Since the pathological changes in OA occur in both bone and cartilage [53,54], contrast-enhanced CT can provide a comprehensive and quantitative analysis of these tissues as well as reveal detailed information on the topology of the cartilage surface. In future studies, the method could be tested for magnetic resonance images, which have been acquired with ultrafast spin-echo sequences that capture signal also from cortical bone.

Slight differences in cartilage topography, such as some roughness on the surface close to the edges, may explain the differences in the contact area and maximum principal stress between the models. For instance, even though the distributions of maximum principal stress were found to be similar (Fig. 2), some subjects displayed relatively high differences in the peak values. These differences may be reduced by improving the surface smoothing algorithm in the semi-automated segmentation method. However, when using the Mimics surface mesh generation, i.e., similar postprocessing of surfaces to that used in manual segmentations, to produce smoother surfaces (modified semi-automatic versus semi-automatic segmentation), the results were not really improved (Figs. 2 and 4). Instead, there were actually greater differences between the simulation results of the models generated from manual and modified semi-automatic segmentations. This again suggests usefulness of our segmentation method for biomechanical modeling purposes.

One limitation of this study is that the number of subjects was rather low which does not allow an extensive comparison of the models. Furthermore, since subject-specific gait information was unavailable, a simplified gait loading was implemented based on the literature [32]. Instead of having hexahedral meshes, tetrahedral meshes were used since they were more suitable for assessing irregular segmented articular cartilage structures. The modeling results were also compared between the models with tetrahedral and hexahedral meshes and they were shown to be similar (Fig. S2 which is available in the Supplemental Materials on the ASME Digital Collection). Only articular cartilages were involved in the models since the present segmentation method is not suitable for handling other tissues, such as menisci. Nonetheless, the effect of menisci was evaluated, and it was observed that differences between the models remained the same whether or not menisci were included (Fig. S3 which is available in the Supplemental Materials on the ASME Digital Collection). Naturally, the lack of menisci influences the absolute values of the analyzed mechanical parameters, such as maximum principal stress. However, it was not necessary to include menisci herein since our primary objective was to compare cartilage mechanics predicted by the FE models generated using different segmentation techniques. Future research could focus on improving the semi-automatic segmentation algorithm to accommodate other soft tissues. Bones were considered as rigid in the models, which is a reasonable assumption, since bone is much stiffer than cartilage. Implementation of bones might slightly change cartilage responses [55]; however, it should not affect the conclusions of this study because the same rigid material assumption for bones was used in all of the models being compared against each other. The motion and forces through the tibiofemoral joint were generated using the known joint contact force and flexion angle of simplified gait [7,31], and, hence, it was not necessary to incorporate ligaments into these models [56].

In conclusion, a novel semi-automatic segmentation method was applied to generate geometries for FE modeling of cartilage biomechanics in the knee. The models generated from manual and semi-automatic segmentations produced similar results. However, semi-automatic segmentation was ten times faster. Thus, the semi-automatic segmentation method described and evaluated herein shows promise for future computational biomechanics studies of the knee and possible clinical application where the model geometry is generated from patient-specific CT arthrography scans.

The research leading to these results has received funding from the Academy of Finland (grant Nos. 286526, 305138, 269315, 324994, 328920, 324529, and 307932), the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 755037), Doctoral Program in Science, Technology, and Computing (SCITECO, University of Eastern Finland), The Research Committee of the Kuopio University Hospital Catchment Area for the State Research Funding (project 5041757) Kuopio, Finland, Sigrid Juselius Foundation, and Finnish Cultural Foundation. CSC—IT Center for Science, Finland, is acknowledged for providing computational resources and modelling software. Mimmi K. Liukkonen (Ph.D.) is acknowledged for the determination of meniscal supporting forces and Santtu Mikkonen (Ph.D.) is acknowledged for expertise in statistical analyses.

• Academy of Finland (Grant Nos. 286526, 305138, 269315, 324994, 328920, 324529, and 307932; Funder ID: 10.13039/501100002341).

• The European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant No. 755037; Funder ID: 10.13039/10001066).

• Doctoral Program in Science, Technology and Computing (SCITECO, University of Eastern Finland).

• The Research Committee of the Kuopio University Hospital Catchment Area for the State Research Funding (project 5041757) Kuopio, Finland (Funder ID: 10.13039/501100004092).

• Sigrid Juselius Foundation.

• Finnish Cultural Foundation.