ABSTRACT

Treatment of bone tumors in the mandible often involves extensive excavation of affected bone, followed by mandibular reconstruction. Prosthetic implants may be needed to restore jaw functionality. The challenges of making prosthetic bone implants include stress shielding and extending the mechanical life of the implant. We have developed a design algorithm to improve the efficiency of prosthesis design. A finite element model of the patient case is constructed from a computer tomography scan, and the computer implements topology optimization techniques to design the prosthesis with limited stress shielding affected by highly biomechanical compatibility. Topology optimization facilitates the design of low weight structures by automatically introducing holes into the structure. This is governed by engineering predetermined constraints to meet certain job specifications. Such a design will be tested for fatigue life before it is ready to be manufactured and used. Topology optimization can be performed as a design process to achieve a final design that takes stress shielding into consideration. The problem of stress shielding is solved by matching the stiffness of the orthopedic implant to the original bone that is being replaced. The material we used was titanium alloy (Ti-6Al-7Nb). Volume fraction of the orthodox implant was used (0.2872 for the studied case) as volume constraints. Compliance of the bulk bone was set as a further constraint to match the stiffness of the bone with the designed structure. Our results show a good life expectancy for the designed parts, with 12% higher life expectancy for stress-based topology optimization than for compliance-based topology optimization.

Keywords:

Topology Optimization, Temporomandibular Joint, Fatigue, Stress Shielding

1. Introduction

Reconstructive surgery of the mandible using 3D printed components is becoming more common, especially with the advancement in the additive manufacturing process and computer simulation techniques. Bone is a connective tissue composed of cells, fibers, and ground substance supported by an extracellular matrix. Mechanically, bone is considered to be a brittle material, with the brittleness varying with age, sex, and conditions [1] . Medical treatment may also cause the bone to become more or less brittle. Tumors in the mandible, such as osteosarcoma and ameloblastoma, need to be removed [2] , and immediate replacement is required to restore mandibular function. The design and choice of materials used in open reduction are based on the surgeon’s “common sense and experience” [3] [4] [5] [6] . Mechanical analysis using finite element analysis can assist surgeons in making optimal choices. Because the bone structure in the head and neck is anatomically complex, facial reconstruction is often lengthy and may require surface scanning of the face. The face is filled with air cavities, soft tissue spots, arteries, and facial nerves; and individual variation also needs to be considered. Computed tomography (CT) imaging is an effective way of constructing three-dimensional (3D) models of the design domain [7] , and its use lowers the probability of malunions and malocclusions. This modeling technique also limits nerve injury within the design domain [8] [9] . Critical locations can be identified and set as conditional design domains. Design using topology optimization sometimes leads to rather complicated design; additive manufacturing could assist with simplifying designs. Automating custom orthopedic implants for use in real-time surgical operations is a distant target at the present time, but with the advance of computer-aided processes with topology optimization, it could be achievable in the future. The domain would be optimized by both the surgeon and the engineers to lower risks before and after surgery, such as infections in elevated spots and soft tissue collateral damage, and to achieve the best available fixing and enclosure. The desired cavities and appropriate loads can be taken into consideration more efficiently [10] . Several cases of optimized plate fixation for fractured mandibles have been reported. The loading condition simulates a real joint in the degree of freedom and actual closure, in addition to muscular tension forces [11] . Bio-tribological characteristics are beyond the scope of this research. Computer-aided design (CAD) provides a good database for greater precision of work piece making. Additive manufacturing controlled by computer-aided manufacturing (CAM) protocols provides control of the mechanical and topological aspects of the implant. The target of this paper was to introduce an optimal design protocol to simulate and test any kind of reconstructive part. Such a design methodology tries to mimic the patient’s actual oral cavity anatomy and bone structure topography. Mimicking the patient’s actual bone topography may help resolve occlusion problems. For example, slight overfilling of a tooth, even by a fraction of a millimeter, can affect normal occlusion. This could lead to a malocclusion that causes chronic headaches or a temporomandibular joint disorder. Stress shielding is a problem that needs to be addressed. Structural behavior matching is a solution presented in this paper, by matching the compliance of the designed part with the original presumed form of the bone. The prosthesis lifespan is studied and incorporated in the design algorithm.

2. CAD-Based Prosthesis Design and Fatigue Life Evaluation

