Hierarchical modelling of in situ elastic deformation of human enamel based on photoelastic and diffraction analysis of stresses and strains

Human enamel is a typical hierarchical mineralized tissue with a two-level composite structure. To date, few studies have focused on how the mechanical behaviour of this tissue is affected by both the rod orientation at the microscale and the preferred orientation of mineral crystallites at the nanoscale. In this study, wide-angle X-ray scattering was used to determine the internal lattice strain response of human enamel samples (with differing rod directions) as a function of in situ uniaxial compressive loading. Quantitative stress distribution evaluation in the birefringent mounting epoxy was performed in parallel using photoelastic techniques. The resulting experimental data was analysed using an advanced multiscale Eshelby inclusion model that takes into account the two-level hierarchical structure of human enamel, and reﬂects the differing rod directions and orientation distributions of hydroxyapatite crystals. The achieved satisfactory agreement between the model and the experimental data, in terms of the values of multidirectional strain components under the action of differently orientated loads, suggests that the multiscale approach captures reasonably successfully the structure–property relationship between the hierarchical architecture of human enamel and its response to the applied forces. This novel and systematic approach can be used to improve the interpretation of the mechanical properties of enamel, as well as of the textured hierarchical biomaterials

Hierarchical modelling of in situ elastic deformation of human enamel based on photoelastic and diffraction analysis of stresses and strains 1

