The mechanical properties of human dentin for 3-D finite element modeling : Numerical and analytical evaluation

Material and methods. Analytical and numerical dentin homogenization according to Luciano and Barbero was performed and E-modulus (E), Poisson’s ratios (v) G-modulus (G) were calculated. The E-modulus of the dentin matrix was 28.0 GPa, Poisson’s ratio (v) was 0.3; finite element models of orthotropic and isotropic dentin were created, loaded and compared using Ansys® 14.5 and CodeAster® 11.2 software.

Dentin is the main mineralized hard dental tissue supporting structure of the tooth.Its mechanical properties are related to microstructural alterations due to ageing, carious lesions, and restorative procedures.
][3][4][5] The FEM enables the modeling of very complex structures (e.g.8][9] However, FE models of the biological structure are usually performed using a simplified approach.The employed method may play different roles depending on the model type, physics-related problems or the investigation's objectives.][5][6][7] Expanding such a method into biomechanics requires a significant and continuous improvement of modeling techniques with respect to the individual anatomic and pathological characteristics of human tissues. The micromechanical approach enables the prediction of the effective elastic properties of the material, based on its microstructure and the material properties of its components; this process is called homogenization.Numerical homogenization is based on the concept of Representative Volume Element (RVE).RVE is a statistical representation of the microstructure on a macro level.In addition, if the microstructure shows a certain periodicity, as in the case of dentin, the periodic boundary conditions of RVE can be imposed. 10One alternative to the numerical method is to use an analytical method in solving the micromechanical problem.Many analytical techniques of homogenization have been based on the equivalent Eigenstrain method, which considers the problem of a single ellipsoidal inclusion embedded in an infinite elastic medium. 11Among such analytical methods, one of the most efficient is the Mori-Tanaka method, which takes into account the interactions between inclusions.The main advantage of the Mori-Tanaka method is the delivery of an explicit formula for effective stiffness of the tensor and for local stresses and strains. 12,13his approach provides the opportunity to simplify the numerical model and, at the same time, takes into consideration the orthotropic nature of dental tissue.Homogenization enables the prediction of effective elastic properties by considering the mechanical properties of the dentin matrix and dentinal tubules fraction.It can take into account the location-dependent variations of dentinal tubules number and diameter.
Despite many years of development and application of the tooth numerical analysis, anisotropy of both dentin and enamel is neglected.Basically, it can be assumed that the simplified impact on the results of strain distributions remains unknown.It can be regarded reasonable to develop the research methodology reflecting the complex structure and nonlinear properties of dentin during FE analyses and identify how tissue anisotropic properties affect strains and displacement distributions in the complex system of materials with different strain-stress characteristics.
The 2 purposes of this study were to evaluate the degree of anisotropy and determine the qualitative and quantitative elastic properties of a location-dependent homogenous material equivalent to dentin and to calculate the differences between the isotropic and layered orthotropic models of dentin.

Material and methods
Two experiments were carried out in order to obtain an equivalent homogenous location-dependent material model of human dentin.The whole dentin, due to its composed structure, was divided into 3 layers: top (superficial external 1/3 dentin close to enamel), middle (middle 1/3 dentin) and bottom (internal 1/3 dentin close to dental pulp).
The elastic properties of the entire dentin and its layers were numerically and analytically determined by the homogenization method, and differences between the isotropic and orthotropic models of the dentin tissue were analyzed.

Analytical method of dentin homogenization
The proposed analytical homogenization of the dentin structure was based on the equivalent Eigenstrain method, and took into account the periodic boundary conditions. 11The dentin structure was assumed to be a transversely isotropic medium with one axis of symmetry, parallel to the direction of the dentinal tubules. 14,15n this case, the micromechanical model used for homogenization was a unidirectional structure with long cylindrical voids embedded in an isotropic matrix.The micromechanical model allowed the strain concentration tensor and corresponding effective stiffness tensor (material properties) of the dentin structure to be determined according to Luciano and Barbero: 10 Where: P is a 4 rank tensor, which describes the geometry of the voids, microstructures, and the distribution of the voids, V f the volume fraction of the voids (lumen tubules), C f , C ms the 4 rank stiffness tensors of the voids (lumen tubules) and matrix (dentin), both isotropic.