Topology optimization is a design method based on a mathematical scheme to introduce nonexistent material within the design space to incorporate design properties such as weight minimization or certain mechanical properties. Michell [12] studied truss design problems to determine the generalized procedure to optimize truss structure for weight minimization, opening the door for design automation. Rozvany et al. [13] studied the optimal layout theory of topology optimization and later introduced perimeter control. Bendsóe and Kikuchi [14] introduced the homogenization method into topology optimization. Bendsøe [15] later developed the solid isotropic microstructure with penalization (SIMP) method. Penalized density is usually considered as a design variable. The direct approach (also known as the artificial density approach), represented by a SIMP method, has been widely used in recent decades. SIMP is a scheme that appliesa discretized design domain to find a smooth optimal structure in which material properties are set to be constant for the discretized domain. Aside of SIMP method, there is the explicit topology optimization approach [16] , which used for solving optimization problem with reduce the need for penalize density [17] .

The solution in the scope of the current scheme faces some challenges, including checkerboards and nonexistence. Nonexistence is a mesh dependency problem that tends to extensively introduce a non-existent element to satisfy the solution with a decreasing value of the objective function. To solve this problem, a relaxation principle is introduced, modifying the density function with a so- called gray region. Gray regions, together with heuristic searching with sufficient constraint, are the remedy. Checkerboards are also byproducts of discretization. The solution tends to be non-applicable because of mesh dependency. With advancement in the design and application of mechanical structures, especially in weight and bulk minimization, static consideration of engineering problems is rendered as unsuited for covering intended functionality. Reparative ergodic loading is a natural behavior of bones. Cyclic deformation [18] [19] [20] and strain life approach (ε-N) [21] is suitable for testing designs and evaluating expected life, especially when considering cyclic hardening. Because of the shape complexity of prostheses designed by topology optimization, rapid prototyping is expected to be used in the machining of such structures. Additive manufacturing using laser 3D printing gives an acceptable surface finish with metal constancy because of the unique characteristics of laser welding. Surface finish is a very important factor determining the lifespan of a printed piece. Microstructurally, residual stresses caused by laser thermal processing are minimized, because of the high heat dissipation of localized heat. Under such assumptions, the surface finish and residual stress will not be highlighted in the fatigue model. Average stress should be used in a 3D case. The signed von Mises stress is a suitable equivalent state for titanium alloys.

3. Stress Shielding Problem Approach

