Finite Element Modeling of Articular Cartilage to Characterize Biomechanical Properties: The Effect of Cartilage Surface Curvature

Degeneration and loss of articular cartilage in the synovial joint have been recognized as the main source of osteoarthritis which leads to pain, swelling and limit the joint mobility. Extensive experimental and computational studies have been performed to study the mechanical behavior and characterize the biomechanical properties of articular cartilage. However, a lack of attention was made on the curvature of the cartilage surface by assuming it was a flat surface. This assumption was inappropriate since the synovial joints possessed curved geometrical shape and may contribute to inaccuracies in characterizing the articular cartilage biomechanical properties. This study aims to examine the effects of the curvature of the cartilage surface in finite element modeling which incorporated with the experiment method to characterize biomechanical properties of articular cartilage. In this study, the biomechanical behavior of contact pressure and pore pressure were investigated at different radius of cartilage surface using the finite element method. The cartilage biomechanical properties of elastic modulus and permeability of the bovine humeral head were then characterized using a combination of indentation test and finite element method. It was found that the cartilage curvature produced a 6% difference in contact pressure and a 39% difference in pore pressure distribution compared to the flat surface cartilage in finite element analysis. Furthermore, significant observation in the characterized biomechanical properties was obtained where the differences of the cartilage curvature reached 33% for elastic modulus and 56% for permeability. Based on the results, the surface curvature of articular cartilage could play an important role in the computational modeling and characterization of its biomechanical properties.


Introduction
Osteoarthritis (OA) is the most common joint disease which normally occurs at the knee, hip, hand and spine. In the early stage of OA, the patient will have difficulties in their movement in daily activities which could lead to chronic joint pain [1,2]. It was reported that nearly 30% of OA moderate patients were led to severe disability [3]. Although OA patients were associated with over 55 years old, an increasing number of patients starts at the age of 30 years old were observed in the recent years [4]. It is well recognized that the common cause of OA is the degeneration of the articular cartilage [4,5].
Articular cartilage is a thin layer tissue that covers the end surface of the articular bone in a synovial joint. It has no blood supply or nerve system and plays a significant role to absorb the load in synovial joints of the human body [6]. The cartilage tissue consists of two distinct phases which are the fluid and solid phases. The fluid phase consists of 60-85% water while the solid phase composed of 15-22% collagen and 4-7% aggrecan by wet weight [7,8].
There are four different layers through the thickness of cartilage tissue where each layer has its own composition and characteristics. The superficial zone, which is the outer 202 Finite Element Modeling of Articular Cartilage to Characterize Biomechanical Properties: The Effect of Cartilage Surface Curvature layer which is 10-20% of the cartilage thickness, has the highest water content [7,8]. The middle zone represents 40-60% of the cartilage volume and contains more proteoglycans to provide swelling pressure and stable hydrated structure [9,10]. The deep zone is the borderline between the cartilage and the subchondral bone represents 30% of the cartilage volume which has a low cell density and provides the greatest resistance to compressive force [11]. This structure and composition contribute to inhomogeneous and uniqueness of the cartilage behavior. There were various constitutive material models have been used to represent cartilage from a single phase to multiphase models. However, the biphasic isotropic model was mostly used to represent the solid and fluid phases [12,13]. The most important biomechanical properties of the cartilage in the biphasic model are the elastic modulus E, and permeability κ [14][15][16]. The permeability indicates the rate of fluid flow in the tissue while the elastic modulus measures the stiffness and reflects the integrity of the type II collagen matrix and the proteoglycan content of cartilage tissue [8]. These biomechanical properties provide reference data for tissue development and regeneration including biomaterial development, functional imaging, functional assessment of engineered and repaired tissues [17,18].
The biphasic biomechanical properties of the cartilage are commonly characterized using a combination of experimental and numerical methods. The most common experimental methods used to characterize the biphasic properties are unconfined compression, confined compression and indentation test. However, the indentation test is often been utilized because of the ease of sample preparation and the test set-up allows the cartilage to be submerged in the fluid during the test. Numerical methods are applied to integrate the experiment data and the biphasic theory to characterize the properties [19,20].
Recent computational advancement has simulated the cartilage using the biphasic constitutive material model in the finite element (FE) method [21,22]. However, the articular cartilage was modeled as an axisymmetric model where the cartilage surface was assumed to be flat in to investigate the biomechanical behavior of the cartilage [19][20][21]. This assumption was inappropriate since the synovial joint possessed curved geometrical shape and this may lead to inaccuracy in characterizing the cartilage biomechanical properties. The curvature of the cartilage surface has been reported to be an important aspect of the computational method to analyze contact stresses in the synovial joint [19,23]. Therefore, this study aims to examine the effects of the curvature of the cartilage surface in finite element modeling which incorporated with the experiment method to characterize biomechanical properties of articular cartilage.