. Introduction
The increasing need to understand and predict the effect of structural alterations on the performance of dental tissues and their artificial replacements arises both in the context of clinical treatment, and the development of novel dental prosthetic materials [1]. The elucidation of the dependence of the mechanical behaviour of human enamel on its complex hierarchical structure remains a challenging task. Two different length scales dominate the structure and the mechanical behaviour of the enamel: At the microscale, $5 lm diameter keyhole-like cross-section aligned prisms (or rods) are oriented towards the crown of the tooth [2]. At the nanoscale level, the rod is known to be a composite made from needle-like biological hydroxyapatite (HAp) crystals ($25-30 nm thick and of length known to be more than 1000 nm, or even spanning the entire thickness of the enamel layer) [3,4], which are held together by the protein matrix between the rods [5][6][7]. The orientation of HAp needle-like crystals within the rod is known to be a gradual variation on the length scale of the rod diameter [8,9].
The focus of most research to date has been on the mechanical properties of the enamel at the macroscale, with microstructural effects rarely taken into account [2,10]. Very few studies have focused on the multiscale analyses required to determine the influence of the nanoscale structure on the macroscopic mechanical response [11]. In fact, both the crystal shape and orientation of the mineral phase nanocrystals have previously been shown to have an effect on the anisotropy of overall stiffness and strength [12]. Therefore, in order to establish a firm understanding of this hierarchical structure-property relationship, further application of advanced nanoscale techniques and the formulation and refinement of systematic models are required.
Synchrotron-based wide-angle X-ray scattering (WAXS) is a non-destructive analytical technique used to quantify the internal strain of atomic lattices on (poly)crystals, both residually stressed and subjected to external loading in situ [13,14]. In addition, WAXS techniques are able to reveal quantitative information about the orientation distribution of crystals (texture) [15]. For example, recent applications of WAXS include the determination of the mechanical behaviour of mineralized biological composites such as bone and bovine teeth [14,[16][17][18] while simultaneously providing insight into the crystallographic parameters and textures of such materials [8]. The strain distribution across the amelo-dentinal junction (ADJ) in bovine teeth was investigated by Almer and Stock [16]. However, this was limited to the strain in the loading direction [16,19]. In practice, the external mastication load is not perfectly aligned with the longitudinal direction due to the complex tooth geometries (e.g. on transverse ridges and cusps in first molars), and also due to the local orientation variation of HAp crystals within the rods. Therefore, the analysis of the relationship between the orientation of the applied load and strain needs to be sought. Moreover, the studies reported earlier did not take into account the nanoparticle distribution effect on the mechanical response. In addition, very few studies devoted to human enamel have ever been published [20]. It is therefore unsurprising that there is a lack of understanding of the effects of different growth histories, species or race on the mineralized tissue morphology, and the corresponding mechanical properties [21].
In order to carry out the in situ mechanical loading experiments in a versatile manner, i.e. allowing different directions of loading with respect to the preferred directions of structural orientation within the tissue, the sample was embedded in a photoelastic epoxy disk. The enamel cubes studied in our preliminary experiments could not withstand high external load applied by direct compression between the platens. Without the protection of the epoxy disk, samples were likely to develop local stress concentrations and/or microcracks that reduce the accuracy of measuring the mechanical behaviour of enamel. Furthermore, due to its birefringent properties, epoxy offers the possibility of deducing information about the stress distribution around the sample using photoelastic techniques. Photoelasticity is a non-destructive, whole-field method widely applied in stress distribution analysis. The fringe pattern arises when the sample is viewed between crossed polarizing plates, with the colour or brightness of the fringe being related to the difference between the principal stresses, which is in turn proportional to the maximum shear stress, or Tresca stress, in the material [22]. We combine the photoelastic analysis technique with WAXS analysis during in situ loading in order to establish the relationship between the external stress distribution and the internal crystal lattice strain.
A number of different models of composite deformation have been previously used to describe the elastic response of mineralized biological tissues that arises through the interaction between different constituent phases [11]. This approach also allowed the unknown properties of the component phases to be determined [23]. However, these models mainly focused on the analysis of deformation in only one direction (loading direction) and therefore were not able to provide adequate consideration of the elastic anisotropy. Recently, a multiscale Eshelby inclusion model has been established and successfully applied for the evaluation of the mechanical response of human dentine [12] and of the enamel subjected to compression along the longitudinal direction of the rods [19] by capturing the relationship between the nanoscale structure and the macroscopic loading. The multiscale modelling approach was shown to capture the micromechanical response reasonably well using the two-level hierarchical description of the structure of dentine and enamel, with each level consisting of an isotropic matrix and a group of anisotropic inclusions.
In this study, in situ photoelasticity and synchrotron WAXS techniques were applied simultaneously to measure the applied stress, the internal HAp crystalline orientation distribution and the elastic lattice strain for three samples of human enamel. The samples were prepared so that the primary rod direction was at a different orientation with respect to the external compressive loading that was applied diametrically to the epoxy disk containing the tissue sample (see Fig. 1). The multiscale Eshelby inclusion model was then applied to the analysis of the results, and the capability of the model to capture the relationship between the nanoscale structure and macroscopic loading was investigated.

Sample preparation
Freshly extracted human third molars with no apparent damage, caries or other dental treatments were used for this study (ethical approval obtained from the National Research Ethics Committee; NHS-REC reference 09.H0405. 33/Consortium R&D No. 1465). An enamel disk 2 mm thick was cut from the each tooth using a low-speed diamond saw (Isomet Buehler Ltd., Lake Bluff, IL, USA) and further prepared into smaller bars. A series of polishing papers were used to refine the final 2 mm Â 2 mm Â 2 mm cube of enamel. The cubes were placed in the centre of a 12 mm diameter cylindrical mould and embedded in epoxy resin (Buehler Epokwick, ITW Test & Measurement GmbH, Dusseldorf, Germany). The disks' surfaces were subsequently polished to expose the enamel surfaces.
In total, the three cubic enamel samples were prepared (designated #6, #7 and #3) with different rod directions with respect to the loading direction (x-direction in Fig. 1a). The predominant direction of rods in sample #6 was parallel to the loading direction with rods lying in the x-y plane, in sample #7 it was perpendicular to the loading direction with rods lying in the x-y plane, and in sample #3 it was perpendicular to the loading direction with rods lying in the y-z plane.

2.2.
In situ X-ray diffraction measurements

Photoelasticity setup
A Sharples S-12 demonstration polariscope was used to collect the in situ photoelastic images. The setup consisted of light source, polarizers, quarter-wave plates and a digital SLR camera as shown in Fig. 1. The quarter-wave plates remained crossed and polarizers were aligned crossed to establish the dark field. A green-light filter was also placed between the light source and the camera lens in order to obtain monochromatic fringes to simplify the analysis. A solid epoxy disk (without the sample in the centre) produced from the same batch of epoxy resin was used as a common calibration specimen.

Mechanical loading setup
A schematic diagram of the experimental setup is shown in Fig. 1a. The epoxy disk which contained the cubic sample of human enamel was slowly deformed along the x-direction in laboratory coordinates (Fig. 1a). Compressive loading was applied along the x-direction at the load levels from 0 to 400 N using a remotely operated and monitored compression rig (Deben, Suffolk, UK) with a 5 kN calibrated load cell. The rig was equipped with custommade jaws, allowing a high-energy transmission X-ray setup to be used. The load was incrementally increased (in 25 N steps and a loading rate of 3.8 N s À1 ) and held constant while the WAXS and photoelastic patterns were collected.

Beamline diffraction setup
The experiment was performed on the B16 test beamline at Diamond Light Source (DLS, Oxford, UK). A monochromatic X-ray beam of 20 keV photon energy (wavelength: k = 0.062 nm) was collimated by slits to a spot size of 0.5 mm Â 0.5 mm. Radiographic images of the samples were initially used to align the samples and determine the position of interest. The incident beam on the sample was perpendicular to the loading direction. Space restrictions of the beamline meant that the sample had to be translated laterally from the WAXS configuration into the photoelastic setup at each consecutive loading increment. WAXS diffraction patterns were recorded using a Photonic Science Image Star 9000 detector (Photonic Science Ltd., UK) which was placed 177.33 mm downstream of the sample. A lightly compacted disk of NIST standard silicon powder was used for precise calibration of the sample-to-detector distance using diffraction pattern analysis [24].

Data analysis
2.3.1. Photoelasticity data analysis [22] The use of a diametrically loaded disk is a standard calibration technique which can be used either to obtain the photoelastic properties of the birefringent material, or to determine the calibration constants if these properties are known. The fringe constant of the material is defined as: where D is the diameter of the epoxy disk (12 mm), N 1 is the number of fringes for the calibration sample and P is the applied load (N). In order to obtain an estimate for f, a series of known loads was applied to an empty epoxy disk and the corresponding fringe numbers recorded. An average value for f was then determined by taking the mean of the f term calculated at each load. The stress distribution in the disk is characterized in terms of the difference between the principal stresses (or Tresca shear stress) using the following expression: where h is the thickness of epoxy (2 mm), and N 2 is the number of fringes for the sample of interest. The deformation state of the plate-shaped epoxy corresponds closely to the plane stress state approximation in two-dimensional photoelasticity. Since the loading was applied along the x-axis, with the epoxy disk lying in the x-y plane, the principal stresses were r 1 < 0, r 2 > 0 and r 3 = 0.

WAXS data analysis
WAXS analysis relies on interpreting the shift of the diffraction peaks obtained from a bundle of HAp crystals, so that the average lattice strain in the crystals can be deduced [13,25]. Typical WAXS patterns of HAp are shown in Fig. 2a-c and only the peaks of interest were selected for interpretation (#6 (0 0 2) peak, #7 (0 0 2) peak and #3 (2 1 0) peak). For samples #6 and #7, the (0 0 2) lattice plane reflection ring from HAp contains information on the orientation of the c-axis of the crystals as well as the fibril orientation. This is due to the fact that the enamel structure has a characteristic strong and distinct parallel orientation of needle-shaped crystals [14]. For sample #3 the c-axes of the crystals lie perpendicular to the measurement plane, so that the (0 0 2) peak is absent or weak. Therefore, another peak (2 1 0) needed to be selected. At each compressive loading step, the X-ray eye detector was used to verify the central position. The Deben compression rig applied the load by displacing the two platens in a balanced way, so as to preserve the centre position approximately stationary. This ensured that a consistent sample position in the centre of the sample was interrogated throughout the WAXS data collection. The elastic lattice The monochromatic X-ray beam was collimated by slits and oriented perpendicular to the sample surface and the loading direction. An X-ray eye detector was used to ensure that the beam was illuminating the cetral position of the enamel specimen. At each detection point, following centring, the X-ray eye detector was translated laterally out of the beam to expose the WAXS detector. WAXS diffraction patterns were recorded at each loading step in three locations on the sample. After each WAXS pattern collection, the compression stage was laterally translated to the photoelastic setup to collect the photoelastic patterns. strain of HAp crystals was calculated from the changes in the interplanar spacing [26]: where d hkl is the deformed d-spacing of the lattice planes with the Miller index {h k l}, and d 0 hkl is the strain-free d-spacing value for the same set of planes.
The data analysis procedure involved pre-processing the 2-D diffraction using Fit2D [27]. Initially the peak of interest on each pattern was ''caked'' (a term used to refer to the selection of a sector in the radial-azimuthal coordinates of each pattern, and binning the data to obtain the equivalent 1-D radial pattern) within a range of 20°. The azimuthal angle is defined as the angle with respect to the x-axis in the x-y plane (perpendicular to the incident beam). For each specimen, a different azimuthal range was selected due to the preferred orientation of crystallites in the sample (for sample #6 between 0°and 45°, for sample #7 between 45°and 90°, and for sample #3 the available range was between 0°and 90°). This is indicated by the red crosses in Fig. 2a-c. The lattice spacing value deduced for the centre direction of each caked sector gives rise to the normal strain value for the corresponding orientation (see Fig. 2a-c). To determine the lattice spacing, the experimental diffraction pattern was converted to integrated 1-D plots of intensity vs. scattering angle. These plots for the three samples with peaks of interest are shown in Fig. 2d-f. The peaks of interest were fitted with a Gaussian function in order to obtain the peak centre. With increasing load, the peak centre position moved with respect to the strain-free reference value. Finally, the calculation of the HAp elastic lattice strain was performed using Eq. (3).
The orientation of the HAp crystals was also carried out using a similar ''caking'' approach. In this case, the 1-D intensity was plotted as a function of azimuthal angle I WAXS (u) $ u WAXS (Fig. 2m and n). Gaussian fitting of the azimuthal centre position of these peaks of plots can then be used to define an associated orientation [8].

Finite-element analysis
The finite-element package ABAQUS Ò v. 6.9 was used to simulate the stress distributions and to obtain the associated photoelastic patterns. Direct comparison of photoelastic patterns between model and experiment provides a way of simultaneously verifying the correctness of the macroscopic properties (of epoxy and enamel) as well as the mechanical model of the stress distribution within the epoxy disk around the enamel samples. In order to visualize the photoelastic fringe patterns the ''Tresca equivalent stress'' r Tresca was selected as the output parameter. Initially a simulation of the solid photoelastic epoxy disk (the calibration specimen) was developed using a model of 12 mm diameter and 2 mm thickness. The elastic constants of the epoxy were identified by comparing r Tresca contour plot with the fringe patterns observed at different loads.
In order to model the enamel sample (2 mm Â 2 mm Â 2 mm), a cubic inclusion was introduced and embedded in the centre of the epoxy disk with the isotropic elastic constants obtained from the calibration. The experimentally applied loads were introduced on the edge of the epoxy model and the stress distribution around the sample (in the x-y plane in global coordinates) was recorded. The tractions present at the boundary of the enamel inclusion obtained from this model were then used for the externally applied stress values in the Eshelby model.

Geometrical assumptions
The geometrical assumptions used in the enamel model have been justified in our previous studies. An outline of these assumptions at this stage is given here [19]. At the first structural level, the geometric model of human enamel involves aligned rods with a keyhole-shaped cross-section. The second-level structural consists of the bundles of HAp needle-shaped crystallites found within each rod. These bundles are roughly aligned with the longitudinal direction of the rods with some minor misorientations. The degree of misorientation can be determined from careful WAXS data interpretation. This analysis assumes that the peak intensity at a given azimuthal angle is proportional to the number of HAp crystals of this orientation [28].
The orientation distributions of the rods in the enamel samples were derived from the analysis of experimental data and are illustrated in Fig. 2g-i. Furthermore, the predicted distributions of needle-shaped HAp crystals (in the measured cross-section, perpendicular to the beam direction) are shown in Fig. 2j-l in which the misorientation distribution is based on the intensity plots in Fig. 2m and n. The principal orientation of sample #3 is aligned so that orientation component is observed perpendicular to the beam, and therefore no equivalent plot is produced.
In the multiscale Eshelby model, both structural levels are considered as non-dilute systems consisting of a number of inhomogeneous inclusions (rods at the first level and HAp crystals at the second level). For simplicity, both rods and HAp crystals are assumed to be of cylindrical shape. Therefore, the classical solution for inclusion in the form of a cylinder is used to simulate both the needle-shaped rods and the individual HAp crystals. In the Eshelby approach, the cylinder is approximated by a prolate spheroid described by the three dimensions, a 1 , a 2 and a 3 , where a 1 = a 2 ( a 3 , i.e. the cross-section of the ellipsoid perpendicular to its longest axis is a circle [29]. The crystal c-axis (corresponding to the (0 0 2) peak) of the needle-shaped HAp crystals is normal to this cross-section of the cylinder, i.e. is aligned with the a 3 direction [28]. In the next section a detailed derivation and discussion of the model formulation and implementation is introduced.

Multiscale Eshelby model
Based on the hierarchical structure of enamel, a two-level Eshelby inclusion model can be used to demonstrate how the HAp crystals respond to external macroscopic loading. The original Eshelby solution gives the analytical expressions for the elastic fields around an ellipsoidal inclusion embedded in an infinite isotropic matrix that has the same elastic properties as the inclusion, but there is a misfit or mismatch strain (eigenstrain) between the inclusion and the matrix [30].
Eshelby introduced a tensor S that only depends on the inclusion shape and the Poisson's ratio of the matrix. For an inhomogeneity with a different stiffness from the matrix, the elastic field can be found using an equivalent inclusion method. This formulation regards the inhomogeneity as equivalent to an inclusion with an appropriate virtual misfit strain which needs to be determined. The function of the virtual misfit strain is to satisfy the equivalence where the elastic field of the inhomogeneity is identical to that of the equivalent inclusion. In this way the classical Eshelby theory can be used to determine the elastic field of any kind of inhomogeneity. This formulation has greater generality and includes the case of rods or HAp needle-like crystals (long cylinders) embedded in the isotropic protein matrix. In the next section, the multiscale model for enamel is briefly introduced, while a detailed derivation of the basic theory is given in the Appendix.

First-level model: multiple aligned rod inclusions within the enamel sample
The purpose of the first-level model for human enamel, which describes multiple aligned rod inclusions within enamel, is to establish the elastic relationship between the externally applied stress r A and the stress in the rod-like inclusions r inclusion = r rod .
According to the Eshelby model derivation [30], the stress in the rods at the first structural level of consideration can be written as: or, expressed more simply: where f 1 is the volume fraction of rods, S 1 is the Eshelby tensor for the long cylinder simulating the rod shape, C matrix and C rod are the stiffnesses of the protein matrix and the rod, and T is the tensor transformation (rotation) matrix. T encapsulates the effect of the three Euler angles that define the orientation of the rods with respect to the fixed laboratory coordinate system. The stiffness of the rod in the simulation is initially unknown and is determined from the second-level model.

Second-level model: partially aligned HAp inclusions within each individual rod
In the second-level model, the relationship between the average strain in the partially aligned HAp crystals (averaged within the considered gauge volume) and the external stress (the rod stress from the first-level model) can be established as follows: matrix r inclusion or, expressed more simply: where hCi HAp and hS 2 i are the average elastic stiffness tensor and the Eshelby tensor for HAp crystals within the gauge volume, respectively, and f 2 is the volume fraction of the HAp crystals with a certain orientation T HAp within the rods. The stiffness of a single HAp has been reported to be transversely isotropic [31]. In Eq. (5), different orientation angles (different T HAp ) would lead to different values of hei HAp aligned .
In order to evaluate the multidirectional normal strain components, each azimuthal angle was taken to correspond to a group of perfectly aligned HAp crystals, as illustrated by Eq. (5), with a certain volume fraction f 2 and orientation T HAp , where f 2 is determined from the WAXS intensity plot (Fig. 2m and n). In other words, the bundled HAp crystals consist of a distribution of multiple groups with different orientations and different volume fractions. Meanwhile, for each group, all the crystals are perfectly aligned, thus hCi HAp and hS 2 i can be represented by the single-crystal values hCi HAp = hCi HAp , hS 2 i = S 2 . Note that a single HAp crystal in the second level can be regarded as a single inhomogeneity with needle-like shape, thus the Eshelby tensor for a single HAp should be the same as the rod inclusion (S 2 = S 1 ).
Further, the stiffness of rod is obtained from the following equation: HAp is the transformed Eshelby tensor and C HAp ¼ T À1 HAp C HAp T ÀT HAp is the transformed stiffness of aligned transversely isotropic HAp crystals. Finally, the relationship between the HAp strain and the overall external loading can be established by the combination of Eqs. 4 and 5 and simplified Eq. (6) to give:

Finite-element modelling of photoelastic patterns
The photoelastic fringe patterns of solid epoxy disks and epoxy with enamel embedded at the disk centre are shown in Fig. 3. The results for Tresca stress distribution simulated in ABAQUS are shown in Fig. 4. The stress variation along the edge of the samples simulated in ABAQUS was not constant, but was concentrated especially at the two corners. Direct comparison and matching of the model to the photoelastic patterns for solid epoxy disks lead to the isotropic elastic constants of epoxy disk being estimated to be E epoxy = 1.2 ± 0.015 GPa, t epoxy = 0.35 [32]. These parameters were used for the samples with embedded enamel to calculate the predicted normal stress distribution around the enamel inclusion illustrated in Fig. 5a that corresponds to the maximum external load of 400 N. At this applied load, the stress values were Fig. 3. Photoelastic images of epoxy disks without (first row) and with embedded enamel (second row) under a series of loading conditions. The number of fringes can be seen to increase with load. The patterns on the first row of solid epoxy disks were used as calibration in order to determine epoxy properties. The patterns on the second row of epoxy with enamel sample were further used for the determination of the stress distribution around the cubic enamel inclusion. The fringe number is marked. averaged over all the boundary nodes of each enamel cube sample edge. This uniform traction was applied in order to eliminate the effect of the singular points or stress oscillations at corners, with the result of À54.1 MPa in the x-direction and 0.27 MPa in the ydirection. Therefore the sample loading condition was found to be close to uniaxial. As expected, the average stress value applied to the sample was also found to be linearly proportional to the value of the externally applied load, as shown in Fig. 5b.  . The textured nature of enamel meant that only limited azimuthal ranges within the rings could be captured in samples #6 and #7. In contrast, in sample #3 the crystals were approximately randomly distributed in the x-y plane, so that complete rings were observed. Fig. 2m and n show the azimuthal intensity variation within the available azimuthal ranges in the measurement plane (x-y plane) in samples #6 and #7, respectively. This provides information about the crystal orienta-tion distribution, as described in Section 3.3. The preferred orientation direction of the basal plane in HAp crystals in samples #6 and #7 is approximately perpendicular to the central bisector of the arc of the peak of interest (see Fig. 2b and c). The precise values are determined from the Gaussian fit in Fig. 2m and n (red line). Since the deformation considered is elastic, the assumption was made that the crystal orientation remains unchanged, i.e. the peak centres found from the plots shown in Fig. 2m and n remain the same at all applied loads.

Nanoscale HAp distribution and mechanical response of enamel
In order to obtain the variation of the average normal elastic lattice strain component of HAp crystals with applied load, the peak shifts at different azimuthal angles were obtained from the WAXS data through caking each pattern into 20°sectors. The available range for each sample was different, as shown by the red marks in Fig. 2a-c: in sample #6 between 0°and 45°, in sample #7 between 45°and 90°, and in sample #3 between 0°and 90°. For each of the samples, the elastic lattice strain as a function of applied stress is plotted as the argument against the applied stress for different azimuthal angles in Fig. 6 (in which 0°represents the loading direction). The experimental results indicate that, as expected, there is an approximately linear relationship in each direction and Fig. 4. ABAQUS simulation images of Tresca stress distribution of the epoxy disks without (first row) and with embedded enamel (second row) under a series of loading conditions. The isotropic elastic constants of the epoxy together with the normal stress distribution around the enamel sample were determined by matching the ABAQUS images with the photoelastic images in Fig. 3. The greyed regions represent the high stress concentration areas and were neglected in the images. that each azimuthal angle corresponds to a different slope. This ratio between the applied uniaxial external stress and the average HAp lattice strain gives the apparent modulus along different directions [33]. At this stage it is also relevant to mention that the initial residual strain in the samples was found to be negligible.

Evaluation of the multiscale Eshelby model
In setting up the multiscale Eshelby model, the average external stress derived from the ABAQUS model was used as the input external stress. The material properties and other parameters were taken from the literature [11,31,34,35], and refined by fitting with the experimental data of HAp crystal lattice strain. The initial volume fraction of HAp crystals was taken from the previous analysis of the self-similar hierarchical two-level model of enamel, i.e. approximately 95% at both the microscale and nanoscale levels [11]. The needle-like HAp crystals in the enamel were also assumed to have a transversely isotropic stiffness with five independent elastic constants [31]. To describe the geometry of rod and HAp crystals, the cylinder Eshelby tensor was used and the length of elliptical axes a 1 and a 2 within the transverse cross-section were expected to be equal to each other, i.e. a 1 /a 2 = 1. In addition, Fig. 6. The experimental results of applied compressive stress vs. elastic lattice strain of HAp crystals for the three samples at different azimuthal angles (every 15°) with respect to the loading direction. (a) Sample #6; due to symmetry, only the available angles from the (0 0 2) ring from 0°to 45°are selected. (b) Sample #7; due to symmetry, only the available angles from the (0 0 2) ring from 45°to 90°are selected. (c) Sample #3; due to the symmetry, only the available angles from the (210) ring from 0°to 90°are selected. Table 1 Preferred orientation and refined structural parameters of the enamel used in the two-level Eshelby model (values reported in the literature are also presented).
Preferred orientation 8°91°- C HAp E xx = 132.1 ± 9 GPa, G xy = G xz = 39. Young's modulus of protein was set to 1 GPa (this does not take into account any effects from viscoelasticity or viscoplasticity [34,35]). The final refined parameters and reference values are summarized in Table 1.
In Fig. 7, the comparison between the model prediction and experimental data of the normal strain variation in the three samples under the maximum external load of 400 N is shown as a plot against the azimuthal angle.

Discussion
The hierarchical nanostructure of human enamel may vary significantly between patients due to personal history, diet, ethnic origin, etc. However, due to ethical legislation in the UK, it is not possible to obtain or publish any further information on the medical history of the patients.
This study was conducted using the penetrating power of synchrotron X-rays to provide a bulk probe for structure and strain analysis. Unlike the vast majority of studies that rely on the surface characterization (scanning electron microscope, atom force microscopy, nanoindentation, Raman, etc., and also the reflection mode of X-ray diffraction), this WAXS study required less preparation effort in terms of cutting and storage and surface conditioning. The use of much thicker samples and gauge volumes ($2 mm) than those for other methods ($0.05 mm) ensured that the bulk response investigated was not affected by minor surface changes. The combination of the data obtained using penetrating radiation (synchrotron X-ray) with model refinement offers the possibility of identifying the nanoscale parameters of bulk enamel. The parameter refinement and validation strategy employed in the model adopted in the present study is particularly helpful in the identification of nanoscale parameters that may be hard to determine in the other experiments. The relationship between the nanoscale structure and the macroscopic mechanical behaviour established in our study improves the understanding of the nanoparticle distribution effects that was lacking in the earlier studies [8,16,20].

Finite-element simulation
The high level of correlation seen between the finite-element analysis result (illustrated in Fig. 4) and the photoelastic patterns (shown in Fig. 3) suggests that the stress values extracted from the finite-element model are a good representation of the tractions applied to the enamel samples. In Fig. 4 the stress concentration region has been shaded grey. As this region is far enough away from the enamel inclusion, it can be neglected. Only the patterns and stress distribution close to the enamel need to be considered. The plot of the normal stress distribution in Fig. 5 is likely to originate from the frictional contact condition established between the enamel sample and epoxy in ABAQUS. However, the use of the averaged value of the stress distribution will minimize the effect of this oscillation. In short, the finite-element results suggest an accurate estimation of the input stress for the two-level Eshelby model.

Refined parameters of the two-level Eshelby model
The input from reference data on material parameters (HAp stiffness, etc.) are only used as the starting guess for the optimization of the elastic constants in order to obtain the best fit to our experimental data. The final parameters used in the multiscale Eshelby model are outlined in Table 1 and were refined iteratively, starting from the values reported in the literature. Of the refinable parameters, the key ones were the volume fraction of inhomogeneities at each level and the elastic constants of a single HAp crystal. Among these parameters, the elastic constants, and the modulus in particular, exert the most profound influence on the result. The three samples studied were cut from the same tooth. Therefore their properties, such as the degree of mineralisation, size of particles, etc., can be assumed to be close and consistent, with the exception of the orientation distribution of nano-HAp crystals which is known to vary with location within the tooth. Model fitting was carried out so as to find the parameters that give the best agreement for all three samples. The single-crystal elastic constants of HAp were assumed to be the same. In the model matching (optimization) process, within each specimen, the volume fraction of the rods in enamel and the HAp crystals in rod were assumed to be the same. After the refinement of the elastic constants (within the range reported in the literature) the volume fraction was further refined. At this stage, further adjustments of other parameters were also attempted and it was found that this only had a minor effect. In order to estimate the error in our evaluation of parameters, we carried out a careful study of error propagation in our model-fitting procedure following the error analysis methods in Ref. [36]. The 95% confidence interval in the fitting results was propagated back into elastic modulus uncertainty. The final result with description has been added to Table 1 Note that since the results of our test are not sensitive to the shear modulus, the uncertainty of the shear modulus was not assessed.
The authors note that environmental and developmental factors may influence the mechanical behaviour of enamel [37]. Following tooth eruption, the interaction of the enamel surface with ions in saliva leads to enamel maturation. This post-eruptive maturation may affect the mechanical behaviour of enamel surface. For example, the incorporation of fluoride into enamel and the formation of fluoroapatite may exert an influence on the superficial layer of dental enamel. However, the transmission X-ray setup we employed ensured that the bulk rather than the surface was characterized for the samples obtained from freshly extracted teeth with no apparent damage, caries or evidence of dental treatments.

Residual strain
During the preparation process, a low-speed diamond saw and polishing papers were used to minimize the effect of induced residual strain. The other origin of residual strain may be the natural growth of teeth. At low load values, there is a transition region in the transversely loading samples (#7 and #3), which might result from the residual tensile strain in the transverse direction. However, as also shown in Fig. 6b and c, the residual strain in the enamel was found to be very small in magnitude and was relieved subsequently at higher loading values. Since the elastic response has been experimentally observed to be dominating after the transition (reflected in the linearity of the experimental stress-strain curve), the presence of this initial residual strain only amounts to an initial offset. This means that there has been no impact from residual strain on the prediction of the elastic properties and that this offset can be neglected in the present analysis.

Normal strain components variation
In this experiment, the main difference between the three different samples was the orientation of the rods with respect to the external loading direction. For sample #6 in Fig. 6b, in which the load was applied along the rod longitudinal direction, the crystal lattice strain at a given stress was larger at 0°(i.e. along the loading direction) than at any other angles (which is a direct result of the Poisson effect). However, for the other two samples (#3 and #7), in which the rod orientation was perpendicular to the loading direction, the lattice strain was found to be larger at 90°than at other angles. These plots also demonstrate that loading along the transverse direction of the rod induces higher levels in stress or strain in the HAp crystals. In addition, the deviation from linearity was found in the three samples (45°for samples #6 and #7, and all the orientation angles for sample #3). Such relatively large experimental errors may result from the relatively weak intensity at the edge of the rings for samples #6 and #7, while for sample #3, the whole ring demonstrates a relatively weak intensity.
The plots of applied stress against lattice strain also have a common trend, in that the linear relationship starts only after a certain applied stress has been reached, especially in samples #3 and #7 (shown in Fig. 6a and c). This initial trend may demonstrate that the HAp crystals in the samples were initially subjected to a hydrostatic stress condition. In any case, this effect appears to have only minor influence on the further response.

The elastic anisotropy of the HAp crystals
The different mechanical responses of HAp crystals with respect to the loading direction shown in Fig. 6 and 7 originate not only from the different rod directions, but also in the strong anisotropic stiffness of HAp crystals. The 3-D directional modulus of a single HAp crystal was calculated by transforming the transversely isotropic stiffness matrix of a single crystal with different rotation matrices. The result is shown in Fig. 8a. The non-spherical shape corresponds to the elastic anisotropy of the HAp crystal. In Fig. 8a, the x-axis represents the longitudinal direction of the needle-shaped HAp crystal. Fig. 8b demonstrates the projections of the resulting modulus on the x-y plane and y-z plane. The projection on the x-y plane is similar to the cases of samples #6 and #7, and the circular projection on the y-z plane is similar to the case of sample #3 in which an isotropic modulus is observed. The orientation with the largest stiffness is found to be along the x-axis, as expected (longitudinal direction of the HAp crystal). However, the most compliant orientation observed in Fig. 8b is not at 90°, but rather around 60°due to the transverse isotropic stiffness of HAp crystals.

Conclusions
In this study, the lattice strain variation related to the nanoscale HAp distribution of human enamel was measured during in situ elastic compression by the combination of synchrotron WAXS and photoelasticity techniques. This study of the nanoscale phenomena and their effect on the macroscopic and microscopic mechanical behaviour provided access to information on both the structural and mechanical aspects of the sample, and allowed progress to be made in understanding the structure-property relationships in the hierarchical biomaterial of human enamel. In addition, as an improvement to an earlier proposed composite model [38,39], a multiscale Eshelby inclusion model was established and modified by introducing the misorientation distribution of HAp crystals at the nanoscale level. This second-level effect cannot be ignored since the mechanical behaviour shows strong dependence on the crystallite orientation distribution. The models were validated by the good agreement observed between the measured and calculated normal lattice strain component variation in different azimuthal directions.
This systematic experimental and modelling approach reported here is able to capture the complete picture of the multiscale structure of human enamel and its evolution under loading. The parameter refinement and validation in the modelling adopted in the present simulation offers the potential identification of nanoscale parameters, which may be hard to determine otherwise. By combining the results with previous studies of human dentine and enamel [12,19], an improved and comprehensive understanding of the multiscale structural-mechanical properties within human dental tissue can be achieved. Such information is essential for developing better prosthetic materials and dental fillings and has the potential to shed light on the mechanical evolution associated with the multiscale structural changes induced in human teeth by disease and treatment. Finally, this approach also enables the characterization of the structure-property relationship of other hard hierarchical biomaterials.

A.3. Enamel first-level model
The first-level model considers multiple aligned rods as the inhomogeneous inclusion and the protein as the surrounding matrix. In an elastic problem, no misfit strain exists in the inhomogeneity. Thus according to Eq. (A5): where C rod is the stiffness of the rods. Further, in terms of Eqs. (A6)-(A8), the rod stress can be determined as shown in Eq. 4, consisting of two parts, namely the contribution from the elastic anisotropy and that from the applied stress.

A.4. Enamel second-level model
The second-level model considers each rod as a composite and the HAp crystals are multiply aligned inhomogeneous inclusions. The obtained rod stress serves as the external load in the second level. Similar to Eq. (A9): where hCi HAp is the average stiffness of the aligned HAp crystals. The average lattice strain in HAp crystals consists of three parts: hei HAp where hei i is the contribution from the elastic anisotropy, hei matrix is due to the finite matrix effect and e A is the contribution from the