2 Materials and methods

2.1 Bone-joint mechanics (humerus project)

▼ 47 (fortgesetzt)

To understand the mechanobiological behavior of a synovial joint the straining of the intact and fractured proximal humerus under physiological-like loading was determined. The model employed in this study was validated by comparison of the bone stiffness in the axial direction and in torsion calculated after numerical analysis to the results obtained from the in vitro analysis realized by Lill and co-workers.

For the analysis of the bone-joint mechanics (large model), the humerus bone was chosen for the following reasons: first, approximately six percent of all fractures occur in the proximal humerus, second: these fractures are critical in osteoporotic patients and remain a clinical challenge, third: new devices have been developed to achieve optimal healing, especially in weak bones, and need to be mechanically analyzed. Lill and co-workers conducted a research project involving all theses premises (Hepp, et al., 2003;Lill, et al., 2002;Lill, et al., 2001;Lill, et al., 2001;Lill, et al., 2003;Lill, et al., 2004). However, mechanical conditions affecting bone, measured in the straining pattern generated after implantation, are still unknown, and were studied in the first part of the present project.

2.1.1 Mechanical in vitro testing of proximal humerus defects

▼ 48 

The procedure described by Lill et al. shall briefly be summarized. 24-paired fresh human humeri were obtained for histomorphometric analysis of the trabecular bone structure (Lill, et al., 2001). The selected specimens were comparable both in age and bone mineral density distribution (BMD) (p>0.05). Subsequently an in vitro test of intact and fractured human proximal humeri was performed. The fractured bones were stabilized with different osteosynthese types and tested under compressive and torsional loads (Lill, et al., 2003). Out of this group seven unpaired humeri were selected for biomechanical testing. The defect was stabilized by a new, internal fixator with locking screws (angular stability) using the standard instrumentation supplied by the producer (Locking Compression Plate for the proximal humerus LCP-PH, Mathys, Bettlach, Switzerland). In contrast to conventional plates the head of these screws is aligned to the body of the implant and provides the screw with angular stability with respect to the implant (Fig. 2.1, center).

Screw length and placement were chosen according to manufacturer’s recommendations. Mechanical testing was performed using an electro-mechanical material testing machine (Type 1455, Zwick, Germany) in displacement control (Lill, et al., 2001). Relative movements of the proximal and distal fracture fragments were determined using an optical measurement system (PC Reflex, Qualysis, Sweden, accuracy < 0.1mm and < 0.1°) with reflective markers (Fig. 2.1, right). The technique and measurement accuracy of this procedure has previously been reported (Chu, et al., 2000). Implant stiffness was calculated from a 0.5 mm axial compressive displacement (5 mm/min) and 4° torsional rotation (50°/min) with an axial preload of 25 N. The specimens were kept moist during testing.

Fig. 2.1: Proximal humerus with defect localization and a LCP-PH implant under compressive loads. Left: Finite element model of the proximal humerus with implant, analyzed in this project. Middle: X-ray of a specimen with implant but prior to the osteotomy. Right: Specimen in material testing machine under pure compression. Reflective markers (distance between marker trees = 60 mm) are attached to each segment to allow determination of the interfragmentary movements. From these movements and loads, the stiffness of the bone implant construct was computed (Maldonado, et al., 2003).

2.1.2 Straining of intact and fractured proximal humerus under physiological loads

2.1.2.1 Simulation of the used human specimens and FEM validation

▼ 49 

Out of the tested group, two representative humeri were selected for this study: one osteoporotic and one normal with reference quality, and scanned using both QCT and DEXA (average DEXA value = 0.26g/cm2 and 0.49g/cm2 respectively). Geometry and density distribution from QCT-scans for these two bones were determined and modeled.

The finite element models were created from the QCT data sets of these two specimens using an automated mesh-generation algorithm to define inhomogeneous material properties based on the CT density information (Seebeck, et al., 2001). A standard phantom (a reference calibrated pattern to estimate bone density in dependence on a gray scale from MRI images) provided a reference to determine apparent wet density. To assign the Young’s modulus as a function of the apparent density, the relationship of Carter and Hayes was used (Carter and Hayes, 1976). According to a method described by Zannoni, a preliminary analysis was performed to define the minimum number of independent elastic modulus of Young’s required to achieve a convergence of the analysis (Zannoni, et al., 1998). Based on this, the material properties of the matching voxel-based QCT data points were averaged for each element, leading to the definition of up to 127 material properties. In each finite element model the Young’s modulus ranged between 22.3 MPa < E < 7837.9 MPa for the reference bone quality (RQ) and 22.3 MPa < E < 7422.5 MPa for the poor bone quality (PQ) specimens. The bones were modeled as isotropic, linear elastic and inhomogeneous material with a Poisson’s ratio of ν = 0.3. The FEM meshes of the RQ and PQ specimens contained 58,048 and 56,682 eight-noded brick elements (52,054 and 50,822 nodes) respectively.