Connective tissues have certain properties that may vary slightly from one patient to another and may change over time. Because it is a composite material, the mechanical properties of osseous tissues are governed by the topology distribution of the materials, such as the microstructure of spongy bone, and the material and its distribution within the composite itself. The mechanical properties of bone stem froman inorganic matrix (consisting mostly of mineral salts) and collagen (with elastic fibers). Because of the nature of live tissues, the mechanical properties change with environmental variation. Wolff’s law is an example of bone adaptation under different loading conditions. In the case of orthopedic implants, the greater stiffness of metallic implants causes stress shielding [21] . There are several engineering techniques used to decrease stress shielding, and these techniques vary with each case. One solution involves the temporal aspect of the implant, by using temporary or biodegradable implants [22] . In some cases, the temporal aspect involves a long-term implant. In such cases, matching of the mechanical properties is a solution; this is achieved by selecting a material with similar properties that is biomechanically compatible. Porous structural materials can be designed from bulk materials such as orthodox bulk implant alloys (e.g. titanium alloys) [23] [24] . Stiffness matching of the designed implant structure for the loading conditions is another solution. Stiffness can be used as a design constraint for the implant, by calculating the stiffness of implant topography of bone as material, such conditions can be introduced into the finite element as its opposite (i.e. compliance $C=F^{\text{T}}U$ . In other words, compliance of the complete bulk design area $\left(\text{Ω}\right)$ of bone as the material is the condition for the real alloying material (Ti alloy) of the topology optimized for $\text{Ω}$ . Stiffness matching of the overall structure was used in this study.

4. Optimization Objective Functions

4.1. Compliance Objective Function

Addressing the generalized Hooke’s law ( $\sigma =E\epsilon$ ), yet still obeying SIMP criteria, is occupying domain Ω. Generally, Ω is considered to be a combination of boundary domains Ω_b, non-design domains Ω_f (whether it is a nonexistent domainor a forced existence domain), and design domains Ω_d. Mechanical properties are proportional to the powered density function $\rho$ such that:

$E\left(\rho \right)={\rho }^{3}E$ $$ (1)

The displacement-based finite element discretization is:

$K\left(\rho \right)U=F$ (2)

Here, K is the stiffness matrix depending on the density function ρ. U is the nodal displacement vector, and F is the nodal force vector. Stiffness is the measure of the ability of the structure to withstand certain loads. It can be introduced in terms of finite element analysis as its opposite (i.e. compliance). Therefore, minimizing compliance will mean maximizing stiffness. The optimization problem based on compliance is:

$\begin{array}{l}\text{minimize}C=F^{\text{T}}U\\ \text{s}\text{.t}\text{.}0<\rho \le 1\\ {\displaystyle \int \rho \text{d}\Omega }\le V_{\max }\end{array}$ (3)

Here $V_{\max }$ is volume fraction constraint.

Increasing the stiffness of the design domain with decreasing mass is the optimal goal for most design problems. Such a target could be achieved by minimizing the compliance [25] . Numerical discretization of the problem is widely used and is a vital aspect in topology optimization such as finite element modeling (FEM). Numerical representation of the initial state, variable behavior, and aggregation problems (checkerboarding) have been studied intensively in recent decades [26] [27] [28] [29] [30] . Compliance provides good results even while neglecting the stress consignation facts. Intense stress spots in a design could easily cause it to fail because of static and dynamic loads. A stress-based optimization method is required to achieve an acceptable design for certain applications.

4.2. Stress-Based Objective Function

The type of stress can be analyzed by considering the theories of elastic failure. In this case, the stress is maximum octahedral shear and we have used the von Mises stress theory. An enveloping stress tensor can be used to point out topological aspects associated with certain constraints. Singularity, localization, and nonlinearity are the most serious challenges for researchers in the topology optimization field [31] . The first step in stress-based design is determining which stress state should be considered for a certain application. Each application requires consideration of the nature of loading and the properties of the material in use. Stresses to be directly addressed are determined by the worst-case scenario, which is failure. The nature of the failure criterion is the scaling condition for stress-based design and the nature of the materials. Maximum distortion energy theory (von Mises ${\sigma }_{vms}$ ) [32] is the failure criterion most commonly referred to in the literature. It has been described asa 3D stress tensor in two categories: von-Mises Equation (4) and signed von Mises Equation (5)

${\sigma }_{vms}=\frac{1}{2}\sqrt{{\left({\sigma }_x-{\sigma }_y\right)}^2+{\left({\sigma }_y-{\sigma }_z\right)}^2+{\left({\sigma }_z-{\sigma }_x\right)}^2+6\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)}$ (4)

${\sigma }_{vms|Sign}=Sign\left(I_1\right)\cdot {\sigma }_{vms}$ (5)

where $I_1$ is the first invariant of the stress tensor.

Using p-norm is adopted in this work. p-norm for Von mises stress take the form

${\Vert {\sigma }_{vms}\Vert }_P=\sqrt[P]{{\displaystyle \sum _{i=1}^N{\left|\frac{{\sigma }_{vms}}{{\sigma }_{\max \_d}}\right|}^P}}$ (6)

Here, ${\sigma }_{\max \_d}$ is the maximum allowable von Mises stress.

Theoretically, efficient optimization could be achieved for a higher value of p-norm power if unlimited computational power was available and if there were no time constraints [16] . The objective function and boundaries would take the form:

$\begin{array}{l}\text{minimize}{\Vert {\sigma }_{vms}\Vert }_P\\ \text{s}\text{.t}\text{.}0<\rho \le 1\end{array}$ (7)

Here, ${\sigma }_{\max \_d}$ is the maximum allowable von Mises stress.

${\Vert {\sigma }_{vms}\Vert }_P$ is widely known as the P-norm function [18] [31] or the KK function as refined by Park and Kikuchi when associated with appropriate constraints. Because of the discretization nature of topology optimization, mesh quality and type play a vital role in the pre-and post-processing of design. As mentioned previously, stress-based topology optimization (SBTO) is affected by FEM accumulative analysis history. High order elements may increase the theoreticalodds of better design, higher resolution design, and increased element geometric density. However, computational and time costs may be a serious problem along with convergence. Living tissue FEM analysis is hardware consuming and time consuming in itself, because of the complex nature of loading and topographical representation. Additional factors are the number of constraints, and the complexity of stresses in nature, which act as objective functions.