Specimen Preparation
The cartilage was harvested from the humeral head of healthy cows aged 2-4 years old. The cows were local breed known as Kedah-Kelantan (KK) with an average weight between 200-250 kg. The humeral heads were obtained within 24 hours after slaughtered.
Excessive tissues and bones of the humeral head joint were removed using an electric saw and a scalpel blade. Throughout the cutting process, the cartilage specimen was constantly moistened with phosphate buffered saline (PBS) solution to avoid dehydration of cartilage tissue. The humeral heads were cut into four sections as shown in Figure 1. The specimens (n =24) were then wrapped with PBS-soaked cotton wool prior to the creep indentation test.   Figure 2. Throughout the test, the cartilage specimen was submerged in the PBS solution to prevent the cartilage from dehydrated. The test was subjected to 0.38 N compression load using a 2 mm radius impermeable spherical indenter which resulted between 10-20% deformation of the cartilage thickness strains [10,14]. The indentation was carried out until the displacement had reached equilibrium. All the data were recorded at every 0.01 seconds using the data acquisition software LabVIEW 8.5.1 (National Instrument Corporation, Austin, TX, USA).

Cartilage Thickness Measurement
The cartilage thickness was measured using a sharp needle indenter with 3.16 N compression load where the load could penetrate the cartilage until it reached the underlying bone [14,24]. The displacement of the needle and the load response was recorded at every 0.001 seconds to enable an accurate measurement of the cartilage thickness. When the needle indenter contacted the surface of cartilage, the load signal was increased until it reached the bone of cartilage where the force and displacement were stabilized. The thickness of the cartilage specimen was obtained by determining the difference between the positions of the needle when it contacted the cartilage surface and the underlying bone as shown in Figure 3.

Cartilage Curve Measurement
The curvature of the cartilage specimen was measured using the SeDoz-G profile projector. The specimen was fixed on the platform to avoid any movement during the measurement. Ten points were marked on the cartilage surface to calculate the radius of the curvature as shown in Figure 4. The surface radius was determined based on the surface points using the software provided by the projector.

Development of Finite Element Model
FE model was developed to simulate the experimental cartilage creep indentation test. An idealized axisymmetric biphasic poroelastic FE model of the cartilage specimen was developed using Abaqus 6.14 (DS Simulia Corp., Providence, RI, USA) software. The thickness of the cartilage was set at 1.20 mm based on the average measured thickness with a 5 mm width as shown in Figure  5. The FE models were then constructed in concave and convex curvatures with the cartilage surface radius ranging from 10 mm and 50 mm to examine the effect of curvature on the contact pressure and pore pressure of the cartilage.   [14,26]. The bottom nodes of the bone were constrained in both horizontal and vertical directions, whilst the nodes on the axis were constrained in the horizontal direction. The spherical indenter was only allowed to move in the vertical direction, while the horizontal and rotational movements were constrained. As can be seen in Figure 5, the fluid flow was prevented at the bottom and the vertical symmetry axis of the cartilage surface whilst the outer edge nodes of the cartilage were maintained at zero pore pressure to allow unrestricted fluid flow. The contact-dependent flow algorithm was applied at the top surface of the cartilage where the flow conditions at the nodes were changed depending on the contact stress between the indenter and cartilage [26].
In this study, creep-deformation was simulated by applying a ramp load of 0.38 N on the spherical indenter for 2 seconds and maintained the load for further 1000 seconds [26]. The 2 seconds ramp period is an instantaneous effect based on an experimental study where the minimum time at which the creep indentation test of the cartilage could be compared reliably was 2 seconds after the load was released [27]. The 2 seconds ramp load was also used in previous computational studies [14,21,26]. In addition, geometric nonlinearities (NLGEOM) parameter was set to allow the geometric nonlinearity deformation occurred within the model. Although the automatic time increments (UTOL) parameter were applied in the model, it was controlled at a maximum pore pressure of 600 kPa in one increment [14,22,26].

