Abstract
Knowledge about the knee cartilage deformation ratio as well as the knee cartilage stress distribution is of particular importance in clinical studies due to the fact that these represent some of the basic indicators of cartilage state and that they also provide information about joint cartilage wear so medical doctors can predict when it is necessary to perform surgery on a patient. In this research, we apply various kinds of sensors such as a system of infrared cameras and reflective markers, three-axis accelerometer, and force plate. The fluorescent marker and accelerometers are placed on the patient’s hip, knee, and ankle, respectively. During a normal walk we are recording the space position of markers, acceleration, and ground reaction force by force plate. Measured data are included in the biomechanical model of the knee joint. Geometry for this model is defined from CT images. This model includes the impact of ground reaction forces, contact force between femur and tibia, patient body weight, ligaments, and muscle forces. The boundary conditions are created for the finite element method in order to noninvasively determine the cartilage stress distribution.
1. Introduction
Sports activities and daily routines such as standing, walking, running, jumping, and other recreational activities impose relatively large loads and movements on the human knee joint. These tasks could cause injuries and degenerations in the joint ligaments, menisci, cartilage, and bones. Thus, knowledge of in vivo joint motion and loading during functional activities is needed to improve our understanding of possible knee joint degeneration and restoration. Internal loadings of knee anatomical structures significantly depend on a lot of factors such as external loads, body weight, ligaments, strengths, and muscles forces.
In the paper [1] the knee implants were used to directly measure loads of participants during daily activities. In vivo knee measurements are very invasive and practically impossible for the case described in our study. There are many techniques for measuring an external variable which can be further used to create biomechanical and mathematical models.
Other studies [2] were based on the registration of fluoroscopic images and computer model of the knee. The most recent method utilized force plate data, CT or MRI skeletal structure data, and motion capture obtained from the infrared position sensor [3, 4]. In another research paper [5] the accelerometer was used to estimate the angle of lower extremities. This is a pretty cheap technique but requires additional processing of collected data and the error of estimated angle is up to six percent. The measure performance can be improved using a combination of accelerometer and gyroscope sensor such as goniometer [6–8].
In the study [9] stress on the knee cartilage during kneeling and standing using finite element models is compared. They used magnetic resonance (MR) images of the flexed knee to build a geometrical model. As a computational tool they used commercial software MIMICS. The results of those studies showed some differences in high-stress regions between kneeling and standing. The conclusion was that the peak von Mises stress and contact pressure on the cartilage were higher in kneeling. The study [10] used computed tomography (CT) images of knee structures during static loading to determine cartilage strains and meniscal movement in a human knee at different time periods of standing and to compare them with the subject-specific 3D finite element (FE) model. The results of these experiments showed that 80% of the maximum strain in cartilage developed immediately and after that cartilage continued to deform slowly. In the study [11] magnetic resonance (MR) images of the right knee of a 27-year-old male subject were used to determine the subsequent alteration in the fluid pressurization in the human knee using a three-dimensional computer model. The results of these studies indicated a redistribution of stresses within the tissue and a relocation of the loading between the tissue matrix and fluid pressure.
The purpose of this study was to estimate stress distribution in the knee cartilage. For that purpose the appropriate system of cameras and force plate platform were used. The deformation of cartilage was measured using marker position data and the 3D model of the lower leg segment. The models were established from computed tomography (CT) images. Simultaneously the deformation is assessed only matching single infrared camera images and CT image using software for image registration technique. In the experimental part the accelerometer sensor was used which potentially can provide more information about the gait. In the computer simulation the FEM analysis was applied with an adaptive change of mechanical parameters of tissue variation in order to match the measured force and deformation.
2. Materials and Methods
2.1. Mechanical Model of the Knee Joint
Knee joint motion represents a complex combination of rotations and translations. The major parts involved in the knee kinematical behavior include femur, tibia and patella. Forces that act at a knee joint are given in Figure 1.
The dominant forces that act at a knee joint are body weight and ground reaction force which are opposite to each other. Muscle forces as well as contact forces between femur and tibia are also included in the model. Equilibrium equations of the knee joint are given below [5]:where is the ground reaction force, is the body weight, and are the two contact point forces, and are the two contact points normal, are the ligament and capsule forces, is the knee joint driver moment along local axis and and are the position vectors where the corresponding forces are being applied [13]. We assumed that femur and tibia could be represented as rigid bodies and that their deformations were very small in contrast to relatively large deformations of cartilage and ligaments. Note that synovial fluid significantly reduces the friction between cartilage surfaces and menisci.
A simplified spring-damper-mass model which was used in this study is shown in Figure 2. It consists of four masses. The upper body was modeled using two masses, one representing its rigid mass, , and the other representing its wobbling masses, . The thigh, leg, and foot of the supporting leg were modeled using two masses, one representing its rigid mass, , and the other representing its wobbling masses, .
The total body mass was obtained from the participant. In the following system of equations (Equation (2)) a dynamics system is described [12] and was later used to calculate the resultant force and moment of the knee cartilage during the stance phase of the gait cycle:In (2), was the lower body rigid mass and was the wobbling mass, was the upper body rigid mass and was the wobbling mass, was the compressive spring and was the damper that connected the upper and lower rigid bodies, was the spring and was the spring-damper unit which connected the lower wobbling mass to the upper and lower rigid bodies, and was the spring and was the spring damper unit which connected the upper wobbling mass to the upper rigid mass [12]. was the vertical contact force which was defined as is contact area and , and are the parameters of the ground reaction model which defined the deformation of the shoes during standing. Parameters for soft and hard shoes are shown in Table 1.
Cartilage was considered as a porous deformable body filled with fluid occupying the whole pore volume. The physical quantities for this analysis were the displacement of solid , relative fluid velocity with respect to the solid (Darcy’s velocity) , fluid pressure , and electrical potential . The governing equations for the coupled problem are described as follows. First, we considered the solid equilibrium equation:where was the stress in the solid phase, was porosity, was the permeability matrix, was the density of the solid, was body force per unit mass, was relative velocity of the fluid, and was acceleration of the solid material. The operator wasThe equilibrium equation of the fluid phase (no electrokinetic coupling) waswhere was pore fluid pressure, was fluid density, and was fluid velocity. This equation is also known as the generalized Darcy’s law. Both equilibrium equations were written per unit volume of the mixture. Combining (3) and (5) we obtainwhere was the total stress which can be expressed in terms of and asand was the mixture density.
Here was a constant vector defined as to indicate that the pressure contributes to the normal stresses only. We also had to take into account the fact that the pressure has a positive sign in compression. Tensional stresses and strains were considered positive as well. In the following analysis we employed the effective stress, , defined aswhich was relevant for the constitutive relations of the solid. Using the definition of relative velocity as the volume of the fluid passing in a unit time through a unit area of the mixture (Darcy’s velocity), we obtainedand transformed (6) intoThe final continuity equation using the elastic constitutive law and fluid incompressibility was given in the form The resulting FE system of equations was solved incrementally [14] with a time step . We imposed the condition that the balance equations are satisfied at the end of each time step . Hence, we derived the following system of equations:where , and were forces in the balance equations for displacement, pressure, fluid velocity, and electrical potential, respectively, and and were mass terms in mass matrix [14].
2.2. Experimental Parts
In this study we used a commercial motion capture system OptiTrack. This system consists of six infrared cameras and four retroreflective markers. The markers are 1.5 cm in diameter and are attached at the precise anatomical locations of the participant’s leg for unilateral gain analysis. These locations were great trochanter region, femoral lateral epicondyle, tuberosity of the tibia, and the center of the anterior region of ankle joint (see Figure 3).
The computerized camera system with accompanying software captures the exact motion of retroreflective markers and thus records their trajectory while the volunteer performs walking over the force plate. The cameras were connected to a computer that collects gait kinematical data. The result of motion tracking is a series of 3D coordinates for each numbered marker. To better understand kinematics and kinetics of gait we used a three-axial accelerometer.
We used Sun SPOT accelerometers (Sun Small Programmable Object Technology) to wirelessly detect the middle of the gait stance phase, that is, the moment when the ground reaction force reaches its maximum value. The Sun SPOT has a sensor board which consists of 2G/6G 3-axis accelerometer, temperature sensor, and light sensor. In the experimental part the Sun SPOT LIS3L02AQ accelerometer is used to measure the orientation or motion in three dimensions, , , and with the sample rate of 100 Hz. Each of these components represents the sum of static accelerations defined by angle of inclination to the corresponding axis and dynamics component stemming from the movement during the walk. The ground reaction force was sampled from the multiaxis AMTI force plate at the rate of 100 Hz. Before the experiment, calibration was conducted to work out the space coordinate system for the camera system field of view. The calibration was performed fully in accordance with the proposed manufacturer’s procedure.
The volunteer performs walking along 2.5 meter distance path away with his own ordinary velocity and attached infrared marker and accelerometers sensor on the left leg. Results for marker coordinates and corresponding accelerations, measured by the three-axial accelerometer, are presented in Figure 4.
According to [15–17] measured values for marker position are influenced by noise due to the wobbling of the participant’s skin. The value of this uncertainty is in range ±2 mm.
The force plate is positioned in the first half of the walking path. During the experiment, the participants were asked to walk along so the force plate records the value of the ground reaction force (Figure 5).
The force value is zero in the beginning of the walk and when the participant stands on the force plate, starting with the heel, the force gradually increases and reaches the maximum and then drops to zero again when the foot is detached by the force plate.
3. Results
One of the main tasks for preparing data for FEM simulation is matching the infrared reflective marker position obtained by infrared camera and landmark position on the CT images when cartilage is unreformed. This procedure is known as image registration and ANTs (Free software for image registration) is used [18].
Main goal of the registration process is obtained from the transformation parameter that maps the marker position in deformed and undeformed knee cartilage. In general this process is finding the optimal transformation that maps every pixel of image with the presence of deformation and another static image when cartilage is undeformed. As similarity measure of image the mutual information is used. The standard algorithm in registration process is elastic deformation procedure. The idea is to model the contours on one of the matching images as elastic object which deformed under the influence of some external forces [19]. In each step of the deformation the images are compared and value of external forces is now proportional to the difference between these cases on the basis of mutual information value. The process is repeated until the difference between the images is greater than some error of convergence. The result of registration between CT and infrared images is shown in Figure 6. The measured difference of the marker positions between deformed and undeformed cartilage is mm. Images showing markers matching these cases are presented in Figures 6(a) and 6(c), respectively. The measured uncertainty of 0.06 mm corresponds to the dimension of a single pixel. The error in estimation of the deformation is significantly greater due to the influence of the distortion of the camera lens, skin movements, and insufficiently accurate registration result.
(a)
(b)
(c)
Simultaneously, we create a 3D model using CT slices. The CT slices are segmented using a threshold value and compared with gray value of the image pixel. The segmented images are submitted by the edge detection operator for boundary extraction. The merging of adjacent boundary is done in each slice and final 3D model is created (Figure 2). This model includes femur, tibia, cartilage, surrounding tissue, and skin.
The anatomical point (tuberosity of tibia and femoral lateral epicondyle) can be easily detected in the model and initial position of the marker can be obtained. Using the measured position of infrared marker (Figure 4) and anatomical position of the marked point on the created 3D model (Figure 7(c)) we can obtain a vertical deformation of cartilage as difference of these values. The observed measured value coincides with the point when the ground reaction force is in the maximum. According to Figure 5 this moment is 1.49 seconds after starting the recording of gait and the deformation is mm. The measured uncertainty of 0.01 mm emerges as a consequence of the limited resolution of the motion capture camera system and resolution of the CT scanner. This methodology for obtained deformation is more precise than registration method but it requires more processor time and memory.
Using the same procedure for image segmentation, a full model of knee joint is created. The model consists of femur, tibia, and cartilage.
We used measured values for the displacement and ground reaction force in order to calculate corresponding matrix elements in the relation (13). Upon model creation we applied boundary conditions: (a) we clamped the distal end of tibia and (b) axially loaded the femur with the maximally measured ground reaction force N (Figure 8(a)).
(a)
(b)
(c)
For modeling the cartilage and meniscus we implemented a finite element formulation where the nodal variables are displacements of solid, ; fluid pressure, ; Darcy’s velocity, ; and electrical potential, with estimated matrix elements. A standard procedure of integration over the element volume was performed and the Gauss’s theorem was employed. An implicit time integration scheme was implemented.
The tetrahedral mesh model had 20537 elements and 4693 nodes (Figure 8(b)). We used PAK solver [20] for the FEM analysis. Total execution time of the analysis was around 30 minutes on the core I7 processor with 12 GB of RAM memory.
The initial mechanical characteristic, Young’s modulus, and Poisson’s ratio, for the femur and tibia, were amounted to GPa and , for the isotropic cartilage MPa and , and for the transversely isotropic menisci MPa and . All these values were taken from the literature [21].
These values were adaptively changed for the purpose of correspondence between the measured deformation and ground reaction force.
The final value of the Young’s modulus of the cartilage is 5.62 MPa with error of estimation of 0.01 MPa. The Young’s modulus for the femur and tibia is adapted to the 18 GPa while the Poison’s ratio is 0.45.
Resultant stress distribution of the FEM analysis for the elements of the knee joint is given in Figure 8(c).
The von Mises stress on the cartilage is presented in Figure 9. As it can be seen the maximal value of stress has MPa magnitude of order and is located on the boundary of the cartilage. This is in compliance with the fact that is cartilage is the weakest part of the knee joint with a tendency to injury and fraying.
The procedure described in this study can offer very useful information for physicians in order to better understand injury formation and improve the process of rehabilitation and, in some perspective, as a support in clinical decision making.
4. Conclusion
The main goal of this study was to introduce a new approach towards a noninvasive effective stress calculation for a specific participant. Input data were provided from the experimental measurements during the participant’s walking test whereupon finite element analysis was performed giving the distribution of the effective stress in the major anatomical elements of the knee joint: femur, tibia, and cartilage. This approach demonstrates a great possibility for preoperative and postoperative surgical planning and treatment of the knee injuries for specific patients.
This study contains some limitations. We neglected the skin movement artifact during the experiment and the influence of ligament presence. However, the current model allows us to investigate the effect of different biomechanical factors and loads on the stress distribution at the knee joint. In this study we used data obtained from infrared cameras and force plate sensors. Besides, in our future work we will try to replace a relatively expensive system of infrared cameras with a much cost effective system of accelerometers so that we will be able to calculate positional data of anatomical points of the lower extremities solely by using their accelerations. The very promising are techniques of the image registration that can be used for assessment of gait parameter using a cheaper mobile cell camera. This will be a consideration of the future work.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This paper is supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia (Projects numbers III41007 and ON 174028).