5. Sensitivity Analysis

Sensitivity analysis plays a major role in achieving converging results while minimizing computational and time input. First order sensitivity analysis is required to be performed for each iteration. The adjoint variable method is used to develop a unified formulation for representing response variation in terms of variation design. Aggregation of objective functions are subjected to constraints with respect to design variable ρ. Compliance sensitivity analysis using the adjoint method gives Equation (8):

$\frac{\text{d}C}{\text{d}\rho }=-U^{\text{T}}\frac{\text{d}K}{\text{d}\rho }U$ (8)

The p-norm derivative is derived to form Equation (9)

$\frac{\text{d}{\Vert {\sigma }_{vms}\Vert }_P}{\text{d}\rho }={\lambda }^{\text{T}}\left(\frac{\text{d}F}{\text{d}\rho }-\frac{\text{d}K}{\text{d}\rho }U\right)+\frac{\text{d}{\Vert {\sigma }_{vms}\Vert }_P}{\text{d}\rho }$ (9)

Here λ is the adjoint variable, which can be calculated from Equation(10)

$K^{\text{T}}\lambda =\frac{\partial {\Vert {\sigma }_{vms}\Vert }_P}{\partial {\sigma }_{vms}}\cdot \frac{\partial {\sigma }_{vms}}{\partial \sigma }\cdot \frac{\partial \sigma }{\partial U}$ (10)

6. Numerical Examples

An in vitro case study was simulated and studied involving a major excavation of the mandible, simulating jaw cancer, as shown in Figure 1.

First order tetrahedron element was used as three-dimensional finite elements for simulating the jaw problem. Each node has three degrees of freedom. Which leads to 12 degrees of freedom for the element. Optistruct solver was used to solve the finite element problem. Temporomandibular (TM) prostheses are commonly used in major reconstructive surgery. An orthodox reconstructive prosthesis based on a design simulation, using a Vitek-Kent prosthesis body with an AO/ASIF TM joint condylar prosthesis, is shown in Figure 2.

The core is the design domain. Bone density varies from 1.5 gm/cm³ to 0.9 gm/cm³ compared with approximately 4.53 gm/cm³ for titanium alloy. Such a difference necessitates coverage of a smaller region of mandibular space to overcome excessive weight. In this case study, the orthodox prosthesis occupied 28.72% of the original replaced bone. The prosthesis material used was titanium alloy (Ti-6Al-7Nb) [33] . To perform finite element analysis, constraints and loads on the mandible need to be defined (Figure 3). Individual variations in biting condition have been ignored. Suitable adaptive loading values are listed in Table 1 and Table 2 based on the work of Koseki et al. [34] .

The design problem that is considered in this work is divided into two regions: the design in which the optimization process will occur; and the non-de- sign domain which is needed to satisfy medical conditions, such as an occlusion mismatch. The implant is made of materials with a high modulus of elasticity, such as titanium alloy. The implant domain is within the boundary domain of the lower modulus of elasticity (bone). This boundary domain will eliminate the need to extract a design model for certain applications as it is being done within the conventional topology optimization design process [35] [36] [37] . The design takes into consideration the original state of the finite element analysis of the mandible. Objective function has been investigated and was compliant and p-norm. The power of p-norm was 50, so it was suitable for the calculation power of the computer we used, creating a reasonably good functional topography for the design norm. In the present computer simulation, a major temporomandibular joint prosthesis reconstruction of the mandible should be done for most of the left-side damaged jawas shown in Figure 4(a). This case simulates resection of the jaw for treatment of bone cancer, necessitatingan immediate large-scale implant. Because of the complexity of the mandibular anatomy, the reconstructive prosthetic part is shaped based on information from the CT scan. The implant consists of a non-design domain and a design domain of titanium alloy. The purpose of the non-design domain is to mimic the actual replaced bone shape to avoid malocclusion. Malocclusion can lead to chronic headaches and serious joint inflammation if it is not treated. Cloning of the actual shape will certainly eliminate malocclusion. As shown in Figure 4(b), the full prosthetic outer shape has been determined by the non-design domain, while the design domain is within the core of the implant.