Characterization of Biomechanical Properties
The biphasic elastic modulus and permeability of the cartilage were characterized using a combination of indentation test and FE method. Axisymmetric biphasic poroelastic FE models were developed for each of the cartilage specimen based on the measured radius and thickness with similar loading and boundary conditions. The elastic modulus and permeability were characterized by fitting the cartilage displacement curve generated from the FE analysis to the indentation test displacement curve. Least-square method was applied to fit the displacement curve using MATLAB R2015b software (MathWorks Inc, MA, USA) as shown in Figure 6. The initial values of equilibrium elastic modulus and permeability of the cartilage were used to initiate the iteration and the final optimized values were obtained when the function achieved the minimum specified convergence criteria or minimal square error occurring between the curves.

Results and Discussion
In the present study, the axisymmetric biphasic FE models were used to observe the effects of surface curvature on the biomechanical behavior of cartilage tissue. Figure 7 and Figure 8 show the contact pressure and pore pressure of the cartilage where the negative and positive values of the radius represent the concave and convex curves respectively.  Based on the results, the cartilage curvature played an important role to examine the contact pressure and pore pressure of the cartilage. The convex curve with 10 mm radius generated the highest contact pressure while the concave curve of the a 10 mm radius produced the lowest contact pressure at 2 seconds and 1000 seconds. The difference between both curves was 7% compared to the flat surface. However, the 10 mm radius of the convex curve produced 39% lowered in pore pressure as compared to the flat surface. This is in agreement with a computational study carried out by Holzapfel and Stadler in 2006 which they reported that the cartilage curvature was crucial and played an important role in the load-bearing characteristics study [30].
The FE analysis clearly indicated that the pore pressure was most affected by the cartilage curvature during the indentation test. This is mainly due to the frictional drag force of interstitial fluid flow being the dominant factor controlling compressive creep behavior. The biphasic model describes the nature of the cartilage in solid and fluid phases where the water content could reach up to 85% [8,9]. Due to the very low permeability of the cartilage, large drag forces were generated from the fluid flow which maintained high fluid pressure over a long period of time [11]. The fluid exudation from the cartilage tissue can be observed with the direction of fluid velocity as shown in Figure 9. At 2 seconds, the direction fluid flow exudation was concentrated at the center of the cartilage and the pore pressure dropped drastically at 1000 seconds as the fluid flow exudation was concentrated towards the edge of the cartilage. Based on the previous study, the cartilage becomes softer and more permeable when the water content is increased [15].
Further assessment of articular cartilage curvature was then carried out to examine the effects on the characterized biomechanical properties of elastic modulus and permeability of the cartilage tissue. The axisymmetric FE models of the cartilage were developed based on the measured thickness and curvature where the radius of the cartilage specimens was 26 mm, 27 mm, 28 mm, 29 mm and 30 mm. It was found that the averaged elastic modulus and permeability of the cartilage were 1.80±1.51 MPa and 0.83±0.50×10 -15 m 4 /Ns which were in a similar range with the previous studies [28,29]. Figure 10 shows the characterized properties of different radius curvature where the percentage difference between the minimum and maximum radius curvature for elastic modulus was 33%. However, a significant difference of 56% was observed on permeability. Similar trends were also found in previous studies where the cartilage permeability was more sensitive compared to other parameters [1,10]. The sensitivity of fluid flow properties of articular cartilage was also observed to be significant in the early phase of OA disease [9,12,13]. Therefore, the cartilage surface curvature could be an important factor to characterize the biomechanical properties of cartilage tissue. Moreover, the biomechanical properties of the cartilage have been reported to be an important aspect in the computational method particularly to examine the contact stress distribution and load-bearing characteristics [19].

Conclusions
The present study examined the effect of cartilage surface curvature on the biomechanical behavior of the cartilage. In FE analysis, a significant effect on pore pressure of the cartilage was observed for smaller curvature as compared to contact pressure. Furthermore, the curvature of the cartilage was observed to affect significantly on the characterized biomechanical properties of elastic modulus and permeability of the articular cartilage. The present results have indicated the importance to consider the curvature of the cartilage surface to study the biomechanical behavior of the cartilage.