Numerical method of dentin homogenization
The prediction of effective material properties of a dentin sample was conducted with the help of the RVE concept, which enabled the modeling of periodic symmetry structures. 10For the numerical homogenization procedure the FEM was used to model the RVE of the dentin structure (Fig. 1).The numerical homogenization procedure was performed with the help of the FEM with additional consideration of the periodic boundary condition on the RVE.The periodic boundary conditions were imposed on a representative FE-model of dentin structure, according to Luciano and Barbero, as followed: 10 Where: u i +j ,u i -j was the imposed displacements of the opposite sides of the RVE cell (side +j and side -j), ∆x k j the distance between the opposite sides of the RVE cell, e ik the mean strain in the assessed RVE cell.
This formulation enabled constant distribution of the strain over the whole RVE to be enforced for each direction.At the same time, the strains in other directions were restrained, which enabled the calculation of average stress.This procedure had to be performed 6 times in order to determine all the independent components of the stiffness tensor of the transversely isotropic material.This meant  ( , , ) ( , , ) that the 6 independent boundary values problems needed to be solved in order to get a full description of the dentin material.The numerical homogenization procedure was performed with the help of the FEM and Ansys ® 14.5 (Ansys, Inc., USA) and CodeAster ® 11.2 (EDF) software.
The same input parameters were used for analytical homogenization (Table 1).Young's modulus (E) was equal to 28.0 GPa, according to Kinney et al. 21Poisson's ratio was equal to 0.3.Three FE-models were prepared: one for each of the top, middle and bottom layers of the dentin structure (Fig. 1).Each model reflected the specific geometry parameters, such as the average diameter of dentin tubules and the volumetric fraction of dentinal tubules in the whole dentin structure.
The results of such analytical and numerical homogenization were used to calculate the anisotropy (%) of the dentin structure as followed: Where: a is anisotropy (relative), E1 the reference value of Young's modulus, Ex = E2 = E3 value of Young's modulus (perpendicular to E1).

Experiment 2: Finite element analysis of the dentin model
The 3-dimensional FE models consisting of a simple cube of bulk dentin were created and analyzed using Code Aster 11.4.The cube had 3 mm sides and was divided into 3 equal parts to represent the top, middle and bottom layers of the dentin structure (Fig. 2).An FE simulation of the evaluated material properties was performed for the orthotropic material model and the isotropic simplified material model (2 models were analyzed).With both models (iso-and orthotropic), 3 load types were analyzed: (1) a shear only with a shear load of 1.0 MPa applied in X direction on the top surface (the bottom surface nodes were fixed in the X-direction, other nodes were fixed on the Y and Z axes), (2) a vertical compression of 1.0 MPa applied towards the Z axis on the top surface (the bottom surface nodes were fixed in all directions), (3) shear with bending with a shear load of 1.0 MPa applied on the top surface, with bottom surface nodes fixed in all directions and other nodes free in the 3 axes.The input data for human dentin were presented in Table 2.The data for the orthotropic model obtained during experiment 1 was used.