The target volume fraction in topology optimization was set to be equal to the volume fraction of the orthodox prosthesis compared with the original bone; in this case, it was 0.2872. To maintain the structural integrity, the pre-design space is skinning the design area as well as a presumable place of fixing pins. As a solution to stress shielding, an additional constraint was added which is compliance constraint to topology optimization based on compliance and SBTO functions. It is necessary to calculate what compliance (C) would be for healthy bone structure. This compliance would be an additional constraint to the upper boundary of the volume fraction. In the case of compliance maximization based objective function, the design scheme follows Equation (3). In the case of minimizing the averaging stress aggregative equation, i.e. p-norm, the design scheme

follows Equation (7). In this case, the compliance condition which is the compliance of the missing bone in case it exists, was measured as 34.04 of 247.27 (the total compliance of the complete healthy mandible).

A fatigue life estimation was performed for the designed prosthesis. The strain-life method was adopted for its accuracy and suitability for the loading conditions. Loading variation was considered as full revisers of 10% of the overall maximum strain, linear variated within the simulated biting time. Modeling of the strain-life is calculated as in Equation (11).

${\epsilon }_a=0.01538157{\left(2N_f\right)}^{-0.095}+0.35{\left(2N_f\right)}^{-0.69}$ (11)

The suggested optimal design procedure started with imaging the case to create a 3D FEM model (Figure 5). In the case of a pre-existing hole, jaw topological similarity can be used to create missing parts within the FEM model.

After constructing the topographical model, the dentist will identify the muscular orientation and move on to then on-design domain associated with some desired design features. Topology optimization will be conducted. Finally, fatigue simulation is under taken as a reference to evaluate the design. Now the design is ready to be 3D printed to use in surgery. The medical and engineering points of view are separated during the whole process to theoretically maximize the design speed. The FEM design and setting of boundary conditions requires a knowledge of anatomy (major aspects of the organ mechanical boundaries), and topological aspects of making holes for blood vessels or marking dangerous places to be avoided (nerve clusters etc.).

7. Results

3D design of the design domain was conditioned to approach bone compliance. Different powers of p-norm showed similar results for the chosen volume

fraction. Results of the higher power (50) were chosen to ensure approaching (Figure 6).

The objective function history for both objective functions is summarized in Figure 7 and Figure 8. Optimization in the case of SBTO was achieved by minimization of the norm function, while in the case of the compliance objective function, the metallic prosthesis stiffness was minimized by maximizing compliance itself to meet bone compliance.

Figure 9 shows the location of the maximum von Mises stress within the designed prosthesis.

Table 3 shows the decrease in life expectancy of the prosthesis implant designed with SBTO compared with the compliance function based topology optimization. Fatigue assessment based on the strain life approach takes into consideration strain hardening within the complicated shape of the produced design

domain as shown in Figure 10.

The SBTO design shows almost 100% improvement in life expectancy for the reconstructive prosthesis compared with the orthodox design. However, the orthodox design of the reconstructive prosthesis usually marginalizes stiffness compatibility. This is because of the lack of individuality in general surgical solutions, so in this case, compliance of the orthodox prosthesis is 107.513, and normal bone compliance (34.04).

8. Conclusion

Stress shielding is a highly anticipated in orthodox reconstructive prostheses, which can be shown in the difference in compliance between the original bone and the implant. Targeting bone compliance as a boundary condition is the method we used in this paper. Stress shielding problems according to stiffness matching is theoretically minimized. Occlusion mismatching is a challenging problem in maxillofacial surgery. By making the design domain match the exact bone topography, high precision machining can give an exact replica of the jaw part to achieve an accurate occlusion. Design using topology optimization gave along-life expectancy in the simulation process. The compliance objective function achieves good results for fatigue life prediction, but stress-based topology optimization achieves better results. The orthodox-based design has a shorter life expectancy than the new optimized designs. Compliance of the orthodox design is almost 210% differ from the original replaced bone. Therefore, stress shielding of orthodox designs is highly desirable. From the design time aspect, maximizing compliance can be considered to be a faster strategy.

Acknowledgements

We thank Helen Jeays, BDSc AE, from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

Cite this paper

Al-Ali, M.A., Al-Ali, M.A., Takezawa, A. and Kitamura, M. (2017) Topology Optimization and Fatigue Analysis of Temporomandibular Joint Prosthesis. World Journal of Mechanics, 7, 323-339. https://doi.org/10.4236/wjm.2017.712025

References