Simulation of a 5mm osteotomy (shaft unimpacted in the Neer classification II (Neer, 1970) or proximal humeri A3 in the AO classification, (Müller, et al., 1989), at the level of the surgical neck was achieved by using low stiffness elements (E = 5 MPa, ν = 0.45) in both finite element models. The defect stabilization was simulated by modeling an angle-stable plate (LCP-PH, Mathys, Bettlach, Switzerland), consisting of 142 eight-noded brick elements and 108 thin 2-noded beam elements representing the screws. The screw length and placement was chosen according to the experimental model described before (Fig. 2.1). The implant materials, screws and the angle-stable plate, were modeled as isotropic, linear elastic and homogeneous with an elastic modulus of 110 GPa and a Poisson’s ratio of ν = 0.3 (commercially used titanium, Mathys, Bettlach, Switzerland).

▼ 50 

To verify the finite element models, the in vitro tests were simulated by applying nodal displacements to the most proximal (anatomical plane) node layer (0.5 mm axial compression and 4° angular rotation).

2.1.2.2 Simulation of physiological loads

To include the implant and muscle attachments the FE models were scaled to a total length of 345 mm by extruding the distal cross-sectional geometry and material properties. The attachments and orientations of 15 muscles that load the shoulder girdle for three different arm positions were taken from the literature and scaled to both humeral models (van der Helm, 1991;van der Helm and Veenbaas, 1991). A single vector represented each muscle. The loading scenarios included the neutral arm position (0°), 90° abduction and 90° forward flexion (Table 2.1). The contact force at the elbow joint was calculated from force vectors of those muscles crossing this joint and the corresponding resultant force vectors. Equilibrium was established by positioning the elbow contact force such that force and moment equilibrium for the whole humerus under physiological-like loading was achieved. The von Misses stress calculation is a common representation of results in continuum mechanics and applies to ductile materials (e.g. metals), which fail under shear loading. In contrast to metallic implants and since in vivo measurements are limited to strains principal strains were selected to represent the load state in the humerus sufficiently. All calculations were performed using the Marc/mentat package (Marc K72/Mentat 3.2; Marc Analysis Research Corp., Palo Alto, CA, USA) on a Unix workstation using linear elastic analyses (MIPS R 10000; Silicon Graphics Inc., Mountain View, CA, USA). The obtained results and clinical implications were reported by (Maldonado, et al., 2003).

2.2 Osteochondral healing (Galileo project)

Given the influence of mechanical conditions on bone healing, and the fact that the finite element method is able to simulate and analyze specific conditions (bone quality and physiological loads) during fracture stabilization, a similar approach should be applied for the study of osteochondral healing: A local model of the joint region necessarily needs to 1. Determine which mechanical conditions have to be selected for analysis, 2. Determine the effect of local conditions such as contact with another joint surface and joint curvature, 3. Consider modeling of specific features such as joint structure (e.g. consideration of the subchondral bone plate) and fluid conditions (interchange of fluids between articular surfaces), 4. Develop a versatile tool that allows combining different mechanical conditions to analyze healing using the proved finite element method.

▼ 51 

Table 2.1: Weight, muscles and joint contact forces in Newton in 90° abduction

Force Vector

Fx
(N)

Fy
(N)

Fz
(N)

Fr
(N)

pectoralis major p. thorak

-29.5

-3.4

-12.4

32.2

pectoralis major p. clav.

-16.9

-0.2

-4.5

17.5

latissimus dorsi

0.0

0.0

0.0

0.0

deltoideus scapularis

-367.7

51.0

53.3

375.0

deltoideus clavicularis

-85.8

2.5

-9.6

86.4

supraspinatus

-25.4

-13.0

9.1

30.0

infraspinatus

-4.4

-5.0

2.0

6.9

subscapularis

-328.7

-140.9

40.3

360.0

teres major

0.0

0.0

0.0

0.0

teres minor

0.0

0.0

0.0

0.0

coracobraquialis

-8.5

1.1

-1.4

8.6

biceps breve

0.0

0.0

0.0

0.0

bicep caput longum p. cran

0.0

0.0

0.0

0.0

bicep caput longum p. cau

24.9

-0.7

1.9

25.0

brachioradialis

0.0

0.0

0.0

0.0

gleno-humeral contact force

346.6

79.9

-159.9

390.0

bone mass

0.0

-16.3

0.0

0.0

elbow contact force

495.4

45.6

81.1

504.1

The first three points concern modeling and simulation and do not impose major planning; the fourth point required the development of a concept, a method to follow (numerical analysis of histological sections, development of an algorithm and its implementation) and, a very important point, the definition of boundary conditions that appropriately represent the observed and measured conditions during the animal experimentation. A detailed description of how this was performed is given in the following:

2.2.1 The tissue differentiation model

2.2.1.1 The concept

Based on acquired knowledge the new differentiation model should contain the following characteristics:

▼ 52 

  1. The model was considered to be poroelastic (biphasic), contrary to previous monophasic analysis, axisymmetric and with a non-linear mechanical behavior. The cartilage properties were defined by permeability, void ratio and elastic modulus of Young’s. In this form cartilage was represented as an incompressible material with 80% water and elastic behavior. This approach then guarantees a more realistic simulation.
  2. The strain limits, which were used as stimuli for differentiation, were determined for the first time from histological analysis (in vivo data).
  3. The model to be developed should be so general that it does not depend on the mesh definition (localization or numeration of the elements).
  4. The values used to define the trilinear curve for each tissue type should be implemented in the algorithms for differentiation to allow evolution from one tissue to another.
  5. A theory describing this simulated healing process should explain how healing occurs and how the observations of histology can be reproduced.
  6. After the validation of the simulated osteochondral repair, the model should be flexible, allowing an easy predictive evaluation of different clinical cases.

2.2.1.2 Numerical analysis of the histological sections

During animal experimentation, histological sections of a surgically created osteochondral defect in the Yucatan minipigs (Fig. 2.2) were stained with Safranin-Orange van Kossa and Safranin-Light Green, and then analyzed for 4, 6, and 12 weeks.

2D contours between different tissues were digitized and defined as B-Splines using Pro/Engineer, a commercial package for solid modeling. With the usage of a developed PYTHON subroutine the B-Splines were imported as composite curves in Marc/mentat, (Fig 2.3). A mesh was created according to the histological sections to determine strains and stresses. The models were loaded with a pressure of 1.35MPa (from gait analysis) acting at the upper boundary of the model. Tissue types during healing were identified from the stained regions of histology. The material properties were taken from literature (Guo, 2001;Lacroix and Prendergast, 2002;Smith and Mansour, 2000) and then assigned in correspondence with these regions.

▼ 53 

Fig. 2.2: Initial osteochondral defect created at the femoral condyle in Yucatan minipigs (Ø = 6mm, depth = 2mm). Left: from: (Bail, et al., 2003), right: from: (Duda, et al., 2005).

A maximal of 6757 quadrilaterals 4-noded elements and 6834 nodes were employed to analyze the histological sections. A convergence test was performed by duplication of the element number and by changing from quads 4 node elements to quads 8 node elements. Since the differences in the strain fields were less than 10%, the models with the minimal number of elements and nodes were considered to be sufficient to evaluate the state of straining during healing.

▼ 54 

Fig. 2.3: Finite element analysis of histological sections to determine factor for growth and resorption. The contours of the different tissues in the histological sections are represented as B-Splines.

After FE analyses maximal and minimal values of the minimum principal strains were used to define the limits of the tri-linear curves given for each tissue type. A maximal strain value for a specific tissue represents in this case the limit for growth and the minimal strain value the limit for resorption. The rates of variation of these values were used to define factors for growth and resorption to simulate differentiation from one tissue type to another one. These new values were employed in the development of the tissue differentiation model described below.

This model is the first using in vivo data for the evaluation and calculation of parameters (limits of the trilinear curve and factors for growth and resorption) to simulate osteochondral healing. The current models generally use parameters to simulate growth and resorption obtained from sensitive parametric analysis using numerical models that is, the parameters are assumed - the model is numerically analyzed using such parameters and compared with histological observations. The parameters are then adjusted upwards to obtain a pattern similar to that observed in histology. To use parameters obtained from in vivo data represents the elimination of supposed parameters, and more importantly, the determination of how tissues respond mechanically under load application. Considering that the load is determined from the gait analysis, and the geometry is reconstructed from histological sections, its numerical analysis represents the in vivo straining of the tissues during healing, which are currently unknown.

2.2.1.3 Implementation of the tissue differentiation model

▼ 55 

To develop this model, a previous study of the physiological loads in human bones allowed understanding of their mechanical behavior, taking into consideration specific characteristics such as geometry and bone quality (Maldonado, et al., 2003). Knowledge of bone strength is important for the understanding of the origin of fractures as well as for optimizing fracture fixations in weak bones (Lauritzen, et al., 1993;Miller, et al., 1996;Saitoh, et al., 1994). Radiological analyses of cadaveric proximal humeri have found an increasing dependency of the rate of bone tissue loss on age, especially in the region beneath the epiphyseal scar, the central zone, and in the greater tuberosity (Hall, 1963). A more recent study has described a strong correlation between radiological measurements, such as high values of bone mineral density (BMD) and mechanical strength of the trabecular bone in the proximal humerus (Lill, et al., 2001). In addition, age related differences in the female and region related differences in both genders were demonstrated for the proximal humerus. These findings provide an insight into the fracture patterns of the proximal humerus and form the basis for a re-design of implants for osteoporotic patients.

Additionally, the humerus project showed that some mechanical conditions influence bone healing stronger than others (bone quality has stronger effects on healing than physiological loads), and therefore the different impact of mechanical conditions in the case of osteochondral healing must be considered as well. Is defect width more critical than defect depth? Is the healing pattern different when the defect thickness is varied? Do the specific mechanical conditions, generated after changes in the joint curvature, have any influence on osteochondral healing? What happens when in the case of usage of grafts to fill the defect, the stiffness is considered to be the same as the surrounding subchondral bone and what changes occur when the graft stiffness is reduced? To find answers to such questions the tissue differentiation model to simulate repair was developed and implemented.

To develop this type of model, some assumptions have been made and were considered: chondrocyte cells shall respond to mechanical signals especially under compressive loads. Due to the fluid phase of the cartilage, which amounts to approximately 80% of the cartilage composition, the effect of this second phase should be taken into account (biphasic material definition). Considering these two assumptions, the project was mainly aimed to develop a tissue differentiation model to analyze the influence of mechanical conditions on osteochondral healing to understand and even predict this healing process.

▼ 56 

Differentiation was simulated by incremental changes in the elastic Young’s modulus. Since chondrocyte activity appears to be principally stimulated by compressive loads (Heiner and Martin, 2004;Li and Herzog, 2004;Wong and Carter, 2003) compressive strains were used as mechanical stimulus to simulate and to maintain differentiation. Factors to simulate growth and resorption for each tissue type were defined after numerical analysis of histological sections. After load application, the stimulus was determined using the schema of a tri-linear curve conventionally used to simulate bone remodeling (Huiskes, 2000;Huiskes, et al., 2000;Ruimerman, et al., 2003;van Rietbergen, et al., 1993). For each tissue type a tri-linear curve was defined (Fig. 2.4; Table 2.2) in the following manner. At the initial state of the defect situation cartilage, cancellous bone, subchondral bone plate and defect (connective tissue) were simulated. In addition, the formation of fibrous tissue and calcified cartilage during healing was expected.

Fig. 2.4: Tri-linear curve used to simulate tissue differentiation. The center point corresponds to the strain value of an intact situation.

Using finite element models (FEM) selected histological sections at 4, 6 and 12 weeks were analyzed. From these numerical analyses two types of parameters to simulate healing were determined for each tissue type: the extreme points of the trilinear curves (F and D in Fig. 2.4) and parameters describing the rates of growth and resorption (Table 2.2). The extreme points were calculated from the average of the maximal and minimal values of the minimum principal strains (mps) measured at each time point (4, 6 and 12 weeks). The parameters describing the rates of growth and resorption were determined from the average of the rate of change in the minimum principal strains between these FE models of histological sections at consecutive time points. The material properties were taken from literature (Guo, 2001;Lacroix and Prendergast, 2002;Smith and Mansour, 2000) (Table 2.3).

▼ 57 

The mechanical stimulus was iteratively calculated as the difference between a current healing state and the state of strains in an intact situation. The new elastic Young’s modulus for each material point (defined as the numerical grid point of the mesh in the model characterized by its momentary tissue type) was calculated following the scheme shown in Fig. 2.5 and using equation 2.1. Finally, the differentiated tissues during the simulated healing in an area (TA) including the defect region and the elements at its basis were quantified, as shown in Fig. 2.6. The percentage (i.e. number of elements) of each differentiated tissue type during healing (after each iteration) was then calculated. The simulated healing was compared qualitatively with histological (location) and quantitatively with histomorphometrical data.

(Equation 2.1)

Where,

If εF < εcurrent < εC then V = Vgrowth,

If εC < εcurrent < εB, then stimulus = 0

If εB < εcurrent < εD, then V= Vresorp

The index A, B, C, D and F refers to points in Fig. 2.4.

▼ 58 

With, Vgrowth:

Factor for growth defined for each tissue type (Table 2.2).

Vresorp:

Factor for resorption defined for each tissue type (Table 2.2).

εA:

Minimum principal strains in the intact model at each material point. Point A in Fig. 2.4.

εB:

Minimal variation found in the mps after numerical analysis of histological sections during healing added to the mps in the intact model at each material point (εA). Point B in Fig. 2.4.

εC:Minimal variation found in the mps after numerical analysis of histological sections during healing subtracted to the mps in the intact model at each material point (εA). Point C in Fig. 2.4.

εcurr:Minimum principal strains in the defect model during healing in the current step at each material point.

▼ 59 

εD:Maximal value of the mps encountered in an intact situation added to the mps in the intact model (εA). Point D in Fig. 2.4.

εF:Maximal value of the mps encountered in an intact situation added to the mps in the intact model (εA). Point F in Fig. 2.4.

When the strain values exceeded the limits defined for each tissue (F and D in Fig. 2.4) the stimulus was set to the limiting value. Thus, a pentacurve was used.

▼ 60 

Table 2.2: Limits used in the tri-linear curves defining the differentiation process for each tissue type where Emin tissue < E elem≤ Emax tissue.

Defect connective tissue

Fibrous tissue

Cartilage

Calcified cartilage

Cancellous bone

Subchondral bone plate

Emin (MPa)

0.2

3.0

8.0

12.0

825

2300

Emax (MPa)

3.0

8.0

12.0

825

2300

22000

Width of the dead zone (CB)

0.1e-4

0.1e-3

0.1e-3

0.1e-3

0

0

Width of the differentiation region (FC; BD)

0.1295

0.1295

0.1295

1.536e-4

0.555e-3

0.555e-3

Factor for growth

1.2

1.15

1.0

1500

400

400

Factor for resorption

1.0

1.0

1.0

900

1500

1500

The factors for growth and resorption, required to simulate differentiation, were defined such that: 1. They were able to reproduce the resorption zone at the basis, observed in histology. The factor for resorption “resorption velocity” in an early state of the healing process should be then significantly higher than the factor for growth at the lateral wall and at the center of the defect. 2. After resorption the increased defect filled with a non-structured, untidy, connective tissue should possess a slow factor for growth. Due to this effect the stiffness of the reabsorbed bone region and the defect should be similar. To simulate these changes each material point was able “to jump” freely from one to another tri-linear curve, in which the factors for growth and resorption, the stimulus and consequently the elastic Young modulus were updated.

In a similar way, pore pressure and fluid velocity were analyzed as possible stimuli to start and to maintain tissue differentiation as well. Pore pressure was set to zero in a perpendicular direction to the model plane at the cartilage joint interface to allow fluid flow. The parameters for growth and resorption to simulate differentiation previously defined remained unchanged. The elastic modulus of Young at each material point was calculated and updated after each iteration. An iteration was defined as a simulating loadmaintained during 100 seconds. Each iteration was subdivided in 6 increments, in which the load was gradually applied. The material properties were iteratively changed at the last increment. The load was only applied after the change of all material points defined in the model. This whole cycle represents a load step. A maximum of 350 steps were required to achieve equilibrium.

▼ 61 

Fig. 2.5: Schematic representation of the tissue differentiation model. Differentiation was simulated by incremental changes in the elastic modulus of Young. Minimum principal strains were used as mechanical stimulus.

When in a load application step, the solution did not converge to the defined accuracy and the load step could, as a result, not be achieved in the initially defined number of increments (6); the increment of time was reduced using a predefined limiting value to assure convergence. The magnitude of this value was calculated using the corresponding equation specified in the ABAUS manual.

To control the numerical stability of the simulated healing, some elements were selected for graphical representation at the interfaces between the defect basis or defect wall with the surrounding tissues. Changes in the elastic Young’s modulus in these previously selected elements during healing were plotted and evaluated until equilibrium was reached. Smooth profiles representing tissue differentiation is an indicator of a normal development of the osteochondral ossification for cancellous bone and the formation of hyaline cartilage from connective tissue.

▼ 62 

In the annex number 1 a brief report about the frequent errors occurring during the run of the present code and how they can be solved is given.

2.2.2 Boundary conditions

To set appropriate boundary conditions for osteochondral healing it is necessary to consider the main results of animal experimentation:

  1. The defect localization influences the loads acting directly on the newly formed tissues. In the majority of the osteochondral healing studies, some authors selected the knee region to create osteochondral defects to analyze its repair process. Whereas at the patella and the lateral femur condyle low levels of compressive loads and high magnitudes of shear stresses are reported, high compressive loads have been measured in the medial femur condyle and the tibia plateau.
  2. The healing process depends strongly on depth and width of the defect (geometrical parameters). Defects grade 1 (ICRS nomenclature) with a surface smaller than 2 cm2 heal without the necessity of special treatments and in the majority of the cases the newly formed cartilage shows an acceptable mechanical stiffness. Defects grade 4, up to the subchondral bone, are able to heal only in a reduced number of cases, but are filled with fibrous cartilage with a low stiffness instead of the normal hyaline cartilage. When the superficial area is higher than 2 cm2 the defect remains unfilled and alternative clinical treatments such as the use of allografts and autografts are more appropriate.
  3. Another important point is the time after the surgery when loading of the affected joint is allowed. Some studies showed that after surgery immediate load application achieves a better healing than when the joint is totally unloaded or when an intermittent load is applied.
  4. Finally the age of the patient is a factor to be considered. Young patients heal faster and the newly formed tissues show a better mechanical quality. The rule of thumb is the younger the patient is and the smaller the defect, the better are the healing chances and the long-term duration of the regenerated tissue.

▼ 63 

All initial material properties used in the tissue differentiation model developed in this project were taken from literature (Guo, 2001;Lacroix and Prendergast, 2002;Smith and Mansour, 2000) Table 2.3. Taking into account that the geometrical limits of the joint were simulated, confined compression (see “cartilage characterization” from the introduction section) behavior was selected to reproduce the mechanical conditions acting at the knee joint. Their boundary conditions were modeled in the finite element model to simulate osteochondral repair. That means fluid exudation was allowed only between the articular surfaces and lateral cartilage deformation was avoided. Pore pressure at the joint curvature interface was set to zero in order to allow the movement of the fluid phase between the cartilage regions in contact. No fluid movement was allowed between the other tissues.

All models were loaded with a displacement applied in an axial direction at the top of the cartilage layer calculated in an intact model under a compressive pressure of 1.35 MPa. This displacement was incrementally applied during the first 100 seconds and maintained constant during healing simulation. Load was applied via contact with an intact joint surface. The other boundary conditions are shown in Fig. 2.6.

Fig. 2.6: Loads and boundary conditions used during healing simulation. TA represents the “total area” in which the percentages (number of elements) of the newly formed tissues were quantified during repair.

▼ 64 

Table 2.3: Initial material properties used in the finite element analysis (Guo, 2001;Lacroix and Prendergast, 2002;Smith and Mansour, 2000).

Defect connective tissue

Cartilage

Cancellous bone

Subchondral bone plate

Young’s module (MPa)

0.2

10

1750

20000

Poisson ratio

0.1627

0.1627

0.30

0.30

Permeability

k(mm4/N Sec)

1.0e-2

5.0e-3

2.019

1.0e-5

Porosity

0.8

0.8

0.8

0.04

2.2.3 Parameter study

After the development of the tissue differentiation model to simulate healing some cases of defect in joints were modeled and analyzed: osteochondral defects with different defect geometries (variations in defect depth and defect width), cartilage with different thickness, osteochondral defects in flat, concave and convex joint surfaces and defect fillings with different stiffness were evaluated. Using the method described above, the quantified tissues were calculated in an area (TA, Fig. 2.6) and compared during healing and when equilibrium was reached. This procedure was applied for all evaluated cases.

2.2.3.1 Analysis of the defect size

In a parameter study, the defect geometry was modified to reflect potential variations in the animal experimental configuration of the defect geometry. Three axisymmetric models were created. The total defect depth was increased from 1.6 to 2.4 mm, the defect diameter from 5.2 to 6.8 mm. Additionally the cartilage thickness was varied from 0.8 to 1.4 mm.

▼ 65 

Convergence tests were performed by increasing the number of the elements and by changing from quads 4-noded finite elements to quads 8-noded finite elements. Differences under 10% in the strain fields indicated that models with 4-noded elements were sufficient to simulate healing.

2.2.3.2 Influence of the local joint curvature

To analyze the influence of the local joint curvature on healing, eight simplified axisymmetric geometries were developed: three models (flat, concave, convex) to represent the intact situation and five models to simulate the injured joint (flat, concave, convex, biomaterial with 100% native bone stiffness, and biomaterial with 50% native bone stiffness). To simulate the joint after the above-mentioned convergence, test models with a flat interface consisting of 1452 4-noded elements and 1564 nodes were created. Cartilage, subchondral bone plate, cancellous bone and an additional defect region (connective tissue) were modeled.

Biphasic behavior was simulated using the ABAQUS soils capabilities (element type CAX4P). The tissue differentiation model was validated qualitatively by comparison with the patterns described histologically. Quantitative validation was performed by comparison of the percentages of different tissue types predicted to form with the histomorphometric analysis at 4, 6 and 12 weeks.

▼ 66 

To simulate osteochondral defects in concave and convex interfaces two models with a defect of the same dimension localized at each side of a curve interface were created. The concave situation was modeled with 1614 4-noded elements and 1721 nodes and a second model with 1603 4-noded elements and 1735 nodes was created to describe the defect on a convex interface. To resemble the tibio-femoral joint geometry of the animals, a radius of curvature of 15 mm was used (Fig. 2.6). The other parameters remained unchanged: material properties, tri-linear curves for each material, and factors for differentiation. Contact was modeled by a frictionless condition at the collinear nodes at the joint interface.

Pore pressure as stimulus for differentiation was analyzed in the convex interface, resembling the local geometry of the joint where the osteochondral defect was surgically created.

2.2.3.3 Analysis of the stiffness of defect fillings

To evaluate the effect of defect filling by means of different implant stiffness values, a cylindrical plug of 6mm diameter in a convex interface with a depth of 1.5mm from the bone plate and two different stiffness values of 100% and 50% of the elastic modulus of the cancellous bone were simulated.

▼ 67 

Generally, all finite element analyses were performed using the Marc/mentat package (Marc Analysis Research Corp.) and ABAQUS (Hibbit, et al., 2003) package on a Unix workstation (Silicon Graphics Inc., MIPS R 10000).

2.2.4 The developed algorithm

The tissue differentiation model was written in FORTRAN. In general the algorithm follows the scheme shown in the Fig. 2.5. The code employs the commercial package ABAQUS as a visualization interface. In this tool the osteochondral defect model is shown as a continuous representation from the initial state of the defect localization to the final state of healing, when the defect was completely filled. As explained above, differentiation is simulated by incremental changes of the elastic Young’s modulus.

The elastic Young’s modulus, in turn, is changed in dependence of a mechanical stimulus (compressive strains). Therefore, the algorithm uses as input data the minimum principal strains after equilibrium at each material point (total number of nodes) in an intact model. These input data are written in form of a formatted external ASCII file. This file contains information about the element number, the node number, the signal magnitude (minimum principal strains, pore pressure or fluid velocity) and the initial material name. These data are then used for defining the middle of the trilinear curve as shown in Fig. 2.4. The equilibrium of the intact situation is altered by the creation of an osteochondral defect. Aim of the algorithm is then to reestablish equilibrium by using a stimulus calculated after comparison of the current minimum principal strain with the minimum principal strains known from the intact situation. As explained above, the algorithm uses the parameters defined in the trilinear curves for each tissue calculated from histological sections to generate changes in the elastic modulus of Young.

▼ 68 

The input data was obtained using the utility subroutine UVARM from ABAQUS. This utility subroutine is called as many times as increments defined in the load step multiplied by the number of nodes and multiplied by the number of desired variables to be visualized. As explained above (description of the tissue differentiation model) one load step consists of 6 increments. The intact situation for a flat interface for example contained 1564 nodes. The variables represent the values that should be known after the FEM analysis. In this case six variables were of interest: three principal strains (maximal, medial and minimum), two fluid velocities (fluid velocity in the X direction and fluid velocity in the Y direction) and pore pressure. In consequence the preanalysis subroutine for this model was called 56304 times (1564 nodes * 6 increments * 6 variables). It is important to know this value in order to set the required variables in the environment files of an ABAQUS user appropriately (see Annex 1). These variables are known as user variables for material information. The variables are calculated, stored in memory after each increment and written in an independent formatted external ASCII file for each variable. In this form after preanalysis of the intact model the corresponding signal information (minimum principal strains) for each node and iteration was recorded. Additionally the element number, node number, increment number and material name were added to this file.

The “user subroutine” (Annex 2) written to simulate healing can now be used since the minimum principal strains and the tissue type for an intact situation are given. After variable definition and dimensioning of necessary arrays, this subroutine contains three principal parts: in the first part, all tissue types (initials and expected to be formed) are defined by limits of maximal and minimal values of Young’s modulus. In the second part all necessary values to define the points conforming the trilinear curves for each tissue and the factors for growth and resorption are given. The third part contains the rhombuses of decisions, in which the current strain values are analyzed and the elastic modulus of Young is calculated and updated after each iteration (schema shown in Fig. 2.5). The last part is related to the manipulation and storage of each variable used in the algorithm.

The finite element method solves the integral equations generated for the osteochondral defect model with the method of the displacements using an algorithm of optimization. This optimization procedure changes the mean bandwidth7 in the solution matrix. Thus, the solution for each node is made following the sequence suggested by the optimization procedure. In this form the ascending numerical sequence defined by the mesh topography could not be used. Therefore, in the algorithm a mechanism was implemented to verify that the current node to be analyzed is the same one as read from the intact situation. To guarantee this, some additional files are created to compare and to write element and node information.

▼ 69 

The elastic Young modulus, which represents the grade of healing at each step, was associated with a state variable (SDV in ABAQUS). A state variable implies that its associated value is stored in a temporal memory field and could be used in a mathematical formula in combination with other numbers or state variables. A state variable can be iteratively saved in memory when its value is associated to a field variable. In this case, its value could be used iteratively in the input model data. Thus, the elastic Young modulus is calculated as a state variable, saved in a field variable and used to show the history of the material stiffness during healing simulation. The input data of the model is thereby updated and will be rewritten after each iteration. However, changes associated with each state variable are performed at the last increment, when the total compressive load applied on the joint during a step was completed. During load application the fluids inside the tissue were exuded to attenuate the compressive load. When a load cycle is completed, new elastic Young’s modulus are calculated for each tissue type. Poisson’s ratio and permeability can then be changed. Thus, the values of these mechanical material properties are adopted in dependence to the tissue type. Tissue differentiation stops when minimum principal strains and tissue types are the same in the defect model as in the intact one. Thus, when a tissue achieves values of the dead zone, the stimulus is set to zero. The last Young modulus is then maintained unchanged during the rest of the healing simulation.

Tissue differentiation is regulated through changes in the elastic modulus of Young principally in accordance to the quantity of stimulus required to achieve the strain field encountered in an equilibrium state. So, this process is started and maintained until all tissues achieve the strain field previously determined in an intact situation. The “center point” of the trilinear curve is the minimum principal strain value for an intact situation in a node with the same coordinates as in the defect model. Obviously the mesh for the intact model and the defect model needs to be identical. Healing can then be described as the intention of convergence of an injured joint to achieve a healthy state represented by the intact model under determined mechanical conditions. As a result, it proved to be feasible to program chondrocyte activity as a biological response to a new local mechanical environment generated when articular damage occurs.

Only for flat a joint, implementation of a treatment of numerical viscosity at the elements conforming the defect wall was necessary. Numerical viscosity reduces possible inconsistencies of Hook’s law at the interface defect wall-host tissues. Such inconsistencies are commonly attributed to singularities in the strain values after load application between two materials in contact, whose elastic Young’s modulus shows a drastic difference (in its values). This numerical instability can create a shock transition zone. In the flat model the instability was formed not only due to a jump in the material properties at the interface between the connective tissue – cancellous bone (0.2 vs. 1750 MPa) and connective tissue – cartilage (0.2 vs. 10MPa), but also due to a high strain concentration at the defect vertex. Numerical viscosity is a well-documented procedure applied in several cases when the discontinuity is detectable. Numerical viscosity was then applied to smooth the mechanical response of the affected region (3% of the total number of the elements of TA, Fig. 2.6).

2.2.4.1 Additional tools developed for analysis of the data

▼ 70 

Some additional tools including mesh generation and treatment of the results were created during the project:

A PYTHON routine was written to convert the IGES data information generated from PRO/Engineer to represent the histological sections as B–Splines to composite curves readable in Marc/mentat. This information was used to perform the finite element mesh in Marc/mentat.

A PYTHON routine was written to convert the mesh information (elements with its node connectivity, nodes coordinates, and material sets) created from Marc/mentat to a readable format to be used as input data in ABAQUS.

▼ 71 

In PYTHON language a routine was developed to write changes of the elastic modulus of Young, calculated after each iteration in each material point, in an external ASCII file.

A PYTHON routine was written to create the steps information automatically to the input data of the defect model to be analyzed.

In MATLAB a short routine was developed to read this data and quantify each tissue type after each iteration.


Fußnoten und Endnoten

7  Mean bandwidth: denoted by B. If a matrix of order N stored in accordance with the skyline format (skymatrix) use S memory locations, the ratio B=S/N is called the mean bandwidth.



© Die inhaltliche Zusammenstellung und Aufmachung dieser Publikation sowie die elektronische Verarbeitung sind urheberrechtlich geschützt. Jede Verwertung, die nicht ausdrücklich vom Urheberrechtsgesetz zugelassen ist, bedarf der vorherigen Zustimmung. Das gilt insbesondere für die Vervielfältigung, die Bearbeitung und Einspeicherung und Verarbeitung in elektronische Systeme.
XDiML DTD Version 4.0Zertifizierter Dokumentenserver
der Humboldt-Universität zu Berlin
HTML-Version erstellt am:
15.11.2007