Experiment 1: Dentin homogenization
The results of analytical and numerical homogenization were presented in Table 3.An example of the results of the FE analysis of the micromechanical model of the dentin structure was shown in Fig. 3.The effective material properties of the dentin structure were determined for the top, bottom and middle dentin layers (from the pulp up to the dentin's exterior surface).A description of the location of dentin layers was presented in Fig. 1.

Experiment 2: Finite element analysis of the dentin model
The displacements and strains recorded on the cross section of the model's surface for the isotropic and orthotropic models at different loading conditions were presented in Fig. 4.An analysis of displacements showed higher values in the isotropic model.The differences reached up to 16% by shear load (Fig. 4a), 37% by compression (Fig. 4c) and 23% in the case of shear with bending (Fig. 4e).Similarly, the strain analysis revealed higher strain values in the isotropic model.Strain values were higher, up to 35% for the shear load (Fig. 4b), 31% for compression (Fig. 4d) and 35% in the case of shear with bending (Fig. 4f).

Discussion
Dentin is a hard tissue that constitutes the bulk of the tooth structure and is composed of approximately 70% inorganic material (hydroxyapatite crystals), 20% organic material (mainly collagen type I), and 10% water by volume.It contains a large number of tubules whose density and diameter decrease from the pulp to the dentinoenamel junction.The dentinal tubules contain the odontoblastic processes, dentinal fluid and nerve fibres, and extend in an S-shaped fashion into the coronal dentin, following a straight course into the radicular dentin.Their wall is formed by very dense, uniformly mineralized peritubular dentin, and they are separated from each other by less densely mineralized intertubular dentin. 22he properties of dentin are described as the response of the tooth to the applied loads, and enable the prediction of tooth strain and fracture resistance.The commonly measured effective properties are Young's modulus, tensile and comprehensive strength and fracture toughness, reflecting the complex interactions of the dentin constituents and the microstructure.The properties depend on tubule density and orientation as well as the local density of the mineral phase caused by the remodeling process.The elastic properties of dentin influencing the tooth strength are usually described as stiffness including Young's modulus, shear modulus and Poisson's ratio. 14he dentin microstructure, in combination with the whole tooth geometry, causes considerable variations in the results of experimental measurements of the dentin properties. 11,14Dental models, like other models from the fields of biology or medicine, have a complex geometry.However, their geometry is relatively easier nowadays thanks to 3D scanners (micro-CT and CT), which offer an easy and fast way to build virtual geometry in any computer-aided engineering (CAE). 6,7On the other hand, the physical description of the material and boundary conditions are often significantly simplified.Some typical FE model inputs may be simplified and/or reduced to the essential parts.However, it is relatively difficult to simplify the material properties without changing the natural behavior of the object.
The properties of dentin are usually derived from experimental measurements, which have their own limitations and yield a wide range of results.The complex in-ner structure of the dentin produces anisotropy, which should be taken into account in the studied samples.All difficulties occur especially when the samples are limited because of the shape and size (such as human teeth).Conducting experiments, therefore, involves the use of more sophisticated methods in order to deliver more reliable results.The load and boundary conditions are complex, as are the contact surfaces of the teeth.Therefore, the loads may be represented in non-uniform pressure distribution applied to the suitable sites (physiologically loaded or reflecting in vitro experiments, depending on the modeling type).During FE simulations the most important parameters to set are the mechanical properties.Both quantitative values and qualitative descriptions (e.g.isotropy, anisotropy, axes orientation) are needed.In most of the FEM studies evaluating dental restoration, the dentin is modeled as an isotropic material.The finite element analysis of the dentin model carried out in our study showed higher strain values in the isotropic model.The significant difference observed between the isotropic and anisotropic dentin model revealed the necessity of an estimation of the dentin anisotropy impact on the global FEM model while the tooth restoration analysis.
The homogenization methods are mainly based on the more or less accurate modeling of the non-homogenous composite microstructure and therefore frequently referred to as micromechanical models.Determination of the material parameters (e.g.Young modulus, Poisson's ratio, thermal expansion coefficient) is performed mathematically by averaging the strain tensor of material properties for individual phases over selected volume (RVE).
Micromechanical models can be divided, with regard to the applied methods, into 3 main groups: empirical, analytical and numerical.The analytical models are particularly noteworthy because of the possibility of easy calculations and a small amount of calibration experiments in the empiric models required.Most of the analytical methods of homogenization are based on the equivalent Eigenstrain method, described mathematically as a singular elliptic inclusion in the elastic infinitely large matrix.Similarly, the numerical homogenization determines the effective elastic parameters for non-homogeneous material based on the properties of material phases and their geometry or internal structure.The numerical homogenization also uses the concept of the RVE and assumes that the material is statistically homogeneous at the macro-size level.Numerical micromechanical models are based mainly on the Finite Element Method (FEM).Furthermore, numerical methods of homogenization allow easy inclusion of various geometric details with complicated shapes.Unfortunately, FE modeling of RVE requires significant effort for modeling and CPU computing time.
An increase of Young's modulus value (maximum on the top layer) alongside a decrease in the volumetric fraction and diameter of tubules was observed.The increase in E1 of 18% and E2 and E3 of 47% on the top layer was also observed.A decrease of the volumetric fraction and diameter of tubules also induced an increase of the shear modulus (G) values.On the top layer of the dentin, the values of G12 and G13 were increased by 34% and G23 by 48% (Table 3).
Based on our data we could suggest that the increase of the tubule diameter (and volumetric fraction of tubules) led to the increase in anisotropy and the decrease in the E-modulus value.
Anisotropy of the dentin varied from min.6.9% in the top layer to max 35.2% at the bottom layer.The mean values of anisotropy obtained analytically and numerically were comparable but they were higher than relatively low anisotropy (~10%) suggested by Kinney et al. 23 Our results underlined the high impact of the dentinal tubules on the elastic properties and anisotropy of the dentin.This was in opposition to Kinney et al., who stated that the intertubular dentin matrix governed the elastic behavior of dentin, and that the tubules did not introduce elastic anisotropy. 21he earlier published dentin tensile strength results confirmed the anisotropic behavior of the dentin and the influence of the tubules direction on the strength. 19,24he high anisotropy value may be also explained by the simplification of geometry when dentinal tubules are modeled as straight tubes.However, the influences of the S-shaped dentinal tubules on the results have not been evaluated in this paper.It should also be noted that the tubules' natural orientation could vary widely across a specimen, so potentially an anisotropy may be undetected in a mechanical test of larger specimens. 25Miura et al., in order to create multi-scale tooth model, measured E and G modules using microindentation. 26They obtained high anisotropy values amounting to E 17.07 GPa, G 1.7 GPa longitudinally and E 5.61 GPa and G 6.0 GPa transversally to the tubules.Unfortunately, the information about measuring locations in the tooth structure was not available.
The dentin anisotropy is also detectable during flexural strength tests; Shinno et al. registered the flexural strength of the anti-plane longitudinal vs transverse group that amounted to be 210 ± 67 MPa and 114 ± 36 MPa, respectively. 27s mentioned above, the decrease in the volumetric fraction and diameter of the tubules induces the increase of the shear modulus (G) and Young's modulus values.Considering the data obtained with the use of the analytical homogenization procedure, we could state that within the same layer of dentin an increase of the tubules' diameter (i.e. increase of volumetric fraction of tubules) resulted in the increase of anisotropy and the decrease of the E-modulus value.Those results, combined with the large deviation of reported mechanical tests results, suggested the necessity for an analysis of the relation be-tween the tubules' diameter (and volumetric fraction of tubules) and mechanical test data.
Shinno et al. reported that a high density of dentinal tubules also decreased the flexural strength of the dentin. 27t the same time, the ageing of the dentin, with the obturation of the tubules, resulting from higher mineralization, could diminish the flexural strength.Panfilov et al. showed higher tensile strength of mature rather than teenage dentin (549.5 ± 26.9 MPa vs 445.2 ± 30.9 MPa, respectively). 28ithin the limitations of this study we can suggest the transversely isotropic behavior of dentin and the anisotropy ranged from 6.9 to 35.2%, depending on the dentin location.The numerical and analytical dentin homogenization presented similar results.The decrease in volume and diameter of tubules resulted in some increase of the shear modulus (G) and Young's modulus (E) values.Both methods give the same level of anisotropy, with relatively small variation.We could assume that the analytical determination of the dentin properties is simpler, more effective and faster, and therefore allows for a recalculation of numerous geometric scenarios.The numerical homogenization requires a new FE model for each new geometric concept leading to a significant increase in the overall cost.However, the analytical model becomes much more complicated (or even unfeasible) for non-parallel or S-shaped tubules, when the numerical homogenization is still applied.
We have attempted to create an anisotropic dentin model imaging its natural nonlinear characteristics and analyze the impact of physiological anisotropy on the dentin strain.The comparison of the isotropic and anisotropic dentin model revealed significant differences in strain, especially in the area with a high density of dentinal tubules.The finite element analysis showed up to over 30% higher strain values in the isotropic model compared to the anisotropic one.It is worth mentioning that Kishen and Vedantam, as in our study, found that the anisotropy of the dentin and changes in the modules E and G (gradients of moduli) during FEM analysis caused significant decrease in displacement and strain values. 8ishen and Vedantam also emphasized the importance of hydrostatic pressure of the dentinal tubule fluid, which was neglected during our study. 8The significant difference observed between the isotropic and anisotropic dentin model pointed to the necessity for an estimation of the impact that the dentin anisotropy in the global model has on the FEM analysis results of the tooth restoration.
The anisotropy of the dentin and changes in modules E and G should be taken into account in future FEM analysis studies.Hence, further study should be directed towards an analysis of the global model of the tooth (whole tooth structure), while applying anisotropy of the dentin and the moduli gradients.This procedure, after model calibration, may be useful in the analysis of the restored teeth and during optimization of restorative process.

Fig. 1 .
Fig. 1.Modeling of dentin structure -the method chain from a single RVE cell (dimensions example for tubule diameter 2.4 μm, axes are visible) (A); dentin layer (B); equivalent orthotropic & homogenous material model (C); global FE model of the tooth (D)

Fig. 2 .Fig. 3 .Fig. 4 .
Fig. 2. FE mesh representing the cube model of the dentin; 3 layers of the dentin; axes of the model and loading directions are visible

Table 1 .
Dentinal tubules diameter and density values gleaned from a literature survey and parameters used in the presented study

Table 2 .
Material properties of dentin tissue used during FEM analysis

Table 3 .
Effective material properties of dentin structure obtained in analytical and numerical homogenization procedure